TD Maple n6 : Méthodes d'intégration numérique

Corrigé

Lundis 28 janvier et 4 février 2008

1 Calculs de primitives et d'intégrales

1. Trivial.
2. On charge la librairie student qui permet de faire des intégrations par parties.
[>with(student) ;

Z x
1
In (x) = dt
0 (1 + t2 )n
Z x
1 + t2 − t2
= dt
0 (1 + t2 )n
Z x Z x
1 t2
= 2 n−1
dt − 2 n
dt
0 (1 + t ) 0 (1 + t )
2t 21 t
Z x
= In−1 (x) − 2 n
dt
0 (1 + t )

Puis on fait une IPP en dérivant u(t) = 21 t. Avec Maple cela donne :

[>Integrale :=Int(tˆ2/(1+tˆ2)ˆn,t=0..x) ;
[>intparts(Integrale,t/2) ;
On obtient alors après quelques réarangements que
" #x
 1  1 1 1
In (x) = 1 + In−1 (x) − t
2(1 − n) 1 − n (1 + t2 )n−1 2
0

2 Intégration numérique

2.1 Méthode des rectangles

1. Voici la procédure permettant de mettre en oeuvre la méthode des rectangles à gauche :

[>rectgauche :=proc(f,a,b,N)
local i,h ;
h :=(b-a)/N ;
evalf(h*sum(f(a+i*h),i=0..N-1)) ;
end ;
2. Pour la méthode des rectangles à droite, la somme va de i = 1 à i = N :
b N
b−aX
Z
f (t)dt ≈ f (ti )
a N i=1

Ce qui donne avec Maple :

1
 [>rectdroite :=proc(f,a,b,N)
local i,h ;
h :=(b-a)/N ;
evalf(h*sum(f(a+i*h),i=1..N)) ;
end ;
3. rectgauche(f,a,b,N) surestime l'intégrale à coup sûr si la fonction f est décroissante sur [a, b].
4. Testons maintenant nos procédures.
[>int(exp,0..1) ;
e−1
[>evalf(%) ;
1.718281828
Puis on regarde ce que donnent nos procédures avec N = 10 intervalles de discrétisation.
[>rectgauche(exp,0,1,10) ;
1.633799400
[>rectdroite(exp,0,1,10) ;
1.805627583
5. Vérication de l'ordre de la méthode : On remplit une liste graphe dans laquelle on met les couples
[ln(1/N),ln(erreur)]. Pour cela, on parcourt la liste L des diérents pas de discrétisation.
[>L :=[10,20,30,40,50,60] ;
[>graphe :=[] ;
[>for N in L do
graphe :=[op(graphe),[ln(1/N),ln(abs(e-1-rectgauche(exp,0,1,N)))]] ;
od :
[>plot(graphe) ;
On obtient bien une droite de pente 1 pour les deux méthodes ce qui veut dire que l'erreur est propor-
tionnelle à la taille h du pas de discrétisation :
Z b
f (t)dt − rectgauche(f,a,b,N) = Ch
a

La constante C est l'exponentielle de l'ordonnée à l'origine de la droite obtenue. Elle ne dépend que de
la fonction f et de l'intervalle [a, b] mais pas de h.

2.2 Méthode du point milieu

1. Voici la procédure mettant en oeuvre la méthode du point milieu :

Cette fois-ci on approche l'intégrale par :
b N −1 N −1
b − a X  ti + ti+1  b − a X 
Z
(2i + 1)h 
f (t)dt ≈ f = f a+
a N i=0 2 N i=0 2

[>pointmilieu :=proc(f,a,b,N)
local i,h ;
h :=(b-a)/N ;
evalf(h*sum(f(a+(2*i+1)*h/2),i=1..N)) ;
end ;
[>pointmilieu(exp,0,1,10) ;
1.717566087
2. On trace le ln de l'erreur en fonction de lnh. On obtient un droite de pente 2. Ceci veut dire que l'erreur
commise est proportionnelle au carré du pas de discrétisation. Donc quand h est petit (N grand), cette
méthode est plus présice que la méthode des rectangles.
3. Si λ est diérent de 21 , on retombe sur une methode d'ordre 1. En eet, on montre que la méthode du
point milieu est la seule méthode à 1 point d'interpolation qui soit d'ordre 2.

2
2.3 Méthode des trapèzes

f (ti ) + f (ti+1 ) f (ti ) + f (ti+1 )

1. L'aire du trapèze délimité par la fonction ane sur [ti , ti+1 ] est (ti+1 −ti ) =h .
2 2
D'où une nouvelle formule d'approximation de l'intégrale de f sur [a, b] :
N −1
!
b
b−a
Z
f (a) + f (b) X
f (t)dt ≈ + f (ti )
a N 2 i=1

2. Voici la procédure qui met en ouvre la méthode des trapèzes :

[>trapeze :=proc(f,a,b,N)
local i,h ;
h :=(b-a)/N ;
evalf(h*((f(a)+f(b))/2+sum(f(a+i*h),i=1..N-1))) ;
end ;
[>trapeze(exp,0,1,10) ;
1.719713492
3. Toujours la même démarche, on obtient un méthode d'odre 2. Il est donc plus ecace d'interpoler la
fonction f par une fonction ane.

2.4 Méthode de Simpson

1. Ce polynôme doit s'annuler en ti+ 12 et ti+1 . Il s'écrit donc sous la forme :

λ(X − ti+ 12 )(X − ti+1 )

On choisit λ de sorte que le polynôme prenne la valeur 1 en ti , donc λ vérie :

1 = λ(ti − ti+ 21 )(ti − ti+1 )

Finalement, le polynôme recherché est :

(X − ti+ 21 )(X − ti+1 )
(ti − ti+ 21 )(ti − ti+1 )

Si on veut que ce polynôme prenne la valeur f (ti ) en ti , il sut de le multiplier par f (ti ).
2. Le polynôme de degré 2 qui coïncide avec f aux points ti , ti+1 et ti+ 12 est la somme des trois polynômes
obtenus comme ci-dessus pour les trois points soit :

(X − ti+ 12 )(X − ti+1 ) (X − ti )(X − ti+1 ) (X − ti )(X − ti+ 21 )

P (X) = f (ti ) + f (ti+ 12 ) + f (ti+1 )
(ti − ti+ 21 )(ti − ti+1 ) (ti+ 21 − ti )(ti+ 12 − ti+1 ) (ti+1 − ti )(ti+1 − ti+ 12 )

3. A l'aide de Maple, on intègre ce polynôme entre ti et ti+1 . Pour simplier les notations, on notera c = ti ,
d = ti+1 . On exprime le polynôme :
[>p :=F(c)*(x-(c+d)/2)*(x-d)/((c-(c+d)/2)*(c-d))
+F((c+d)/2)*(x-c)*(x-d)/(((c+d)/2-c)*((c+d)/2-d))
+F(d)*(x-c)*(x-(c+d)/2)/((d-c)*(d-(c+d)/2)) ;

Puis on le développe :
[>expand(%) ;

Puis on calcule l'aire sous la parabole :

[>int(p,x=c..d) ;
[>simplify(%) ;

3
 [>factor(%) ;

Finalement, on voit que l'aire sous la parabole est donnée par :

f (ti ) + 4f (ti+ 21 ) + f (ti+1 ) f (ti ) + 4f (ti+ 12 ) + f (ti+1 )
(ti+1 − ti ) =h
6 6
4. Ce qui donne pour le calcul approché de l'intégrale :
Z b N −1 N −1    
h h X 2h X (2i + 1)h 2h h
f (t)dt ≈ (f (a) + f (b)) + f (ti ) + f a+ + f a+
a 6 3 i=1 3 i=1 2 3 2

5. Voici la procédure qui met en oeuvre la méthode de Simpson :

[>simpson :=proc(f,a,b,N)
local i,h ;
h :=(b-a)/N ;
evalf(h/6*(f(a)+f(b)+2*sum(f(a+i*h)+2*f(a+(2*i+1)*h/2),i=1..N-1)+4*f(a+h/2))) ;
end ;
[>simpson(exp,0,1,10) ;
1.718282781
6. On trace le ln de l'erreur en fonction de lnh. On obtient un droite de pente 4. Ceci veut dire que l'erreur
commise est proportionnelle à h4 . Donc quand h est petit (N grand), cette méthode à trois points
d'interpolation est plus précise que les méhodes à un point et deux points précédentes.