In [1]:
from IPython.core.display import HTML, display
css_file = './custom.css'
HTML(open(css_file, "r").read())
Out[1]:

2021 CR : le sujet comporte trois exercices indépendents.

Exercice 1 : étude d'un schéma

On étudie le schéma de résolution du problème de Cauchy $y'(t)=\varphi(t,y(t))$ et $y(t_0)=y_0$ défini par la formule de récurrence

$$ \begin{cases} u_0 = y_0 \\ K_1 = \varphi\left(t_n,u_n\right)\\ K_2 = \varphi\left(t_{n}+\frac{1}{2}h,u_n+\frac{1}{2}hK_1\right)\\ K_3 = \varphi\left(t_{n}+h,u_n+hK_1\right)\\ u_{n+1} = u_n + h\left(aK_1+bK_2+cK_3\right) & n=0,1,\dots N-1 \end{cases} $$

où $a,b,c \in [0;1]$.

Q1 [3 points]
Pour quelles valeurs du triplet $(a,b,c)$ retrouve-t-on

  • la méthode d'Euler explicite (EE) ?
  • la méthode d'Euler modifiée (EM) ?
  • la méthode de Heun (H) ?
  • Si $(a,b,c)=(1,0,0)$ on retrouve la méthode d'Euler explicite (EE): $$ u_{n+1} = u_n + h K_1 = u_n + h \varphi\left(t_n,u_n\right) $$
  • Si $(a,b,c)=(0,1,0)$ on retrouve la méthode d'Euler modifiée (EM): $$ u_{n+1} = u_n + h K_2 = u_n + h \varphi\left(t_{n}+\frac{1}{2}h,u_n+\frac{1}{2}h\varphi\left(t_n,u_n\right)\right) $$
  • Si $(a,b,c)=(\frac{1}{2},0,\frac{1}{2})$ on retrouve la méthode de Heun (H): $$ u_{n+1} = u_n + \frac{h}{2} (K_1+K_3) = u_n + \frac{h}{2} \left(\varphi\left(t_n,u_n\right)+\varphi\left(t_{n}+h,u_n+h\varphi\left(t_n,u_n\right)\right)\right) $$

Q2 [2 points]
Le schéma peut être vu comme une méthode de type Runge-Kutta. Écrire la matrice de Butcher associée.
Le schéma est explicite, semi-implicite ou implicite ?

$$ \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \dfrac{1}{2} & \dfrac{1}{2} & 0 & 0 \\ 1 & 1 & 0 & 0 \\ \hline & a & b & c \end{array} $$

Le schéma est explicite (la matrice $\mathbb{A}$ est triangulaire inférieure à diagonale nulle).

Q3 [1 points]
Sans faire de calculs, indiquer quel est l'ordre maximale du schéma.

Soit $\omega$ l'ordre de la méthode et $s=3$ le nombre d'étages. Étant un schéma explicite, $\omega\le s=3$.

Q4 [3 points]
Étudier théoriquement l'ordre du schéma en fonction de $a,b,c$.

  • 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 $\omega\ge1$ et $\displaystyle\sum_{j=1}^s b_j c_j=\frac{1}{2}$ alors $\omega\ge2$

  • Si $\omega\ge2$ et $\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$

Il est donc

  • au moins d'ordre 1 ssi $a=1-b-c$
  • il est d'ordre 2 ssi, de plus, $\frac{b}{2}+c=\frac{1}{2}$ donc ssi $\begin{cases}a=c\\b=1-2c\end{cases}$
  • il ne peut pas être d'ordre 3 car la deuxième condition n'est jamais satisfaite.

On peut vérifier nos calculs comme ci-dessous:

In [2]:
import sympy as sym
sym.init_printing()

from IPython.display import display, Math

sym.var('a,b,c')

cc=[0,sym.S(1)/2,1]
bb=[a,b,c]
A=[[cc[0],0,0],[cc[1],0,0],[cc[2],0,0]]
s=len(cc)
print('Étude de la consistance (ordre 1)')
calc=[sum(bb)]
display(Math(r'\sum_{j=1}^s b_j='+sym.latex(calc[0]) +'\\text{ doit être egale à 1}'    ))
for i in range(s):
    calc.append(sum(A[i])-cc[i]) 
    display(Math(r'i={}'+str(i)+'\quad \sum_{j=1}^s a_{ij}-c_i='+sym.latex(calc[-1])+'\\text{ doit être egale à 0}'))

print("\nÉtude de l'Ordre 2\n")
calc=sum([bb[i]*cc[i] for i in range(s)])  
display(Math(r'\sum_{j=1}^s b_j c_j='+sym.latex(calc )+'\\text{ doit être egale à }\\frac{1}{2}')) 

print("\nÉtude de l'Ordre 3\n")   
calc3_1=sum([bb[i]*cc[i]**2 for i in range(s)])
display(Math(r'\sum_{j=1}^s b_jc_j^2='+sym.latex(calc3_1)+'\\text{ doit être egale à }\\frac{1}{3}'     ))
calc3_2=sum([bb[i]*A[i][j]*cc[j] for i in range(s) for j in range(s)])
display(Math(r'\sum_{i=1}^s\sum_{j=1}^s b_ia_{ij}c_j='+sym.latex(calc3_2)+'\\text{ doit être egale à }\\frac{1}{6}'     ))

Étude de la consistance (ordre 1)
$\displaystyle \sum_{j=1}^s b_j=a + b + c\text{ doit être egale à 1}$
$\displaystyle i={}0\quad \sum_{j=1}^s a_{ij}-c_i=0\text{ doit être egale à 0}$
$\displaystyle i={}1\quad \sum_{j=1}^s a_{ij}-c_i=0\text{ doit être egale à 0}$
$\displaystyle i={}2\quad \sum_{j=1}^s a_{ij}-c_i=0\text{ doit être egale à 0}$ 
Étude de l'Ordre 2
$\displaystyle \sum_{j=1}^s b_j c_j=\frac{b}{2} + c\text{ doit être egale à }\frac{1}{2}$ 
Étude de l'Ordre 3
$\displaystyle \sum_{j=1}^s b_jc_j^2=\frac{b}{4} + c\text{ doit être egale à }\frac{1}{3}$
$\displaystyle \sum_{i=1}^s\sum_{j=1}^s b_ia_{ij}c_j=0\text{ doit être egale à }\frac{1}{6}$

Q5 [2 points]
Étudier théoriquement la A-stabilité du schéma pour les valeurs du triplet $(a,b,c)$ qui donnent l'ordre maximal.

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 u_n\\ K_2 = -\beta \left(u_n+\frac{1}{2}hK_1\right)=-\beta \left(u_n-\frac{1}{2}\beta h u_n\right)\\ K_3 = -\beta \left(u_n+hK_1\right)=-\beta \left(u_n-\beta h u_n\right)\\ u_{n+1} = u_n + h\left(aK_1+bK_2+cK_3\right) & n=0,1,\dots N-1 \end{cases} $$ d'où $$ u_{n+1} = \left( 1-(a+b+c)\beta h + \left(\frac{b}{2}+c\right)(\beta h)^2 \right)u_n $$

L'ordre maximal (i.e. l'ordre 2) est obtenu pour $a+b+c=1$ et $b=1-2c$, on trouve alors $$u_{n+1} = \left(\dfrac{1}{2}(\beta h)^2 -(\beta h)+1\right)u_n$$ Vérifions ces calculs:

In [3]:
import sympy as sym
sym.init_printing()
from IPython.display import display, Math

# pour ne pas effacer l'affectation de "h", ici je vais l'appeler "dt"
sym.var('u_n,u_np1,beta,dt,x, a,b,c')

K1 = -beta*u_n
display(Math('K_1='+sym.latex(K1)))
K2 = -beta*(u_n+sym.S(1)/2*dt*K1).factor()
display(Math('K_2='+sym.latex(K2)))
K3 = -beta*(u_n+dt*K1).factor()
display(Math('K_3='+sym.latex(K3)))
u_np1 = (u_n+dt*(a*K1+b*K2+c*K3)).factor()
display(Math('u_{n+1}='+sym.latex(u_np1)))
print("Avec les paramètres qui donnent l'ordre 2 on obtient")
display(Math('u_{n+1}='+sym.latex(u_np1.subs(a,1-b-c).subs(b,1-2*c).simplify())))
$\displaystyle K_1=- \beta u_{n}$
$\displaystyle K_2=\frac{\beta u_{n} \left(\beta dt - 2\right)}{2}$
$\displaystyle K_3=\beta u_{n} \left(\beta dt - 1\right)$
$\displaystyle u_{n+1}=- \frac{u_{n} \left(2 a \beta dt - b \beta^{2} dt^{2} + 2 b \beta dt - 2 \beta^{2} c dt^{2} + 2 \beta c dt - 2\right)}{2}$ 
Avec les paramètres qui donnent l'ordre 2 on obtient
$\displaystyle u_{n+1}=\frac{u_{n} \left(\beta^{2} dt^{2} - 2 \beta dt + 2\right)}{2}$

On note $x=\beta h$ et on étudie la fonction $q\colon \mathbb{R}^+\to\mathbb{R}$ définie par $q(x)=\frac{1}{2}x^2-x+1$ car le schéma est A-stable ssi $|q(x)|<1$.
Il s'agit d'une parabole convexe de sommet en $x_s=1$ et $q(x_s)=\frac{1}{2}>-1$ donc $q(x)<1$ pour $0<x<2$ et $q(x)>0$ pour tout $x$. On conclut que le schéma est A-stable pour tout $h<\dfrac{2}{\beta}$.

On peut vérifier les calculs avec sympy:

In [4]:
import sympy as sym
sym.init_printing()
from IPython.display import display, Math
sym.var('x',nonnegative=True)

q=sym.Rational(1,2)*x**2-x+1
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("q est une parabole dont le sommet est")
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)))
sym.plot(q,1,-1,xlim=[-1,15],ylim=[-1.5,1.5]);
print("Conclusion: -1<q(x)<1 pour tout x<2. Vérifions ce calcul:")
print("q(x)<1")
display(sym.solve(q<1))
print("q(x)>-1")
display(sym.solve(q>-1))
print("Notons que q(x)>0 pour tout 0<x<2")
#display(sym.solve(q>0))
$\displaystyle q(x)=\frac{x^{2}}{2} - x + 1$
$\displaystyle q(0)=1$
$\displaystyle \lim_{x \to \infty}\left(\frac{x^{2}}{2} - x + 1\right)=\infty$ 
On cherche les points stationnaires dans R^+
$\displaystyle q'(x)=x - 1$ 
q est une parabole dont le sommet est
$\displaystyle q'(x)=0 \iff x\in\left[ 1\right]$
$\displaystyle x=1\text{ est un minimum et}$
$\displaystyle q(1)=\frac{1}{2}$ 
Conclusion: -1<q(x)<1 pour tout x<2. Vérifions ce calcul:
q(x)<1
$\displaystyle 0 < x \wedge x < 2$ 
q(x)>-1
$\displaystyle \text{True}$ 
Notons que q(x)>0 pour tout 0<x<2

Pour la question 6 on va fixer $(a,b,c)$. Pour cela, écrivez votre nom et prénom dans la cellule ci-dessous et vous obtiendrez un choix pour les paramètres.

In [5]:
from IPython.display import display, Math
import sympy 
sympy.init_printing()

nom="Faccanoni"
prenom="Gloria"

DEN=list(range(3,11))
idx=sum([ord(c) for c in nom+prenom])%len(DEN)
#for idx in range(len(DEN)):
my_c=sympy.S(1)/DEN[idx]
my_b=1-2*my_c
display(Math(r'a='+sympy.latex(my_c)+',\\quad b='+sympy.latex(my_c)+',\\quad  c='+sympy.latex(my_b)))
$\displaystyle a=\frac{1}{3},\quad b=\frac{1}{3},\quad c=\frac{1}{3}$

Q6 [2 points]
Implémenter le schéma avec les valeurs $a,b,c$ trouvées avec vos nom et prénom et approcher la solution du problème de Cauchy suivant avec $N+1$ points (quelles valers de $N$ peut-on choisir? justifier la réponse) $$ \begin{cases} y'(t)=-50y(t), &t\in[0;1]\\ y(0)=1 \end{cases} $$ Vérifier ensuite empiriquement l'ordre de convergence.

Pour la A-stabilité il faut que $h$ soit $<\frac{2}{50}$ donc $N>\frac{1-0}{h}=\frac{50}{2}$

In [6]:
%reset -f

from matplotlib.pylab import *
from scipy.optimize import fsolve

# def RK(tt,phi,y0): # comme suite geometrique
#     x=beta*h
#     q=x**2/2-x+1
#     uu = [y0]
#     for i in range(len(tt)-1):
#         uu.append( q*uu[i] )
#     return uu


a,b,c = 1/3, 1/3, 1/3


def RK(tt,phi,y0):
    uu = [y0]
    h = tt[1]-tt[0]
    for i in range(len(tt)-1):
        K1=phi(tt[i],uu[i])
        K2=phi(tt[i]+h/2,uu[i]+h/2*K1)
        K3=phi(tt[i+1],uu[i]+h*K1)
        uu.append( uu[i]+h*(a*K1+b*K2+c*K3) )
    return uu


t0, y0, tfinal = 0, 1, 1
beta=50

phi = lambda t,y : -beta*y
sol_exacte = lambda t : exp(-beta*t)

figure(figsize=(10,5))
N=int(50/2)+1
tt = linspace(t0,tfinal,N+1)
h = tt[1]-tt[0]

if h>=2/beta:
    print('Attention: h > limite de stabilité')
    
yy = [sol_exacte(t) for t in tt]    
# uu = RK(tt)
uu = RK(tt,phi,y0)
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)


H = []
err = []
N = 50
for k in range(7):
    N+=100
    tt = linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = [sol_exacte(t) for t in tt]
    uu = RK(tt,phi,y0)
    H.append(h)
    err.append(max([abs(uu[i] - yy[i]) for i in range(len(yy))]))
    
omega=polyfit(log(H),log(err), 1)[0]
figure(figsize=(8,5))
plot(log(H), log(err), 'c-o')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
title(r'Ordre $\omega$={:1.2f}'.format(omega))
grid(True);

Exercice 2 : étude d'un schéma multipas

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$.

Soit la méthode multipas

$$ u_{n+1}=u_{n-1} +h\left(b_{-1} \varphi_{n+1} +b_0\varphi_{n}+b_{1} \varphi_{n-1}\right). $$

Q1 [2 points]
Pour quelles valeurs des paramètres $b_0, b_1, b_2$ est-elle zéro-stable?

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=(r-1)(r+1) $$ dont les racines sont $$ r_0=1,\quad r_1=-1. $$ La méthode est donc zéro-stable pour tout choix des paramètres car $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1. \end{cases} $$

Q2 [3 points]
Pour quelles valeurs des paramètres la méthode est convergente d'ordre 4?

  • Une méthode consistante d'ordre $\omega\ge1$ et zéro-stable est convergente d'ordre $\omega$.
  • On a $p=1$: c'est une méthode à $q=p+1=2$ pas avec $a_0=0$ et $a_1=1$.
  • La méthode est implicite si $b_{-1}\neq0$.

La première barrière de Dahlquist affirme qu'un schéma implicite à $q=2$ étapes consistante et zéro-stable peut être au mieux d'ordre $\omega=q+2=4$.

Pour que la méthode soit consistante d'ordre 4 il faut que $$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=a_0+a_1=1, \\ \displaystyle\sum_{j=0}^p (-j)a_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 \end{cases} $$ Dans notre cas cela reviens à imposer $$ \begin{cases} 0+1=1 \\ {}-1+\big(b_{-1}+b_0+b_1\big)=1 \\ 1-2\big(-b_{-1}+b_1\big)=1 \\ {}-1+3\big(b_{-1}+b_1\big)=1 \\ 1-4\big(-b_{-1}+b_1\big)=1 \end{cases} $$ On cherche donc $b_{-1}$, $b_0$ et $b_1$ tels que $$ \begin{pmatrix} 1&1&1\\ 2&0&-2\\ 3&0&3\\ 4&0&-4 \end{pmatrix} \begin{pmatrix} b_{-1}\\b_0\\b_1 \end{pmatrix} {}= \begin{pmatrix} 2\\0\\2\\0 \end{pmatrix} $$ Il s'agit d'un système surdeterminé qui admet l'unique solution $$ b_{-1}=b_1=\frac{1}{3},\qquad b_0=\frac{4}{3} $$

In [7]:
import sympy as sym
sym.init_printing()

sym.var('b_m1,b_0,b_1')

eq2 = sym.Eq(b_m1 + b_0 + b_1, 2)
eq3 = sym.Eq(2 * b_m1 - 2 * b_1, 0)
eq4 = sym.Eq(3 * b_m1 + 3 * b_1, 2)
eq5 = sym.Eq(4 * b_m1 - 4 * b_1, 0)
sym.solve([eq2, eq3, eq4, eq5])
Out[7]:
$\displaystyle \left\{ b_{0} : \frac{4}{3}, \ b_{1} : \frac{1}{3}, \ b_{m1} : \frac{1}{3}\right\}$

Q3 [3 points]
Vérifier empiriquement l'ordre de convergence sur le problème de Cauchy $$ \begin{cases} y'(t) = y(t)\ln(1+y(t)), &\forall t \in I=[0,1],\\ y(0) = 2, \end{cases} $$ après en avoir calculé la solution exacte.

  • On définit la solution exacte (calculée en utilisant le module sympy) pour estimer les erreurs:
In [8]:
%reset -f
%matplotlib inline

from matplotlib.pylab import *

# variables globales
t0     = 0
tfinal = 1
y0     = 2

phi = lambda t,y:y*log(1/y)
In [9]:
import sympy as sym
sym.init_printing()

t  = sym.Symbol('t')
y  = sym.Function('y')
edo= sym.Eq( sym.diff(y(t),t) , y(t)*sym.log(1/y(t)) )
display(edo)
solgen = sym.dsolve(edo)
display(solgen)

consts = sym.solve( sym.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
display(consts)
solpar=solgen.subs(consts).simplify()
display(solpar)

sol_exacte = sym.lambdify(t,solpar.rhs,'numpy')
$\displaystyle \frac{d}{d t} y{\left(t \right)} = y{\left(t \right)} \log{\left(\frac{1}{y{\left(t \right)}} \right)}$
$\displaystyle y{\left(t \right)} = e^{C_{1} e^{- t}}$
$\displaystyle \left\{ C_{1} : \log{\left(2 \right)}\right\}$
$\displaystyle y{\left(t \right)} = e^{e^{- t} \log{\left(2 \right)}}$
  • On calcule la solution approchée pour différentes valeurs de $N$:
In [10]:
from scipy.optimize import fsolve

def multipas(phi, tt, y0):
    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

H = []
err_mp = []
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_mp = multipas(phi, tt, y0)
    H.append(h)
    err_mp.append(max([abs(uu_mp[i] - yy[i]) for i in range(len(yy))]))

print ('Multipas ordre=%1.2f' %(polyfit(log(H),log(err_mp), 1)[0]))

figure(figsize=(8,5))
plot(log(H), log(err_mp), 'r-o')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
grid(True)
show()

Multipas ordre=4.04
In [ ]: