Étude de la convergence¶
Considérons le problème de Cauchy
trouver la fonction $y \colon I\subset \mathbb{R} \to \mathbb{R}$ définie sur l'intervalle $I=[0,1]$ telle que $$ \begin{cases} y'(t) = y(t), &\forall t \in I=[0,1],\\ y(0) = 1 \end{cases} $$ dont la solution est $y(t)=e^{t}$.
- On se propose d'estimer l'ordre de convergence de la méthode d'Euler explicite.
On calcule d'abord la solution approchée avec différentes valeurs de $h_k=1/N_k$, à savoir $1/2$, $1/4$, $1/8$, \dots, $1/1024$ (ce qui correspond à différentes valeurs de $N_k$, à savoir $2$, $2^2$, $2^3$, ... $2^{10}$).
Pour chaque valeur de $h_k$, on calcule le maximum de la valeur absolue de l'erreur et on sauvegarde toutes ces erreurs dans le vecteur
err_epde sort queerr_ep[k]contient $e_k=\max_{i=0,\dots,N_k}|y(t_i)-u_{i}|$ avec $N_k=2^{k+1}$. Pour estimer l'ordre de convergence on applique alors la formule $$ p_k=\frac{\ln\left| \frac{e_k}{e_{k-1}} \right|}{\ln\left| \frac{h_k}{h_{k-1}} \right|}, \qquad k=1,\dots,10. $$ On trouve ainsi la pente de la droite qui relie l'erreur au pas $k$ à l'erreur au pas $k+1$. Pour estimer la pente globale on pourra utiliser la fonctionpolyfitbasée sur la régression linéaire. - Même exercice pour estimer l'ordre de convergence des méthodes d'Euler modifié, de Heun, AB$_2$, AB$_3$, AB$_4$, AB$_5$, N$_2$, N$_3$, N$_4$, RK$_4$.
- Même exercice pour les méthodes predictor-corrector AM$_4$-AB$_1$, AM$_4$-AB$_2$, AM$_4$-AB$_3$.
- Puisque la fonction $\varphi(t,y)=y$, toute méthode implicite peut être rendu explicite par un calcul élémentaire. On peut alors estimer de la même manière l'ordre de convergence des méthodes d'Euler implicite, de Crank-Nicolson, AM$_1$, AM$_2$, AM$_3$, AM$_4$, BDF$_2$, BDF$_3$ en explicitant directement pour chaque schéma l'expression de $u_{n+1}$. Cependant, nous pouvons utiliser le le module
SciPysans modifier l'implémentation des schémas.
Attention: les schémas multistep ont besoin d'initialiser plusieurs pas de la suite définie pas récurrence pour pouvoir démarrer. Dans cette étude, au lieu d'utiliser un schéma d'ordre inférieur pour initialiser la suite, on utilisera la solution exacte (en effet, l'utilisation d'un schéma d'ordre inférieur dégrade l'ordre de précision).
On commence par importer les modules math et matplotlib :
import math
from matplotlib.pylab import *
%matplotlib inline
Pour résoudre les équations implicites présentes dans les schémas implicites, on peut par exemple importer la fonction fsolve du module scipy.optimize :
from scipy.optimize import fsolve
On initialise le problème de Cauchy
t0, y0, tfinal = 0., 1., 3.
On définit la solution exacte:
def sol_exacte(t):
return math.exp(t)
#return math.exp(-t)
#return math.sqrt(2.*t+1.)
#return math.sqrt(t**2+1.)
#return 1./math.sqrt(1.-2.*t)
On définit l'équation différentielle : phi est une fonction python qui contient la fonction mathématique $\varphi(t, y)=2ty$ dépendant des variables $t$ et $y$.
def phi(t,y):
return y
#return -y
#return 1./y
#return t/y
#return y**3
On écrit les schémas numériques :
les nœuds d'intégration $[t_0,t_1,\dots,t_{N}]$ sont contenus dans le vecteur tt (qui change en fonction de h) et les valeurs $[u_0,u_1,\dots,u_{N}]$ pour chaque méthode sont contenues dans le vecteur uu.
Schémas explicites¶
Schéma d'Euler progressif : $$ \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi(t_n,u_n)& n=0,1,2,\dots N-1 \end{cases} $$
def euler_progressif(phi,tt):
uu = [y0]
for i in range(N):
uu.append(uu[i]+h*phi(tt[i],uu[i]))
return uu
Schéma de Adam-Bashforth d'ordre 2: $$ \begin{cases} u_0=y_0,\\ u_{1}=u_0+h\varphi(t_0,u_0),\\ u_{n+1}=u_n+\frac{h}{2}\Bigl(3\varphi(t_n,u_n)-\varphi(t_{n-1},u_{n-1})\Bigr)& n=1,2,3,4,5,\dots N-1 \end{cases} $$
def AB2(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
for i in range(1,N):
k1 = h * phi( tt[i], uu[i] )
k2 = h * phi( tt[i-1], uu[i-1] )
uu.append( uu[i] + (3.*k1-k2) /2.0 )
return uu
Schéma de Adam-Bashforth d'ordre 3: $$ \begin{cases} u_0=y_0,\\ u_{1}=u_0+h\varphi(t_0,u_0),\\ u_{2}=u_1+\frac{h}{2}\Bigl(3\varphi(t_1,u_1)-\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_n+\frac{h}{12}\Bigl(23\varphi(t_n,u_n)-16\varphi(t_{n-1},u_{n-1})+5\varphi(t_{n-2},u_{n-2})\Bigr)& n=2,3,4,5,\dots N-1 \end{cases} $$
def AB3(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
uu.append(sol_exacte(tt[2]))
for i in range(2,N):
k1 = h * phi( tt[i], uu[i] )
k2 = h * phi( tt[i-1], uu[i-1] )
k3 = h * phi( tt[i-2], uu[i-2] )
uu.append( uu[i] + (23.*k1-16.*k2+5.*k3) /12.0 )
return uu
Schéma de Adam-Bashforth d'ordre 4: $$ \begin{cases} u_0=y_0,\\ u_{1}=u_0+h\varphi(t_0,u_0),\\ u_{2}=u_1+\frac{h}{2}\Bigl(3\varphi(t_1,u_1)-\varphi(t_{0},u_{0})\Bigr),\\ u_{3}=u_2+\frac{h}{12}\Bigl(23\varphi(t_2,u_2)-16\varphi(t_{1},u_{1})+5\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_n+\frac{h}{24}\Bigl(55\varphi(t_n,u_n)-59\varphi(t_{n-1},u_{n-1})+37\varphi(t_{n-2},u_{n-2})-9\varphi(t_{n-3},u_{n-3})\Bigr)& n=3,4,5,\dots N-1 \end{cases} $$
def AB4(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
uu.append(sol_exacte(tt[2]))
uu.append(sol_exacte(tt[3]))
for i in range(3,N):
k1 = h * phi( tt[i], uu[i] )
k2 = h * phi( tt[i-1], uu[i-1] )
k3 = h * phi( tt[i-2], uu[i-2] )
k4 = h * phi( tt[i-3], uu[i-3] )
uu.append( uu[i] + (55.*k1-59.*k2+37.*k3-9.*k4) /24.0 )
return uu
Schéma de Adam-Bashforth d'ordre 5: $$ \begin{cases} u_0=y_0,\\ u_{1}=u_0+h\varphi(t_0,u_0),\\ u_{2}=u_1+\frac{h}{2}\Bigl(3\varphi(t_1,u_1)-\varphi(t_{0},u_{0})\Bigr),\\ u_{3}=u_2+\frac{h}{12}\Bigl(23\varphi(t_2,u_2)-16\varphi(t_{1},u_{1})+5\varphi(t_{0},u_{0})\Bigr),\\ u_{4}=u_3+\frac{h}{24}\Bigl(55\varphi(t_3,u_3)-59\varphi(t_{2},u_{2})+37\varphi(t_{1},u_{1})-9\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_n+\frac{h}{720}\Bigl(1901\varphi(t_n,u_n)-2774\varphi(t_{n-1},u_{n-1})+2616\varphi(t_{n-2},u_{n-2})-1274\varphi(t_{n-3},u_{n-3})+251\varphi(t_{n-4},u_{n-4})\Bigr)& n=4,5,\dots N-1 \end{cases} $$
def AB5(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
uu.append(sol_exacte(tt[2]))
uu.append(sol_exacte(tt[3]))
uu.append(sol_exacte(tt[4]))
for i in range(4,N):
k1 = h * phi( tt[i], uu[i] )
k2 = h * phi( tt[i-1], uu[i-1] )
k3 = h * phi( tt[i-2], uu[i-2] )
k4 = h * phi( tt[i-3], uu[i-3] )
k5 = h * phi( tt[i-4], uu[i-4] )
uu.append( uu[i] + (1901.*k1-2774.*k2+2616.*k3-1274.*k4+251*k5) /720.0 )
return uu
Schéma de Nylström d'ordre 2: $$ \begin{cases} u_0=y_0,\\ u_{1}=u_0+h\varphi(t_0,u_0),\\ u_{n+1}=u_{n-1}+2h\varphi(t_{n},u_{n})& n=1,2,3,4,5,\dots N-1 \end{cases} $$
def N2(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
for i in range(1,N):
k1 = h * phi( tt[i], uu[i] )
uu.append( uu[i-1] + 2.0*k1 )
return uu
Schéma de Nylström d'ordre 3: $$ \begin{cases} u_0=y_0,\\ u_{1}=u_0+h\varphi(t_0,u_0),\\ u_{2}=u_{0}+2h\varphi(t_{1},u_{1}),\\ u_{n+1}=u_{n-1}+\frac{h}{3}\Bigl(7\varphi(t_{n},u_{n})-2\varphi(t_{n-1},u_{n-1})+\varphi(t_{n-2},u_{n-2})\Bigr)& n=2,3,4,5,\dots N-1 \end{cases} $$
def N3(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
uu.append(sol_exacte(tt[2]))
for i in range(2,N):
k1 = h * phi( tt[i], uu[i] )
k2 = h * phi( tt[i-1], uu[i-1] )
k3 = h * phi( tt[i-2], uu[i-2] )
uu.append( uu[i-1] + (7.*k1-2.*k2+k3)/3. )
return uu
Schéma de Nylström d'ordre 4: $$ \begin{cases} u_0=y_0,\\ u_{1}=u_0+h\varphi(t_0,u_0),\\ u_{2}=u_{0}+2h\varphi(t_{1},u_{1}),\\ u_{3}=u_{1}+\frac{h}{3}\Bigl(7\varphi(t_{2},u_{2})-2\varphi(t_{1},u_{1})+\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_{n-1}+\frac{h}{3}\Bigl(8\varphi(t_{n},u_{n})-5\varphi(t_{n-1},u_{n-1})+4\varphi(t_{n-2},u_{n-2})-\varphi(t_{n-3},u_{n-3})\Bigr)& n=3,4,5,\dots N-1 \end{cases} $$
def N4(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
uu.append(sol_exacte(tt[2]))
uu.append(sol_exacte(tt[3]))
for i in range(3,N):
k1 = h * phi( tt[i], uu[i] )
k2 = h * phi( tt[i-1], uu[i-1] )
k3 = h * phi( tt[i-2], uu[i-2] )
k4 = h * phi( tt[i-3], uu[i-3] )
uu.append( uu[i-1] + (8.*k1-5.*k2+4.*k3-k4)/3. )
return uu
Schéma de Runke-Kutta RK4: $$\begin{cases} u_0=y(t_0)=y_0,\\ \tilde u_{n+1/2}=u_n+\frac{h}{2} \varphi(t_{n},u_{n}),\\ \check u_{n+1/2}=u_n+\frac{h}{2} \varphi(t_{n}+\frac{h}{2},\tilde u_{n+1/2}),\\ \hat u_{n+1}=u_n+h\varphi(t_{n+1},\check u_{n+1/2}),\\ u_{n+1}=u_n+\frac{h}{6}\left(\varphi(t_{n},u_{n})+2\varphi(t_{n}+\frac{h}{2},\tilde u_{n+1/2} )+2\varphi(t_{n}+\frac{h}{2}, \check u_{n+1/2})+\varphi(t_{n+1},\hat u_{n+1} \right)& n=0,1,2,\dots N-1 \end{cases} $$
def RK4(phi,tt):
uu = [y0]
for i in range(N):
k1 = h * phi( tt[i], uu[i] )
k2 = h * phi( tt[i]+h/2. , uu[i]+k1/2. )
k3 = h * phi( tt[i]+h/2. , uu[i]+k2/2. )
k4 = h * phi( tt[i+1] , uu[i]+k3 )
uu.append( uu[i] + (k1+2.0*k2+2.0*k3+k4) /6.0 )
return uu
Schémas implicites¶
Schéma d'Euler régressif : $$ \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+h\varphi(t_{n+1},u_{n+1})& n=0,1,2,\dots N-1 \end{cases} $$
Attention : $u_{n+1}$ est solution de l'équation $x=u_n+h\varphi(t_{n+1},x)$, c'est-à-dire un zéro de la fonction (en générale non linéaire) $\lambda \colon x\mapsto -x+u_n+h\varphi(t_{n+1},x)$.
def euler_regressif(phi,tt):
uu = [y0]
for i in range(N):
temp = fsolve(lambda x: -x+uu[i]+h*phi(tt[i+1],x), uu[i])
uu.append(temp)
return uu
Schéma de Crank-Nicolson : $$ \begin{cases} u_0=y_0,\\ u_{n+1}=u_n+\frac{h}{2}\Bigl(\varphi(t_n,u_n)+\varphi(t_{n+1},u_{n+1})\Bigr)& n=0,1,2,\dots N-1 \end{cases} $$
Attention : $u_{n+1}$ est solution de l'équation $x=u_n+\frac{h}{2}(\varphi(t_n,u_n)+\varphi(t_{n+1},x))$, c'est-à-dire un zéro de la fonction (en générale non linéaire) $\lambda \colon x\mapsto -x+u_n+\frac{h}{2}(\varphi(t_n,u_n)+\varphi(t_{n+1},x))$.
def CN(phi,tt):
uu = [y0]
for i in range(len(tt)-1):
temp = fsolve(lambda x: -x+uu[i]+0.5*h*( phi(tt[i+1],x)+phi(tt[i],uu[i]) ), uu[i])
uu.append(temp)
return uu
Schéma de AM-2 : $$ \begin{cases} u_0=y_0,\\ u_1=u_0+\frac{h}{2}\Bigl(\varphi(t_1,u_1)+\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_n+h\Bigl(5\varphi(t_{n+1},u_{n+1})+8\varphi(t_n,u_n)-\varphi(t_{n-1},u_{n-1})\Bigr)& n=1,2,\dots N-1 \end{cases} $$
def AM2(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
for i in range(1,N):
temp = fsolve(lambda x: -x+uu[i]+h*( 5.*phi(tt[i+1],x)+8.*phi(tt[i],uu[i])-phi(tt[i-1],uu[i-1]) )/12., uu[i])
uu.append(temp)
return uu
Schéma de AM-3 : $$ \begin{cases} u_0=y_0,\\ u_1=u_0+\frac{h}{2}\Bigl(\varphi(t_1,u_1)+\varphi(t_{0},u_{0})\Bigr),\\ u_{2}=u_n+h\Bigl(5\varphi(t_{2},u_{2})+8\varphi(t_1,u_1)-\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_n+h\Bigl(9\varphi(t_{n+1},u_{n+1})+19\varphi(t_n,u_n)-5\varphi(t_{n-1},u_{n-1})+\varphi(t_{n-2},u_{n-2})\Bigr)& n=2,3,\dots N-1 \end{cases} $$
def AM3(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
uu.append(sol_exacte(tt[2]))
for i in range(2,N):
temp = fsolve(lambda x: -x+uu[i]+h*( 9.*phi(tt[i+1],x)+19.*phi(tt[i],uu[i])-5.*phi(tt[i-1],uu[i-1])+phi(tt[i-2],uu[i-2]) )/24., uu[i])
uu.append(temp)
return uu
Schéma de AM-4 : $$ \begin{cases} u_0=y_0,\\ u_1=u_0+\frac{h}{2}\Bigl(\varphi(t_1,u_1)+\varphi(t_{0},u_{0})\Bigr),\\ u_{2}=u_1+h\Bigl(5\varphi(t_{2},u_{2})+8\varphi(t_1,u_1)-\varphi(t_{0},u_{0})\Bigr),\\ u_{3}=u_2+h\Bigl(9\varphi(t_{3},u_{3})+19\varphi(t_2,u_2)-5\varphi(t_{1},u_{1})+\varphi(t_{0},u_{0})\Bigr),\\ u_{n+1}=u_n+h\Bigl(251\varphi(t_{n+1},u_{n+1})+646\varphi(t_n,u_n)-264\varphi(t_{n-1},u_{n-1})+106\varphi(t_{n-2},u_{n-2})-19\varphi(t_{n-3},u_{n-3})\Bigr)& n=3,4,\dots N-1 \end{cases} $$
def AM4(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
uu.append(sol_exacte(tt[2]))
uu.append(sol_exacte(tt[3]))
for i in range(3,N):
temp = fsolve(lambda x: -x+uu[i]+h*( 251.*phi(tt[i+1],x)+646.*phi(tt[i],uu[i])-264.*phi(tt[i-1],uu[i-1])+106.*phi(tt[i-2],uu[i-2])-19.*phi(tt[i-3],uu[i-3]) )/720., uu[i])
uu.append(temp)
return uu
def AM5(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
uu.append(sol_exacte(tt[2]))
uu.append(sol_exacte(tt[3]))
uu.append(sol_exacte(tt[4]))
for i in range(4,N):
temp = fsolve(lambda x: -x+uu[i]+h*( 475.*phi(tt[i+1],x)+1427.*phi(tt[i],uu[i])-798.*phi(tt[i-1],uu[i-1])+482.*phi(tt[i-2],uu[i-2])-173.*phi(tt[i-3],uu[i-3])+27.*phi(tt[i-4],uu[i-4]))/1440., uu[i])
uu.append(temp)
return uu
Schéma BDF2 : $$ \begin{cases} u_0=y_0,\\ u_1=u_0+h\varphi(t_1,u_1),\\ u_{n+1}=\frac{4}{3}u_n-\frac{1}{3}u_{n-1}+\frac{2}{3}\varphi(t_{n+1},u_{n+1})& n=1,2,\dots N-1 \end{cases} $$
def BDF2(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
for i in range(1,N):
temp = fsolve(lambda x: -x+4./3.*uu[i]-1./3.*uu[i-1] + 2./3.*h*phi(tt[i+1],x) , uu[i])
uu.append(temp)
return uu
Schéma BDF3 : $$ \begin{cases} u_0=y_0,\\ u_1=u_0+h\varphi(t_1,u_1),\\ u_{2}=\frac{4}{3}u_1-\frac{1}{3}u_{0}+\frac{2}{3}\varphi(t_{2},u_{2}),\\ u_{n+1}=\frac{4}{3}u_n-\frac{1}{3}u_{n-1}+\frac{2}{3}\varphi(t_{n+1},u_{n+1})& n=2,3,\dots N-1 \end{cases} $$
def BDF3(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
uu.append(sol_exacte(tt[2]))
for i in range(2,N):
temp = fsolve(lambda x: -x+18./11.*uu[i]-9./11.*uu[i-1] + 2./11.*uu[i-2]+6./11.*h*phi(tt[i+1],x) , uu[i])
uu.append(temp)
return uu
Schémas predicteur-correcteur¶
Schéma d'Euler modifié: $$ \begin{cases} u_0=y_0,\\ \tilde u = u_n+\frac{h}{2}\varphi(t_n,u_n),\\ u_{n+1}=u_n+h\varphi\left(t_n+\frac{h}{2},\tilde u\right)& n=0,1,2,\dots N-1 \end{cases} $$
def euler_modifie(phi,tt):
uu = [y0]
for i in range(N):
k1 = h * phi( tt[i], uu[i] )
uu.append( uu[i]+h*phi(tt[i]+h/2.,uu[i]+k1/2.) )
return uu
Schéma de Heun: $$ \begin{cases} u_0=y_0,\\ \tilde u = u_n+h\varphi(t_n,u_n)\\ u_{n+1}=u_n+\frac{h}{2}\Bigl(\varphi(t_n,u_n)+\varphi(t_{n+1},\tilde u)\Bigr)& n=0,1,2,\dots N-1 \end{cases} $$
def heun(phi,tt):
uu = [y0]
for i in range(N):
k1 = h * phi( tt[i], uu[i] )
k2 = h * phi( tt[i+1], uu[i] + k1 )
uu.append( uu[i] + (k1+k2) /2.0 )
return uu
Schéma AM-2 AB-1 : $$ \begin{cases} u_0=y_0,\\ u_1=u_0+h\varphi(t_0,u_0),\\ \tilde u=u_n+h\varphi(t_n,u_n),\\ u_{n+1}=u_n+\frac{h}{12}\Bigl(5\varphi(t_{n+1},\tilde u)+8\varphi(t_n,u_n)-\varphi(t_{n-1},u_{n-1})\Bigr)& n=1,2,\dots N-1 \end{cases} $$
def AM2AB1(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
for i in range(1,N):
pred = uu[i] + h*phi(tt[i],uu[i])
uu.append(uu[i]+h*(5.*phi(tt[i+1],pred)+8.*phi(tt[i],uu[i])-phi(tt[i-1],uu[i-1]))/12.)
return uu
Schéma AM-3 AB-2 : $$ \begin{cases} u_0=y_0,\\ u_1=u_0+h\varphi(t_0,u_0),\\ u_2=u_0+\frac{h}{2}(3\varphi(t_1,u_1)-\varphi(t_{0},u_{0})),\\ \tilde u=u_n+\frac{h}{2}(3\varphi(t_n,u_n)-\varphi(t_{n-1},u_{n-1})),\\ u_{n+1}=u_n+\frac{h}{24}\Bigl(9\varphi(t_{n+1},\tilde u)+19\varphi(t_n,u_n)-5\varphi(t_{n-1},u_{n-1})+\varphi(t_{n-2},u_{n-2})\Bigr)& n=2,3,\dots N-1 \end{cases} $$
def AM3AB2(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
uu.append(sol_exacte(tt[2]))
for i in range(2,N):
k1 = h * phi( tt[i], uu[i] )
k2 = h * phi( tt[i-1], uu[i-1] )
pred = uu[i] + (3.*k1-k2) /2.0
uu.append(uu[i]+h*(5.*phi(tt[i+1],pred)+8.*phi(tt[i],uu[i])-phi(tt[i-1],uu[i-1]))/12.)
return uu
def AM4AB2(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
uu.append(sol_exacte(tt[2]))
uu.append(sol_exacte(tt[3]))
for i in range(3,N):
k1 = h * phi( tt[i], uu[i] )
k2 = h * phi( tt[i-1], uu[i-1] )
pred = uu[i] + (3.*k1-k2) /2.0
uu.append(uu[i]+h*( 251.*phi(tt[i+1],pred)+646.*phi(tt[i],uu[i])-264.*phi(tt[i-1],uu[i-1])+106.*phi(tt[i-2],uu[i-2])-19.*phi(tt[i-3],uu[i-3]) )/720.)
return uu
def AM4AB3(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
uu.append(sol_exacte(tt[2]))
uu.append(sol_exacte(tt[3]))
for i in range(3,N):
k1 = h * phi( tt[i], uu[i] )
k2 = h * phi( tt[i-1], uu[i-1] )
k3 = h * phi( tt[i-2], uu[i-2] )
pred = uu[i] + (23.*k1-16.*k2+5.*k3) /12.0
uu.append(uu[i]+h*( 251.*phi(tt[i+1],pred)+646.*phi(tt[i],uu[i])-264.*phi(tt[i-1],uu[i-1])+106.*phi(tt[i-2],uu[i-2])-19.*phi(tt[i-3],uu[i-3]) )/720.)
return uu
def AM4AB4(phi,tt):
uu = [y0]
uu.append(sol_exacte(tt[1]))
uu.append(sol_exacte(tt[2]))
uu.append(sol_exacte(tt[3]))
for i in range(3,N):
k1 = h * phi( tt[i], uu[i] )
k2 = h * phi( tt[i-1], uu[i-1] )
k3 = h * phi( tt[i-2], uu[i-2] )
k4 = h * phi( tt[i-3], uu[i-3] )
pred = uu[i] + (55.*k1-59.*k2+37.*k3-9.*k4) /24.0
uu.append(uu[i]+h*( 251.*phi(tt[i+1],pred)+646.*phi(tt[i],uu[i])-264.*phi(tt[i-1],uu[i-1])+106.*phi(tt[i-2],uu[i-2])-19.*phi(tt[i-3],uu[i-3]) )/720.)
return uu
Calculs¶
Pour chaque schéma, on calcule la solution approchée avec différentes valeurs de $h_k=1/N_k$, à savoir $1/2$, $1/4$, $1/8$, \dots, $1/1024$ (ce qui correspond à différentes valeurs de $N_k=2^{k+1}$, à savoir $2$, $2^2$, $2^3$, ... $2^{10}$). On sauvegarde les valeurs de $h_k$ dans le vecteur H.
Pour chaque valeur de $h_k$, on calcule le maximum de la valeur absolue de l'erreur et on sauvegarde toutes ces erreurs dans le vecteur err_schema de sort que err_schema[k] contient $e_k=\max_{i=0,\dots,N_k}|y(t_i)-u_{i}|$ avec $N_k=2^{k+1}$.
H = []
err_ep = []
err_AB2 = []
err_AB3 = []
err_AB4 = []
err_AB5 = []
err_N2 = []
err_N3 = []
err_N4 = []
err_RK4 = []
err_er = []
err_CN = []
err_AM2 = []
err_AM3 = []
err_AM4 = []
err_AM5 = []
err_BDF2 = []
err_BDF3 = []
err_em = []
err_heun = []
err_AM4AB2 = []
err_AM4AB3 = []
err_AM4AB4 = []
for k in range(8):
N = 2**(k+3)
h = (tfinal-t0)/N
H.append(h)
tt = [t0+i*h for i in range(N+1)]
yy = [sol_exacte(t) for t in tt]
# schemas explicites
uu_ep = euler_progressif(phi,tt)
uu_AB2 = AB2(phi,tt)
uu_AB3 = AB3(phi,tt)
uu_AB4 = AB4(phi,tt)
uu_AB5 = AB5(phi,tt)
uu_N2 = N2(phi,tt)
uu_N3 = N3(phi,tt)
uu_N4 = N4(phi,tt)
uu_RK4 = RK4(phi,tt)
# schemas implicites
uu_er = euler_regressif(phi,tt)
uu_CN = CN(phi,tt)
uu_AM2 = AM2(phi,tt)
uu_AM3 = AM3(phi,tt)
uu_AM4 = AM4(phi,tt)
uu_AM5 = AM5(phi,tt)
uu_BDF2 = BDF2(phi,tt)
uu_BDF3 = BDF3(phi,tt)
# schemas predictor-corrector
uu_em = euler_modifie(phi,tt)
uu_heun = heun(phi,tt)
uu_AM4AB2 = AM4AB2(phi,tt)
uu_AM4AB3 = AM4AB3(phi,tt)
uu_AM4AB4 = AM4AB4(phi,tt)
err_ep.append(max([abs(uu_ep[i]-yy[i]) for i in range(N+1)]))
err_AB2.append(max([abs(uu_AB2[i]-yy[i]) for i in range(N+1)]))
err_AB3.append(max([abs(uu_AB3[i]-yy[i]) for i in range(N+1)]))
err_AB4.append(max([abs(uu_AB4[i]-yy[i]) for i in range(N+1)]))
err_AB5.append(max([abs(uu_AB5[i]-yy[i]) for i in range(N+1)]))
err_N2.append(max([abs(uu_N2[i]-yy[i]) for i in range(N+1)]))
err_N3.append(max([abs(uu_N3[i]-yy[i]) for i in range(N+1)]))
err_N4.append(max([abs(uu_N4[i]-yy[i]) for i in range(N+1)]))
err_RK4.append(max([abs(uu_RK4[i]-yy[i]) for i in range(N+1)]))
err_er.append(max([abs(uu_er[i]-yy[i]) for i in range(N+1)]))
err_CN.append(max([abs(uu_CN[i]-yy[i]) for i in range(N+1)]))
err_AM2.append(max([abs(uu_AM2[i]-yy[i]) for i in range(N+1)]))
err_AM3.append(max([abs(uu_AM3[i]-yy[i]) for i in range(N+1)]))
err_AM4.append(max([abs(uu_AM4[i]-yy[i]) for i in range(N+1)]))
err_AM5.append(max([abs(uu_AM5[i]-yy[i]) for i in range(N+1)]))
err_BDF2.append(max([abs(uu_BDF2[i]-yy[i]) for i in range(N+1)]))
err_BDF3.append(max([abs(uu_BDF3[i]-yy[i]) for i in range(N+1)]))
err_em.append(max([abs(uu_em[i]-yy[i]) for i in range(N+1)]))
err_heun.append(max([abs(uu_heun[i]-yy[i]) for i in range(N+1)]))
err_AM4AB2.append(max([abs(uu_AM4AB2[i]-yy[i]) for i in range(N+1)]))
err_AM4AB3.append(max([abs(uu_AM4AB3[i]-yy[i]) for i in range(N+1)]))
err_AM4AB4.append(max([abs(uu_AM4AB4[i]-yy[i]) for i in range(N+1)]))
Pour estimer l'ordre de convergence on applique alors la formule
$$
p_k=\frac{\ln\left| \frac{e_k}{e_{k-1}} \right|}{\ln\left| \frac{h_k}{h_{k-1}} \right|},
\qquad k=1,\dots,10.
$$
On trouve ainsi la pente de la droite qui relie l'erreur au pas $k$ à l'erreur au pas $k+1$.
Pour estimer la pente globale on utilise la fonction polyfit basée sur la régression linéaire.
print ('Euler progressif %1.2f' %(polyfit(log(H),log(err_ep), 1)[0]))
print ('AB2 %1.2f' %(polyfit(log(H),log(err_AB2), 1)[0]))
print ('AB3 %1.2f' %(polyfit(log(H),log(err_AB3), 1)[0]))
print ('AB4 %1.2f' %(polyfit(log(H),log(err_AB4), 1)[0]))
print ('AB5 %1.2f' %(polyfit(log(H),log(err_AB5), 1)[0]))
print ('N2 %1.2f' %(polyfit(log(H),log(err_N2), 1)[0]))
print ('N3 %1.2f' %(polyfit(log(H),log(err_N3), 1)[0]))
print ('N4 %1.2f' %(polyfit(log(H),log(err_N4), 1)[0]))
print ('RK4 %1.2f' %(polyfit(log(H),log(err_RK4), 1)[0]))
print('\n')
print ('Euler retrograde %1.2f' %(polyfit(log(H),log(err_er), 1)[0]))
print ('CN %1.2f' %(polyfit(log(H),log(err_CN), 1)[0]))
print ('AM2 %1.2f' %(polyfit(log(H),log(err_AM2), 1)[0]))
print ('AM3 %1.2f' %(polyfit(log(H),log(err_AM3), 1)[0]))
print ('AM4 %1.2f' %(polyfit(log(H),log(err_AM4), 1)[0]))
print ('AM5 %1.2f' %(polyfit(log(H),log(err_AM5), 1)[0]))
print ('BDF2 %1.2f' %(polyfit(log(H),log(err_BDF2), 1)[0]))
print ('BDF3 %1.2f' %(polyfit(log(H),log(err_BDF3), 1)[0]))
print('\n')
print ('Euler modifie %1.2f' %(polyfit(log(H),log(err_em), 1)[0]))
print ('Heun %1.2f' %(polyfit(log(H),log(err_heun), 1)[0]))
print ('AM4AB2 %1.2f' %(polyfit(log(H),log(err_AM4AB2), 1)[0]))
print ('AM4AB3 %1.2f' %(polyfit(log(H),log(err_AM4AB3), 1)[0]))
print ('AM4AB4 %1.2f' %(polyfit(log(H),log(err_AM4AB4), 1)[0]))
Pour afficher l'ordre de convergence on utilise une échelle logarithmique : on représente $\ln(h)$ sur l'axe des abscisses et $\ln(\text{err})$ sur l'axe des ordonnées. Le but de cette représentation est clair: si $\text{err}=Ch^p$ alors $\ln(\text{err})=\ln(C)+p\ln(h)$. En échelle logarithmique, $p$ représente donc la pente de la ligne droite $\ln(\text{err})$.
figure(1)
loglog(H,err_ep, 'r-o',label='Euler Explicite')
loglog(H,err_AB2, 'g-+',label='AB-2')
loglog(H,err_AB3, 'c-D',label='AB-3')
loglog(H,err_AB4, 'y-*',label='AB-4')
loglog(H,err_AB5, 'r-.',label='AB-5')
loglog(H,err_N2, 'y-*',label='N-2')
loglog(H,err_N3, 'r-<',label='N-3')
loglog(H,err_N4, 'b-+',label='N-4')
loglog(H,err_RK4, 'g-o',label='RK4')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
title("Schemas explicites")
legend(bbox_to_anchor=(1.04,1),loc='upper left')
grid(True)
figure(2)
loglog(H,err_er, 'r-o',label='Euler Implicite')
loglog(H,err_CN, 'b-v',label='Crank-Nicolson')
loglog(H,err_AM2, 'm->',label='AM-2')
loglog(H,err_AM3, 'c-D',label='AM-3')
loglog(H,err_AM4, 'g-+',label='AM-4')
loglog(H,err_AM5, 'b->',label='AM-5')
loglog(H,err_BDF2, 'y-*',label='BDF-2')
loglog(H,err_BDF3, 'r-<',label='BDF-3')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
title("Schemas implicites")
legend(bbox_to_anchor=(1.04,1),loc='upper left')
grid(True)
figure(3)
loglog(H,err_em, 'c-o',label='Euler modifie')
loglog(H,err_heun, 'y->',label='Heun')
loglog(H,err_AM4AB2, 'r-<',label='AM4-AB2')
loglog(H,err_AM4AB3, 'b-v',label='AM4-AB3')
loglog(H,err_AM4AB4, 'm-*',label='AM4-AB4')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
title("Schemas predicteur-correcteur")
legend(bbox_to_anchor=(1.04,1),loc='upper left')
grid(True)
show()
