Interpolation d'une fonction pubic.ipynb

Kernel: SageMath 7.6

%display latex

On commence avec la fonction suivante : $$f(x) = \sin(3x) + e^x$$

x = var('x')
f = sin(3*x) + e^x;
p1 = plot(f,(-2,2)); p1

Ensuite on trouvera $5$ pointes pour notre fonction: $$f(0),\quad f\left(\frac{1}{4}\right),\quad f\left(\frac{1}{2}\right),\quad f\left(\frac{3}{4}\right), \quad f(1)$$

for i in range(5):
    print "f(x_"+str(i)+") = f("+str(float(i/4))+") = " + str(f(float(i/4)))

f(x_0) = f(0.0) = 1
f(x_1) = f(0.25) = 1.96566417671108
f(x_2) = f(0.5) = 2.64621625730418
f(x_3) = f(0.75) = 2.89507321350060
f(x_4) = f(1.0) = 2.85940183651891

p2 = scatter_plot([(i/4,f(i/4)) for i in range(5)]); p2 + p1

On commence avec nos équations :

$$\begin{align*} a_0 + a_1x_0 + a_2x_0^2+ a_3x_0^3 + a_4x_0^4 &= f(x_0)\\ a_0 + a_1x_1 + a_2x_1^2+ a_3x_1^3 + a_4x_1^4 &= f(x_1)\\ a_0 + a_1x_2 + a_2x_2^2+ a_3x_2^3 + a_4x_2^4 &= f(x_2)\\ a_0 + a_1x_3 + a_2x_3^2+ a_3x_3^3 + a_4x_3^4 &= f(x_3)\\ a_0 + a_1x_4 + a_2x_4^2+ a_3x_4^3 + a_4x_4^4 &= f(x_4) \end{align*}$$

On remplace les $x_i$ par les vrais nombres :

$$\begin{align*} a_0 &= 1\\ a_0 + a_1(0,25) + a_2(0,25)^2+ a_3(0,25)^3 + a_4(0,25)^4 &= 1.96566417671108\\ a_0 + a_1(0,5) + a_2(0,5)^2+ a_3(0,5)^3 + a_4(0,5)^4 &= 2.64621625730418\\ a_0 + a_1(0,75) + a_2(0,75)^2+ a_3(0,75)^3 + a_4(0,75)^4 &= 2.89507321350060\\ a_0 + a_1 + a_2+ a_3 + a_4 &= 2.85940183651891 \end{align*}$$

En forme matricielle $AX = B$, on a :

$$\begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 1 & (0,25) & (0,25)^2 & (0,25)^3 & (0,25)^4 \\ 1 & (0,5) & (0,5)^2 & (0,5)^3 & (0,5)^4 \\ 1 & (0,75) & (0,75)^2 & (0,75)^3 & (0,75)^4 \\ 1 & 1 & 1 & 1 & 1 \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ a_4 \end{bmatrix} = \begin{bmatrix} 1 \\ 1.96566417671108 \\ 2.64621625730418 \\ 2.89507321350060 \\ 2.85940183651891 \end{bmatrix}$$

A = matrix([[(i/4)^j for j in range(5)] for i in range(5)])
B = matrix([[f(round(float(i/4),4))] for i in range(5)])
AB = block_matrix([[A,B]], subdivid=False)
AB

r = AB.rref()
r

P = r[0][5] + r[1][5] * x + r[2][5] * x^2 + r[3][5] * x^3 + r[4][5] * x^4
P

p3 = plot(P,(-1,2),color='green',linestyle='--'); p3

x = p1 + p2 + p3; x

x.xmin(0); x.xmax(1); x.ymin(1); x.ymax(3); x