Implémentation des schémas

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) = 2ty(t), &\forall t \in I=[0,1],\\ y(0) = 1 \end{cases} $$ dont la solution est $y(t)=e^{t^2}$.

On se propose de

  1. calculer la solution approchée obtenue avec la méthode d'Euler explicite avec $h=1/N$ et $N=8$ (pour bien visualiser les erreurs);
  2. même exercice pour les méthodes AB$_2$, AB$_3$, AB$_4$, AB$_5$, N$_2$, N$_3$, N$_4$, RK$_4$;
  3. même exercice pour les méthodes d'Euler implicite, de Crank-Nicolson, AM$_1$, AM$_2$, AM$_3$, AM$_4$, BDF$_2$, BDF$_3$ (pour résoudre à chaque pas de temps l'équation non linéaire on peut implémenter une méthode numérique comme la dichotomie ou la méthode de Newton ou, plus simplement, utiliser le module SciPy);
  4. même exercice pour les méthodes predictor-corrector Euler modifié, Heun, AM$_2$-AB$_1$, AM$_3$-AB$_2$.

On commence par importer les modules math et matplotlib :

In [1]:
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 :

In [2]:
from scipy.optimize import fsolve

On initialise le problème de Cauchy

In [3]:
t0 = 0.
tfinal = 1.
y0 = 1.

On définit la solution exacte:

In [4]:
def sol_exacte(t):
	return y0*math.exp(t**2)

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$.

In [5]:
def phi(t,y):
	return 2.*y*t

On introduit la discrétisation: les nœuds d'intégration $[t_0,t_1,\dots,t_{N}]$ sont contenus dans le vecteur tt

In [6]:
N = 8 
h = (tfinal-t0)/N
tt = [ t0+i*h for i in range(N+1) ]

On écrit les schémas numériques : les valeurs $[u_0,u_1,\dots,u_{N}]$ pour chaque méthode sont contenues dans le vecteur uu.

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} $$

In [7]:
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} $$

In [8]:
def AB2(phi,tt):
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    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} $$

In [9]:
def AB3(phi,tt):
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    uu.append(uu[1]+h*(3.*phi(tt[1],uu[1])-phi(tt[0],uu[0]))/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} $$

In [10]:
def AB4(phi,tt):
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    uu.append(uu[1]+h*(3.*phi(tt[1],uu[1])-phi(tt[0],uu[0]))/2.)
    uu.append(uu[2]+h*(23*phi(tt[2],uu[2])-16*phi(tt[1],uu[1])+5*phi(tt[0],uu[0]))/12.)
    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} $$

In [11]:
def AB5(phi,tt):
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    uu.append(uu[1]+h*(3.*phi(tt[1],uu[1])-phi(tt[0],uu[0]))/2.)
    uu.append(uu[2]+h*(23*phi(tt[2],uu[2])-16*phi(tt[1],uu[1])+5*phi(tt[0],uu[0]))/12.)
    uu.append(uu[3]+h*(55*phi(tt[3],uu[3])-59*phi(tt[2],uu[2])+37*phi(tt[1],uu[1])-9*phi(tt[0],uu[0]))/24.)
    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} $$

In [12]:
def N2(phi,tt):
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    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} $$

In [13]:
def N3(phi,tt):
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    uu.append(uu[0]+2*h*phi(tt[1],uu[1]))
    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} $$

In [14]:
def N4(phi,tt):
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    uu.append(uu[0]+2*h*phi(tt[1],uu[1]))
    uu.append(uu[1]+(7*h*phi(tt[2],uu[2])-2.*h*phi(tt[1],uu[1])+h*phi(tt[0],uu[0]))/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} $$

In [15]:
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)$.

In [16]:
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))$.

In [17]:
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} $$

In [18]:
def AM2(phi,tt):
    uu = [y0]
    temp = fsolve(lambda x: -x+uu[0]+0.5*h*( phi(tt[1],x)+phi(tt[0],uu[0]) ), uu[0])
    uu.append(temp)
    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} $$

In [19]:
def AM3(phi,tt):
	uu = [y0]
	temp = fsolve(lambda x: -x+uu[0]+0.5*h*( phi(tt[1],x)+phi(tt[0],uu[0]) ), uu[0])
	uu.append(temp)
	temp = fsolve(lambda x: -x+uu[1]+h*( 5.*phi(tt[2],x)+8.*phi(tt[1],uu[1])-phi(tt[0],uu[0]) )/12., uu[1])
	uu.append(temp)
	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} $$

In [20]:
def AM4(phi,tt):
	uu = [y0]
	temp = fsolve(lambda x: -x+uu[0]+0.5*h*( phi(tt[1],x)+phi(tt[0],uu[0]) ), uu[0])
	uu.append(temp)
	temp = fsolve(lambda x: -x+uu[1]+h*( 5.*phi(tt[2],x)+8.*phi(tt[1],uu[1])-phi(tt[0],uu[0]) )/12., uu[1])
	uu.append(temp)
	temp = fsolve(lambda x: -x+uu[2]+h*( 9.*phi(tt[3],x)+19.*phi(tt[2],uu[2])-5.*phi(tt[1],uu[1])+phi(tt[0],uu[0]) )/24., uu[1])
	uu.append(temp)
	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

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} $$

In [21]:
def BDF2(phi,tt):
	uu = [y0]
	temp = fsolve(lambda x: -x+uu[0]+h*phi(tt[1],x), uu[0])
	uu.append(temp)
	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} $$

In [22]:
def BDF3(phi,tt):
	uu = [y0]
	temp = fsolve(lambda x: -x+uu[0]+h*phi(tt[1],x), uu[0])
	uu.append(temp)
	temp = fsolve(lambda x: -x+4./3.*uu[1]-1./3.*uu[0] + 2./3.*h*phi(tt[2],x), uu[1])
	uu.append(temp)
	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} $$

In [23]:
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} $$

In [24]:
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} $$

In [25]:
def AM2AB1(phi,tt):
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    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} $$

In [26]:
def AM3AB2(phi,tt):
    uu = [y0]
    uu.append(uu[0]+h*phi(tt[0],uu[0]))
    uu.append(uu[0]+0.5*h*(3.*phi(tt[1],uu[1])-phi(tt[0],uu[0])))
    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

Calculs

On calcule les solutions exacte et approchées:

In [27]:
yy = [sol_exacte(t) for t in tt]

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)

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_BDF2 = BDF2(phi,tt)
uu_BDF3 = BDF3(phi,tt)

uu_em   = euler_modifie(phi,tt)
uu_heun = heun(phi,tt)
uu_AM2AB1 = AM2AB1(phi,tt)
uu_AM3AB2 = AM3AB2(phi,tt)

On compare les graphes des solutions exacte et approchées et on affiche le maximum de l'erreur:

In [28]:
figure()
plot(tt,yy,'b-',tt,uu_ep,'r-D')
title('Euler explicite - max(|erreur|)=%1.10f'%(max([abs(uu_ep[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_AB2,'r-D')
title('AB2 - max(|erreur|)=%1.10f'%(max([abs(uu_AB2[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_AB3,'r-D')
title('AB3 - max(|erreur|)=%1.10f'%(max([abs(uu_AB3[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_AB4,'r-D')
title('AB4 - max(|erreur|)=%1.10f'%(max([abs(uu_AB4[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_AB5,'r-D')
title('AB5 - max(|erreur|)=%1.10f'%(max([abs(uu_AB5[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_N2,'r-D')
title('N2 - max(|erreur|)=%1.10f'%(max([abs(uu_N2[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_N3,'r-D')
title('N3 - max(|erreur|)=%1.10f'%(max([abs(uu_N3[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_N4,'r-D')
title('N4 - max(|erreur|)=%1.10f'%(max([abs(uu_N4[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_RK4,'r-D')
title('RK4 - max(|erreur|)=%1.10f'%(max([abs(uu_RK4[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_er,'r-D')
title('Euler implicite - max(|erreur|)=%1.10f'%(max([abs(uu_er[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_CN,'r-D')
title('Crank Nicolson - max(|erreur|)=%1.10f'%(max([abs(uu_CN[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_AM2,'r-D')
title('AM2 - max(|erreur|)=%1.10f'%(max([abs(uu_AM2[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_AM3,'r-D')
title('AM3 - max(|erreur|)=%1.10f'%(max([abs(uu_AM3[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_AM4,'r-D')
title('AM4 - max(|erreur|)=%1.10f'%(max([abs(uu_AM4[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_BDF2,'r-D')
title('BDF2 - max(|erreur|)=%1.10f'%(max([abs(uu_BDF2[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_BDF3,'r-D')
title('BDF3 - max(|erreur|)=%1.10f'%(max([abs(uu_BDF3[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_em,'r-D')
title('Euler modifie - max(|erreur|)=%1.10f'%(max([abs(uu_em[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_heun,'r-D')
title('Heun - max(|erreur|)=%1.10f'%(max([abs(uu_heun[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_AM2AB1,'r-D')
title('AM2AB1 - max(|erreur|)=%1.10f'%(max([abs(uu_AM2AB1[i]-yy[i]) for i in range(N+1)])))

figure()
plot(tt,yy,'b-',tt,uu_AM3AB2,'r-D')
title('AM3AB2 - max(|erreur|)=%1.10f'%(max([abs(uu_AM3AB2[i]-yy[i]) for i in range(N+1)])))

show()