from IPython.core.display import HTML, display
css_file = './custom.css'
HTML(open(css_file, "r").read())
M62 DM 2020¶
Exercice : implémentation et comparaison de schémas¶
Le déplacement $x(t)$ d’un système oscillant composé d’une masse et d’un ressort, soumis à une force de frottement proportionnelle à la vitesse, est décrit par l’équation différentielle du second ordre $$ x''(t)+ 5x'(t)+6x(t)=0, $$ avec $x(0) = 1$ et $x'(0) = 0$, pour $t \in [0, 5]$.
Q1 [1 point]
Écrire le système sous forme matricielle. Calculer ensuite la solution exacte avecsympy. Afficher $t\mapsto x(t)$, $t\mapsto x'(t)$ et $x\mapsto x'(x)$.
Si on note $y(t)=x(t)$ et $z(t)=x'(t)$ on a $y(0)=1$, $z(0)=0$ et $$ \begin{cases} y'(t)=z(t),\\ z'(t)=-5z(t)-6y(t) \end{cases} \quad\text{i.e.}\quad \begin{pmatrix}y\\z\end{pmatrix}'(t) {}= \begin{pmatrix}0&1\\-6&-5\end{pmatrix} \begin{pmatrix}y\\z\end{pmatrix}(t) $$
%reset -f
%matplotlib inline
from matplotlib.pylab import *
# VARIABLES GLOBALES
t0 = 0
tfinal = 5
N = 50 #300
tt = linspace(t0,tfinal,N+1)
y0 = 1
z0 = 0
phi1 = lambda t,y,z : z
phi2 = lambda t,y,z : -6*y-5*z
def affichage(tt,yy,zz,uu,ww,s):
subplot(1,2,1)
plot(tt,uu,'--',label=r'$u(t) \approx y(t)$',color="red")
plot(tt,ww,'--',label=r'$w(t) \approx z(t)$',color="blue")
plot(tt,yy,label='$y(t)$',color="red")
plot(tt,zz,label='$z(t)$',color="blue")
xlabel('t')
legend()
title(f'{s} - $y(t)$ et $z(t)$')
grid()
subplot(1,2,2)
plot(uu,ww,'--',label='Approchée')
plot(yy,zz,label='Exacte')
xlabel('y')
ylabel('z')
legend()
title(f'{s} - z(y)')
grid()
axis('equal');
La solution exacte est
\begin{align*}
x(t)=y(t)&=3e^{-2t}−2e^{-3t},\\
x'(t)=z(t)&=-6e^{-2t}+6e^{-3t}.
\end{align*}
Vérifions-le avec le module sympy:
import sympy as sym
sym.init_printing()
t = sym.Symbol('t')
y = sym.Function('y')
z = sym.Function('z')
edo1 = sym.Eq( sym.diff(y(t),t) , phi1(t,y(t),z(t)) )
edo2 = sym.Eq( sym.diff(z(t),t) , phi2(t,y(t),z(t)) )
display(edo1)
display(edo2)
solgen = sym.dsolve([edo1,edo2],[y(t),z(t)])
# display(solgen)
consts = sym.solve( [ sym.Eq( y0, solgen[0].rhs.subs(t,t0)) , sym.Eq( z0, solgen[1].rhs.subs(t,t0)) ] , dict=True)[0]
# display(consts)
solpar_1 = solgen[0].subs(consts)
solpar_2 = solgen[1].subs(consts)
display(solpar_1)
display(solpar_2)
func_1 = sym.lambdify(t,solpar_1.rhs,'numpy')
func_2 = sym.lambdify(t,solpar_2.rhs,'numpy')
from matplotlib.pylab import *
figure(figsize=(17,7))
yy = func_1(tt)
zz = func_2(tt)
subplot(1,2,1)
plot(tt,yy,label=r'$t\mapsto y(t)$',color="red")
plot(tt,zz,label=r'$t\mapsto z(t)$',color="blue")
legend()
grid()
subplot(1,2,2)
plot(yy,zz)
xlabel(r'$y$')
ylabel(r'$z$')
grid()
axis('equal');
Q2 [2 points]
Calculer la solution approchée obtenue par la méthode d'Euler Progressif avec $301$ points. Afficher $t\mapsto x(t)$, $t\mapsto x'(t)$ et $x\mapsto x'(x)$ en comparant solution exacte et solution approchée.
On notera $u_n\approx y_n=x(t_n)$ et $w_n\approx z(t_n)$.
Euler explicite $$ \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+h\varphi_1(t_n,u_n,w_n),\\ w_{n+1}=w_n+h\varphi_2(t_n,u_n,w_n). \end{cases} $$
def euler_progressif(phi1,phi2,tt,y0,z0):
uu = [y0]
ww = [z0]
h = tt[1]-tt[0]
for i in range(len(tt)-1):
uu.append(uu[i]+h*phi1(tt[i],uu[i],ww[i]))
ww.append(ww[i]+h*phi2(tt[i],uu[i],ww[i]))
return [uu,ww]
[uu, ww] = euler_progressif(phi1,phi2,tt,y0,z0)
figure(figsize=(18,7))
affichage(tt,yy,zz,uu,ww,"Euler explicite")
Q3 [2 points]
Calculer la solution approchée obtenue par la méthode d'Euler Regressif avec $301$ points. Afficher $t\mapsto x(t)$, $t\mapsto x'(t)$ et $x\mapsto x'(x)$ en comparant solution exacte et solution approchée.
Euler implicite $$ \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+h\varphi_1(t_{n+1},u_{n+1},w_{n+1}),\\ w_{n+1}=w_n+h\varphi_2(t_{n+1},u_{n+1},w_{n+1}), \end{cases} $$ qu'on peut rendre explicite par un calcul élementaire: $$ w_{n+1} = w_n + h(-6u_{n+1}-5w_{n+1}) = w_n + h(-6(u_{n}+hw_{n+1})-5w_{n+1}) \quad\implies\quad (1+5h+6h^2)w_{n+1} = w_n-6hu_{n}$$ ainsi $$ \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+h\dfrac{w_n-6hu_{n}}{1+5h+6h^2},\\ w_{n+1}=\dfrac{w_n-6hu_{n}}{1+5h+6h^2}. \end{cases} $$
# VERSION explicite
def euler_regressif(phi1,phi2,tt,y0,z0):
uu = [y0]
ww = [z0]
h = tt[1]-tt[0]
for i in range(len(tt)-1):
uu.append( uu[i]+h*(ww[i]-6*h*uu[i])/(1+5*h+6*h**2) )
ww.append( (ww[i]-6*h*uu[i])/(1+5*h+6*h**2) )
return [uu,ww]
# VERSION AVEC fsolve
# from scipy.optimize import fsolve
# def euler_regressif(phi1,phi2,tt,y0,z0):
# uu = [y0]
# ww = [z0]
# h = tt[1]-tt[0]
# for i in range(len(tt)-1):
# systeme = lambda z : [ -z[0]+uu[i]+h*phi1(tt[i+1],z[0],z[1]) , -z[1]+ww[i]+h*phi2(tt[i+1],z[0],z[1]) ]
# utemp,wtemp = fsolve( systeme , (uu[i],ww[i]) )
# uu.append(utemp)
# ww.append(wtemp)
# return [uu,ww]
[uu, ww] = euler_regressif(phi1,phi2,tt,y0,z0)
figure(figsize=(18,7))
affichage(tt,yy,zz,uu,ww,"Euler implicite")
Q4 [1 point]
Calculer la solution approchée obtenue par la fonctionodeintdu modulescipy. Afficher $t\mapsto x(t)$, $t\mapsto x'(t)$ et $x\mapsto x'(x)$ en comparant solution exacte et solution approchée.
from scipy.integrate import odeint
# Bien noter l'ordre des variables dans phi
pphi = lambda t,yy : [ phi1(t,yy[0],yy[1]), phi2(t,yy[0],yy[1]) ]
yy0 = [y0,z0]
uu,ww = odeint(pphi,yy0,tt, tfirst=True).T
figure(figsize=(18,7))
affichage(tt,yy,zz,uu,ww,"Odeint")
Exercice : étude d'un schéma RK¶
Soit le schéma de Runge-Kutta dont la matrice de Butcher est $$ \begin{array}{c|cc} \dfrac{1}{3} & \dfrac{1}{3} & 0\\ 1 & 1 & 0\\ \hline & \dfrac{3}{4} & \dfrac{1}{4} \end{array} $$
Q5 [1 point] Écrire le schéma sous la forme d'une suite définie par récurrence.
Le schéma associé à cette matrice est semi-implicite ($K_1$ dépend de lui même) et permet de calculer $u_{n+1}$ à partir de $u_n$ par la formule de récurrence $$\begin{cases} u_0 = y_0 \\ K_1 = \varphi\left(t_n+\frac{h}{3},u_n+\frac{h}{3}K_1\right)\\ K_2 = \varphi\left(t_{n+1},u_n+hK_1\right)\\ u_{n+1} = u_n + \frac{h}{4}\left(3K_1+K_2\right) & n=0,1,\dots N-1 \end{cases}$$
Q6 [2 points] Étudier théoriquement l'ordre du schéma.
Soit $\omega$ l'ordre de la méthode.
C'est un schéma semi-implicite à $2$ étages ($s=2$) donc $\omega\le2s=4$
Si $\begin{cases} \displaystyle\sum_{j=1}^s b_{i}=1& \\ \displaystyle c_i=\sum_{j=1}^s a_{ij}&i=1,\dots,s \end{cases}$ alors $\omega\ge1$
Si $\displaystyle\sum_{j=1}^s b_j c_j=\frac{1}{2}$ alors $\omega\ge2$
- Si, de plus, $\begin{cases} \displaystyle\sum_{j=1}^s b_j c_j^2=\frac{1}{3}\\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i a_{ij} c_j=\frac{1}{6} \end{cases}$ alors $\omega\ge3$
- Si, de plus,
$\begin{cases} \displaystyle\sum_{j=1}^s b_j c_j^3=\frac{1}{4}& \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i c_i a_{ij} c_j=\frac{1}{8} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i a_{ij} c_j^2=\frac{1}{12} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s\sum_{k=1}^s b_i a_{ij}a_{jk} c_k=\frac{1}{24} \end{cases}$ alors $\omega\ge4$
Calculons donc toutes ces sommes:
import sympy as sym
sym.init_printing()
from IPython.display import display, Math
c=[sym.Rational(1,3),1]
b=[sym.Rational(3,4),sym.Rational(1,4)]
A=[[sym.Rational(1,3),0],[1,0]]
s=len(c)
print('Ordre 1 (=consistance): si la première somme vaut 1 et les sommes suivantes 0')
display(Math(r'\sum_{j=1}^s b_j='+sym.latex(sum(b)) ))
for i in range(s):
display(Math(r'i={}'+str(i)+'\quad \sum_{j=1}^s a_{ij}-c_i='+sym.latex(sum(A[i])-c[i]) ))
print('Ordre 2 si la somme suivante est égale à 1/2')
display(Math(r'\sum_{j=1}^s b_j c_j='+sym.latex(sum([b[i]*c[i] for i in range(s)])) ))
print('Ordre 3 si les sommes suivantes sont respectivement égales à 1/3 et 1/6')
display(Math(r'\sum_{j=1}^s b_j c_j^2='+sym.latex(sum([b[i]*c[i]**2 for i in range(s)]).simplify()) ))
display(Math(r'\sum_{i,j=1}^s b_i a_{ij} c_j='+sym.latex(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])) ))
print('Ordre 4 si les sommes suivantes sont respectivement égales à 1/4, 1/8, 1/12 et 1/24')
display(Math(r'\sum_{j=1}^s b_j c_j^3='+sym.latex(sum([b[i]*c[i]**3 for i in range(s)]).simplify()) ))
display(Math(r'\sum_{i,j=1}^s b_i c_i a_{ij} c_j='+sym.latex(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])) ))
display(Math(r'\sum_{i,j=1}^s b_i a_{ij} c_j^2='+sym.latex(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)])) ))
display(Math(r'\sum_{i,j,k=1}^s b_i a_{ij}a_{jk} c_k='+sym.latex(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)])) ))
D'après ces résultats, le schéma est d'ordre $3$.
# import sympy as sym
# sym.init_printing()
# from IPython.display import display, Math
# c=[sym.Rational(1,3),1]
# b=[sym.Rational(3,4),sym.Rational(1,4)]
# A=[[sym.Rational(1,3),0],[1,0]]
# s=len(c)
# B = lambda p : sum([sum([b[i]*c[i]**(q-1) for i in range(s)]) != 1/q for q in range(1, p+1)])==0
# C = lambda e : sum([sum([A[i][j]*c[j]**(q-1) for j in range(s)]) != c[i]**q/q for q in range(1, e+1) for i in range(s)])==0
# D = lambda z : sum([sum([b[i]*c[i]**(q-1)*A[i][j] for i in range(s)]) != (b[j]/q) * (1-c[j]**q) for q in range(1, z+1) for j in range(s)])==0
# omega=0
# for z in range(1, 2*s):
# for e in range(1, 2*s):
# if C(e) and D(z):
# p=1
# while B(p) and p<=min(z+e+1,2*e+2,2*s):
# p+=1
# omega=max(omega,p-1)
# display(Math(f"\omega = {omega}"))
Q7 [2 points] Étudier théoriquement la A-stabilité du schéma.
Soit $\beta>0$ un nombre réel positif et considérons le problème de Cauchy $$\begin{cases} y'(t)=-\beta y(t), &\text{pour }t>0,\\ y(0)=1. \end{cases}$$ Sa solution est $y(t)=e^{-\beta t}$ donc $$\lim_{t\to+\infty}y(t)=0.$$
Le schéma appliqué à ce problème de Cauchy s'écrit $$\begin{cases} u_0 = y_0 \\ K_1 = -\beta\left(u_n+\frac{h}{3}K_1\right)\\ K_2 = -\beta\left(u_n+hK_1\right)\\ u_{n+1} = u_n + \frac{h}{4}\left(3K_1+K_2\right) & n=0,1,\dots N-1 \end{cases}$$
L'équation $K_1 = -\beta\left(u_n+\frac{h}{3}K_1\right)$ donne $K_1=\frac{-3\beta}{\beta h+3}u_n$ ainsi $K_2=\beta\frac{2\beta h-3}{\beta h+3}u_n$ et enfin $$ u_{n+1} = u_n + \frac{h}{4}\left( \frac{-9\beta}{\beta h+3}+\beta\frac{2\beta h-3}{\beta h+3}\right)u_ n= \frac{(\beta h)^2-4\beta h+6}{2(\beta h+3)}u_n $$
Vérifions ce calcul:
import sympy as sym
sym.init_printing()
from IPython.display import display, Math
# pour ne pas effacer l'affectation de "h", ici on vas l'appeler "dt" mais on affiche "h"
dt = sym.Symbol('h')
u_np1 = sym.Symbol('u_{n+1}')
beta = sym.Symbol(r'\beta')
sym.var('u_n,x')
K1 = sym.solve(-x-beta*(u_n+dt/3*x),x)[0]
display(Math('K_1='+sym.latex(K1)))
K2 = -beta*(u_n+dt*K1).factor()
display(Math('K_2='+sym.latex(K2)))
u_np1 = (u_n+dt/4*(3*K1+K2)).factor()
display(Math('u_{n+1}='+sym.latex(u_np1)))
On note $x=\beta h$ et on étudie la fonction $q\colon \mathbb{R}^+\to\mathbb{R}$ définie par $q(x)=\dfrac{x^2-4x+6}{2(x+3)}=\dfrac{1}{2}\left(x-7+\dfrac{27}{x+3}\right)$. Le schéma est A-stable ssi $|q(x)|<1$. D'après les calculs ci-dessous on conclut que le schéma est A-stable ssi $h<\dfrac{6}{\beta}$.
import sympy as sym
sym.init_printing()
from IPython.display import display, Math
sym.var('x',nonnegative=True)
q=((x**2-4*x+6)/(2*x+6)).apart()
display(Math('q(x)='+sym.latex(q)))
display(Math('q(0)='+sym.latex(q.subs(x,0))))
lim=sym.Limit(q,x,sym.oo)
display(Math(sym.latex(lim)+'='+sym.latex(lim.doit())))
print('On cherche les points stationnaires dans R^+')
dq=(q.diff(x)).factor()
display(Math("q'(x)="+sym.latex(dq)))
print("Le signe de q' est le signe du numérateur qui est une parabole. On cherche où elle s'annulle pour x>0")
sol=sym.solve(dq)
display(Math("q'(x)=0 \iff x\in"+sym.latex(sol)))
minimum=sol[0]
display(Math("x="+sym.latex(minimum)+"\\text{ est un minimum et}"))
qmin=q.subs(x,sol[0])
display(Math("q("+sym.latex(minimum)+")="+sym.latex(qmin)+"\\approx"+sym.latex(sym.N(qmin))))
print("Conclusion: q(x)>0 pour tout x>0. Vérifions ce calcul:")
print("q(x)>0")
display(sym.solve(q>0))
print("On resout maintenant q(x)=1:")
display(sym.solve(q-1))
print('On trouve q(0)=q(6)=1')
print("Conclusion: q(x)<1 pour x<6. Vérifions ce calcul:")
print("q(x)<1")
display(sym.solve(q<1))
sym.plot(q,1,-1,xlim=[-1,15],ylim=[-1.5,1.5]);
Q8 [2 points] Implémenter le schéma et approcher la solution du problème de Cauchy suivant avec $N+1$ points et $N=8$. $$ \begin{cases} y'(t)=-6y(t), &t\in[0;1]\\ y(0)=1 \end{cases} $$
Dans chaque point $t_i$, il faut approcher $(K_1)_i$ en resolvant une équation. Si on utilise la fonction fsolve du module scipy.optimize, il faut initialiser fsolve avec une approximation de $(K_1)_i$. On choisira donc $\varphi(t_i,u_i)$:
%reset -f
from matplotlib.pylab import *
from scipy.optimize import fsolve
def RK(phi,tt):
uu = [y0]
for i in range(len(tt)-1):
k1 = fsolve( lambda x : -x+phi( tt[i]+h/3, uu[i]+h/3*x ) , phi(tt[i],uu[i]) )[0]
k2 = phi( tt[i+1], uu[i]+h*k1 )
uu.append( uu[i] + h/4*(3*k1+k2 ))
return uu
t0, y0, tfinal = 0, 1, 3
sol_exacte = lambda t : exp(-6*t)
phi = lambda t,y : -6*y
print('A-stable ssi h <',6/6)
figure(figsize=(10,5))
N=8
tt = linspace(t0,tfinal,N+1)
h = tt[1]-tt[0]
yy = [sol_exacte(t) for t in tt]
uu = RK(phi,tt)
plot(tt,yy,'b-',label=("Exacte"))
plot(tt,uu,'ro',label=("RK"))
title(r' $N$={}, $h$={}'.format(N,h))
xlabel('$t$')
ylabel('$u$')
legend()
grid(True)
Exercice : étude d'un schéma Predictor-Corrector multipas¶
On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Soit les méthodes multipas \begin{align*} u_{n+1}&=u_{n-3}+\frac{4h}{3}\left(2\varphi_{n}-\varphi_{n-1}+2\varphi_{n-2}\right)&\text{[P]} \\ u_{n+1}&=u_{n-1}+\frac{h}{3}\left(\varphi_{n+1}+4\varphi_{n}+\varphi_{n-1}\right)&\text{[C]} \end{align*}
Q9 [3 point] Étudier consistance, ordre et zéro-stabilité de la méthode P.
P est une méthode à $q=4$ pas explicite :
- $p=3$
- $b_{-1}=0$
- $a_0=a_1=a_2=0$ et $a_3=1$
- $b_0=\frac{8}{3}$, $b_1=\frac{-4}{3}$, $b_2=\frac{8}{3}$ et $b_3=0$
- La méthode est consistante d'ordre $\omega=4$ car $$ \begin{cases} \xi(0)=\displaystyle\sum_{j=0}^p a_j=1, \\ \xi(1)=\displaystyle\sum_{j=0}^p (-j)^{1}a_j+1\sum_{j=-1}^p (-j)^{0}b_j=1 \\ \xi(2)=\displaystyle\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j=1 \\ \xi(3)=\displaystyle\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j=1 \\ \xi(4)=\displaystyle\sum_{j=0}^p (-j)^{4}a_j+4\sum_{j=-1}^p (-j)^{3}b_j=1 \\ \xi(5)=\displaystyle\sum_{j=0}^p (-j)^{5}a_j+5\sum_{j=-1}^p (-j)^{4}b_j\neq1 \end{cases} $$ Attention $(-0)^0=1$ dans les formules précédentes.
Vérifions ces calculs ci-dessous:
import sympy as sym
sym.init_printing()
from IPython.display import display, Math
p=3
bm1=0
a0=0
a1=0
a2=0
a3=1
b0=sym.Rational(8,3)
b1=sym.Rational(-4,3)
b2=sym.Rational(8,3)
b3=0
aa=[a0,a1,a2,a3]
bb=[b0,b1,b2,b3]
cond=[]
cond.append(sum(aa))
cond.append(sum([-j*aa[j]+bb[j] for j in range(len(aa))])+bm1)
if cond==[1,1]:
display(Math("\omega\ge1"))
for n in range(2,9):
cond.append(sum( [ (-j)**n*aa[j]+n*(-j)**(n-1)*bb[j] for j in range(len(aa)) ])+n*bm1)
if cond[n]==1:
display(Math("\omega\ge"+str(n)))
else:
display(Math("\omega<"+str(n)))
break
#print(cond)
Le premier polynôme caractéristique est $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^4-a_0r^3-a_1r^2-a_2r-a_3=r^4-1 $$ dont les racines sont $$ r_0=1,\ r_1=-1,\ r_2=i,\ r_3=-i $$ La méthode est zéro-stable ssi $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1, \end{cases} $$ ce qui est bien vérifié.
La méthode est convergente car consistante et zéro-stable.
Q10 [3 point] Étudier consistance, ordre et zéro-stabilité de la méthode C.
C est une méthode à $q=2$ pas implicite :
- $p=1$
- $b_{-1}=\frac{1}{3}$
- $a_0=0$ et $a_1=1$
- $b_0=\frac{4}{3}$ et $b_1=\frac{1}{3}$
- La méthode est consistante d'ordre 4 car $$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=1, \\ \displaystyle\sum_{j=0}^p -ja_j+\sum_{j=-1}^p b_j=1 \\ \displaystyle\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j=1 \\ \displaystyle\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j=1 \\ \displaystyle\sum_{j=0}^p (-j)^{4}a_j+4\sum_{j=-1}^p (-j)^{3}b_j=1 \\ \displaystyle\sum_{j=0}^p (-j)^{5}a_j+5\sum_{j=-1}^p (-j)^{4}b_j\neq1 \end{cases} $$ Vérifions ces calculs ci-dessous:
import sympy as sym
sym.init_printing()
from IPython.display import display, Math
p=1
bm1=sym.Rational(1,3)
a0=0
a1=1
b0=sym.Rational(4,3)
b1=sym.Rational(1,3)
aa=[a0,a1]
bb=[b0,b1]
cond=[]
cond.append(sum(aa))
cond.append(sum([-j*aa[j]+bb[j] for j in range(len(aa))])+bm1)
if cond==[1,1]:
display(Math("\omega\ge1"))
for n in range(2,9):
cond.append(sum( [ (-j)**n*aa[j]+n*(-j)**(n-1)*bb[j] for j in range(len(aa)) ])+n*bm1)
if cond[n]==1:
display(Math("\omega\ge"+str(n)))
else:
display(Math("\omega<"+str(n)))
break
#print(cond)
Le premier polynôme caractéristique est $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^2-a_0r-a_1=r^2-1 $$ dont les racines sont $$ r_0=1,\ r_1=-1 $$ La méthode est zéro-stable ssi $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1, \end{cases} $$ ce qui est bien vérifié.
La méthode est convergente car consistante et zéro-stable.
Q11 [2 point] Calculer empiriquement l'ordre de convergence des méthodes P et C ainsi que de la méthode Predictor-Corrector associée aux deux méthodes P et C sur le problème de Cauchy $$ \begin{cases} y'(t) = 1+\big(t-y(t)\big)^2, &\forall t \in I=[2,3],\\ y(2) = 1, \end{cases} $$ dont la solution exacte est $y(t)=t+\dfrac{1}{1-t}$.
Il est d'ordre $4$ car le prédictor est d'ordre $4$ et le corrector d'ordre $4$.
On définit la solution exacte pour estimer les erreurs et on calcule la solution approchée pour différentes valeurs de $N$ pour les trois schémas:
%reset -f
%matplotlib inline
from matplotlib.pylab import *
# variables globales
t0 = 2
tfinal = 3
y0 = 1
phi = lambda t,y : 1+(t-y)**2
sol_exacte = lambda t : t+1/(1-t)
from scipy.optimize import fsolve
def multipas_P(phi, tt):
h = tt[1] - tt[0]
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,len(tt) - 1):
uu.append( uu[i-3]+4*h/3* ( 2*phi(tt[i],uu[i])-phi(tt[i-1],uu[i-1])+2*phi(tt[i-2],uu[i-2]) ) )
return uu
def multipas_C(phi, tt):
h = tt[1] - tt[0]
uu = [y0]
uu.append(sol_exacte(tt[1]))
for i in range(1,len(tt) - 1):
temp = fsolve ( lambda x : -x+uu[i-1]+h/3* ( phi(tt[i+1],x) + 4*phi(tt[i],uu[i])+phi(tt[i-1],uu[i-1]) ), uu[i])
uu.append(temp)
return uu
def multipas_PC(phi, tt):
h = tt[1] - tt[0]
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,len(tt) - 1):
u_tilde = uu[i-3]+4*h/3* ( 2*phi(tt[i],uu[i])-phi(tt[i-1],uu[i-1])+2*phi(tt[i-2],uu[i-2]) )
uu.append( uu[i-1]+h/3* ( phi(tt[i+1],u_tilde) + 4*phi(tt[i],uu[i])+phi(tt[i-1],uu[i-1]) ) )
return uu
H = []
err_p = []
err_c = []
err_pc = []
N = 10
for k in range(7):
N+=20
tt = linspace(t0, tfinal, N + 1)
h = tt[1] - tt[0]
yy = [sol_exacte(t) for t in tt]
uu_p = multipas_P(phi, tt)
uu_c = multipas_C(phi, tt)
uu_pc = multipas_PC(phi, tt)
H.append(h)
err_p.append(max([abs(uu_p[i] - yy[i]) for i in range(len(yy))]))
err_c.append(max([abs(uu_c[i] - yy[i]) for i in range(len(yy))]))
err_pc.append(max([abs(uu_pc[i] - yy[i]) for i in range(len(yy))]))
print ('Multipas P ordre=%1.2f' %(polyfit(log(H),log(err_p), 1)[0]))
print ('Multipas C ordre=%1.2f' %(polyfit(log(H),log(err_c), 1)[0]))
print ('Multipas PC ordre=%1.2f' %(polyfit(log(H),log(err_pc), 1)[0]))
figure(figsize=(8,5))
plot(log(H), log(err_p), 'r-o', label='Multipas P')
plot(log(H), log(err_c), 'b-o', label='Multipas C')
plot(log(H), log(err_pc), 'c-o', label='Multipas PC')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
legend(bbox_to_anchor=(1.04, 1), loc='upper left')
grid(True);
On a bien
- ordre du corrector $\omega_{[C]}=4$
- ordre du predictor $\omega_{[P]}=4$
et donc
ordre du predictor-corrector $\omega_{[PC]}=\min\{\omega_{[C]},\omega_{[P]}+1\}=4$