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
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 d'une méthode numérique.
H
. 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}$. Remarque: puisque la fonction $\varphi(t,y)=y$ est linéaire, toute méthode implicite peut être rendue explicite par un calcul élémentaire en explicitant directement pour chaque schéma l'expression de $u_{n+1}$. Cependant, nous pouvons utiliser le le module SciPy
sans modifier l'implémentation des schémas (mais on payera l'ordre de convergence de fsolve
).
Estimer l'ordre de convergence des méthodes:
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 écrit les schémas numériques :
tt
(qui change en fonction de h
) uu
.def euler_progressif(phi,tt):
uu = [y0]
for i in range(N):
uu.append(uu[i]+h*phi(tt[i],uu[i]))
return uu
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
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
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
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
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
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
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
def RK4(phi,tt):
uu = [y0]
for i in range(N):
k1 = phi( tt[i], uu[i] )
k2 = phi( tt[i]+h/2 , uu[i]+h*k1/2 )
k3 = phi( tt[i]+h/2 , uu[i]+h*k2/2 )
k4 = phi( tt[i+1] , uu[i]+h*k3 )
uu.append( uu[i] + h*(k1+2*k2+2*k3+k4) /6 )
return uu
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
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
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
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
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
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
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
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
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
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
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
On initialise le problème de Cauchy
t0, y0, tfinal = 0, 1, 3
On définit la solution exacte:
def sol_exacte(t):
return exp(t)
#return exp(-t)
#return sqrt(2.*t+1.)
#return sqrt(t**2+1.)
#return 1./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)$ dépendant des variables $t$ et $y$.
def phi(t,y):
return y
#return -y
#return 1./y
#return t/y
#return y**3
Pour chaque schéma, on calcule la solution approchée avec différentes valeurs de $h_k=1/N_k$, à savoir $N_k=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,...,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)
# erreurs
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 estime la pente de la droite qui relie l'erreur au pas $k$ à l'erreur au pas $k+1$ en echelle logarithmique en utilisant la fonction polyfit
basée sur la régression linéaire.
print ('EE\t %1.2f' %(polyfit(log(H),log(err_ep), 1)[0]))
print ('AB2\t %1.2f' %(polyfit(log(H),log(err_AB2), 1)[0]))
print ('AB3\t %1.2f' %(polyfit(log(H),log(err_AB3), 1)[0]))
print ('AB4\t %1.2f' %(polyfit(log(H),log(err_AB4), 1)[0]))
print ('AB5\t %1.2f' %(polyfit(log(H),log(err_AB5), 1)[0]))
print ('N2\t %1.2f' %(polyfit(log(H),log(err_N2), 1)[0]))
print ('N3\t %1.2f' %(polyfit(log(H),log(err_N3), 1)[0]))
print ('N4\t %1.2f' %(polyfit(log(H),log(err_N4), 1)[0]))
print ('RK4\t %1.2f' %(polyfit(log(H),log(err_RK4), 1)[0]))
print('\n')
print ('EI\t %1.2f' %(polyfit(log(H),log(err_er), 1)[0]))
print ('CN\t %1.2f' %(polyfit(log(H),log(err_CN), 1)[0]))
print ('AM2\t %1.2f' %(polyfit(log(H),log(err_AM2), 1)[0]))
print ('AM3\t %1.2f' %(polyfit(log(H),log(err_AM3), 1)[0]))
print ('AM4\t %1.2f' %(polyfit(log(H),log(err_AM4), 1)[0]))
print ('AM5\t %1.2f' %(polyfit(log(H),log(err_AM5), 1)[0]))
print ('BDF2\t %1.2f' %(polyfit(log(H),log(err_BDF2), 1)[0]))
print ('BDF3\t %1.2f' %(polyfit(log(H),log(err_BDF3), 1)[0]))
print('\n')
print ('EM\t %1.2f' %(polyfit(log(H),log(err_em), 1)[0]))
print ('Heun\t %1.2f' %(polyfit(log(H),log(err_heun), 1)[0]))
print ('AM4AB2\t %1.2f' %(polyfit(log(H),log(err_AM4AB2), 1)[0]))
print ('AM4AB3\t %1.2f' %(polyfit(log(H),log(err_AM4AB3), 1)[0]))
print ('AM4AB4\t %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,figsize=(20,5))
subplot(1,3,1)
loglog(H,err_ep, 'r-o',label='AB-1=EE')
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='RK41')
xlabel('$h$')
ylabel('$e$')
title("Schemas explicites")
legend()
grid(True)
subplot(1,3,2)
loglog(H,err_er, 'r-o',label='AM-0=EI')
loglog(H,err_CN, 'b-v',label='AM-1=CN')
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('$h$')
ylabel('$e$')
title("Schemas implicites")
legend()
grid(True)
subplot(1,3,3)
loglog(H,err_em, 'c-o',label='EM')
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('$h$')
ylabel('$e$')
title("Schemas predicteur-correcteur")
legend()
grid(True)
show()