Revision history [back]
 | 1 | initial version |
Here is my solution more involved than the one of dsejas, but perhaps it can help (Thanks to Sebastien)
# Polynomial interpolation
# Passing the matrix
A = Matrix([[0,0],[0.5,0.25],[0.6,0.36],[2,4]])
var('a0, a1, a2, a3')
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 for coefficients
sol=solve([ssr_a0==0,ssr_a1==0,ssr_a2==0,ssr_a3==0,ssr_a4==0], a0, a1, a2, a3, solution_dict=True)
f(x)= a0 + a1*x+ a2*x^2 + a3*x^3
# g is a the desired function
g(x)=f(x).subs(sol[0])
# derivative of g
gp=g.diff(x)
pl=plot(g(x),(x,A[0][0],A[3][0]))
pl+=points((A[i] for i in range(4)),color='red',size=20)
pl+=plot(gp(x),(x,A[0][0],A[3][0]),color="green")
pl+=text("$f(x)$",(1.2,1),color="blue", fontsize='small', rotation=0)
pl+=text("$f'(x)$",(1.2,2.8),color="green", fontsize='small', rotation=0)
show(pl)
| | 2 | No.2 Revision |
Here is my solution more involved than the one of dsejas, but perhaps it can help (Thanks to Sebastien)
# Polynomial interpolation
# Passing the matrix
A = Matrix([[0,0],[0.5,0.25],[0.6,0.36],[2,4]])
var('a0, a1, a2, a3')
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 for coefficients
sol=solve([ssr_a0==0,ssr_a1==0,ssr_a2==0,ssr_a3==0,ssr_a4==0], a0, a1, a2, a3, solution_dict=True)
f(x)= a0 + a1*x+ a2*x^2 + a3*x^3
# g is a the desired function
g(x)=f(x).subs(sol[0])
# derivative of g
gp=g.diff(x)
pl=plot(g(x),(x,A[0][0],A[3][0]))
pl+=points((A[i] for i in range(4)),color='red',size=20)
pl+=plot(gp(x),(x,A[0][0],A[3][0]),color="green")
pl+=text("$f(x)$",(1.2,1),color="blue", fontsize='small', rotation=0)
pl+=text("$f'(x)$",(1.2,2.8),color="green", fontsize='small', rotation=0)
show(pl)
