Définir les fonctions \(I=]-1,1[\to\R\) suivantes :
et calculer leurs valeurs exactes \(R_1,\dots R_5\).
f0(x):=cos(x);
f1(x):=1/(1+x^2);
f2(x):=exp(-x^2);
f3(x):=sqrt(1+x^2);
f4(x):=x<0?0:1;
f5(x):=cos(x)*exp(sin(x));
F := [f0,f1,f2,f3,f4,f5];
R := seq(int(F[j](x),x,-1,1),j=0..5);
Écrire une fonction rectangles(f,a,b,n) qui calcule la somme des
rectangles de \(f\) sur \(n\) points équirépartis sur \([a,b]\).
rectangles(f,a,b,n) : = {
1/n*sum(f(a+k*(b-a)/n),k,0,n-1)
}
Représenter graphiquement la convergence de la méthode des rectangles (diagramme du nombre de chiffres corrects en fonction du nombre de subdivisions) pour les fonctions ci-dessus. On fera des calculs en précision de 60 chiffres.
Digits:=80;
for(k:=0;k<=5;k++) {
plotlist(seq([n,-logb(R[k]-rectangles(F[k],-1,1,n),10)],n=1..10),color=k,epaisseur=2)
}
%plotlist(seq([n,-log(Pi/4-rectangles(f1,-1,1,n))],n=1..10))
%plotlist(seq([n,[seq(R[k]-rectangles(F[k],-1,1,n),k=0..5)]],n=1..25))
Écrire la fonction trapezes(f,a,b,n) qui implante la méthode des
trapèzes, comparer sa convergence.
trapezes(f,a,b,n) := {
1/n*sum(f(a+k*(b-a)/n),k,0,n-1)
}
plotlist(seq([n,-logb(R[k]-rectangles(F[k],-1,1,n),10)],n=1..10),color=k,epaisseur=2)
La méthode de Newton-Cotes d’ordre \(k\) consiste à remplacer sur chaque subdivision une fonction par son polynôme interpolateur de Lagrange \(L_f\) d’ordre \(k\) en des points régulièrement espacés (comprenant les extrémités du segment).
où \(x_j=a+\frac{(b-a)j}k\), \(w_i=\int_0^1 L_i(x)\d x\) et
Par exemple, la méthode des trapèzes est la méthode de Newton-Cotes d’ordre 1 (la fonction est remplacée par le segment affine qui joint les deux extrémités du segment).
Calculer les coefficients \(w_i\) pour les méthodes d’ordre \(3\) (Simpson) et
\(5\) (Boole), en vérifiant qu’ils sont indépendants de l’intervalle \([a,b]\). On
pourra utiliser la fonction lagrange.
int(lagrange(seq(j/2,j,0,2),seq(X^j,j,0,2)),x=0..1)
int(lagrange(seq(j/4,j,0,4),seq(X^j,j,0,4)),x=0..1)
Écrire une fonction coeffNewtonCotes(k) qui renvoie la liste des
coefficients de la méthode pour \(k\) points. Que se produit-il à partir du
rang 8 ?
coeffNewtonCotes := proc(k)
I := int(lagrange(seq(j/k,j,0,k),seq(X^j,j,0,k)),x=0..1)
return coeffs(I,X)
end proc;
Écrire une fonction NewtonCotes(f,k,a,b,n) qui calcule l’intégrale de
\(f\) sur \([a,b]\) avec \(n\) points par la méthode de Newton-Cotes d’ordre \(k\)
(effectuer \(\floor{n/k}\) subdivisions).
NewtonCotes:=proc(f,k,a,b,n)
wk := coeffNewtonCotes(k);
delta := (b-a)/(n*k);
xj := a;
res := 0;
for(l:=0;l<n;l++) {
for(j:=0;j<=k;j++) {
if(j!=0) xj += delta;
res += wk[j]*f(xj);
}
}
return res*(b-a)/n;
end proc;
Faire le graphique de la convergence de la méthode de Boole (Newton-Cotes d’ordre 5).
for(k:=0;k<=5;k++) {
plotlist(seq([n,-logb(R[k]-NewtonCotes(F[k],5,-1,1,n),10)],n=1..10),color=k,epaisseur=2);
plotlist(seq([n,-logb(R[k]-rectangles(F[k],-1,1,n),10)],n=2..100),color=k,epaisseur=2)
}
On va appliquer un procédé d’accélération de convergence à la formule des trapèzes. Pour ce faire, il faut avoir connaissance de l’ordre de convergence de la méthode des trapèzes de pas \(h=\frac{b-a}n\). On démontre (par exemple via la formule d’Euler-MacLaurin) que, en notant \(T(f,h)\) la formule des trapèzes de pas \(h\),
Le principe de l’accélération de convergence est de faire sauter ce terme \(\alpha_f h^2\) à l’aide d’un système [1] de deux relations prises en \(h\) et \(h'<h\). On obtient alors une convergence en \(h^4\), et ainsi de suite.
| [1] | La formule d’Euler-MacLaurin permet d’exprimer explicitement la constante \(\alpha_f\) en terme de valeurs de \(f'\), mais l’astuce et la généralité de l’accélération de convergence consiste à se débarasser de \(\alpha_f h^2\) sans le calculer. |
Pour la méthode de Romberg, on commence par calculer \(T(f,h)\) et \(T(f,h/2)\). Alors
On note
ce premier convergent, puis on itère le procédé.
La théorie de la formule des trapèzes montre qu’il existe \(\beta_f\) telle que \(R_{0,1}(f,h)=\int_a^b f+\beta_f h^4+O(h^6)\). On évalue \(R_{2,1}(f,h)=R_{1,1}(f,h/2)\) et on calcule
qui vérifie
On définit ainsi pour tout \(m\geq n\) les convergents \(R_{m,n}\) (avec les valeurs \(R_{m,0}(f,h)=T(f,h/2^m)\)). La méthode d’intégration de Romberg consiste à calculer les convergents \(R_{m,m}\). Pour simplifier, on prend \(h=(b-a)\) dans le calcul de \(R_{0,0}\), de sorte que \(h=\frac{b-a}{2^m}\) dans \(R_{m,0}\).
La méthode de Gauss-Lengendre (de même que ses variantes Gauss-, Gauss-Hermitte…) est aussi une méthode d’interpolation. Cependant, par un choix plus judicieux des points d’interpolation, elle est exacte à un ordre beaucoup plus élevé (avec \(n+1\) points, la méthode est exacte à l’ordre \(2n-1\)).