2019 Contrôle Terminale session 1

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from scipy.optimize import fsolve
from IPython.display import display, Markdown, Math

EXERCICE 1 : choix d’un schéma

Soit un problème de Cauchy mathématiquement bien posé mais numériquement mal posé (i.e. sensible aux perturbations sur la donnée initiale). Parmi les schémas EE (Euler Explicite), EI (Euler Implicite) et H (Heun), lequel faut-il choisir? Justifier la réponse.

Toute méthode numérique induisant une erreur sur la solution, cette erreur va s’amplifier avec les itérations successives puisque le problème est numériquement mal posé. La seule issue pour traiter ce genre de problème est d’utiliser des schémas d’ordre élevé (ici donc le schéma de Heun qui est d’ordre 2 tandis que les schémas d’Euler sont d’ordre 1) avec un pas \(h\) petit (donc ce n’est pas la condition de A-stabilité qui limite le pas).

EXERCICE 2 : choix d’un schéma

Soit le problème de Cauchy \( \begin{cases} y'(t)=-10^5 y(t), &t\in[0;100]\\ y(0)=1 \end{cases} \) Parmi les schémas EE (Euler Explicite), EI (Euler Implicite) et H (Heun), lequel faut-il choisir? Justifier la réponse.

Le problème est bien posé mathématiquement et numériquement.
On ne peut pas utiliser ni le schéma EE ni le schéma H car la condition de stabilité est trop contraignante : dans les deux cas il faut un pas \(h<\frac{2}{10^5}\), i.e. une discrétisation avec \(N>\frac{100}{h}= \frac{10^7}{2}\) points.

Le schéma EI, en revanche, est inconditionnellement A-stable: on peut donc choisir le pas \(h\) sans contrainte de stabilité.

EXERCICE 3 : schéma predictor-corrector

Soit le problème de Cauchy \( \begin{cases} y'(t)=1+\vartheta\left(t-y(t)\right), &t\in[0;2]\\ y(0)=1 \end{cases} \)

1. Calculer la solution exacte de ce problème de Cauchy.
2. Soit les deux schémas \( \begin{aligned} &\begin{cases} u_0=y_0,\\ u_1=u_0+h\varphi(t_1,u_1),\\ u_{n+1}=\frac{4}{3}u_n-\frac{1}{3}u_{n-1}+\frac{2}{3}h\varphi(t_{n+1},u_{n+1})& n=1,2,3,\dots N-1 \end{cases} \\ &\begin{cases} u_0=y_0,\\ u_{1}=u_0+h\varphi(t_0,u_0),\\ u_{n+1}=u_n+\frac{h}{2}\Bigl(3\varphi(t_n,u_n)-\varphi(t_{n-1},u_{n-1})\Bigr)& n=1,2,3,4,5,\dots N-1 \end{cases} \end{aligned} \) Écrire le schéma predictor-corrector associé à ce couple de schémas sous la forme d’une suite définie par récurrence.
3. Implémenter le schéma predictor-corrector de la question précédente et le tester avec le problème de Cauchy. On choisira \(N=1+5\vartheta\) points de discrétisation et on affichera la solution approchée et la solution exacte sur le même repère.
4. Vérifier empiriquement l’ordre de convergence du schéma predictor-corrector de la question précédente en affichant la courbe de convergence et en estimant sa pente avec le problème de Cauchy.

Chaque sujet correspond à une valeur de \(\vartheta\) parmi les suivantes:

THETA = list(range(10,54,5))
THETA

[10, 15, 20, 25, 30, 35, 40, 45, 50]

Correction

On peut calculer la solution à la main: \( a(t)y'(t)+b(t)y(t)=g(t) \) avec \(a(t)=1\), \(b(t)=\vartheta\) et \(g(t)=1+\vartheta t\) donc + \(A(t)=\int \frac{b(t)}{a(t)}\mathrm{d}t=\int \vartheta\mathrm{d}t= \vartheta t\), + \(B(t)=\int \frac{g(t)}{a(t)} e^{A(t)}\mathrm{d}t=\int (1+\vartheta t) e^{\vartheta t}\mathrm{d}t= t e^{\vartheta t}\),

d’où \(y(t)=c e^{-\vartheta t}+ t\).

En imposant la condition \(1=y(0)\) on trouve l’unique solution du problème de Cauchy donné: \( \boxed{y(t)=e^{-\vartheta t}+t.} \)

Sinon, on peut utiliser le module de calcul symbolique.

t = sp.symbols('t')
y = sp.Function('y')

f = 1 + sp.symbols('theta')*(t-y(t))
edo = sp.Eq(sp.diff(y(t),t),f)
display(edo)

t0 = 0
y0 = 1
IC = sp.Eq(y(t0),y0)
display(IC)

sol = sp.dsolve(edo,y(t),ics={y(t0):y0})
display(sol)

sol_exacte = sp.lambdify((t, sp.symbols('theta')), sol.rhs, 'numpy')

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = \theta \left(t - y{\left(t \right)}\right) + 1\)

\(\displaystyle y{\left(0 \right)} = 1\)

\(\displaystyle y{\left(t \right)} = t + e^{- t \theta}\)

\( \begin{cases} u_0=y_0,\\ \tilde u_1=u_0+h\varphi(t_0,u_0),\\ u_1=u_0+h\varphi(t_1,\tilde u_1),\\ \tilde u_{n+1} = u_n+\frac{h}{2}\Bigl(3\varphi(t_n,u_n)-\varphi(t_{n-1},u_{n-1})\Bigr)& \\ u_{n+1}=\frac{4}{3}u_n-\frac{1}{3}u_{n-1}+\frac{2}{3}h\varphi(t_{n+1},\tilde u_{n+1})& n=1,2,3,\dots N-1 \end{cases} \)

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

t0, y0, tfinal = 0, 1, 2

i = 0
for theta in THETA:
    N = 1+5*theta
    tt = np.linspace(t0,tfinal,N+1)
    phi = lambda t,y : 1+theta*(t-y)
    yy = sol_exacte(tt,theta)
    uu = PC(phi,tt,y0)
    i += 1
    plt.figure(i,figsize=(20,5))
    plt.plot(tt,yy,'r-',label=('Exacte'))
    plt.plot(tt,uu,'bo',label=('PC'))
    plt.title('$\\theta$={}'.format(theta))
    plt.xlabel('$t$')
    plt.ylabel('$y$')
    plt.legend() 
    plt.grid(True)

On calcule la solution approchée avec différentes valeurs de \(h_k=2/N_k\), à savoir \(N_k=(1+5\vartheta)2^k\), \(k=0,...7\). On sauvegarde les valeurs de \(h_k\) dans le vecteur H.

Pour chaque valeur de \(h_k\), on calcule le maximum de la valeur absolue de l’erreur et on sauvegarde toutes ces erreurs dans le vecteur err de sort que err[k] contient \(e_k=\max_{i=0,...,N_k}|y(t_i)-u_{i}|\).

Pour afficher l’ordre de convergence on utilise une échelle logarithmique : on représente \(\ln(h)\) sur l’axe des abscisses et \(\ln(\text{err})\) sur l’axe des ordonnées ainsi, si \(\text{err}=Ch^p\), alors \(\ln(\text{err})=\ln(C)+p\ln(h)\). En échelle logarithmique, \(p\) représente donc la pente de la ligne droite \(\ln(\text{err})\).

Pour estimer l’ordre de convergence on estime la pente de la droite qui relie l’erreur au pas \(k\) à l’erreur au pas \(k+1\) en echelle logarithmique en utilisant la fonction polyfit basée sur la régression linéaire.

t0, y0, tfinal = 0, 1, 2

i = 0
for theta in THETA:
    H = []
    err = []
    for k in range(8):
        N = (1+5*theta)*2**k  
        tt = np.linspace(t0,tfinal,N+1)
        h = tt[1]-tt[0]
        H.append(h)
        phi = lambda t,y : 1+theta*(t-y)
        yy = sol_exacte(tt,theta)
        uu = PC(phi,tt,y0)
        err.append(np.linalg.norm(np.array(uu)-np.array(yy),np.inf))
    i += 1
    plt.figure(i,figsize=(20,5))
    plt.loglog(H,err, 'r-o')
    pente = np.polyfit(np.log(H), np.log(err), 1)[0]
    plt.title('$\\theta=${}, Ordre = {:1.4f}'.format(theta,pente))
    plt.xlabel('$h$')
    plt.ylabel('$e$')
    plt.grid(True)

EXERCICE 4 : Schéma Runge-Kutta

Soit le schéma de Runge-Kutta dont la matrice de Butcher est \( \begin{array}{c|cc} 1 & 1 & 0\\ \alpha & \alpha & 0\\ \hline & (1-\gamma) & \gamma \end{array} \)

1. Écrire le schéma sous la forme d’une suite définie par récurrence.
2. Étudier théoriquement l’ordre du schéma.
3. Implémenter le schéma et le tester avec le même problème de Cauchy qu’à l’exercice précedent. On choisira \(N=1+2\vartheta\) points de discrétisation et on affichera la solution approchée et la solution exacte sur le même repère. Attention : le schéma est semi-implicite, la fonction fsolve du module scipy.optimize est votre amie…

Chaque sujet correspond à un couple \((\alpha,\gamma)\) choisi parmi les suivants :

zeta = 1
GAMMA = [sp.Rational(13,24),sp.Rational(14,24),sp.Rational(15,24),sp.Rational(16,24),sp.Rational(17,24),
         sp.Rational(18,24),sp.Rational(19,24),sp.Rational(20,24),sp.Rational(21,24),sp.Rational(22,24),sp.Rational(23,24)]
for gamma in GAMMA:
    alpha = sp.Rational(sp.Rational(1,2)+zeta*(gamma-1),gamma)
    display(Markdown(rf"\(\alpha = {sp.latex(alpha)},\qquad \gamma = {sp.latex(gamma)}\)"))

\(\alpha = \frac{1}{13},\qquad \gamma = \frac{13}{24}\)

\(\alpha = \frac{1}{7},\qquad \gamma = \frac{7}{12}\)

\(\alpha = \frac{1}{5},\qquad \gamma = \frac{5}{8}\)

\(\alpha = \frac{1}{4},\qquad \gamma = \frac{2}{3}\)

\(\alpha = \frac{5}{17},\qquad \gamma = \frac{17}{24}\)

\(\alpha = \frac{1}{3},\qquad \gamma = \frac{3}{4}\)

\(\alpha = \frac{7}{19},\qquad \gamma = \frac{19}{24}\)

\(\alpha = \frac{2}{5},\qquad \gamma = \frac{5}{6}\)

\(\alpha = \frac{3}{7},\qquad \gamma = \frac{7}{8}\)

\(\alpha = \frac{5}{11},\qquad \gamma = \frac{11}{12}\)

\(\alpha = \frac{11}{23},\qquad \gamma = \frac{23}{24}\)

Correction

Considérons le problème de Cauchy \(\begin{cases} y'(t) = \varphi(t,y(t)), &\forall t \in I=[t_0,T],\\ y(t_0) = y_0, \end{cases}\) avec \(y_0\) une valeur donnée et supposons que l’on ait montré l’existence et l’unicité d’une solution \(y\) pour \(t\in I\).

Pour \(h>0\) soit \(t_n\stackrel{\text{déf.}}{=} t_0+nh\) avec \(n=0,1,2,\dots,N\) une suite de \(N+1\) nœuds de \(I\) induisant une discrétisation de \(I\) en \(N\) sous-intervalles \(I_n=[t_n;t_{n+1}]\) chacun de longueur \(h>0\) (appelé le pas de discrétisation).

Pour chaque nœud \(t_n\), on cherche la valeur inconnue \(u_n\) qui approche la valeur exacte \(y_n\equiv y(t_n)\).
- L’ensemble de \(N+1\) valeurs \(\{t_0, t_1=t_0+h,... , t_{N}=T \}\) représente les points de la discrétisation.
- L’ensemble de \(N+1\) valeurs \(\{y_0, y_1,... , y_{N} \}\) représente la solution exacte.
- L’ensemble de \(N+1\) valeurs \(\{u_0 = y_0, u_1,... , u_{N} \}\) représente la solution numérique. Cette solution approchée sera obtenue en construisant une suite récurrente.

Le schéma qu’on va construire permet de calculer implicitement \(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+h,u_n+hK_1\right)\\ K_2 = \varphi\left(t_n+\alpha h,u_n+\alpha hK_1\right)\\ u_{n+1} = u_n + h \left((1-\gamma)K_1+\gamma K_2\right) & n=0,1,\dots N-1 \end{cases}\)

• C’est un schéma semi-implicite à deux étages: l’ordre est au plus 4.
• \(c_1=a_{11}+a_{12}\), \(c_2=a_{21}+a_{22}\) et \(1=b_1+b_2\) \(\implies\) il est consistante
• \(c_1 b_1 + c_2b_2=\frac{1}{2}\) \(\implies\) il est d’ordre \(2\)
• \(b_1c_1^2+b_2c_2^2=\frac{1}{2\gamma}\) et \(b_1a_{11}c_1+b_2a_{12}c_2=\frac{1}{2}\) \(\implies\) il n’est pas d’ordre \(3\)

On peut vérifier ces calculs avec sympy:

gamma = sp.Symbol(r'\gamma')
alpha = (gamma-sp.S(1)/2)/gamma
c = [1,alpha]
b = [1-gamma,gamma]
A = [[1,0],[alpha,0]]
s = len(c)


display(Markdown(  r'\(\sum_{j=1}^s b_j='+sp.latex(sum(b))+'\)'    ))

for i in range(s):
    display(Markdown(
        rf"\(i = {i} \quad \sum_{{j=1}}^{{{s}}} a_{{ij}} - c_i = {sp.latex(sum(A[i]) - c[i])}\)"
    ))

display(Markdown(
    rf"\(\sum_{{j=1}}^{{{s}}} b_j c_j = {sp.latex(sum(b[i]*c[i] for i in range(s)))}\)"
))
display(Markdown(
    rf"\(\sum_{{j=1}}^{{{s}}} b_j c_j^2 = {sp.latex(sum(b[i]*c[i]**2 for i in range(s)).simplify())}\)"
))
display(Markdown(
    rf"\(\sum_{{i,j=1}}^{{{s}}} b_i a_{{ij}} c_j = {sp.latex(sum(b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)))}\)"
))

\(\sum_{j=1}^s b_j=1\)

\(i = 0 \quad \sum_{j=1}^{2} a_{ij} - c_i = 0\)

\(i = 1 \quad \sum_{j=1}^{2} a_{ij} - c_i = 0\)

\(\sum_{j=1}^{2} b_j c_j = \frac{1}{2}\)

\(\sum_{j=1}^{2} b_j c_j^2 = \frac{1}{4 \gamma}\)

\(\sum_{i,j=1}^{2} b_i a_{ij} c_j = \frac{1}{2}\)

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

def RK(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = np.empty_like(tt)
    uu[0] = y0
    for n in range(len(tt)-1):
        k1 = fsolve( lambda x : x - phi( tt[n+1], uu[n]+h*x ) , phi(tt[n],uu[n]) )[0]
        k2 = phi( tt[n]+alpha*h , uu[n]+alpha*h*k1 )
        uu[n+1] = uu[n] + h*((1-gamma)*k1+gamma*k2)
    return uu

Comparons solution exacte et solution approchée (pour trouver votre sujet vous devez chercher les valeurs de \(\vartheta\) et \(\beta\) correspondants):

t0, y0, tfinal = 0, 1, 2

for theta in THETA:

    N = 1+2*theta
    tt = np.linspace(t0,tfinal,N+1)
    sol_exacte = lambda t : t+np.exp(-theta*t)
    phi = lambda t,y : 1+theta*(t-y)
    yy = sol_exacte(tt)
    fig = plt.figure(theta,figsize=(20,5))
    fig.suptitle(rf'Cas $\theta$={theta}', fontsize=16)
    
    g = len(GAMMA)
    for i, gamma in enumerate(GAMMA, start=1):
        alpha = float(sp.Rational(sp.Rational(1,2)+zeta*(gamma-1),gamma))
        uu = RK(phi,tt,y0)
        plt.subplot(1,g,i)
        plt.plot(tt,yy,'b-',label=("Exacte"))
        plt.plot(tt,uu,'ro',label=("RK"))
        plt.title(rf'$\alpha={alpha:.3f}$'+'\n'+rf'$\gamma={float(gamma):.3f}$')
        plt.xlabel('$t$')
        plt.ylabel('$u$')
        plt.legend() 
        plt.grid(True)

    fig.tight_layout(rect=[0, 0, 1, 0.92])