None 2020-DM-corrige
In [38]:
from IPython.core.display import HTML, display
css_file = './custom.css'
HTML(open(css_file, "r").read())
Out[38]:

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 avec sympy. 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) $$

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

In [40]:
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');
$\displaystyle \frac{d}{d t} y{\left(t \right)} = z{\left(t \right)}$
$\displaystyle \frac{d}{d t} z{\left(t \right)} = - 6 y{\left(t \right)} - 5 z{\left(t \right)}$
$\displaystyle \left[ y{\left(t \right)} = - \frac{C_{1} e^{- 3 t}}{3} - \frac{C_{2} e^{- 2 t}}{2}, \ z{\left(t \right)} = C_{1} e^{- 3 t} + C_{2} e^{- 2 t}\right]$
$\displaystyle \left\{ C_{1} : 6, \ C_{2} : -6\right\}$
$\displaystyle y{\left(t \right)} = 3 e^{- 2 t} - 2 e^{- 3 t}$
$\displaystyle z{\left(t \right)} = - 6 e^{- 2 t} + 6 e^{- 3 t}$

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

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

In [42]:
# 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 fonction odeint du module scipy. Afficher $t\mapsto x(t)$, $t\mapsto x'(t)$ et $x\mapsto x'(x)$ en comparant solution exacte et solution approchée.

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

In [18]:
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)]))     ))
Ordre 1 (=consistance): si la première somme vaut 1 et les sommes suivantes 0
$\displaystyle \sum_{j=1}^s b_j=1$
$\displaystyle i={}0\quad \sum_{j=1}^s a_{ij}-c_i=0$
$\displaystyle i={}1\quad \sum_{j=1}^s a_{ij}-c_i=\mathtt{\text{0}}$
Ordre 2 si la somme suivante est égale à 1/2
$\displaystyle \sum_{j=1}^s b_j c_j=\frac{1}{2}$
Ordre 3 si les sommes suivantes sont respectivement égales à 1/3 et 1/6
$\displaystyle \sum_{j=1}^s b_j c_j^2=\frac{1}{3}$
$\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=\frac{1}{6}$
Ordre 4 si les sommes suivantes sont respectivement égales à 1/4, 1/8, 1/12 et 1/24
$\displaystyle \sum_{j=1}^s b_j c_j^3=\frac{5}{18}$
$\displaystyle \sum_{i,j=1}^s b_i c_i a_{ij} c_j=\frac{1}{9}$
$\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j^2=\frac{1}{18}$
$\displaystyle \sum_{i,j,k=1}^s b_i a_{ij}a_{jk} c_k=\frac{1}{18}$

D'après ces résultats, le schéma est d'ordre $3$.

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

In [20]:
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)))
$\displaystyle K_1=- \frac{3 \beta u_{n}}{\beta h + 3}$
$\displaystyle K_2=\frac{\beta u_{n} \left(2 \beta h - 3\right)}{\beta h + 3}$
$\displaystyle u_{n+1}=\frac{u_{n} \left(\beta^{2} h^{2} - 4 \beta h + 6\right)}{2 \left(\beta h + 3\right)}$

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

In [21]:
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]);
$\displaystyle q(x)=\frac{x}{2} - \frac{7}{2} + \frac{27}{2 \left(x + 3\right)}$
$\displaystyle q(0)=1$
$\displaystyle \lim_{x \to \infty}\left(\frac{x}{2} - \frac{7}{2} + \frac{27}{2 \left(x + 3\right)}\right)=\infty$
On cherche les points stationnaires dans R^+
$\displaystyle q'(x)=\frac{x^{2} + 6 x - 18}{2 \left(x + 3\right)^{2}}$
Le signe de q' est le signe du numérateur qui est une parabole. On cherche où elle s'annulle pour x>0
$\displaystyle q'(x)=0 \iff x\in\left[ -3 + 3 \sqrt{3}\right]$
$\displaystyle x=-3 + 3 \sqrt{3}\text{ est un minimum et}$
$\displaystyle q(-3 + 3 \sqrt{3})=-5 + 3 \sqrt{3}\approx0.196152422706632$
Conclusion: q(x)>0 pour tout x>0. Vérifions ce calcul:
q(x)>0
$\displaystyle \text{True}$
On resout maintenant q(x)=1:
$\displaystyle \left[ 0, \ 6\right]$
On trouve q(0)=q(6)=1
Conclusion: q(x)<1 pour x<6. Vérifions ce calcul:
q(x)<1
$\displaystyle 0 < x \wedge x < 6$

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)$:

In [22]:
%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)
A-stable ssi h < 1.0

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:

In [23]:
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)            
$\displaystyle \omega\ge1$
$\displaystyle \omega\ge2$
$\displaystyle \omega\ge3$
$\displaystyle \omega\ge4$
$\displaystyle \omega<5$
  • 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:
In [24]:
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)            
$\displaystyle \omega\ge1$
$\displaystyle \omega\ge2$
$\displaystyle \omega\ge3$
$\displaystyle \omega\ge4$
$\displaystyle \omega<5$
  • 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:

In [25]:
%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);
Multipas P ordre=4.26
Multipas C ordre=4.04
Multipas PC ordre=4.12

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$