Feuille Principale.ipynb

Kernel: SageMath (stable)

$$Accélération$$ $$de$$ $$convergence$$

Vitesse de Convergence

1.Présentation du problème

Présentation et Contexte $$\\$$ On dispose d'une suite réelle $(u_n)$ qui converge avec un certain ordre vers $l$ et on souhaiterait construire, à partir de $(u_n)$, par un procédé général, une suite $(v_n)$ qui converge encore vers $l$ mais avec un ordre supérieur. On s'intéresse essentiellemnt à deux procédé : $$\\$$ Le $\Delta^2$ d'Aitken et l'extrapolation de Richardson.

2. Définitions

Accélération de la convergence :

Soit $(v_n)$ une suite convergente vers $l$ avec $v_n\neq l$ $\forall n \in \mathbb{N}$, on dit que la convergence de cette suite vers $l$ est plus rapide que celle de $(u_n)$ si : $$\lim\limits_{n \to \infty}\frac{v_n-l}{u_n-l}=0$$ Accélérer la convergence d'une suite consiste à construire à partir de cette dernière une autre suite qui converge plus vite vers la même limite. $$\\$$

Notion d'ordre et Vitesse de convergence:

Soient $l$ la limite de la suite (u$_n$) et q un entier appelé ordre de convergence, on a : $$\lim\limits_{n \to \infty}\frac{|u_{n+1}-l|}{|u_n-l|^q}=k$$ avec $0<k<1$ Ici $k$ désigne le taux de convergence, plus précisemment la vitesse à laquelle (u$_n$) converge.Si $q=1$ alors la convergence de la suite est dite linéaire, en revanche si $q=2$ celle-ci est dite quadratique. $$\\$$ La vitesse de convergence est un facteur important de la qualité des algorithmes. Si la vitesse de convergence est élevée, l'algorithme converge rapidement et le temps de calcul est moindre. Ces préoccupations de rapidité de convergence ont conduit à diversifier les modes de convergenceet à chercher des processus optimaux.

Exemple :

Étudions la convergence des deux suites $u_n=1+\frac{1}{n}$ et $v_n=1+\frac{1}{n^5}$. Il est clair que ces deux suites convergent vers 1. Observons leur vitesse de convergence:

u(n)=1+1/n
v(n)=1+1/(n^5)
for i in range (1,10):
    print i,RR(u(i)),RR(v(i))

1 2.00000000000000 2.00000000000000
2 1.50000000000000 1.03125000000000
3 1.33333333333333 1.00411522633745
4 1.25000000000000 1.00097656250000
5 1.20000000000000 1.00032000000000
6 1.16666666666667 1.00012860082305
7 1.14285714285714 1.00005949901827
8 1.12500000000000 1.00003051757812
9 1.11111111111111 1.00001693508781

On remarque de $v_n$ convergerait plus vite que $u_n$. Vérifions que $$\lim\limits_{n \to \infty}\frac{v_n-1}{u_n-1}=0$$

print(limit((v(n)-1)/(u(n)-1),n=oo))

0

3.Aitken

Le $\Delta^2$ Aitken: C'est un procédé d’accélération de la convergence de suites en analyse numérique popularisé par le mathématicien Alexander Aitken en 1926. Lorsqu'on applique le procédé d'accélération de la convergence d'Aitken, cela consiste de créer à partir $(u_n)$ une nouvelle suite $(v_n)$ ayant la même limite que celle-ci; convergent plus rapidement vers $l$. On définit alors: $$v_n=u_n-\frac{(u_{n+1}-u_n)^2}{u_n-2u_{n+1}+u_{n+2}}=u_n-\frac{(\Delta u_n)^2}{\Delta^2u_n}$$ La suite définie par $u_{n+1}=g(u_n)$, à l'aide de ce procéder on vient de montrer qu'à partir d'un certain rang, $v_n$ est une meilleure approximation de $l$ (la limite) que $u_n$. Il est naturel de poursuivre les itérations en remplaçant $u_n$ par $v_n$.

Exemples d'application: Méthodes itératives (point fixe: $u_{n+1}=f(u_n)$)$$\\$$ Sommes de séries : $\frac{\pi}{4}=\sum_{k=0}^{\infty}\frac{(-1)^k}{2k+1}$ $$\\$$

Hypothèses: Le procédé d'Aitken s'applique aux suites ayant une convergence linéaire comme par exemple la méthode du point fixe ou alors les suites géométriques.$$\\$$ De plus $(u_n) \neq (v_n)$ $$\\$$ Il accélère toutes les suites dont le rapport de deux écarts consécutifs converge vers une limite non nulle comrpise entre -1 et 1. Aussi, (u$_n$)est mal défini si $\Delta^2u_n$ contient un "0", c'est-à-dire, les premières différences se répètent donc il faut le restreindre à un indice $n>0$. $$\\$$

Instabilité de la méthode: Le $\Delta^2 Aitken$ est très instable numériquement car le numérateur et le dénominateur sont proches de 0, il faut alors calculer beaucoup de décimales. Il vaut mieux alors utiliser une autre écriture équivalente à l'initiale, mais beaucoup plus stable telle que :

$$v_n=u_{n+1}+\frac{1}{\frac{1}{u_{n+2}-u_{n+1}}-\frac{1}{u_{n+1}-u_n}}$$

Démonstration de l'accélération: Soit $(u_n)$ une suite qui converge vers l et $\lim\limits_{n \to \infty}\frac{u_{n+1}-l}{u_n-l}= k$ avec $0 < k< 1$ $$\\$$ Posons S$_n$= $\frac{u_{n+1}-u_n}{u_n-u_{n-1}}$ , alors la suite $v_n=\frac{u_{n+1}-s_nu_n}{1-s_n}$ converge vers $l$ plus vite que $(u_n)$ $$\\$$ Soit $l=0$, $$\begin{align*}v_n &=\frac{u_n^2-u_{n-1}u_{n+1}}{2u_n-u_{n-1}-u_{n+1}}\\ &=u_{n-1}-\frac{(u_n-u_{n-1})^2}{u_{n+1}-2u_n+u_{n-1}}\\&= u_{n-1}-\frac{(\Delta u_n)^2}{(\Delta)^2u_n}\end{align*}$$

D'où $$\begin{align*}v_n &=\frac{v_n}{u_n}\\ &=\frac{1}{u_n}v_n\\ &=\frac{u_nu_{n+1}}{u_nu_{n+1}}\frac{(\frac{u_n}{u_{n+1}}-\frac{u_{n-1}}{u_n})}{(2\frac{u_n}{u_{n+1}}-\frac{u_{n-1}}{u_{n+1}}-1)}\\ &=\frac{\frac{u_n}{u_{n+1}}-\frac{u_{n-1}}{u_n}}{2\frac{u_n}{u_{n+1}}-\frac{u_{n-1}}{u_{n+1}}-1}\end{align*}$$ Le dénominateur tends vers $2k-k^2-1=-(k-1)^2\neq0$, et le numérateur tend vers $k-k=0$. Par conséquent, (v$_n$) tend "plus vite" vers 0 que (u$_n$). $$\\$$

4. Codage Aitken et Steffensen

Aitken:

def delta(u,n):
    return u(n+1)-u(n)
def delta2(u,n):
    return u(n)-2*u(n+1)+u(n+2)

def Aitken(u):
    var('x')
    return u(x-1)-(delta(u,x-1)**2)/delta2(u,x-1)

def AitkenV(u,n):
    return [RR(u(k)-(delta(u,k)**2)/delta2(u,k)) for k in range(n+1)]
R=RealField(100)

Avec le point fixe

Soit l'équation $e^x=2\ln(1+x)$, on sait qu'il existe une unique solution dans l'intervalle $[1;5]$ notée $r$. Soit $w(x)=ln(1+x)+ln(2)$ vérifiant $w(r)=r$ et $|w'(r)|<1$.

Étudions la fonction $w(n)=\ln(1+n)+\ln(2)$. Nous allons donc essayer plusieurs "méthodes" du point fixe : Version Classique; Aitken; Steffensen.

Méthode du point fixe:

w(n)=ln(1+n)+ln(2)
x0=1
p=w(x0)-x0
L=[]
LL=[]
for i in range(20):
    z=x0
    x0=w(x0)
    t=log(p/(x0-z))
    L.append((i,t))
    LL.append(R(x0))
    print R(x0)
pp=point((i for i in L),ymin=-5,ymax=10,color='red')
pp.axes_labels(['Iterations','Decimales Won'])
pp

print'--------------'
print'Ordre de convergence'
print
for i in range (9):
    print (math.log10((abs(LL[i+3]-LL[i+2]))/(abs(LL[i+2]-LL[i+1])))/(math.log10((abs(LL[i+2]-LL[i+1]))/abs((LL[i+1]-LL[i])))))

1.3862943611198906188344642429
1.5628888667518891735995210268
1.6342822663899545545473495528
1.6617579407780844965756466083
1.6721339650603005203749444976
1.6760245716019100257735064889
1.6774795049967419284155390661
1.6780230493407211374206254570
1.6782260347152329265494128200
1.6783018285670211840314179215
1.6783301281871213108396160315
1.6783406943863035067686494819
1.6783446394482765826026323419
1.6783461123972094453350854755
1.6783466623445560810577112230
1.6783468676754831145709669317
1.6783469443387926292167378820
1.6783469729621607757033850637
1.6783469836491138062058034466
1.6783469876392439562818574938
--------------
Ordre de convergence

1.05438478868
1.0197913856
1.00732125318
1.00272416998
1.00101581736
1.00037909071
1.00014151405
1.00005283285
1.00001972543

Appliquons Aitken au point fixe $$\\$$ Aitken:

LA=[]
w(n)=ln(1+n)+ln(2)
x0=1
wA=Aitken(w)
p=wA(x0)-x0
n=10
for i in range(n):
    z=x0
    x0=R(wA(x0))
    t=log(p/(x0-z))
    LA.append((i,t))
    print x0

q=point((i for i in LA),ymin=-5,ymax=10,color='blue')
q.axes_labels(['Iterations','Decimales Won'])
q.save('test.png')
q
#for i in LA:
    #print i[0],float(real_part(i[1]))

2.3632304423479432064700760983
2.8984305361316223563194660132
3.0481005868806815015414405654
3.0861918387170173352874499185
3.0956562795278705025157770662
3.0979939049107127015279846477
3.0985704259810465475319950255
3.0987125598861069586226081838
3.0987475980465441807998182821
3.0987562352931055209097977237

Maintenant, comparons la méthode du point fixe et Aikten :

q+pp

Calculons la pente des courbes:

Aitken: On prend deux points par lecture graphique A(4;5) et B(6;8). En calculant le coefficient directeur on a CoeffA=$\frac{8-5}{6-4}=\frac{3}{2}$

Point fixe: On a deux points C(0;0) et D(1;1). CoeffPF=$\frac{1-0}{1-0}=1$

Or on remarque que CoeffA>CoeffPF, alors la fonction $w(n)$ converge plus rapidement avec le procédé d'Aitken qu'en l'appliant avec la méthode du point fixe.

On remarque que la converge est beaucoup plus rapide avec le procédé d'Aitken. $$\\$$ Essayons avec la méthode de Steffensen. $$\\$$

5. Steffensen

http://www.lamfa.u-picardie.fr/chehab/Ens/cours_zeros.pdf

Méthode de Steffensen: Dans l'analyse numérique, la méthode de Steffensen est un procédé de recherche de racines, similaire à la méthode de Newton, nommée d'après Johan Frederik Steffensen. La méthode de Steffensen atteint aussi la convergence quadratique, mais sans utiliser de dérivées comme la méthode de Newton. Ce procédé aussi accélère la vitesse de convergence d'une suite, mais plus rapidement que Aitken. $$\\$$ En partant de la méthode du point fixe et de newton. Il existe une fontion $g$ définie comme pour la méthode de Newton, mais n'étant pas la dérivée de $f$. $$\\$$ A savoir : Méthode de Newton = $$u_{n+1} = u_n-\frac{f(u_n)}{f'(u_n)}$$

définition du procédé: Soit une méthode itérative engendrée par g. $$\\$$ Lorsqu'on applique le procédé d'accélération de la convergence d'Aitken à la suite définie par $u_{n+1}=g(u_n)$, on vient à partir d'un certain rang, $v_n$ est une meilleure approximation de $l$ que $u_n$. Il est donc naturel de poursuivre les itérations en remplaçant $u_n$ par $v_n$, puis de calculer $g(v_n)$ et $g[g(v_n)]$ pour appliquer l'accélération d'Aitken. $$\\$$ Ce qui revient à calculer successivement: $$\begin{align*}y_n &=v_n\end{align*}$$ Et$$\begin{align*}y_{n+1} &=g(y_n)\end{align*}$$ Et$$\begin{align*}y_{n+2} &=g(y_{n+1})\end{align*}$$

donc: $$v_{n+1}=y_n-\frac{(y_{n+1}-y_n)^2}{y_{n+2}-2y_{n+1}+y_n}$$ Or, on peut définir $v_{n+1}$ en fonction de $v_n$ à partir de g: $$v_{n+1}=v_n-\frac{(g(v_n)-v_n)^2}{g(g(v_n))-2g(v_n)+v_n}$$ C'est la méthode de Steffensen que l'on peut définir formellement par la fonction G suivante qui donne le processus itératif (avec $u_n=x$ (par le point fixe)): $$G(x)=\frac{xg(g(x))-(g(x))^2}{g(g(x))-2g(x)+x}$$ En effet, supposons $l$ une racine de l'équaton $x=g(x)$ Si la méthode engendré par $g$ est d'ordre 1, alors la méthode de Steffensen est d'ordre 2 au moins.

Ordre de la méthode de Steffensen: $$\\$$ La méthode de Steffensen est d'ordre 2, ainsi G'(x)=0. $$\\$$ Calculons G'(x): $$G(x)=\frac{xg(g(x))-(g(x))^2}{g(g(x))-2g(x)+x}$$ $G(x)=\frac{U(x)}{V(x)}$ avec $U(x)= xg(g(x))-(g(x))^2$ et $V(x)= g(g(x))-2g(x)+x$ $$\\$$ Soit $U'(x)=g(g(x))+xg'(g(x))g'(x)-2g(x)$ et $V'(x)=g'(g(x))-2g'(x)+1$ $$\\$$ Alors $G'(x)=\frac{U'(x)V(x)-V'(x)U(x)}{V(x)^2}$ $$\\$$ Donc $$\begin{align*}G'(x) &=\frac{(g(g(x)+xg'(g(x))g'(x)-2g(x))(g(g(x))-2g(x)+x)-(g'(x)-2g'(x)+1)(xg(x)-(g(x))^2)}{(V(x))^2}\end{align*}$$ Et en prenant le dénominateur $$\begin{align*}&=(-g(x)+xg'(x)^2)(-g(x)+x)-(-g'(x)+1)(xg(x)-(g(x))^2)\\ &=g(x)^2-xg(x)-xg'(x)^2g(x)+x^2g'(x)^2-(-xg(x)g'(x)+g(x)^2g'(x)+xg(x)-g(x)^2)\\ &=g(x)^2-xg(x)-xg'(x)^2g(x)+x^2g'(x)^2+xg(x)g'(x)-g(x)^2g'(x)-xg(x)+g(x)^2\\ &=2g(x)^2-2xg(x)-xg'(x)^2(g(x)-x)+g(x)g'(x)(x-g(x))\end{align*}$$ $$\\$$ En ayant supposer une racine $l$ de l'équaton $x=g(x)$, il nous manque qu'à remplacer dans l'expression $$\\$$$$\begin{align*}&=2x^2-2x^2-xg'(x)^2(x-x)+xg'(x)(x-x)\\ G'(x) &=0\end{align*}$$

def Steffensen(u):
    var('x')
    return u(x)+1/(1/(u(u(x))-u(x))-1/(u(x)-x))
LS=[]
w(n)=ln(1+n)+ln(2)
wS=Steffensen(w)
RR=RealField(500)
x0=1
p=x0
for i in range (6):
    z=x0
    x0=RR(wS(z))
    t=log(p/(x0-z))
    t=real_part(t)
    LS.append((i,t))

    
qq=point((i for i in LS),ymin=-5,ymax=103,color='blue')
qq.axes_labels(['Iterations','Decimales Won'])
qq

Appliquons maintenant $Aitken$ à la fameuse série de Leibniz : $$s_n=\sum_{k=0}^{\infty}\frac{(-1)^k}{2k+1}=\frac{\pi}{4}$$ On a $s_n=\frac{\pi}{4}$. Donc $4 s_n=\pi$ (notons $S_n=4 s_n$). $$\\$$ Etudions la convergence de $S_n$ avec le procédé d'Aitken (en utilisant la série $sA_n$) et sans le procédé(avec $s_n$).

from time import time
def s(x):
    var('k')
    return 4*sum((-1)^k/(2*k+1),k,0,x)
def sA(x):
    return AitkenV(s,x)
eps=10e-4 #Si on prend un précision plus grande le temps de calcul est long donc on se contente de 10e-4.
i=0
j=0
#Avec le procéde d'Aitken
debut1=time()
while(abs(sA(i)[i]-pi)>eps):
    i=i+1
fin1=time()-debut1
#Sans le procédé d'Aitken
debut2=time()
while(abs(s(j)-pi)>eps):
    j=j+1
fin2=time()-debut2

    
print i,'itérations en',fin1,'sec'
print j,'itérations',fin2,'sec'
print 'Aitken est',floor(fin2/fin1),'fois plus rapide ! Trop super!!'
# En testant avec eps=10e-5 on obtient i=12 et j=9999

5 itérations en 0.0976889133453 sec
999 itérations 12.1485378742 sec
Aitken est 124 fois plus rapide ! Trop super!!

6. Richardson

En analyse numérique, le procédé d'extrapolation de Richardson est une technique d'accélération de la convergence.

Définition: L'extrapolation est le calcul d'un point d'une courbe dont on ne dispose pas d'équation, à partir d'autres points, lorsque l'abscisse du point à calculer est au-dessus du maximum ou en dessous du minimum des points connus.

Théorème: Soit $(u_n)$ une suite de nombres réels qui converge vers $l$. On suppose qu'on a un développement asymptotique de la forme $$\\$$ $u_n=l+\lambda k^n+O(k'^n)$ avec $\lambda \in \mathbb{R}$ et $|k'|<|k|<1$, de sorte que la convergence de $(u_n)$ vers $l$ est géométrique de rapport $k$. On pose $$\\$$ $v_n=\frac{u_{n+1}-ku_n}{1-k}$ $$\\$$ Alors, $v_n-l$ est un $O(k'^n)$ et donc $v_n$ converge ver $l$ avec une convergence ui est (au moins) géométrique de rapport $k'$.

Hypothèses: Cette méthode s'applique lorsqu'on a une suite $(u_n)$ qui converge vers $l$ avec une convergence géométrique de rapport $k<1$ et qu'on connaît le rapport $k$. $$\\$$ Il faut connaître le rapport $k$, indispensable pour calculer $v_n$. $$\\$$ Il faut ne pas connaître le coefficient $\lambda$ (sinon il suffit de retrancher $\lambda k^n$ pour avoir aussitôt une suite qui converge comme un $O(k'^n)$).

Convergence géométrique: On dit qu'une suite $(u_n)$ converge (au moins) géométriquement vers $l$ si elle est dominée par une suite géométrique $(k^n)$ avec $0<k<1$.

Démonstration: On va commencer par expliquer d'où sort la suite $v_n$. Il s'agit d'éliminer le terme en $k^n$ dans l'expression de $u_n$.L'idée, comme on ne connaît pas $\lambda$, est de regarder $u_{n+1}$, puis, en calculant $u_{n+1}-ku_n$, d'éliminer les $k^n$ entre les deux termes. Maintenant , on note $u_{n+1}-ku_n$ converge vers $(1-k)l$ et, en divisant par $1-k$, on trouve une suite $v_n$ qui converge vers $l$. $$\\$$ Précisément, posons $u_n=l+\lambda k^n+wn$, de sorte que $w_n$ est un $O(k^n)$. Un calcul immédiat donne $$\\$$ $v_n-l=\frac{w_{n+1}-kw_n}{1-k}$ $$\\$$ Si on a $|w_n|\leq A|k'|^n$, on a alors $$\\$$ $$\begin{align*}|v_n-l| &\leq \frac{|w_{n+1}|+|k||w_n|}{1-k}\\ &\leq \frac{A(|k'|+|k|)}{1-k} |k'^n|\\ &= M|k'^n|\end{align*}$$ CQFD. $$\\$$

Si on prend $$\begin{align*}\frac{v_n-l}{u_n-l} &=\frac{M|k'^n|}{\lambda |k|^n}\\ &Et\lim\limits_{n \to \infty}\frac{v_n-l}{u_n-l}\\ &=\lim\limits_{n \to \infty}\frac{M}{\lambda}|\frac{k'}{k}|^n\\\ &=0\end{align*}$$ car $|k'|<|k|<1$ et $|\frac{k'}{k}|<1$ Alors $\lim\limits_{n \to \infty}|\frac{k'}{k}|^n=0$. Ce qui confirme que $(v_n)$ converge vers $l$ plus vite que $(u_n)$.

Remarque: Cette méthode s'applique notamment à la méthode des trapèzes pour calculer la valeur approchée d'une intégrale. Elle prend alors le nom de Romberg-Richardson, et est effectivement employée dans les logiciels comme Maple ou les calculatrices.

Applications:$$\\$$ -Approximation de $\pi$ par la convergence de la suite d'Archimède $A_n=2^nsin(\frac{\pi}{2^n})$. $$\\$$ -Approximation de $e$ par la convergence de la suite $u_n=\exp(n\ln(1+\frac{1}{n}))$.

Prenons l'exemple de l'approximation de $\pi$:

On utilise la methode de Richardson (procédé d’accélérations de convergence successives) pour accelerer la convergence de la suite d'archimède (sans connaitre la valeur de $\pi$)

$$A_{n+1}=2^{n+1}\sin(\frac{\pi}{2^n})$$

qui tend vers $\pi$ quand n tend vers l'infini .

Rappelons le developpement limité de $\sin(x)$ à l'ordre $2k+1$ est : $$\sin(x)=x-\frac{x^3}{3!}+\frac{x^5}{5!}+...+(-1)^k\frac{x^{2k+1}}{(2k+1)!}+o(x^{2k+3})$$

On peut alors écrire pour tout entier naturel $n$:

$$\begin{align*}A_n&=2^{n+1}\left(\frac{\pi}{2^{n+1}}-\frac{1}{3!}\Bigg(\frac{\pi}{2^{n+1}}\Bigg)^3+\frac{1}{5!}\Bigg(\frac{\pi}{2^{n+1}}\Bigg)^5-...+(-1)^k\frac{1}{(2k+1)!}\Bigg(\frac{\pi}{2^{n+1}}\Bigg)^{2k+1}\right)+o(...)\\ &=\pi-\frac{1}{3!}\frac{\pi^3}{2^{2n+2}}+\frac{1}{5!}\frac{\pi^5}{2^{4n+4}}-...+(-1)^k\frac{1}{(2k+1)!}\frac{\pi^{2k+1}}{2^{2kn+2k}}+o(...)\\ &=\pi-\frac{\pi^3}{3!2^2}\frac{1}{4^n}+\frac{\pi^5}{5!2^4}\frac{1}{4^{2n}}-...+\frac{(-1)^k\pi^{2k+1}}{(2k+1)!2^{2k}}\frac{1}{4^{kn}}+o(...)\\ &=\pi+\frac{e_1}{4^n}+\frac{e_2}{4^{2n}}+...+\frac{e_k}{4^{kn}}+o\left(\frac{1}{4^{kn}}\right)~forme~recherchée\end{align*}$$

avec $e_k=\frac{(-1)^k\pi^{2k+1}}{(2k+1)!2^{2k}}$ les termes d'erreurs.

La suite ($A_n$) converge vers $\pi$ avec une vitesse de convergence égal à $\frac{1}{4}$. Nous pouvons le vérifier avec le code suivant :

#convergence de An
A(n)=2^n*sin(pi/2^n)
R=RealField(100)
for i in range(10):
    print RR(abs(R(A(i+1)-pi)/R(A(i)-pi)))
    

0.363380227632419
0.274323356811065
0.255855728965272
0.251450263739790
0.250361717852688
0.250090376554618
0.250022590833378
0.250005647501788
0.250001411862538
0.250000352964828

L’idée de l'accélération de convergence de Richardson consiste à faire disparaître les termes d’erreur en $\frac{1}{4^{kn}}a_k$ les un après les autres et fabriquer un nouvelle suite $\big(A_{n}^{1}\big)$ avec une convergence vers $\pi$ plus rapide. Cette élimination peut se fait algébriquement : $$\left\{ \begin{matrix}\begin{align*} A_n&=\pi-\frac{1}{4^n}a_1+o(\frac{1}{4^n})\\ A_{n+1}&=\pi-\frac{1}{4^{n+1}}a_1+o(\frac{1}{4^{n+1}})\end{align*}\end{matrix}\right .$$

On multiplie par 4 la seconde ligne:

$$\left\{ \begin{matrix}\begin{align*} A_n&=\pi-\frac{1}{4^n}a_1+o(\frac{1}{4^n})\\ 4A_{n+1}&=4\pi-\frac{1}{4^n}a_1+o(\frac{1}{4^n})\end{align*}\end{matrix}\right .$$

On a ensuite

$$\begin{align*} 4A_{n+1}-A_n&=3\pi+o\left(\frac{1}{4^n}\right)\\ \Rightarrow~\frac{4A_{n+1}-A_n}{3}&=\pi+o(\frac{1}{4^n})\end{align*}$$

On considère la suite $\big(A_{n}^{1}\big)$ définie par $$A_{n}^{1}=\frac{4A_{n+1}-A_n}{3}$$

le terme d'erreur en $\frac{1}{4^n}$ disparait car il est négligeable devant $\frac{1}{4^n}$ et son terme d'erreur est en $o\Big(\frac{1}{4^{2n}}\Big)$, sa vitesse de convergence est $\frac{1}{16}$ (voir code ci-dessous). $\big(A_{n}^{1}\big)$ converge vers $\pi$ plus vite que $A_n$ nous le montrerons après dans un tableau car ce n'est pas fini.

A1(n)=(4*A(n+1)-A(n))/3
for i in range (10):
    print RR(abs(R(A1(i+1)-pi)/R(A1(i)-pi)))

0.0779556287666338
0.0660420039391321
0.0633667112218643
0.0627155233646191
0.0625538089965189
0.0625134477636382
0.0625033616606413
0.0625008403976447
0.0625002100983165
0.0625000525245107

On peut donc programmer Richardson :

def Richardson(s,n):
    L=[[RR(s(k)) for k in range (1,n)]]
    LL=[]
    for i in range (1,n-1):
        for j in range (n-i-1):
            LL.append((4^(i)*L[i-1][j+1]-L[i-1][j])/(4^i-1))
        L.append(LL)
        LL=[]
    return L
u(n)=2^n*sin(pi/2^n)
t(n)=exp(n*ln(1+1/n))
v(n)=t(2^n)#On applique Romberg pour aller plus vite
print("ATTENTION LE RESULTAT SE LIT DE DROITE A GAUCHE ET DE HAUT EN BAS")
print("Approximation de pi")
L=Richardson(u,9)
for i in range (len(L)):
    print list(reversed(L[i]))
L=Richardson(v,9)
print('-------------')
print("Approximation de e")
for i in range (len(L)):
    print list(reversed(L[i]))

ATTENTION LE RESULTAT SE LIT DE DROITE A GAUCHE ET DE HAUT EN BAS
Approximation de pi
[3.14151380114430, 3.14127725093277, 3.14033115695475, 3.13654849054594, 3.12144515225805, 3.06146745892072, 2.82842712474619, 2.00000000000000]
[3.14159265121481, 3.14159261559211, 3.14159204575769, 3.14158293664190, 3.14143771670383, 3.13914757031223, 3.10456949966159]
[3.14159265358966, 3.14159265358107, 3.14159265303208, 3.14159261797111, 3.14159039312994, 3.14145277502227]
[3.14159265358979, 3.14159265358979, 3.14159265358860, 3.14159265328605, 3.14159257754434]
[3.14159265358979, 3.14159265358979, 3.14159265358979, 3.14159265358307]
[3.14159265358979, 3.14159265358979, 3.14159265358979]
[3.14159265358979, 3.14159265358979]
[3.14159265358979]
-------------
Approximation de e
[2.71299162425343, 2.70773901968802, 2.69734495256510, 2.67699012937818, 2.63792849736660, 2.56578451395035, 2.44140625000000, 2.25000000000000]
[2.71474249244191, 2.71120370872899, 2.70412989362740, 2.69001067338204, 2.66197649183868, 2.60724393526713, 2.50520833333333]
[2.71497841135610, 2.71167529640243, 2.70507117497710, 2.69187961881827, 2.66562532894345, 2.61404630872938]
[2.71503084175219, 2.71178012372665, 2.70528056475739, 2.69229635357819, 2.66644404355003]
[2.71504358966602, 2.71180561219319, 2.70533148323261, 2.69239773518614]
[2.71504675484439, 2.71181194076500, 2.70534412619257]
[2.71504754478677, 2.71181352020690]
[2.71504774218802]

7.Méthode de Romberg:

A)La méthode des trapèzes composites

La méthode de Romberg utilise le procédé d’extrapolation de Richardson appliqué à la méthode des trapèzes (calcul d'intrégrale) et permet d'améliorer l'ordre de convergence.

Soit $f$ une fonction de classe $C^{\infty}$ sur le segment [$a,b$] avec $a<b$, on cherche a approximer l'intégrale $$\int_{a}^{b}f(t)dt$$ que l'on note $F(t)$. Par la méthode des trapèzes-composites, avec un découpage de [$a,b$] en $n$ segments réguliers, on obtient la formule d'approximation de l’intégrale $I$ de $f$ suivante:

$$\begin{align*}I(h)=h\Big(\frac{f(a)+f(b)}{2}+\sum_{i=1}^{n-1}f(a+ih)\Big)~~~~ avec ~~~~ h=\frac{b-a}{n}.\end{align*}$$

La méthode est dite d’ordre 1 elle donne un résultat exact pour un polynôme de degré 1.

def compositetrapeze(f,a,b,n):
    h=(b-a)/n
    var('i')
    I=h*(((f(a)+f(b))/2)+sum(f(a+i*h),i,1,n-1))
    return RR(I)

#test
f(x)=1/x
print compositetrapeze(f,1,2,1000)

F=integral(f,x,1,2)
float(F)

0.693147243059937

0.6931471805599453

B) Formule d'Euler-Maclaurin :

La formule d'Euler-Maclaurin est une relation entre sommes discrètes et intégrales, découverte par Euler en 1735 pour accéléré le calcul de limites de séries lentement convergentes et par Colin Maclaurin pour calculer des valeurs approchées d'intégrales.

Elle nous prouve l'existence d'un développement limité à tout ordre de l’erreur commise dans la méthode des trapèzes, de la forme

$$\begin{align*}E(h)=a_1h^2+a_2h^4+...+a_{k-1}h^{2k-2}+o(h^{2k})~ou~&E(h)=F(t)-I(h)\end{align*}$$

D'après la formule d'Euler-Maclaurin, le développement limité de l'erreur commise ne contient que des puissances paires de h. L'idée est d'annuler un à un les terme du développement E en effectuant une combinaison linéaire entre $I(h)$ et $I(h/2)$:

$$\begin{align*}I(h)&=F(t)+a_1h^2+a_2h^4+...+a_{k-1}h^{2k-2}+o(h^{2k})\\ I(\frac{h}{2})&=F(t)+a_1\frac{h^2}{2^2}+a_2\frac{h^4}{2^4}+...+a_{k-1}\frac{h^{2k-2}}{2^{2k-2}}+o(h^{2k})\\ \hspace{1.5cm}--------&------------------------\\ (*) ~~~ \frac{4I(\frac{h}{2})-I(h)}{3}&=F(t) \hspace{1.5cm} +\underbrace{b_2h^4+...+b_{k-1}h^{2k-2}+o(h^{2k})}\\ &\hspace{5.6cm}o(h^4)\end{align*}$$

On remarque ainsi, que l'ordre de l'erreur est en $h^4$ donc en $\frac{1}{n^4}$ au lieu de $\frac{1}{n^2}$ (car $b-a$ et est constante et $n$ on le choisit mais tend vers l'infini pour plus de précision).

• Démonstration : l'expression $(*)$ obtenue correspond à la méthode de Simpson (erreur en $\frac{1}{n^4}$ une interpolation parabolique):

En developpant la dernière expression on obtient

$$\begin{align*}F(t)&=\frac{1}{3}\left(2h\Big(\frac{f(a)+f(b)}{2}+\sum_{i=1}^{n-1}f(a+ih)+\sum_{i=1}^{n-1}f(a+(2i-1))\frac{h}{2})\Big)-h\Big((\frac{f(a)+f(b)}{2})+\sum_{i=1}^{n-1}f(a+ih)\Big)\right)\\ &=\frac{1}{3}\left(h\Big(\frac{f(a)+f(b)}{2}+\sum_{i=1}^{n-1}f(a+ih)\Big)+2h\sum_{i=1}^{n-1}f(a+(2i-1)\frac{h}{2})\right)\\ &=\frac{h}{3}\left(\frac{f(a)+f(b)}{2}+\sum_{i=1}^{n-1}f(a+ih)+2\sum_{i=1}^{n-1}f(a+(2i-1)\frac{h}{2})\right)\\ &=\frac{1}{2}\frac{h}{3}\left(f(a)+f(b)+2\sum_{i=1}^{n-1}f(a+ih)+4\sum_{i=1}^{n-1}f(a+(2i-1)\frac{h}{2})\right)\end{align*}$$

C)Extrapolation de Richardson - Formule récurrente de Romberg $A_{n,k}$:

On utilise l'extrapolation de Richardson à l'expresion $(*)$ ce qui va nous permettre d'avoir une approximation de plus en plus précise et voir ensuite la formule recurrente de Romberg.

$$\begin{align*} F(t)&=I(h)+a_{2k}h^{2k}+o(h^{2k}) \hspace{1cm} (**)\\ F(t)&=I(\frac{h}{2})+a_{2k}(\frac{h}{2})^{2k}+o(h^{2k})\\ F(t)&=I(\frac{h}{2})+a_{2k}\frac{h^{2k}}{4^k}+o(h^{2k})\\ 4^kF(t)&=4^kI(\frac{h}{2})+a_{2k}h^{2k}+o(h^{2k}) \hspace{1cm} (***)\\ (***)-(**) \hspace{2cm} 4^kF(t)-F(t)&=4^kI(\frac{h}{2})-I(h)+o(h^{2k})\\ F(t)&=\frac{4^kI(\frac{h}{2})-I(h)}{4^k-1}+o(h^{2k})\end{align*}$$

et on note $A(h)=\frac{4^kI(\frac{h}{2})-I(h)}{4^k-1}$ et en choisissant A(h) comme approximation de F(t), l'erreur est en $h^{2k+2}$.

Maintenant on prend $h_n=\frac{(b-a)}{2^n}$ un nombre pair d'intervalles de la subdivision.

• Notons désormais A(n,0) les précédents nombres I(h_n) résultats de la méthode des trapèzes. $A_{n,0}=I(h_n)$

pour n=0, $h_0=\frac{(b-a)}{2^0}=b-a$ donc $A_{0,0}=I(h_0)$ $$\begin{align*}A_{0,0} &=\frac{1}{2}(b-a)(f(a)+f(b))\\ A_{1,0} &=\frac{b-a}{4}(f(a)+f(b)+2f(a+\frac{b-a}{2}))\\ A_{n,0} &=I(h_n)=\frac{1}{2}A_{n-1,0}+h_n\sum_{k=0}^{2^{n+1}-1}f(a+(2k+1)\frac{b-a}{2^n})\end{align*}$$

Ces valeurs seront placées dans la première colonne, elle correspond (par définition) à la méthode du trapèze.

• Notons A(n,1) obtenues par les extrapolations de Richardson de la forme $A(h_n)=\frac{4^kI(\frac{h_n}{2})-I(h_n)}{4^k-1}=\frac{4^kI(h_{n+1})-I(h_n)}{4^k-1}$ a partir $I(h_n)$ et $I(h_n/2)$ avec k$\ge$1 et n$\ge$1.

pour n=1, on a $A(1,1)=\frac{4I(h_1)-I(h_0)}{3}=\frac{4A_{1,0}-A_{0,0}}{3}$ car $\frac{h_0}{2}=h_1$

et on a $A_{n,1}=\frac{4A_{n,0}-A_{n-1,0}}{3}$

Ces valeurs seront placées dans la deuxième colonne, c'est la méthode de Simpson. En troisième colonne on aura la méthode de Boole qui est la formule de Newton-Cotes.

Nous pouvons désormais accélérer la convergence au seul moyen des $A_{n,k}$ issus de l'extrapolation de Richardson est obtenir la formule récurrente de Romberg $A_{n,k}=\frac{4^kA_{n,k-1}-A_{n-1,k-1}}{4^k-1}$ qui est démontrer en dessous, c'est une approximation de l'intégrale qui à l'étape $A_{n,n}$ est en $o(h_n^{2n+2})$ ceci à condition que la fonction soit de classe $C^{2n+2}$.Le calcul de $A_{n,n}$ est exact pour les polynômes de degré 2n+1 : elle est d’ordre 2n+1.

La méthode de W.Romberg utilise l'extrapolation de Richardson à partir de 2^n applications de la méthode des trapèzes. Soit $A_{n,0}$ les évaluations de l'intégrale par la méthode des trapèzes $$\begin{align*}A_{0,0} &=\frac{b-a}{2}(f(a)+f(b))\end{align*}$$ $$\begin{align*}A_{1,0} &=\frac{b-a}{4}(f(a)+f(b)+2f(a+\frac{b-a}{2}))\end{align*}$$ $$\begin{align*}A_{n,0} &=\frac{1}{2}A_{n-1,0}+\frac{(b-a)}{2^n}\sum_{k=0}^{2^{n+1}-1}f(a+(2k+1)\frac{b-a}{2^n})\end{align*}$$ Si la dérivée seconde de f est continue bornée sur [a;b], la suite $A_{n,0}$ converge vers la valeur exacte de l'intégrale. Pour accélérer la vitese de la convergence, on aplique l'extrapolation de Richardson, au couple $A_{n,0},A_{n-1,0}$ pour définir $A_{n,1}$ qui converge vers la valeur de l'intégrale si la dérivée quatrième de f est continue bornée. $$\\$$ $A_{n,1}=\frac{4A_{n,0}-A{n-1,0}}{3}$ $$\\$$ De proche en proche, on définit ainsi les valeurs extrapolées $A_{n,l}=\frac{4^lA_{n,l-1}-A_{n-1,l-1}}{4^l-1}$ $$\\$$ Lorsque n tend vers l'infini, on alors $A_{n,l}=\int_{a}^{b}f(x)dx + O(4^{-n(l+1)})$

'C:'

$$\begin{align*} F(t)&=I(h)+a_{2k}h^{2k}+o(h^{2k}) \hspace{1cm} (***)\\ F(t)&=I(\frac{h}{2})+a_{2k}(\frac{h}{2})^{2k}+o(h^{2k})\\ F(t)&=I(\frac{h}{2})+a_{2k}\frac{h^{2k}}{4^k}+o(h^{2k})\\ 4^kF(t)&=4^kI(\frac{h}{2})+a_{2k}h^{2k}+o(h^{2k}) \hspace{1cm} (**)\\ (***)-(**) \hspace{2cm} 4^kF(t)-F(t)&=4^kI(\frac{h}{2})-I(h)+o(h^{2k})\\ F(t)&=\frac{4^kI(\frac{h}{2})-I(h)}{4^k-1}+o(h^{2k})\end{align*}$$

'D'

Notons $A_{n,0}$ les précédents nombres $I(h_n)$ résultats de la méthode des trapèzes qui seront placées dans la première colonne. $A_{n,0}=I(h_n)$

Notons $A_{n,1}$ obtenues par les extrapolations de Richardson de la forme $A(h_n)=\frac{4^kI(\frac{h_n}{2})-I(h_n)}{4^k-1}=\frac{4^kI(h_{n+1})-I(h_n)}{4^k-1}$ a partir $I(h_n)$ et $I(h_n/2)$ avec k$\ge$1 et n$\ge$1 correspondant à la méthode de Simpson placé dans la deuxieme colonne.

Nous pouvons désormais accélérer la convergence au seul moyen des $A_{n,k}$ issus de l'extrapolation de Richardson est obtenir la formule récurrente de Romberg $A_{n,k}=\frac{4^kA_{n,k-1}-A_{n-1,k-1}}{4^k-1}$, c'est une approximation de l'intégrale qui à l'étape $A_{n,n}$ est en $o(h_n^{2n+2})$ ceci à condition que la fonction soit de classe $C^{2n+2}$.Le calcul de $A_{n,n}$ est exact pour les polynômes de degré 2n+1 : elle est d’ordre 2n+1.