>

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*(56*cos(3/2*pi)*pi-20*pi+72*cos(1/2*pi)*pi-12*cos(2*pi)*pi-96*cos(pi)*pi)*t^3/pi^4+1/6*(-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*(56*cos(3/2*pi)*pi-20*pi+72*cos(1/2*pi)*pi-12*cos(2*pi)*pi-96*cos(pi)*pi)*t^3/pi^4+1/6*(-5...

> 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]

>

Interpolation de lagrange.

Ecrire une procédure donnant le polynome d'interpolation de Lagrange de n points donnés.

> Lagrange:=proc(X,Y)
local P,n,L,k;

n:=nops(X); P:=0;

for k from 1 to n do

P:=P+Y[k]*product((x-X[j])/(X[k]-X[j]),j=1..k-1)*product((x-X[j])/(X[k]-X[j]),j=k+1..n);

od;

unapply(collect(P,x),x);

end proc;

Lagrange := proc (X, Y) local P, n, L, k; n := nops(X); P := 0; for k to n do P := P+Y[k]*product((x-X[j])/(X[k]-X[j]), j = (1 .. k-1))*product((x-X[j])/(X[k]-X[j]), j = (k+1 .. n)) end do; unapply(co...Lagrange := proc (X, Y) local P, n, L, k; n := nops(X); P := 0; for k to n do P := P+Y[k]*product((x-X[j])/(X[k]-X[j]), j = (1 .. k-1))*product((x-X[j])/(X[k]-X[j]), j = (k+1 .. n)) end do; unapply(co...Lagrange := proc (X, Y) local P, n, L, k; n := nops(X); P := 0; for k to n do P := P+Y[k]*product((x-X[j])/(X[k]-X[j]), j = (1 .. k-1))*product((x-X[j])/(X[k]-X[j]), j = (k+1 .. n)) end do; unapply(co...Lagrange := proc (X, Y) local P, n, L, k; n := nops(X); P := 0; for k to n do P := P+Y[k]*product((x-X[j])/(X[k]-X[j]), j = (1 .. k-1))*product((x-X[j])/(X[k]-X[j]), j = (k+1 .. n)) end do; unapply(co...Lagrange := proc (X, Y) local P, n, L, k; n := nops(X); P := 0; for k to n do P := P+Y[k]*product((x-X[j])/(X[k]-X[j]), j = (1 .. k-1))*product((x-X[j])/(X[k]-X[j]), j = (k+1 .. n)) end do; unapply(co...Lagrange := proc (X, Y) local P, n, L, k; n := nops(X); P := 0; for k to n do P := P+Y[k]*product((x-X[j])/(X[k]-X[j]), j = (1 .. k-1))*product((x-X[j])/(X[k]-X[j]), j = (k+1 .. n)) end do; unapply(co...Lagrange := proc (X, Y) local P, n, L, k; n := nops(X); P := 0; for k to n do P := P+Y[k]*product((x-X[j])/(X[k]-X[j]), j = (1 .. k-1))*product((x-X[j])/(X[k]-X[j]), j = (k+1 .. n)) end do; unapply(co...

> Lagrange([1,2,3,4],[1,-1,0,4]); interp([1,2,3,4],[1,-1,0,4],x);

proc (x) options operator, arrow; 6+3/2*x^2-13/2*x end proc

6+3/2*x^2-13/2*x

> P:=Lagrange([1,3,4,2],[4,2,5,1]);

P := proc (x) options operator, arrow; 13-1/3*x^3+4*x^2-38/3*x end proc

Sur le modèle déjà vu plus haut, on crèe une procédure permettant d'interpoler une fonction par la méthode de Lagrange que l'on vient de programmer puis on la compare avec la méthode présente.

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

Y:=map(f,X);

P:=Lagrange(X,Y);

end;

IntpLAGR := proc (f, X) local Y, P; Y := map(f, X); P := Lagrange(X, Y) end proc

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

[Plot]

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

[Plot]

>

Interpolation : méthode de Newton

Donner une procédure calculant le polynome de Newton interpolant n points donnés. Le transformer, là encore, en une procédure permettant d'interpoler une fonction f quelconque. On pourra commencer par implémenter l'opérateur de différences divisées.

On peut très bien donner un algorithme récursif (très simple) qui "paraphrase" la définition des différences divisées. On va plutôt donner l'algorithme qui calcule la table des différences divisées (et donc nous donner tous les termes).

> ListeDiffDiv:=proc(X,Y)
local L,i,j,n;

n:=nops(X); L:=Y;

for i from 1 to n-1 do

for j from n to i+1 by -1 do

L[j]:=(L[j]-L[j-1])/(X[j]-X[j-i]);

od;

od;

return L;

end;

ListeDiffDiv := proc (X, Y) local L, i, j, n; n := nops(X); L := Y; for i to n-1 do for j from n by -1 to i+1 do L[j] := (L[j]-L[j-1])/(X[j]-X[j-i]) end do end do; return L end procListeDiffDiv := proc (X, Y) local L, i, j, n; n := nops(X); L := Y; for i to n-1 do for j from n by -1 to i+1 do L[j] := (L[j]-L[j-1])/(X[j]-X[j-i]) end do end do; return L end procListeDiffDiv := proc (X, Y) local L, i, j, n; n := nops(X); L := Y; for i to n-1 do for j from n by -1 to i+1 do L[j] := (L[j]-L[j-1])/(X[j]-X[j-i]) end do end do; return L end procListeDiffDiv := proc (X, Y) local L, i, j, n; n := nops(X); L := Y; for i to n-1 do for j from n by -1 to i+1 do L[j] := (L[j]-L[j-1])/(X[j]-X[j-i]) end do end do; return L end procListeDiffDiv := proc (X, Y) local L, i, j, n; n := nops(X); L := Y; for i to n-1 do for j from n by -1 to i+1 do L[j] := (L[j]-L[j-1])/(X[j]-X[j-i]) end do end do; return L end procListeDiffDiv := proc (X, Y) local L, i, j, n; n := nops(X); L := Y; for i to n-1 do for j from n by -1 to i+1 do L[j] := (L[j]-L[j-1])/(X[j]-X[j-i]) end do end do; return L end procListeDiffDiv := proc (X, Y) local L, i, j, n; n := nops(X); L := Y; for i to n-1 do for j from n by -1 to i+1 do L[j] := (L[j]-L[j-1])/(X[j]-X[j-i]) end do end do; return L end proc

> ListeDiffDiv([x0,x1],[y0,y1]);

[y0, (y1-y0)/(x1-x0)]

> ListeDiffDiv([x0,x1,x2,x3],[y0,y1,y2,y3]);

[y0, (y1-y0)/(x1-x0), ((y2-y1)/(x2-x1)-(y1-y0)/(x1-x0))/(x2-x0), (((y3-y2)/(x3-x2)-(y2-y1)/(x2-x1))/(x3-x1)-((y2-y1)/(x2-x1)-(y1-y0)/(x1-x0))/(x2-x0))/(x3-x0)]

> ListeDiffDiv([1,3,4,2],[4,2,5,1]);

[4, -1, 4/3, (-1)/3]

> PolynomeNewton:=proc(X,Y)
local P,i,j,coeff;

P:=0;

coeff:=ListeDiffDiv(X,Y);

for i from 1 to nops(X) do

P:=P+coeff[i]*product(x-X[j],j=1..i-1);

od;

unapply(P,x);

end;

PolynomeNewton := proc (X, Y) local P, i, j, coeff; P := 0; coeff := ListeDiffDiv(X, Y); for i to nops(X) do P := P+coeff[i]*product(x-X[j], j = (1 .. i-1)) end do; unapply(P, x) end procPolynomeNewton := proc (X, Y) local P, i, j, coeff; P := 0; coeff := ListeDiffDiv(X, Y); for i to nops(X) do P := P+coeff[i]*product(x-X[j], j = (1 .. i-1)) end do; unapply(P, x) end procPolynomeNewton := proc (X, Y) local P, i, j, coeff; P := 0; coeff := ListeDiffDiv(X, Y); for i to nops(X) do P := P+coeff[i]*product(x-X[j], j = (1 .. i-1)) end do; unapply(P, x) end procPolynomeNewton := proc (X, Y) local P, i, j, coeff; P := 0; coeff := ListeDiffDiv(X, Y); for i to nops(X) do P := P+coeff[i]*product(x-X[j], j = (1 .. i-1)) end do; unapply(P, x) end procPolynomeNewton := proc (X, Y) local P, i, j, coeff; P := 0; coeff := ListeDiffDiv(X, Y); for i to nops(X) do P := P+coeff[i]*product(x-X[j], j = (1 .. i-1)) end do; unapply(P, x) end procPolynomeNewton := proc (X, Y) local P, i, j, coeff; P := 0; coeff := ListeDiffDiv(X, Y); for i to nops(X) do P := P+coeff[i]*product(x-X[j], j = (1 .. i-1)) end do; unapply(P, x) end procPolynomeNewton := proc (X, Y) local P, i, j, coeff; P := 0; coeff := ListeDiffDiv(X, Y); for i to nops(X) do P := P+coeff[i]*product(x-X[j], j = (1 .. i-1)) end do; unapply(P, x) end proc

> P:=PolynomeNewton([1,3,4,2],[4,2,5,1]);

P := proc (x) options operator, arrow; 5-x+4/3*(x-1)*(x-3)-1/3*(x-1)*(x-3)*(x-4) end proc

> plot(P);

[Plot]

>

On l'adapte pour interpoler une fonction.

> IntpNewton:=proc(f,X)
PolynomeNewton(X,map(f,X));

end proc;

IntpNewton := proc (f, X) PolynomeNewton(X, map(f, X)) end proc

> IntpNewton(exp,[1,2,3]);

proc (x) options operator, arrow; exp(1)+(exp(2)-exp(1))*(x-1)+(1/2*exp(3)-exp(2)+1/2*exp(1))*(x-1)*(x-2) end proc

>

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

[Plot]

>