M62 : Correction Examen 2019

1 M62 : Correction Examen 2019

1.1 Choix d’un schéma

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

Correction Q2
Toute méthode numérique induisant une erreur sur la solution, cette erreur va s’amplifier avec les itérations succéssives 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 élévé (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).

1.2 Choix d’un schéma

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

Correction Q3
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é.

1.3 Solution exacte

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

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]

Q4 Calculer la solution exacte de ce problème de Cauchy.

Correction Q4
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 mais attention, pour ne pas avoir des mauvaises intéractions entre le module matplotlib et le module sympy, on utilisera l’import de sympy avec un alias:

import sympy as sym

import sympy as sym
sym.init_printing()

sym.var('t,C1,theta')
u=sym.Function('u')
f=1+theta*(t-u(t))
edo=sym.Eq(sym.diff(u(t),t),f)
print("EDO:")
display(edo)

solgen=sym.dsolve(edo,u(t)).simplify()
print("Solution générale:")
display(solgen)

t0=0
u0=1
consts = sym.solve( [solgen.rhs.subs(t,t0)-u0 ],C1, dict=True)[0]
print("Prise en compte de la condition initiale u({})={}:".format(t0,u0))
display(consts)

solpar=solgen.subs(consts)
print("Solution du problème de Cauchy:")
display(solpar)

EDO:
Solution générale:
Prise en compte de la condition initiale u(0)=1:
Solution du problème de Cauchy:

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

\(\displaystyle u{\left(t \right)} = C_{1} e^{- t \theta} + t\)

\(\displaystyle \left\{ C_{1} : 1\right\}\)

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

1.4 Schéma Predictor-corrector

Soit les deux schémas \(\begin{align} &\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{align}\)

Q5 Écrire le schéma predictor-corrector associé à ce couple de schémas sous la forme d’une suite définie par récurrence.

Correction Q5 \( \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} \)

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

Correction Q6

from matplotlib.pylab import *

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 = linspace(t0,tfinal,N+1)
    sol_exacte = lambda t : t+exp(-theta*t)
    phi = lambda t,y : 1+theta*(t-y)
    yy = [sol_exacte(t) for t in tt]
    uu = PC(phi,tt,y0)
    i+=1
    figure(i,figsize=(20,5))
    plot(tt,yy,'r-',label=('Exacte'))
    plot(tt,uu,'bo',label=('PC'))
    title('$\\theta$={}'.format(theta))
    xlabel('$t$')
    ylabel('$y$')
    legend() 
    grid(True)

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

Correction Q7

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 = linspace(t0,tfinal,N+1)
        h = tt[1]-tt[0]
        H.append(h)
        sol_exacte = lambda t : t+exp(-theta*t)
        phi = lambda t,y : 1+theta*(t-y)
        yy = [sol_exacte(t) for t in tt]
        uu = PC(phi,tt,y0)
        err.append(max([abs(uu[i]-yy[i]) for i in range(len(uu))]))
    i+=1
    figure(i,figsize=(20,5))
    loglog(H,err, 'r-o')
    title('$\\theta=${}, Ordre = {:1.4f}'.format(theta,polyfit(log(H),log(err), 1)[0]))
    xlabel('$h$')
    ylabel('$e$')
    #legend() 
    grid(True)

1.5 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} \)

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

from IPython.display import display, Math
from fractions import Fraction
zeta = 1
GAMMA = [Fraction(13,24),Fraction(14,24),Fraction(15,24),Fraction(16,24),Fraction(17,24),
         Fraction(18,24),Fraction(19,24),Fraction(20,24),Fraction(21,24),Fraction(22,24),Fraction(23,24)]
for gamma in GAMMA:
    alpha = Fraction(Fraction(1,2)+zeta*(gamma-1),gamma)
    display(Math(rf"\alpha = {</span>alpha<span class="sc">},\qquad \gamma = {</span>gamma<span class="sc">}"))

\(\displaystyle \alpha = 1/13,\qquad \gamma = 13/24\)

\(\displaystyle \alpha = 1/7,\qquad \gamma = 7/12\)

\(\displaystyle \alpha = 1/5,\qquad \gamma = 5/8\)

\(\displaystyle \alpha = 1/4,\qquad \gamma = 2/3\)

\(\displaystyle \alpha = 5/17,\qquad \gamma = 17/24\)

\(\displaystyle \alpha = 1/3,\qquad \gamma = 3/4\)

\(\displaystyle \alpha = 7/19,\qquad \gamma = 19/24\)

\(\displaystyle \alpha = 2/5,\qquad \gamma = 5/6\)

\(\displaystyle \alpha = 3/7,\qquad \gamma = 7/8\)

\(\displaystyle \alpha = 5/11,\qquad \gamma = 11/12\)

\(\displaystyle \alpha = 11/23,\qquad \gamma = 23/24\)

Q8 Écrire le schéma sous la forme d’une suite définie par récurrence.

Correction Q8

Considérons le problème de Cauchy >trouver une fonction \(y \colon I\subset \mathbb{R} \to \mathbb{R}\) définie sur un intervalle \(I=[t_0,T]\) telle que \(\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}\)

Q9 Étudier théoriquement l’ordre du schéma.

Correction Q9

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

import sympy as sym
sym.init_printing()

from IPython.display import display, Math

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

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

display(Math(r'\sum_{j=1}^s b_j c_j='+sym.latex(sum([b[i]*c[i] for i in range(s)]))     ))

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

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

\(\displaystyle i={}0\quad \sum_{j=1}^s a_{ij}-c_i=\mathtt{\text{0}}\)

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

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

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

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

Q10 Implémenter le schéma de la question précédente et le tester avec le problème de Cauchy.
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…

Correction Q10

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

from scipy.optimize import fsolve

def RK(phi,tt,y0):
    h = tt[1]-tt[0]
    uu = [y0]
    for i in range(len(tt)-1):
        k1 = fsolve( lambda x : x - phi( tt[i+1], uu[i]+h*x ) , phi(tt[i],uu[i]) )[0]
        k2 = phi( tt[i]+alpha*h , uu[i]+alpha*h*k1 )
        uu.append( uu[i] + 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):

from matplotlib.pylab import *

t0, y0, tfinal = 0, 1, 2

sol_exacte = lambda t : t+exp(-theta*t)
phi = lambda t,y : 1+theta*(t-y)


for theta in THETA:
    N = 1+2*theta
    tt = linspace(t0,tfinal,N+1)
    yy = [sol_exacte(t) for t in tt]
    figure(theta,figsize=(20,5))
    g = len(GAMMA)
    i = 0
    for gamma in GAMMA:
        alpha=Fraction(Fraction(1,2)+zeta*(gamma-1),gamma)
        uu = RK(phi,tt,y0)
        i+=1
        subplot(1,g,i)
        plot(tt,yy,'b-',label=("Exacte"))
        plot(tt,uu,'ro',label=("RK"))
        title(r' $\theta$={}, $\alpha$={}'.format(theta,alpha))
        xlabel('$t$')
        ylabel('$u$')
        legend() 
        grid(True)