mait1.html

Université de Lille 1

Maitrise de mathématiques

Analyse Numérique des EDP

J.-P. Chehab
2000 / 2001

Feuille illustrative 1

Projections dans

Approximation au sens des moindres carrés

Le problème

Soit f une fonction continue sur l'intervalle [a,b]. On désire déterminer le polynôme P de degré inférieur ou égal à n tel que la quantité

soit minimale. On montre qu'il existe un unique polynôme P ayant cette propriété et que P est caracterisé

par

(R) =0 pour tout q, polynôme de degré inférieur ou égal à n

Maintenant cherchons à déterminer P. P peut s'écrire sous la forme

P(x)=

Les monomes ei(x)= forment une base de l'espace vectoriel des polynômes de degré inférieur ou égal à n, si bien que la relation (R) est équivalente à

(R') = = pour tout j=0,...n .

Les sont donc solution du système linéaire

H.a = F

où a=( ,..., ), F=( ,..., ) et où H est la matrice (symétrique) de coefficients

Si l'intervalle considéré est [0,1] alors

= =

H est la matrice de Hilbert. On montre que cond(H) est équivalent à lorsque n tend vers + , ce qui est tres mauvais.

Une Illustration

Construisons la matrice de Hilbert pour une valeur donnée de n.

> restart;

> with(linalg):

> n := 4;

Warning, new definition for norm

Warning, new definition for trace

> f:=(i,j)->1/(i+j-1);

> Hil:=matrix(n,n,f);

Calculons le nombre de condition et l'inverse de Hil.

> evalc(cond(Hil));

> inverse(Hil);

Remarque : en fait, on peut manier la matrice de Hilbert en Maple sans avoir à la construire. Pour cela on peut utiliser la commande hilbert (n) , où n est la dimension.

Application.

Détérminons le polynôme de degré inferieur ou égal à 6 qui approche la fonction

f(x)= sur l'intervalle [0,1] au sens des moindres carrés.

Construction du second membre

> Fi:=i->int(exp(y)*y^(i-1),y=0..1);

> F:=vector(n,Fi);

Résolution du système

> Ai:=linsolve(hilbert(n),F);

Comparons le polynôme obtenu avec la fonction de départ.

> Pol:=sum('Ai[k]*y^(k-1)','k'=1..n);

> plot([Pol,exp(y)],y=0..1);

Le polynôme a été déterminé en intégrant exactement les termes du second membre du système

(le vecteur F) ; cela est possible pour cet exemple. A présent examinons ce qui se passe si le second membre est intégré par une méthode numérique.

Cela peut se réaliser en Maple à l'aide de l'instruction

evalf(int(fonction_a_integrer,y=0..1,ndigit);

qui retournera une valeur de l'intégrale avec ndigit chiffres significatifs. En choisissant différentes valeurs de ndigit , on peut voir une traduction du mauvais conditionnement de la matrice de Hilbert en comparant le polynôme obtenu avec la fonction de depart. A partir de quelle valeur de ndigit les courbes semblent-elles confondues ?

Comme auparavant, on construit d'abord le second membre :

> Gi:=i->evalf(int(exp(y)*y^(i-1),y=0..1),3);

Resolvons le (nouveau) système

> G:=vector(n,Gi);

> Bi:=linsolve(hilbert(n),G);

> Pol2:=sum('Bi[k]*y^(k-1)','k'=1..n);

> plot([Pol2,exp(y)],y=0..1);

Packages utilisés : linalg

Remarque

Plus généralement

On peut énoncer le problème d'approximation au sens des moindres carrés par :

Déterminer le polynôme P de degré inférieur ou égal à n tel que la quantité

(*)

soit minimale. Ici est une fonction intégrable, positive p.p. Dans le cas précédent nous avions pris . On appelle "la fonction poids".

Les polynômes de Tchebytcheff

Prenons et [a,b]=[-1,1].

Ces polynômes sont orthogonaux par rapport au produit scalaire à poids (*).

> restart;

> with(orthopoly):

> pdt_scal_poids:=proc(p,q)
local pdt;
pdt:=int(p*q/sqrt(1-x^2),x=-1..1);
RETURN(pdt);
end;

Considérons deux polynômes de Tchebytcheff :

> pdt_scal_poids(T(4,x),T(9,x));

Reprenons l'exercice précédent en cherchant à minimiser la norme à poids et en exprimant le polynôme dans la base de Tchebytcheff :

> with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> construction_matrice:=proc(n)#version brutale sans exploiter l'orthogonalité
local mat,i,j;
mat:=matrix(n+1,n+1,0):
for i from 0 to n do
for j from 0 to n do
mat[i+1,j+1]:=pdt_scal_poids(T(i,x),T(j,x));
#mat[i+1,j+1]:=evalf(pdt_scal_poids(T(i,x),T(j,x)),3);
od;
od;
RETURN(evalm(mat));
end;

> n:=4;

> A:=construction_matrice(n);

Calculons le second membre

> Fi:=i->pdt_scal_poids(exp(x),T(i,x));

> F:=vector(n+1,Fi);

> Ti:=linsolve(A,F);

> PolT:=sum('Ti[k]*T(k-1,x)','k'=1..n);

> plot([PolT,exp(y)],y=0..1);

Représentation de fonctions

Pour construitre et représenter graphiquement des fonctions de pour n=1 et n=2, on utilise la commande plot pour n=1 et plot3d pour n=2.

Représentation graphique de fonctions de

n=1

> restart;

> f:=x->exp(1/(x^2-1))*piecewise( -1 <=x and x <=1,1,0);

> plot(f(x),x=-2..2);

n=2

> epsilon:=10^(-6);

> g:=(x,y)->f(x)*f(y);

> g(x,y);

> plot3d(g(x,y),x=-1+epsilon..1-epsilon,y=-1+epsilon..1-epsilon);

> h:=(x,y)->exp(1/(x^2+y^2-1))*piecewise(x^2+y^2<=1,1,0);

> plot3d(h(x,y),x=-1+epsilon..1-epsilon,y=-1+epsilon..1-epsilon);

Un autre façon de construire des fonction de D(R^n) est de composer une des deux fonctions présentées ci-dessus avec une fonction régulière s'annulant en 0. Exemple

> plot3d(sin(6*Pi*h(x,y)),x=-1+epsilon..1-epsilon,y=-1+epsilon..1-epsilon);

>

Représentation graphique de fonctions de non continues

Exemple 1

> restart;

> f:=(x,y)->log(abs(log(x^2+y^2)));

On a vu en cours ou en TD que cette fonction est dans H1(Omega) lorsque Omega est le disque ouvert de centre 0 et de rayon 1/2 mais qu'elle n'admet pas de représentant continu, elle est en effet non bornée en l'origine. Ici, on représente graphiquement cette fonction et on calcule sa norme H1.

Passons en polaires.

Calculons la norme L2 au carré de cette fonction

> nl2:=2*Pi*int(subs(x^2+y^2=rho,simplify(f(x,y)^2))*rho,rho=0..1/2);