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

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

Exercice 1 : étude d'un schéma RK

Soit le schéma de Runge-Kutta dont la matrice de Butcher est $$ \begin{array}{c|cc} 0 & 0 & 0\\ \dfrac{\eta}{4} & \dfrac{\eta}{4} & 0\\ \hline & \dfrac{3-\mu}{3} & \dfrac{\mu}{3} \end{array} $$

Q1
Écrire le schéma sous la forme d'une suite définie par récurrence.
Le schéma est explicite, semi-implicite ou implicite ?

Le schéma associé à cette matrice est explicite (car la matrice $\mathbb{A}$ est triangulaire inférieure à diagonale nulle).
On peut 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,u_n\right)\\ K_2 = \varphi\left(t_{n}+\dfrac{\eta}{4}h,u_n+\dfrac{\eta}{4}hK_1\right)\\ u_{n+1} = u_n + \dfrac{h}{3}\left((3-\mu) K_1+\mu K_2\right) & n=0,1,\dots N-1 \end{cases}$$

Q2
Sans faire de calculs, indiquer quel est l'ordre maximale du schéma.

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

Q3
Étudier théoriquement l'ordre du schéma en fonction de $\eta$ et $\mu$.

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

D'après les calculs ci-dessous il est toujours au moins d'ordre 1 et il est d'ordre 2 ssi $\eta\mu=6$.

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

from IPython.display import display, Math

sym.var('eta,mu')

c=[0,sym.S(eta)/4]
b=[sym.S(3-mu)/3,sym.S(mu)/3]
A=[[0,0],[sym.S(eta)/4,0]]
s=len(c)

print('Étude de la consistance (ordre 1)\nOn a')
calc=[sum(b)]
display(Math(r'\sum_{j=1}^s b_j='+sym.latex(calc[0])     ))
for i in range(s):
    calc.append(sum(A[i])-c[i]) 
    display(Math(r'i={}'+str(i)+'\quad \sum_{j=1}^s a_{ij}-c_i='+sym.latex(calc[-1])))

if calc[0]!=1 or sum([calc[i+1]==0 for i in range(s)])<s:
    print(f"donc le schéma n'est pas consistant\n")
    
else:
    print(f"alors le schéma est consistant (donc au moins d'ordre 1)\n")
    print("Étude de l'Ordre 2\non a")
    calc=sum([b[i]*c[i] for i in range(s)])  
    display(Math(r'\sum_{j=1}^s b_j c_j='+sym.latex(calc )))
    if calc==sym.S(1)/2:
        print(f"donc le schéma est d'ordre 2\n")
    else:
        print(f"alors le schéma est d'ordre 2 ssi\n")
        display(sym.solve(2*calc-1)[0])

Étude de la consistance (ordre 1)
On a
$\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=0$ 
alors le schéma est consistant (donc au moins d'ordre 1)

Étude de l'Ordre 2
on a
$\displaystyle \sum_{j=1}^s b_j c_j=\frac{\eta \mu}{12}$ 
alors le schéma est d'ordre 2 ssi
$\displaystyle \left\{ \eta : \frac{6}{\mu}\right\}$

Q4
Étudier théoriquement la A-stabilité du schéma pour les valeurs de $\eta$ et $\mu$ 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+\dfrac{\eta}{4}hK_1\right)\\ u_{n+1} = u_n + \dfrac{h}{3}\left((3-\mu)K_1+\mu K_2\right) & n=0,1,\dots N-1 \end{cases}$$

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, eta,mu')

K1 = -beta*u_n
display(Math('K_1='+sym.latex(K1)))
K2 = -beta*(u_n+sym.S(eta)/4*dt*K1).factor()
display(Math('K_2='+sym.latex(K2)))
u_np1 = (u_n+dt*(sym.S(3-mu)/3*K1+sym.S(mu)/3*K2)).factor()
display(Math('u_{n+1}='+sym.latex(u_np1)))
$\displaystyle K_1=- \beta u_{n}$
$\displaystyle K_2=\frac{\beta u_{n} \left(\beta dt \eta - 4\right)}{4}$
$\displaystyle u_{n+1}=\frac{u_{n} \left(\beta^{2} dt^{2} \eta \mu - 12 \beta dt + 12\right)}{12}$

Donc, étant donné que pour avoir l'ordre deux on a $\eta\mu=6$, alors $$u_{n+1} = \left(\dfrac{1}{2}(\beta h)^2 -(\beta h)+1\right)u_n$$

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 5 on va fixer $\eta$ et $\mu$. Pour cela, écrire 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"

L=[1,2,3,6]
idx=sum([ord(c) for c in nom+prenom])%len(L)
display(Math(r'\eta='+sympy.latex(L[idx])))
display(Math(r'\mu='+sympy.latex(sympy.S(6)/L[idx])))
$\displaystyle \eta=1$
$\displaystyle \mu=6$

Q5
Implémenter le schéma avec les valeurs de $\eta$ et $\mu$ 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.

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


eta, mu = 1, 6


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]+eta/4*h,uu[i]+eta/4*h*K1)
        uu.append( uu[i]+h/3*((3-mu)*K1+mu*K2) )
    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=50
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 : choix d'un schéma

Soit le problème de Cauchy $$ \begin{cases} y'(t)=-10^2 \Big(y(t)-\cos(t)\Big), &t\in[0;2]\\ y(0)=0 \end{cases} $$

  1. Calculer la solution exacte.

  2. Choisir le schéma le plus adapté parmi les suivants (en justifiant la réponse)

    • EE (Euler Explicite),
    • EI (Euler Implicite),
    • EM (Euler Modifié),
    • H (Heun).
  3. Approcher ensuite la solution du problème donné avec $N=40$ avec le schéma choisi (on affichera la solution exacte et la solution approchée sur le même repère).
  • On définit la solution exacte (calculée en utilisant le module sympy):
In [7]:
%reset -f
%matplotlib inline

import sympy as sym
sym.init_printing()

t  = sym.Symbol('t')
y  = sym.Function('y')

t0     = 0
tfinal = 2
y0     = 0

# NB il faut utiliser le cos du module sympy dans l'EDO
edo= sym.Eq( sym.diff(y(t),t) , -50*(y(t)-sym.cos(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)
$\displaystyle \frac{d}{d t} y{\left(t \right)} = - 50 y{\left(t \right)} + 50 \cos{\left(t \right)}$
$\displaystyle y{\left(t \right)} = \left(C_{1} + \frac{50 e^{50 t} \sin{\left(t \right)}}{2501} + \frac{2500 e^{50 t} \cos{\left(t \right)}}{2501}\right) e^{- 50 t}$
$\displaystyle \left\{ C_{1} : - \frac{2500}{2501}\right\}$
$\displaystyle y{\left(t \right)} = \frac{50 \sin{\left(t \right)}}{2501} + \frac{2500 \cos{\left(t \right)}}{2501} - \frac{2500 e^{- 50 t}}{2501}$
  • Il s'agit d'une EDO stiff (le terme en exponentielle décroît très rapidement): on ne peut utiliser ni les schémas EE et EM ni le schéma H car la condition de stabilité est trop contraignante. Le schéma EI, en revanche, est inconditionnellement A-stable: on peut donc choisir le pas $h$ sans contrainte de stabilité.
In [8]:
from matplotlib.pylab import *

N = 40
tt = linspace(t0, tfinal, N + 1)

sol_exacte = sym.lambdify(t,solpar.rhs,'numpy')
yy = sol_exacte(tt); # numpy donc la fct est vectorisée

plot(tt,yy);
  • On calcule la solution approchée
In [9]:
from scipy.optimize import fsolve

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


phi = lambda t,y : -50*(y-cos(t))

h = tt[1] - tt[0]
print(f'h={h}')
uu = EI(phi, tt, y0)
err=(max([abs(uu[i] - yy[i]) for i in range(len(yy))]))

figure(figsize=(8,5))
plot(tt, uu, 'o',tt,yy,'r-')
xlabel('$t$')
ylabel('$y$')
grid(True)
show()

h=0.05

Exercice 3 : choix d'un schéma

Soit le problème de Cauchy $$ \begin{cases} y'(t)=\alpha (y(t)-t), &t\in[0;8]\\ y(0)=\dfrac{1}{\alpha}+\varepsilon \end{cases} $$

Pour fixer $\alpha$, écrire votre nom et prénom dans la cellule ci-dessous et vous obtiendrez un choix pour ce paramètre.

In [10]:
%reset -f
%matplotlib inline

from IPython.display import display, Math
import sympy
sympy.init_printing()

nom="Faccanoni"
prenom="Gloria"

L=[2,3,4,5]
idx=sum([ord(c) for c in nom+prenom])%len(L)
alpha=L[idx]
display(Math(r'\alpha='+sympy.latex(alpha)))

alpha=5
$\displaystyle \alpha=2$

>

  1. Calculer la solution exacte.

  2. Soit $\varepsilon=0$. Parmi les schémas suivants, choisir le plus adapté (plusieurs réponses possibles, choix à justifier):

    • EE (Euler Explicite),
    • EI (Euler Implicite),
    • EM (Euler Modifié),
    • H (Heun).
  3. Soit $\varepsilon=0$. Approcher la solution du problème donné avec un des schémas choisis (on affichera la solution exacte et la solution approchée sur le même repère).
  • On calcule la solution exacte (calculée en utilisant le module sympy) :
In [11]:
import sympy as sym
sym.init_printing()

t       = sym.Symbol('t')
epsilon = sym.Symbol('epsilon')
y       = sym.Function('y')

t0     = 0
y0     = sym.Rational(1,alpha)+epsilon

edo= sym.Eq( sym.diff(y(t),t) , alpha*(y(t)-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)


from matplotlib.pylab import *
tt=linspace(0,8,101)
for val in [1.e-7,0,-1.e-7]:
    sol_exacte = sym.lambdify(t,solpar.rhs.subs(epsilon,val),'numpy')
    plot(tt,sol_exacte(tt))
$\displaystyle \frac{d}{d t} y{\left(t \right)} = - 5 t + 5 y{\left(t \right)}$
$\displaystyle y{\left(t \right)} = \left(C_{1} + \frac{\left(5 t + 1\right) e^{- 5 t}}{5}\right) e^{5 t}$
$\displaystyle \left\{ C_{1} : \epsilon\right\}$
$\displaystyle y{\left(t \right)} = \epsilon e^{5 t} + t + \frac{1}{5}$
  • On remarque que c'est un problème de Cauchy numériquement mal posé (si $\varepsilon=0$ la fonction est linéaire, avec $\varepsilon>0$ elle croît comme une exponentielle et avec $\varepsilon<0$ elle decroît comme une exponentielle): on doit donc choisir le schéma le plus précis parmi les schémas proposés, à savoir le schéma de Heun ou le schéma d'Euler Modifié qui sont d'ordre 2.
In [12]:
sol_exacte = sym.lambdify(t,solpar.rhs.subs(epsilon,0),'numpy')

t0     = 0
tfinal = 8
y0     = 1/alpha

N = 40

phi = lambda t,y : alpha*(y-t)

def EM(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        uu.append( uu[i]+h*phi(tt[i]+h/2,uu[i]+k1*h/2) )
    return uu

def heun(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        k1 = phi( tt[i], uu[i] )
        k2 = phi( tt[i+1], uu[i] + h*k1  )
        uu.append( uu[i] + (k1+k2)*h/2 )
    return uu


tt = linspace(t0, tfinal, N + 1)
yy = [sol_exacte(t) for t in tt]
uu_EM = EM(phi, tt, y0)
uu_H  = heun(phi, tt, y0)

figure(figsize=(16,5))
subplot(1,2,1)
plot(tt,uu_EM,'o', tt,yy,'r-')
xlabel('$t$')
ylabel('$y$')
grid(True)
title('EM')
subplot(1,2,2)
plot(tt,uu_H,'o', tt,yy,'r-')
xlabel('$t$')
ylabel('$y$')
grid(True)
title('H')
show()

En effet, si on utilise les schémas d'ordre 1 proposés avec le même nombre de points de discrétisation on obtient

In [13]:
from scipy.optimize import fsolve

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

def EI(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        uu.append( fsolve(lambda x: -x+uu[i]+h*phi(tt[i+1],x) , uu[i] ))
    return uu


uu_EE = EE(phi, tt, y0)
uu_EI = EI(phi, tt, y0)

figure(figsize=(16,5))
subplot(1,2,1)
plot(tt,uu_EE,'o', tt,yy,'r-')
xlabel('$t$')
ylabel('$y$')
grid(True)
title('EE')
subplot(1,2,2)
plot(tt,uu_EI,'o', tt,yy,'r-')
xlabel('$t$')
ylabel('$y$')
grid(True)
title('EI')
show()

/usr/lib/python3/dist-packages/scipy/optimize/minpack.py:162: RuntimeWarning: The iteration is not making good progress, as measured by the 
  improvement from the last ten iterations.
  warnings.warn(msg, RuntimeWarning)
In [ ]: