# Error in polynomial interpolation

I have edited this question to insist on what goes wrong

The following code ends with an indetermined parameter.

# Polynomial interpolation
# Passing the matrix
A = Matrix([[0,0],[0.5,0.25],[0.6,0.36],[2,4]])
var('a0, a1, a2, a3, a4')
B =[[A[i][1], a3*A[i][0]^3+ a2*A[i][0]^2+ a1*A[i][0]+ a0] for i in range(0,4)]
C=[[(B[i][0]-B[i][1])^2] for i in range(0,4)]
# sum of the square of errors
ssr=sum(C[i][0] for i in range(4))
# First order condition of optimality
ssr_a0=ssr.diff(a0)
ssr_a1=ssr.diff(a1)
ssr_a2=ssr.diff(a2)
ssr_a3=ssr.diff(a3)
ssr_a4=ssr.diff(a4)
# solution
sol=solve([ssr_a0==0,ssr_a1==0,ssr_a2==0,ssr_a3==0,ssr_a4==0], a0, a1, a2, a3, a4, solution_dict=True)
sol


The proposed solution being [{a0: 0, a1: 0, a2: 1, a3: 0, a4: r3}].

But the ssr_ai $i = 0,\ldots, 3$ are all linear. So sol should gives the same result that the following code

C=matrix([[8, 6.2, 9.22, 16.682], [6.2, 9.22,16.68,32.3842],[9.22,16.682,32.3842,64.21802],[16.68,32.3842,64.21802,128.124562]])
b=vector([9.22,16.68,32.3842,64.21902])
C.solve_right(b)


So why it is not the case ? (There is one and only one solution)

This lead to a second question : for the ssr_ai how can I select the parameters to construct the C matrix and the b vector not by hand as I have been obliged to do ?

edit retag close merge delete

1

( 2020-04-30 13:02:04 +0100 )edit

I should obtain a value for a0, a1, a3, and a4. So what is this result a4 : r3 ? What means r3 ? I already encounter such a result but I have forgoten how to resolve it.

( 2020-04-30 17:21:47 +0100 )edit

r3 is a parameter in the space of solution. In the documentation of solve, you can see this example: "If there is a parameter in the answer, that will show up as a new variable. In the following example, "r1" is an arbitrary constant (because of the "r")":

sage: forget()
sage: x, y = var('x,y')
sage: solve([x+y == 3, 2*x+2*y == 6],x,y)
[[x == -r1 + 3, y == r1]]

( 2020-04-30 18:03:52 +0100 )edit

Thanks Sebastien for your comment. But this is weird. My system is linear so either there is no solution either there is one. I have tried to substitute find_root() to solve but I suspect that it is only programmed to solve one variable equation. The solution would be to write the system in matrix notation but I do not know how to select the parameters to wtrite a matrix.

( 2020-05-01 07:00:28 +0100 )edit

"But this is weird. My system is linear so either there is no solution either there is one." This is not true as a system of linear equations can have 0, 1 or an infinity of solutions. See the example in my previous comment.

( 2020-05-01 09:34:10 +0100 )edit

Sort by » oldest newest most voted

If you look at the system of linear equations, you can see that it really has 4 equations with 4 unknowns (the 5th equation is 0=0 and the fifth unknown a4 just doesn't appear)

sage: [ssr_a0==0,ssr_a1==0,ssr_a2==0,ssr_a3==0,ssr_a4==0]
[8*a0 + 6.20*a1 + 9.220000*a2 + 16.68200*a3 - 9.2200000 == 0,
6.20000*a0 + 9.2200000*a1 + 16.6820000*a2 + 32.3842000*a3 - 16.682000 == 0,
9.2200*a0 + 16.6820000*a1 + 32.384200*a2 + 64.21802000*a3 - 32.38420000 == 0,
16.68200*a0 + 32.384200*a1 + 64.2180200*a2 + 128.12456200*a3 - 64.2180200 == 0,
0 == 0]


That system has a solution involving only the first four unknown, which explains what you obtain:

[{a0: 0, a1: 0, a2: 1, a3: 0, a4: r3}]


Indeed, we can see that the coefficient of a2 matches the constant terms. Alternatively, you can just remove a4 in the call to solve since that variable does not appear in the system:

sage: solve([ssr_a0==0,ssr_a1==0,ssr_a2==0,ssr_a3==0,ssr_a4==0], a0, a1, a2, a3, solution_dict=True)
[{a0: 0, a1: 0, a2: 1, a3: 0}]


EDIT (answering a question in the comments): To transform the system of linear equations to matrices and solve the system, you may do:

sage: SSR = [ssr_a0,ssr_a1,ssr_a2,ssr_a3,ssr_a4]
sage: A = matrix([[ssr_a.coefficient(a) for a in [a0,a1,a2,a3]] for ssr_a in SSR])
sage: b = vector([-ssr_a.subs(a0=0,a1=0,a2=0,a3=0) for ssr_a in SSR])
sage: A \ b
(-4.81096644000283e-16, 8.08282642120849e-14, 0.999999999999817, 7.15728811560799e-14)
sage: A.change_ring(QQ) \ b.change_ring(QQ)  # or in the ring of rationals
(0, 0, 1, 0)

more