from IPython.display import display, Latex
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
Vous pouvez commencer à exécuter un serveur de notebook en écrivant dans un terminal :
$$\text{`jupyter notebook &`}$$
Cela imprimera des informations dans votre console et ouvrira un navigateur Web.
La page d'arrivée est le tableau de bord qui affiche les notebook actuellement disponibles (par défaut, le répertoire à partir duquel le serveur du bloc-notes a été démarré).
Vous pouvez
Nouveau
Pour ce TP, on importera tous les modules necessaires comme suit:
%reset -f
%matplotlib inline
%autosave 300
from math import *
from matplotlib.pylab import *
from scipy.optimize import fsolve
On écrit les schémas numériques :
tt
uu
.Schéma d'Euler progressif = de Adam-Bashforth à 1 pas: $$ \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):
h = tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
uu.append(uu[i]+h*phi(tt[i],uu[i]))
return uu
Schéma de Adam-Bashforth à 2 pas: $$ \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):
h = tt[1]-tt[0]
uu = [y0]
uu.append(uu[0]+h*phi(tt[0],uu[0]))
for i in range(1,len(tt)-1):
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 à 3 pas: $$ \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):
h = tt[1]-tt[0]
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,len(tt)-1):
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 à 4 pas: $$ \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):
h = tt[1]-tt[0]
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,len(tt)-1):
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 à 5 pas: \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):
h = tt[1]-tt[0]
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,len(tt)-1):
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 à 2 pas: $$ \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):
h = tt[1]-tt[0]
uu = [y0]
uu.append(uu[0]+h*phi(tt[0],uu[0]))
for i in range(1,len(tt)-1):
k1 = h * phi( tt[i], uu[i] )
uu.append( uu[i-1] + 2.0*k1 )
return uu
Schéma de Nylström à 3 pas: $$ \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):
h = tt[1]-tt[0]
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,len(tt)-1):
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 à 4 pas: $$ \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):
h = tt[1]-tt[0]
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,len(tt)-1):
k1 = phi( tt[i], uu[i] )
k2 = phi( tt[i-1], uu[i-1] )
k3 = phi( tt[i-2], uu[i-2] )
k4 = phi( tt[i-3], uu[i-3] )
uu.append( uu[i-1] + (8.*k1-5.*k2+4.*k3-k4)*h/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):
h = tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
k1 = phi( tt[i], uu[i] )
k2 = phi( tt[i]+h/2. , uu[i]+k1/2. )
k3 = phi( tt[i]+h/2. , uu[i]+k2/2. )
k4 = phi( tt[i+1] , uu[i]+k3 )
uu.append( uu[i] + (k1+2.0*k2+2.0*k3+k4)*h/6 )
return uu
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):
h = tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
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):
h = tt[1]-tt[0]
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):
h = tt[1]-tt[0]
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,len(tt)-1):
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):
h = tt[1]-tt[0]
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,len(tt)-1):
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):
h = tt[1]-tt[0]
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,len(tt)-1):
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} $$
def BDF2(phi,tt):
h = tt[1]-tt[0]
uu = [y0]
temp = fsolve(lambda x: -x+uu[0]+h*phi(tt[1],x), uu[0])
uu.append(temp)
for i in range(1,len(tt)-1):
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):
h = tt[1]-tt[0]
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,len(tt)-1):
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é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):
h = tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
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):
h = tt[1]-tt[0]
uu = [y0]
for i in range(len(tt)-1):
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):
h = tt[1]-tt[0]
uu = [y0]
uu.append(uu[0]+h*phi(tt[0],uu[0]))
for i in range(1,len(tt)-1):
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):
h = tt[1]-tt[0]
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,len(tt)-1):
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
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
On définit la solution exacte:
sol_exacte = lambda t,y0 : y0*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$.
phi = lambda t,y : 2.*y*t
On initialise le problème de Cauchy
t0 = 0.
tfinal = 1.
y0 = 1.
On introduit la discrétisation: les nœuds d'intégration $[t_0,t_1,\dots,t_{N}]$ sont contenus dans le vecteur tt
N = 8
h = (tfinal-t0)/N
tt = linspace(t0,tfinal,N+1)
On calcule les solutions exacte et approchées:
yy = [sol_exacte(t,y0) 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:
figure(1,figsize=(20,16))
subplot(4,5,1)
plot(tt,yy,'b-',tt,uu_ep,'r-D')
title('AB1=EE - max(|err|)=%1.3f'%(max([abs(uu_ep[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,2)
plot(tt,yy,'b-',tt,uu_AB2,'r-D')
title('AB2 - max(|err|)=%1.3f'%(max([abs(uu_AB2[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,3)
plot(tt,yy,'b-',tt,uu_AB3,'r-D')
title('AB3 - max(|err|)=%1.3f'%(max([abs(uu_AB3[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,4)
plot(tt,yy,'b-',tt,uu_AB4,'r-D')
title('AB4 - max(|err|)=%1.3f'%(max([abs(uu_AB4[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,5)
plot(tt,yy,'b-',tt,uu_AB5,'r-D')
title('AB5 - max(|err|)=%1.3f'%(max([abs(uu_AB5[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,6)
plot(tt,yy,'b-',tt,uu_N2,'r-D')
title('N2 - max(|err|)=%1.3f'%(max([abs(uu_N2[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,7)
plot(tt,yy,'b-',tt,uu_N3,'r-D')
title('N3 - max(|err|)=%1.3f'%(max([abs(uu_N3[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,8)
plot(tt,yy,'b-',tt,uu_N4,'r-D')
title('N4 - max(|err|)=%1.3f'%(max([abs(uu_N4[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,9)
plot(tt,yy,'b-',tt,uu_RK4,'r-D')
title('RK4 - max(|err|)=%1.3f'%(max([abs(uu_RK4[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,10)
plot(tt,yy,'b-',tt,uu_er,'r-D')
title('AM0=EI - max(|err|)=%1.3f'%(max([abs(uu_er[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,11)
plot(tt,yy,'b-',tt,uu_CN,'r-D')
title('AM1=CN - max(|err|)=%1.3f'%(max([abs(uu_CN[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,12)
plot(tt,yy,'b-',tt,uu_AM2,'r-D')
title('AM2 - max(|err|)=%1.3f'%(max([abs(uu_AM2[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,13)
plot(tt,yy,'b-',tt,uu_AM3,'r-D')
title('AM3 - max(|err|)=%1.3f'%(max([abs(uu_AM3[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,14)
plot(tt,yy,'b-',tt,uu_AM4,'r-D')
title('AM4 - max(|err|)=%1.3f'%(max([abs(uu_AM4[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,15)
plot(tt,yy,'b-',tt,uu_BDF2,'r-D')
title('BDF2 - max(|err|)=%1.3f'%(max([abs(uu_BDF2[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,16)
plot(tt,yy,'b-',tt,uu_BDF3,'r-D')
title('BDF3 - max(|err|)=%1.3f'%(max([abs(uu_BDF3[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,17)
plot(tt,yy,'b-',tt,uu_em,'r-D')
title('Euler modifie - max(|err|)=%1.3f'%(max([abs(uu_em[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,18)
plot(tt,yy,'b-',tt,uu_heun,'r-D')
title('Heun - max(|err|)=%1.3f'%(max([abs(uu_heun[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,19)
plot(tt,yy,'b-',tt,uu_AM2AB1,'r-D')
title('AM2AB1 - max(|err|)=%1.3f'%(max([abs(uu_AM2AB1[i]-yy[i]) for i in range(N+1)])))
subplot(4,5,20)
plot(tt,yy,'b-',tt,uu_AM3AB2,'r-D')
title('AM3AB2 - max(|err|)=%1.3f'%(max([abs(uu_AM3AB2[i]-yy[i]) for i in range(N+1)])));