>

Interpolation polynomiale

> X:=[1,2,3,4,5]: Y:=[2,-2,0,4,5];

Y := [2, -2, 0, 4, 5]

> P:=interp(X,Y,t);

P := -1/24*t^4-1/4*t^3+133/24*t^2-73/4*t+15

> eval(P,t=1);

2

> P:=unapply(P,t);

P := proc (t) options operator, arrow; -1/24*t^4-1/4*t^3+133/24*t^2-73/4*t+15 end proc

> map(P,X);

[2, -2, 0, 4, 5]

> map(P,[1.5,2.5]);

[-.96093750, -1.52343750]

> plot(P,-1..6);

[Plot]

> A:=zip((x,y)->[x,y],X,Y);

A := [[1, 2], [2, -2], [3, 0], [4, 4], [5, 5]]

> plot(A, style=point, symbol=circle);

[Plot]

On rassemble les deux (par la fonction plots)

> plots[display](plot(A, style=point, symbol=circle),plot(P,-1..6));

[Plot]

>

Automatiser le tout dans une procédure (qui prend deux vecteurs X et Y et trace leur polynôme d'interpolation ainsi que le nuage de points).

> dessin:=proc(X,Y)
local A,x1,x2,P,pas;

P:=unapply(interp(X,Y,t),t);

x1:=min(op(X));x2:=max(op(X));pas:=(x2-x1)/10;

A:=zip((x,y)->[x,y],X,Y);

plots[display](plot(A, style=point, symbol=circle),plot(P,x1-pas..x2+pas));

end proc;

dessin := proc (X, Y) local A, x1, x2, P, pas; P := unapply(interp(X, Y, t), t); x1 := min(op(X)); x2 := max(op(X)); pas := 1/10*x2-1/10*x1; A := zip(proc (x, y) options operator, arrow; [x, y] end pr...dessin := proc (X, Y) local A, x1, x2, P, pas; P := unapply(interp(X, Y, t), t); x1 := min(op(X)); x2 := max(op(X)); pas := 1/10*x2-1/10*x1; A := zip(proc (x, y) options operator, arrow; [x, y] end pr...dessin := proc (X, Y) local A, x1, x2, P, pas; P := unapply(interp(X, Y, t), t); x1 := min(op(X)); x2 := max(op(X)); pas := 1/10*x2-1/10*x1; A := zip(proc (x, y) options operator, arrow; [x, y] end pr...dessin := proc (X, Y) local A, x1, x2, P, pas; P := unapply(interp(X, Y, t), t); x1 := min(op(X)); x2 := max(op(X)); pas := 1/10*x2-1/10*x1; A := zip(proc (x, y) options operator, arrow; [x, y] end pr...dessin := proc (X, Y) local A, x1, x2, P, pas; P := unapply(interp(X, Y, t), t); x1 := min(op(X)); x2 := max(op(X)); pas := 1/10*x2-1/10*x1; A := zip(proc (x, y) options operator, arrow; [x, y] end pr...dessin := proc (X, Y) local A, x1, x2, P, pas; P := unapply(interp(X, Y, t), t); x1 := min(op(X)); x2 := max(op(X)); pas := 1/10*x2-1/10*x1; A := zip(proc (x, y) options operator, arrow; [x, y] end pr...dessin := proc (X, Y) local A, x1, x2, P, pas; P := unapply(interp(X, Y, t), t); x1 := min(op(X)); x2 := max(op(X)); pas := 1/10*x2-1/10*x1; A := zip(proc (x, y) options operator, arrow; [x, y] end pr...dessin := proc (X, Y) local A, x1, x2, P, pas; P := unapply(interp(X, Y, t), t); x1 := min(op(X)); x2 := max(op(X)); pas := 1/10*x2-1/10*x1; A := zip(proc (x, y) options operator, arrow; [x, y] end pr...dessin := proc (X, Y) local A, x1, x2, P, pas; P := unapply(interp(X, Y, t), t); x1 := min(op(X)); x2 := max(op(X)); pas := 1/10*x2-1/10*x1; A := zip(proc (x, y) options operator, arrow; [x, y] end pr...

> dessin([1,3,4,5],[-2,-6,2,6]);

[Plot]

>

Interpolation d'une fonction.

On interpole une fonction f de la façon suivante. Soit [a,b] un  intervalle donné, on prend n points dans cet intervalle (ce nombre n fait partie des données), x1, ..., xn et on calcule f(x1), ..., f(xn). Puis, on calcule le polynome d'interpolation en ces points. Le mieux est de prendre les points équidistants dans l'intervalle de la façon suivante:

> Equidist:=(a,b,n)->[seq(a+k*(b-a)/(n-1),k=0..n-1)];

Equidist := proc (a, b, n) options operator, arrow; [seq(a+k*(b-a)/(n-1), k = (0 .. n-1))] end proc

> Equidist(0,2,5);

[0, 1/2, 1, 3/2, 2]

> Intp:=proc(f,X)
local t,Y,P;

Y:=map(f,X);

P:=unapply(interp(X,Y,t),t);

end;

Intp := proc (f, X) local t, Y, P; Y := map(f, X); P := unapply(interp(X, Y, t), t) end proc

> Intp(cos,Equidist(0,2*pi,5));

proc (t) options operator, arrow; 1+1/6*(4+24*cos(pi)-16*cos(1/2*pi)+4*cos(2*pi)-16*cos(3/2*pi))*t^4/pi^4+1/6*(-20*pi-12*cos(2*pi)*pi+56*cos(3/2*pi)*pi+72*cos(1/2*pi)*pi-96*cos(pi)*pi)*t^3/pi^4+1/6*(1...proc (t) options operator, arrow; 1+1/6*(4+24*cos(pi)-16*cos(1/2*pi)+4*cos(2*pi)-16*cos(3/2*pi))*t^4/pi^4+1/6*(-20*pi-12*cos(2*pi)*pi+56*cos(3/2*pi)*pi+72*cos(1/2*pi)*pi-96*cos(pi)*pi)*t^3/pi^4+1/6*(1...

> plot({cos,Intp(cos,Equidist(0,7,5))},0..7);

[Plot]

Attention, dans cet exemple, on peut voir une certaine imprécision aux bornes (-1 et 1) de la fonction d'interpolation (même si évidemment elle est fidèle aux points d'interpolation. Le deuxième exemple accentue le phénomène (essayez avec d'autres valeurs)
.

> f:=t->1/(1+8*t^2): plot({f,Intp(f,Equidist(-1,1,6))},-1..1);

[Plot]

> f:=t->1/(1+8*t^2): plot({f,Intp(f,Equidist(-1,1,12))},-1..1);

[Plot]

>