2021 - Contrôle Terminal session 1

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

1. Écrire le schéma sous la forme d’une suite définie par récurrence.
   Le schéma est explicite, semi-implicite ou implicite ?
2. Sans faire de calculs, indiquer quel est l’ordre maximale du schéma.
3. Étudier théoriquement l’ordre du schéma en fonction de $\eta$ et $\mu$.
4. Étudier théoriquement la A-stabilité du schéma pour les valeurs de $\eta$ et $\mu$ qui donnent l’ordre maximal.
5. 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.

%reset -f
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

from IPython.display import display, Math, Markdown
prlat= lambda *args: display(Math(''.join( sympy.latex(arg) for arg in args )))

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

def ordre_RK(s, A=None, b=None, c=None):

    j = sympy.symbols('j')
    
    if A is None:
        A = sympy.MatrixSymbol('a', s, s)
    else:
        A = sympy.Matrix(A)
    if c is None:
        c = sympy.symbols(f'c_0:{</span>s<span class="sc">}')    
    if b is None:
        b = sympy.symbols(f'b_0:{</span>s<span class="sc">}')
    
    
    display(Markdown("**Matrice de Butcher**"))
    But = matrice_Butcher(s, A, b, c)
    
    Eqs = []

    display(Markdown(f"**On a {</span>s<span class="op">+</span><span class="dv">1</span><span class="sc">} conditions pour avoir consistance = pour être d'ordre 1**"))
    ordre_1(s, A, b, c, j)
    
    display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
    Eq_2 = ordre_2(s, A, b, c, j)
    Eqs.append(Eq_2)

    return Eqs


def matrice_Butcher(s, A, b, c):
    But = sympy.Matrix(A)
    But = But.col_insert(0, sympy.Matrix(c))   
    last = [sympy.Symbol(" ")]
    last.extend(b)
    But = But.row_insert(s,sympy.Matrix(last).T)
    display(Math(sympy.latex(sympy.Matrix(But))))
    return But

def ordre_1(s, A, b, c, j):
    texte = "\sum_{j=1}^s b_j =" + f"{</span><span class="bu">sum</span>(b)<span class="sc">.</span>simplify()<span class="sc">}"
    texte += r"\text{ doit être égale à }1"
    display(Math(texte))
    for i in range(s):
        somma = sympy.summation(A[i,j],(j,0,s-1)).simplify()
        texte = r'\sum_{j=1}^s a_{'</span><span class="op">+</span><span class="bu">str</span>(i)<span class="op">+</span><span class="vs">r'j}=' + sympy.latex( somma )
        texte += r"\text{ doit être égale à }"+sympy.latex(c[i])
        display(Math( texte ))
    return None    

def ordre_2(s, A, b, c, j):
    texte = '\sum_{j=1}^s b_j c_j=' +sympy.latex(sum([b[i]*c[i] for i in range(s)]).simplify())
    texte += r"\text{ doit être égale à }\frac{1}{2}"
    display(Math(texte))
    Eq_2 = sympy.Eq( sum([b[i]*c[i] for i in range(s)]).simplify(), sympy.Rational(1,2) )
    return Eq_2




##################################################
# Conditions avec s étages
##################################################

eta = sympy.Symbol(r'\eta')
mu = sympy.Symbol(r'\mu')

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


Eqs = ordre_RK(s, A, b, c)

display(Markdown("**Solutions**"))

sol = sympy.solve(Eqs)[0]
sol
# for eq in Eqs:
#     display(eq)

Matrice de Butcher

$\displaystyle \left[\begin{matrix}0 & 0 & 0\\\frac{\eta}{4} & \frac{\eta}{4} & 0\\ & 1 - \frac{\mu}{3} & \frac{\mu}{3}\end{matrix}\right]$

On a 3 conditions pour avoir consistance = pour être d’ordre 1

$\displaystyle \sum_{j=1}^s b_j =1\text{ doit être égale à }1$

$\displaystyle \sum_{j=1}^s a_{0j}=0\text{ doit être égale à }0$

$\displaystyle \sum_{j=1}^s a_{1j}=\frac{\eta}{4}\text{ doit être égale à }\frac{\eta}{4}$

On doit ajouter 1 condition pour être d’ordre 2

$\displaystyle \sum_{j=1}^s b_j c_j=\frac{\eta \mu}{12}\text{ doit être égale à }\frac{1}{2}$

Solutions

$\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 $y'(t)=-\beta y(t)$ pour $t>0$ et $y(0)=1$. 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}$

Dans le cas optimal on a $\eta\mu=6$, ainsi
$u_{n+1} = \left(\dfrac{1}{2}(\beta h)^2 -(\beta h)+1\right)u_n$

u_n = sympy.Symbol('u_n')
dt = sympy.Symbol('h') # je change le nom de la variable "h" en "dt" pour ne pas effacer l'affectation de "h"
beta = sympy.Symbol(r'\beta')

K1 = -beta*u_n
display(Math('K_1='+sympy.latex(K1)))
K2 = -beta*(u_n+sympy.S(eta)/4*dt*K1).factor()
display(Math('K_2='+sympy.latex(K2)))
u_np1 = (u_n+dt*(sympy.S(3-mu)/3*K1+sympy.S(mu)/3*K2)).factor()
display(Math('u_{n+1}='+sympy.latex(u_np1)))

display(Markdown("Sachant que"))
display(sol)
display(Markdown("on a"))
u_np1 = u_np1.subs(sol)
display(Math('u_{n+1}='+sympy.latex(u_np1)))

$\displaystyle K_1=- \beta u_{n}$

$\displaystyle K_2=\frac{\beta u_{n} \left(\beta \eta h - 4\right)}{4}$

$\displaystyle u_{n+1}=\frac{u_{n} \left(\beta^{2} \eta \mu h^{2} - 12 \beta h + 12\right)}{12}$

Sachant que

$\displaystyle \left\{ \eta : \frac{6}{\mu}\right\}$

on a

$\displaystyle u_{n+1}=\frac{u_{n} \left(6 \beta^{2} h^{2} - 12 \beta h + 12\right)}{12}$

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:

x = sympy.Symbol('x',nonnegative=True)

q = (u_np1/u_n).subs(beta*dt,x)
display(Math('q(x)='+sympy.latex(q)))

display(Markdown("$|q(x)|<1$ pour"))
display(sympy.solve(abs(q)<1))

sympy.plot(q,xlim=[0,3],ylim=[-1.5,1.5]);


# display(Math('q(0)='+sympy.latex(q.subs(x,0))))
# lim=sympy.Limit(q,x,sympy.oo)
# display(Math(sympy.latex(lim)+'='+sympy.latex(lim.doit())))

# print('On cherche les points stationnaires dans R^+')
# dq=(q.diff(x)).factor()
# display(Math("q'(x)="+sympy.latex(dq)))
# print("q est une parabole dont le sommet est")
# sol=sympy.solve(dq)
# display(Math("q'(x)=0 \iff x\in"+sympy.latex(sol)))
# minimum=sol[0]
# display(Math("x="+sympy.latex(minimum)+"\\text{ est un minimum et}"))
# qmin=q.subs(x,sol[0])
# display(Math("q("+sympy.latex(minimum)+")="+sympy.latex(qmin)))
# sympy.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(sympy.solve(q<1))
# print("q(x)>-1")
# display(sympy.solve(q>-1))
# print("Notons que q(x)>0 pour tout 0<x<2")
# #display(sympy.solve(q>0))

$\displaystyle q(x)=\frac{x^{2}}{2} - x + 1$

$|q(x)|<1$ pour

$\displaystyle 0 < x \wedge 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.

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])))

eta, mu = L[idx], 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.

# ==================================================
# SCHEMA DE RUNGE-KUTTA
# ==================================================

# def RK(tt,phi,y0): # comme suite geometrique
#     x = beta*h
#     q = x**2/2-x+1
#     uu = y0*q**np.arange(len(tt))
#     return uu

def RK(tt,phi,y0):
    uu = np.zeros_like(tt)
    uu[0] = y0
    h = tt[1]-tt[0]
    for n in range(len(tt)-1):
        K1 = phi(tt[n],uu[n])
        K2 = phi(tt[n]+eta/4*h,uu[n]+eta/4*h*K1)
        uu[n+1] = uu[n]+h/3*((3-mu)*K1+mu*K2) 
    return uu


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

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


# ==================================================
# COMPARAISON POUR N DONNE
# ==================================================
plt.figure(figsize=(20,5))

plt.subplot(1,2,1)

N = 50
tt = np.linspace(t0,tfinal,N+1)
h = tt[1]-tt[0]

if h>=2/beta:
    print('Attention: h > limite de stabilité')
    
yy = sol_exacte(tt)   
# uu = RK(tt)
uu = RK(tt,phi,y0)
plt.plot(tt,yy,'b-',label=("Exacte"))
plt.plot(tt,uu,'ro',label=("RK"))
plt.title(f'${</span>N<span class="op">=</span><span class="sc">}$, ${</span>h<span class="op">=</span><span class="sc">} < {</span><span class="dv">2</span><span class="op">/</span>beta<span class="sc">}$')
plt.xlabel('$t$')
plt.ylabel('$u$')
plt.legend() 
plt.grid(True)



# ==================================================
# ORDRE DE CONVERGENCE
# ==================================================

plt.subplot(1,2,2)

H = []
err = []
N = 50
for k in range(7):
    N += 100
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu = RK(tt,phi,y0)
    H.append(h)
    err.append(np.linalg.norm(yy-uu, np.inf))
    
omega = np.polyfit(np.log(H),np.log(err), 1)[0]
plt.plot(np.log(H), np.log(err), 'c-o')
plt.xlabel('$\ln(h)$')
plt.ylabel('$\ln(e)$')
plt.title(rf'Ordre $\omega$={</span>omega<span class="sc">:1.2f}')
plt.grid(True);

2 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).

%reset -f
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

from IPython.display import display, Math, Markdown
prlat= lambda *args: display(Math(''.join( sympy.latex(arg) for arg in args )))

• On définit la solution exacte (calculée en utilisant le module sympy):

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

t0     = 0
tfinal = 2
y0     = 0

# NB il faut utiliser le cos du module sympy et non de numpy
edo = sympy.Eq( sympy.diff(y(t),t) , -50*(y(t)-sympy.cos(t)) )
display(edo)
solpar = sympy.dsolve(edo, y(t), ics={y(t0):y0})
display(solpar)

sympy.plot(solpar.rhs, (t, t0, tfinal), ylabel='$y(t)$', xlabel='$t$');

sol_exacte = sympy.lambdify(t,solpar.rhs,'numpy')

$\displaystyle \frac{d}{d t} y{\left(t \right)} = - 50 y{\left(t \right)} + 50 \cos{\left(t \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é.

• On calcule la solution approchée

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


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

N = 40
tt = np.linspace(t0, tfinal, N + 1)
h = tt[1] - tt[0]
yy = sol_exacte(tt)

# uu = EE(phi, tt, y0) # pour montrer les pb de stabilité avec le même nombre de points
uu = EI(phi, tt, y0)
err = np.linalg.norm(yy-uu, np.inf)

plt.figure(figsize=(8,5))
plt.plot(tt, uu, 'o', label='EI')
plt.plot(tt, yy, 'r-', label='Exacte')
plt.xlabel('$t$')
plt.ylabel('$y$')
plt.grid(True)
plt.legend()
plt.title(f'{</span>h <span class="op">=</span> <span class="sc">}')
plt.show();

3 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$, écrivez votre nom et prénom dans la cellule ci-dessous et vous obtiendrez un choix pour ce paramètre.

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).

%reset -f
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve

import sympy 
sympy.init_printing()

from IPython.display import display, Math, Markdown
prlat= lambda *args: display(Math(''.join( sympy.latex(arg) for arg in args )))

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$

• On calcule la solution exacte (calculée en utilisant le module sympy) :

t       = sympy.Symbol('t')
epsilon = sympy.Symbol(r'\varepsilon')
y       = sympy.Function('y')

t0 = 0
y0 = sympy.Rational(1,alpha)+epsilon
display(Math('y(0)='+sympy.latex(y0)))

edo = sympy.Eq( sympy.diff(y(t),t) , alpha*(y(t)-t) )
display(edo)
solpar = sympy.dsolve(edo, y(t), ics={y(t0):y0})
display(solpar)

tt = np.linspace(0,8,101)
plt.grid()
for val in [1.e-7,0,-1.e-7]:
    sol_exacte = sympy.lambdify(t,solpar.rhs.subs(epsilon,val),'numpy')
    plt.plot(tt,sol_exacte(tt), label=f'$\epsilon={</span>val<span class="sc">}$')
plt.legend();

$\displaystyle y(0)=\varepsilon + \frac{1}{2}$

$\displaystyle \frac{d}{d t} y{\left(t \right)} = - 2 t + 2 y{\left(t \right)}$

$\displaystyle y{\left(t \right)} = \varepsilon e^{2 t} + t + \frac{1}{2}$

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

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

t0     = 0
tfinal = 8
y0     = 1/alpha

N = 5
tt = np.linspace(t0, tfinal, N + 1)

sol_exacte = sympy.lambdify(t,solpar.rhs.subs(epsilon,0),'numpy')
yy = sol_exacte(tt)

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

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


fig,ax = plt.subplots(nrows=1,ncols=2,figsize=(16,5))
ax[0].plot(tt, yy, 'r-', label='Exacte')
ax[1].plot(tt, yy, 'r-', label='Exacte')
for i,schema in enumerate([EM, heun]):
    uu = schema(phi, tt, y0)
    ax[i].plot(tt, uu, 'o', label=schema.__name__)
    ax[i].set_xlabel('$t$')
    ax[i].set_ylabel('$y$')
    ax[i].grid(True)
    ax[i].legend();

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

def EE(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
        k1 = phi( tt[n], uu[n] )
        uu[n+1] = uu[n]+h*k1
    return uu

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


fig,ax = plt.subplots(nrows=1,ncols=2,figsize=(16,5))
ax[0].plot(tt, yy, 'r-', label='Exacte')
ax[1].plot(tt, yy, 'r-', label='Exacte')
for i,schema in enumerate([EE, EI]):
    uu = schema(phi, tt, y0)
    ax[i].plot(tt, uu, 'o', label=schema.__name__)
    ax[i].set_xlabel('$t$')
    ax[i].set_ylabel('$y$')
    ax[i].grid(True)
    ax[i].legend();

Remarque : dans les anciennes versions de python, les calculs avec des flottants donnaient des résultats différents.