2019 Contrôle Continu

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp

from IPython.display import display, Markdown, Math

EXERCICE : schéma de Runge-Kutta

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

Remarque : chaque sujet correspond à une valeur de \(\alpha\) choisie parmi les suivantes : \( \left\{\frac{1}{6},\frac{1}{4},\frac{1}{3},\frac{5}{12},\frac{7}{12},\frac{2}{3},\frac{3}{4},\frac{5}{6},\frac{11}{12}\right\} \)

1. Écrire le schéma sous la forme d’une suite définie par récurrence.
2. Étudier théoriquement la A-stabilité.
3. Implémenter le schéma et le tester avec le problème de Cauchy \(y'(t)=\beta y(t)\), \(y(0)=1\) sur l’intervalle \([0;1]\). On choisira le nombre minimal de points pour avoir A-stabilité. Remarque : chaque sujet correspond à une valeur de \(\beta\) choisie parmi les suivantes : \(\left\{50,60,70,80,90,100\right\}\)
4. Étudier théoriquement l’ordre du schéma.
5. Calculer la solution exacte du problème de Cauchy \( \begin{cases} y'(t)=-2ty^2(t), & t\in[0;3]\\ y(0)=2 \end{cases} \)
6. Vérifier empiriquement l’ordre de convergence en affichant la courbe de convergence et en estimant la pente sur ce problème de Cauchy.

ALPHA = [sp.Rational(2,12), sp.Rational(3,12), sp.Rational(4,12), sp.Rational(5,12),
         sp.Rational(7,12), sp.Rational(8,12), sp.Rational(9,12), sp.Rational(10,12), sp.Rational(11,12)]
for alpha in ALPHA:
    display(Markdown(f"\(c_2 = a_{{21}} = {sp.latex(alpha)} \\qquad b_1={sp.latex(1-sp.Rational(1,2*alpha))} \\qquad b_2={sp.latex(sp.Rational(1,2*alpha))}\)"))

\(c_2 = a_{21} = \frac{1}{6} \qquad b_1=-2 \qquad b_2=3\)

\(c_2 = a_{21} = \frac{1}{4} \qquad b_1=-1 \qquad b_2=2\)

\(c_2 = a_{21} = \frac{1}{3} \qquad b_1=- \frac{1}{2} \qquad b_2=\frac{3}{2}\)

\(c_2 = a_{21} = \frac{5}{12} \qquad b_1=- \frac{1}{5} \qquad b_2=\frac{6}{5}\)

\(c_2 = a_{21} = \frac{7}{12} \qquad b_1=\frac{1}{7} \qquad b_2=\frac{6}{7}\)

\(c_2 = a_{21} = \frac{2}{3} \qquad b_1=\frac{1}{4} \qquad b_2=\frac{3}{4}\)

\(c_2 = a_{21} = \frac{3}{4} \qquad b_1=\frac{1}{3} \qquad b_2=\frac{2}{3}\)

\(c_2 = a_{21} = \frac{5}{6} \qquad b_1=\frac{2}{5} \qquad b_2=\frac{3}{5}\)

\(c_2 = a_{21} = \frac{11}{12} \qquad b_1=\frac{5}{11} \qquad b_2=\frac{6}{11}\)

Correction : écriture schéma

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 explicitement \(u_{n+1}\) à partir de \(u_n\) et il est donc possible de calculer successivement \(u_1\), \(u_2\),…, en partant de \(u_0\) 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+\alpha h,u_n+\alpha hK_1\right)\\ u_{n+1} = u_n + \frac{h}{2\alpha} \left((2\alpha-1)K_1+K_2\right) & n=0,1,\dots N-1 \end{cases}\)

Correction : étude de la A-stabilité

Définition de A-stabilité

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)=y_0 \end{cases}\) où \(y_0\neq0\) est une valeur donnée. Sa solution est \(y(t)=y_0e^{-\beta t}\) donc \(\lim_{t\to+\infty}y(t)=0.\)

Si, sous d’éventuelles conditions sur \(h\), on a \( \lim_{n\to+\infty} u_n =0, \) alors on dit que le schéma est A-stable.

Le schéma donné appliqué à cette équation devient

h = sp.symbols('h')
Alpha = sp.symbols(r'\Alpha')
beta = sp.symbols(r'\beta')
u_n = sp.symbols('u_n')
u_np1 = sp.symbols('u_{n+1}')

Phi = lambda y : -beta*y

K_1 = Phi(u_n)
display(sp.Eq(sp.Symbol('K_1'), K_1))

K_2 = Phi(u_n+Alpha*h*K_1)
display(sp.Eq(sp.Symbol('K_2'), K_2.expand().factor(u_n)))

display(sp.Eq(u_np1,(u_n+h/(2*Alpha)*((2*Alpha-1)*K_1+K_2  )).factor()))

\(\displaystyle K_{1} = - \beta u_{n}\)

\(\displaystyle K_{2} = \beta u_{n} \left(\Alpha \beta h - 1\right)\)

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

c’est-à-dire \(\begin{cases} u_0 = y_0 \\ u_{n+1} = \dfrac{2 - 2(\beta h) + (\beta h)^2}{2}u_n & n=0,1,\dots N-1 \end{cases}\) Par induction on obtient \( u_{n}=\left(1 - (\beta h) + \frac{1}{2}(\beta h)^2 \right)^nu_0. \) Notons que \(u_n\) ne dépend pas du choix de \(\alpha\).

Étant une suite géométrique, \(\lim\limits_{n\to+\infty}u_n=0\) si et seulement si \( \left|1 - (\beta h) + \frac{1}{2}(\beta h)^2\right|<1. \) Notons \(x\) le produit \(\beta h\) et \(q\) le polynôme \(q(x)=\frac{1}{2}x^2-x+1\).

# from sympy.plotting import plot

x = sp.symbols('x', positive=True)
q = 1 - 1*x + x**2/2
display(Markdown(f"\(q(x)={sp.latex(q)}\)"))

display(Markdown("Calculons les valeurs de $x$ pour lesquelles $|q(x)|<1$:"))
display(sp.solve(sp.Abs(q)<1))

display(Markdown("Calculons la dérivée de $q$ pour trouver les extremums locaux:"))
dq = sp.diff(q,x)
display(sp.Eq(sp.Symbol("q'(x)"), dq))

display(Markdown("Calculons les valeurs de $x$ pour lesquelles $q'(x)=0$:"))
x_sommet = sp.solve(dq)[0]
display(x_sommet)

display(Markdown("Calculons la valeur de $q$ au sommet:"))
y_sommet=q.subs(x,x_sommet)
display(y_sommet)

display(Markdown("Comme $y_{sommet}=q(x_{sommet})>0$, on en déduit que $q(x)>0$ pour tout $x$ et que $q(x)<1$ ssi $x<2x_{sommet}$:"))
display(Markdown(f"ainsi $q(x)>0$ pour tout $x$ et $q(x)<1$ ssi $x<{sp.latex(2*x_sommet)}$"))
sp.plot(q,1,-1,(x,0,x_sommet*2));

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

Calculons les valeurs de \(x\) pour lesquelles \(|q(x)|<1\):

\(\displaystyle x < 2\)

Calculons la dérivée de \(q\) pour trouver les extremums locaux:

\(\displaystyle q'(x) = x - 1\)

Calculons les valeurs de \(x\) pour lesquelles \(q'(x)=0\):

\(\displaystyle 1\)

Calculons la valeur de \(q\) au sommet:

\(\displaystyle \frac{1}{2}\)

Comme \(y_{sommet}=q(x_{sommet})>0\), on en déduit que \(q(x)>0\) pour tout \(x\) et que \(q(x)<1\) ssi \(x<2x_{sommet}\):

ainsi \(q(x)>0\) pour tout \(x\) et \(q(x)<1\) ssi \(x<2\)

Il s’agit de la parabole convexe de sommet \(\left(1,\frac{1}{2}\right)\).
Nous avons \(|q(x)|<1\) si et seulement si \(x<2\). La relation \(\lim\limits_{n\to+\infty}u_n=0\) est donc satisfaite si et seulement si \( h <\frac{2}{\beta}. \) Cette condition de stabilité limite le pas \(h\) d’avance en \(t\) lorsqu’on utilise le schéma.

Nous avons \(q(x)>0\) pour tout \(x\) donc cette convergence est monotone.

Correction : implémentation et test

\(h<\frac{2}{\beta}\) sur l’intervalle \([0;1]\) impose \(N>\frac{\beta}{2}\).

def RKalpha(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]+h*alpha , uu[i]+h*alpha*k1 )
        uu.append( uu[i] + h*((2*alpha-1)*k1+k2 )/(2*alpha))
    return uu

BETA = np.arange(50,101,10)

t0, y0, tfinal = 0, 1, 1

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

i = 0


for beta in BETA:
    N=(int(beta/2)+1)
    tt = np.linspace(t0,tfinal,N+1)
    yy = [sol_exacte(t) for t in tt]
    alpha=ALPHA[0]
    uu = RKalpha(phi,tt,y0) 
    i+=1
    plt.figure(i,figsize=(20,5))
    plt.plot(tt,yy,lw=2,label=("Exacte"))
    plt.plot(tt,uu,'r*',label=("N="+str(N)+" beta="+str(beta)))
    plt.xlabel('$t$')
    plt.ylabel('$u$')
    plt.legend() 
    plt.grid(True)

Correction : étude de l’ordre

• C’est un schéma explicite à deux étage: l’ordre est au plus 2.
• \(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\)

for alpha in ALPHA:
    print("alpha = ",alpha)
    c = [0,alpha]
    b = [1-1/(2*alpha),1/(2*alpha)]
    A = [[0,0],[alpha,0]]
    s = len(c)
    ordre = []
    ordre.append(sum(b)==1)
    for i in range(s):
        ordre.append(sum(A[i])==c[i])
    ordre.append(sum([b[i]*c[i] for i in range(s)])==sp.Rational(1,2))
    print(set(ordre))
    print("---------------------")

alpha =  1/6
{True}
---------------------
alpha =  1/4
{True}
---------------------
alpha =  1/3
{True}
---------------------
alpha =  5/12
{True}
---------------------
alpha =  7/12
{True}
---------------------
alpha =  2/3
{True}
---------------------
alpha =  3/4
{True}
---------------------
alpha =  5/6
{True}
---------------------
alpha =  11/12
{True}
---------------------

Correction : solution exacte

On peut soit calculer la solution “à la main” soit utiliser sympy.

Il s’agit d’une EDO à variables séparables. La fonction \(y(t)=0\) pour tout \(t\) est solution de l’EDO mais elle ne vérifie pas la CI. Toute autre solution de l’EDO sera non nulle et se trouve formellement comme suit: \( y'(t)=-2ty^2(t) \quad\implies\quad \frac{y'(t)}{y^2(t)}=-2t \quad\implies\quad \int y^{-2}\mathrm{d}y=-2\int t\mathrm{d}t \quad\implies\quad y(t)=\frac{1}{t^2+c}, \ c\in\mathbb{R}. \) En imposant la CI on obtient \(2=1/C\) d’où l’unique solution du problème de Cauchy: \(y(t)=\frac{2}{2t^2+1}\).

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

edo = sp.Eq(sp.diff(u(t),t),-2*t*(u(t))**2)
display(edo)

t0 = 0
u0 = 2
IC = sp.Eq(u(t0),u0)
display(IC)

sol = sp.dsolve(edo,u(t),ics={u(t0):u0})
print('Solution :')
display(sol)

sp.plot(sol.rhs,(t,0,3))

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

\(\displaystyle \frac{d}{d t} u{\left(t \right)} = - 2 t u^{2}{\left(t \right)}\)

\(\displaystyle u{\left(0 \right)} = 2\)

Solution :

\(\displaystyle u{\left(t \right)} = \frac{1}{t^{2} + \frac{1}{2}}\)

Correction : étude empirique de l’ordre de convergence

On calcule la solution approchée avec différentes valeurs de \(h_k=1/N_k\), à savoir \(N_k=2^3\), … \(2^{10}\). 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}|\) avec \(N_k=2^{k+1}\).

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, 2, 3
phi = lambda t,y : -2*t*y**2

i = 0
for alpha in ALPHA:
    alpha = float(alpha)
    HH = []
    err = []
    for k in range(8):
        N = 2**(k+3)  
        tt = np.linspace(t0,tfinal,N+1)
        HH.append( tt[1]-tt[0] )
        yy = sol_exacte(tt)
        uu = RKalpha(phi,tt,y0)
        err.append(np.linalg.norm(uu-yy, np.inf))
    
    i += 1
    plt.figure(i,figsize=(20,5))
    plt.loglog(HH,err, 'r-o')
    pente = np.polyfit(np.log(HH), np.log(err), 1)[0]
    plt.title(f'Alpha={alpha:g}, Ordre = {pente:.4f}')
    plt.xlabel('$h$')
    plt.ylabel('$e$')
    #legend() 
    plt.grid(True)