M64 - Data Science & Scientific Computing 6 - M62 Devoir Maison 2020

M62 Devoir Maison 2020

Auteur·rice

Gloria FACCANONI

Date de publication

20 mars 2024

1 M62 Devoir Maison 2020

1.1 Exercice : implémentation et comparaison de schémas

Le déplacement \(x(t)\) d’un système oscillant composé d’une masse et d’un ressort, soumis à une force de frottement proportionnelle à la vitesse, est décrit par l’équation différentielle du second ordre \( x''(t)+ 5x'(t)+6x(t)=0, \) avec \(x(0) = 1\) et \(x'(0) = 0\), pour \(t \in [0, 5]\).

  1. Écrire le système sous forme matricielle. Calculer ensuite la solution exacte avec sympy. Afficher \(t\mapsto x(t)\), \(t\mapsto x'(t)\) et \(x\mapsto x'(x)\).
  2. Calculer la solution approchée obtenue par la méthode d’Euler Progressif avec \(301\) points. Afficher \(t\mapsto x(t)\), \(t\mapsto x'(t)\) et \(x\mapsto x'(x)\) en comparant solution exacte et solution approchée.
  3. Même exercice avec la méthode d’Euler Regressif.
  4. Même exercice avec la fonction odeint du module scipy.
Code 
%reset -f

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})

import sympy as sym
sym.init_printing()

Q1 [1 point]
Écrire le système sous forme matricielle. Calculer ensuite la solution exacte avec sympy. Afficher \(t\mapsto x(t)\), \(t\mapsto x'(t)\) et \(x\mapsto x'(x)\).

Si on note \(y(t)=x(t)\) et \(z(t)=x'(t)\) on a \(y(0)=1\), \(z(0)=0\) et \( \begin{cases} y'(t)=z(t),\\ z'(t)=-6y(t)-5z(t) \end{cases} \quad\text{i.e.}\quad \begin{pmatrix}y\\z\end{pmatrix}'(t) {}= \begin{pmatrix}0&1\\-6&-5\end{pmatrix} \begin{pmatrix}y\\z\end{pmatrix}(t) \)

Code 
# VARIABLES GLOBALES
t0 = 0
tfinal = 5
N = 300

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

y0 = 1
z0 = 0

phi1 = lambda t,y,z : z
phi2 = lambda t,y,z : -6*y-5*z

# FONCTION UTILE

def affichage(tt, yy, zz, uu=None, ww=None, s=""):
    plt.figure(figsize=(17,7)) 
   
    plt.subplot(1,2,1)
    if uu is not None and ww is not None:
        plt.plot(tt,uu,'--',label=r'$u(t) \approx y(t)$',color="red")
        plt.plot(tt,ww,'--',label=r'$w(t) \approx z(t)$',color="blue")
    plt.plot(tt,yy,label='$y(t)$',color="red")
    plt.plot(tt,zz,label='$z(t)$',color="blue")
    plt.xlabel('t')
    plt.legend()
    plt.title(f'{</span>s<span class="sc">} - $y(t)$ et $z(t)$') 
    plt.grid()

    plt.subplot(1,2,2)
    if uu is not None and ww is not None:
        plt.plot(uu,ww,'--',label='Approchée')
    plt.plot(yy,zz,label='Exacte')
    plt.xlabel('y')
    plt.ylabel('z')
    plt.legend()
    if s!="":
        plt.title(f'{</span>s<span class="sc">} - z(y)')
    plt.grid()
    plt.axis('equal');

La solution exacte est \(\begin{align*} x(t)=y(t)&=3e^{-2t}−2e^{-3t},\\ x'(t)=z(t)&=-6e^{-2t}+6e^{-3t}. \end{align*}\) Vérifions-le avec le module sympy:

Code 
t = sym.Symbol('t')
y = sym.Function('y')
z = sym.Function('z')
edo1 = sym.Eq( sym.diff(y(t),t) , phi1(t,y(t),z(t)) )
edo2 = sym.Eq( sym.diff(z(t),t) , phi2(t,y(t),z(t)) )
display(edo1)
display(edo2)
solpar_1, solpar_2 = sym.dsolve([edo1,edo2],[y(t),z(t)], ics={y(t0):y0, z(t0):z0})
display(solpar_1)
display(solpar_2)

func_1 = sym.lambdify(t,solpar_1.rhs,'numpy')
func_2 = sym.lambdify(t,solpar_2.rhs,'numpy')

yy = func_1(tt)
zz = func_2(tt)

affichage(tt,yy,zz,s="Solution exacte")

\(\displaystyle \frac{d}{d t} y{\left(t \right)} = z{\left(t \right)}\)

\(\displaystyle \frac{d}{d t} z{\left(t \right)} = - 6 y{\left(t \right)} - 5 z{\left(t \right)}\)

\(\displaystyle y{\left(t \right)} = 3 e^{- 2 t} - 2 e^{- 3 t}\)

\(\displaystyle z{\left(t \right)} = - 6 e^{- 2 t} + 6 e^{- 3 t}\)

Q2 [2 points]
Calculer la solution approchée obtenue par la méthode d’Euler Progressif avec \(301\) points. Afficher \(t\mapsto x(t)\), \(t\mapsto x'(t)\) et \(x\mapsto x'(x)\) en comparant solution exacte et solution approchée.

On notera \(u_n\approx y_n=x(t_n)\) et \(w_n\approx z(t_n)\).

Euler explicite \( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+h\varphi_1(t_n,u_n,w_n),\\ w_{n+1}=w_n+h\varphi_2(t_n,u_n,w_n). \end{cases} \)

Code 
def euler_progressif(phi1,phi2,tt,y0,z0):
    h = tt[1]-tt[0]
    uu, ww = np.zeros_like(tt), np.zeros_like(tt)
    uu[0], ww[0] = y0, z0
    for n in range(len(tt)-1):
        uu[n+1] = uu[n] + h*phi1(tt[n],uu[n],ww[n])
        ww[n+1] = ww[n] + h*phi2(tt[n],uu[n],ww[n])
    return [uu,ww]

[uu, ww] = euler_progressif(phi1,phi2,tt,y0,z0)
affichage(tt,yy,zz,uu,ww,"Euler explicite");

Q3 [2 points]
Calculer la solution approchée obtenue par la méthode d’Euler Regressif avec \(301\) points. Afficher \(t\mapsto x(t)\), \(t\mapsto x'(t)\) et \(x\mapsto x'(x)\) en comparant solution exacte et solution approchée.

Euler implicite \( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+h\varphi_1(t_{n+1},u_{n+1},w_{n+1}),\\ w_{n+1}=w_n+h\varphi_2(t_{n+1},u_{n+1},w_{n+1}), \end{cases} \) qu’on peut rendre explicite par un calcul élementaire: \( w_{n+1} = w_n + h(-6u_{n+1}-5w_{n+1}) = w_n + h(-6(u_{n}+hw_{n+1})-5w_{n+1}) \quad\implies\quad (1+5h+6h^2)w_{n+1} = w_n-6hu_{n}\) ainsi \( \begin{cases} u_0=1,\\ w_0=0,\\ u_{n+1}=u_n+h\dfrac{w_n-6hu_{n}}{1+5h+6h^2},\\ w_{n+1}=\dfrac{w_n-6hu_{n}}{1+5h+6h^2}. \end{cases} \)

Code 
# VERSION explicite
# def euler_regressif(phi1,phi2,tt,y0,z0):
#   h = tt[1]-tt[0]
#   uu, ww = np.zeros(len(tt)), np.zeros(len(tt))
#   uu[0], ww[0] = y0, z0
#   for n in range(len(tt)-1):
#       uu[n+1] = uu[n] + h*(ww[n]-6*h*uu[n])/(1+5*h+6*h**2)
#       ww[n+1] = (ww[n]-6*h*uu[n])/(1+5*h+6*h**2)  
#   return [uu,ww]


#  VERSION AVEC fsolve
from scipy.optimize import fsolve
def euler_regressif(phi1,phi2,tt,y0,z0):
    h = tt[1]-tt[0]
    uu, ww = np.zeros_like(tt), np.zeros_like(tt)
    uu[0], ww[0] = y0, z0
    for n in range(len(tt)-1):
        systeme = lambda k : [ -k[0]+uu[n]+h*phi1(tt[n+1],k[0],k[1]) , 
                               -k[1]+ww[n]+h*phi2(tt[n+1],k[0],k[1]) ]
        uu[n+1], ww[n+1] = fsolve( systeme , (uu[n],ww[n]) ) 
    return [uu,ww]

[uu, ww] = euler_regressif(phi1,phi2,tt,y0,z0)
affichage(tt,yy,zz,uu,ww,"Euler implicite")

Q4 [1 point]
Calculer la solution approchée obtenue par la fonction odeint du module scipy. Afficher \(t\mapsto x(t)\), \(t\mapsto x'(t)\) et \(x\mapsto x'(x)\) en comparant solution exacte et solution approchée.

Code 
pphi = lambda t,yy : [ phi1(t,yy[0],yy[1]), phi2(t,yy[0],yy[1])  ]
yy0  = [y0,z0]

# OLD
from scipy.integrate import odeint
uu,ww = odeint(pphi,yy0,tt, tfirst=True).T

# NEW
# from scipy.integrate import solve_ivp
# uu,ww = solve_ivp(pphi,[t0,tfinal],yy0, t_eval=tt, 
#                 #   method='RK45'
#                   method='BDF'
#                   ).y


affichage(tt,yy,zz,uu,ww,"Scipy")

1.2 Exercice : étude d’un schéma RK

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

  1. Écrire le schéma sous la forme d’une suite définie par récurrence.
  2. Étudier théoriquement l’ordre de ce schéma.
  3. Étudier théoriquement la A-stabilité de ce schéma.
  4. Implémenter ce schéma et le tester sur le problème de Cauchy suivant avec \(N+1\) points et \(N=8\) : \( \begin{cases} y'(t)=-6y(t), &t\in[0;1]\\ y(0)=1 \end{cases} \)
Code 
%reset -f

import sympy 
sympy.init_printing()

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

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

Q5 [1 point] Écrire le schéma sous la forme d’une suite définie par récurrence.

Le schéma associé à cette matrice est semi-implicite (\(K_0\) dépend de lui même) et permet de calculer \(u_{n+1}\) à partir de \(u_n\) par la formule de récurrence \(\begin{cases} u_0 = y_0 \\ K_0 = \varphi\left(t_n+\frac{h}{3},u_n+\frac{h}{3}K_0\right)\\ K_1 = \varphi\left(t_{n+1},u_n+hK_0\right)\\ u_{n+1} = u_n + \frac{h}{4}\left(3K_0+K_1\right) & n=0,1,\dots N-1 \end{cases}\)

Q6 [2 points] Étudier théoriquement l’ordre du schéma.

Soit \(\omega\) l’ordre de la méthode.

  • C’est un schéma semi-implicite à \(2\) étages (\(s=2\)) donc \(\omega\le2s=4\)

  • \(\omega\ge1\) si \(\begin{cases} b_1+b_2=1 \\ c_1=a_{11}+a_{12} \\ c_2=a_{21}+a_{22} \end{cases}\)

  • \(\omega\ge2\) si, de plus, \(b_1c_1+c_2b_2=\frac{1}{2}\)

  • \(\omega\ge3\) si, de plus, \(\begin{cases} b_1c_1^2+b_2c_2^2=\frac{1}{3}\\ \left(b_1 a_{11} c_1 +b_1 a_{12} c_2\right) + \left(b_2 a_{21} c_1+b_2 a_{22} c_2\right) =\frac{1}{6} \end{cases}\)

  • \(\omega\ge4\) si, de plus,
    \(\begin{cases} b_1c_1^3+b_2c_2^3=\frac{1}{4}& \\ \left(b_1 c_1 a_{11} c_1+b_1 c_1 a_{12} c_2\right)+\left(b_2 c_2 a_{21} c_1+b_2 c_2 a_{22} c_2\right) =\frac{1}{8} \\ \left(b_1 a_{11} c_1^2 +b_1 a_{12} c_2^2\right) + \left(b_2 a_{21} c_1^2+b_2 a_{22} c_2^2\right)=\frac{1}{12} \\ b_1 a_{11}a_{11} c_1 + b_1 a_{11}a_{22} c_2 + b_1 a_{12}a_{21} c_1+b_1 a_{12}a_{22} c_2 + b_2 a_{21}a_{11} c_1 + b_2 a_{21}a_{22} c_2 + b_2 a_{22}a_{21} c_1+b_2 a_{22}a_{22} c_2 =\frac{1}{24} \end{cases}\)

Calculons donc toutes ces sommes (on pose directement \(a_{12}=a_{22}=0\) et \(a_{21}=1\)):

  • \(\begin{cases} b_1+b_2=\frac{3}{4}+\frac{1}{4}=1 \\ c_1=\frac{1}{3}=a_{11} \\ c_2=1=a_{21} \end{cases}\)

  • \(b_1c_1+c_2b_2=\frac{3}{4}\frac{1}{3}+1\frac{1}{4}=\frac{1}{2}\)

  • \(\begin{cases} b_1c_1^2+b_2c_2^2=\frac{3}{4}\frac{1}{3^2}+\frac{1}{4}1^2=\frac{1}{3}\\ \left(b_1 a_{11} c_1 \right) + \left(b_2 c_1\right) = \frac{3}{4}\frac{1}{3}\frac{1}{3}+ \frac{1}{4}\frac{1}{3}=\frac{1}{6} \end{cases}\)

  • \(b_1c_1^3+b_2c_2^3=\frac{3}{4}\frac{1}{3^3}+\frac{1}{4}1^3 =\frac{5}{18} \neq \frac{1}{4}\)

Code 
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**"))
    Eq_1 = ordre_1(s, A, b, c, j)
    Eqs.append(Eq_1)
    
    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])

    display(Markdown("**On doit ajouter 2 conditions pour être d'ordre 3**"))
    Eq_3 = ordre_3(s, A, b, c, j)  
    Eqs.append(Eq_3)  
    
    display(Markdown("**On doit ajouter 4 conditions pour être d'ordre 4**"))
    Eq_4 = ordre_4(s, A, b, c, j)
    Eqs.append(Eq_4)

    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))
    Eq_1 = [sympy.Eq( sum(b).simplify(), 1 )]
    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 ))
        Eq_1.append(sympy.Eq( somma, c[i] ))
    return Eq_1    


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


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


def ordre_4(s, A, b, c, j):
    texte = '\sum_{j=1}^s b_j c_j^3='
    texte += sympy.latex( sum([b[i]*c[i]**3 for i in range(s)]).simplify() )
    texte += r"\text{ doit être égale à }\frac{1}{4}"
    display(Math(texte))
    Eq_4_1 = sympy.Eq( sum([b[i]*c[i]**3 for i in range(s)]).simplify(), sympy.Rational(1,4) )
    texte = r'\sum_{i,j=1}^s b_i c_i a_{ij} c_j='
    somma = sum([b[i]*c[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify()
    texte = texte + sympy.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{8}"
    display(Math(texte))
    Eq_4_2 = sympy.Eq( sum([b[i]*c[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify(), sympy.Rational(1,8) )
    texte = r'\sum_{i,j=1}^s b_i a_{ij} c_j^2='
    somma = sum([b[i]*A[i,j]*c[j]**2 for j in range(s) for i in range(s)]).simplify()
    texte = texte + sympy.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{12}"
    display(Math(texte))
    Eq_4_3 = sympy.Eq( sum([b[i]*A[i,j]*c[j]**2 for j in range(s) for i in range(s)]).simplify(), sympy.Rational(1,12) )
    texte = r'\sum_{i,j,k=1}^s b_i a_{ij} a_{jk} c_k='
    somma = sum([b[i]*A[i,j]*A[j,k]*c[k] for k in range(s) for j in range(s) for i in range(s)]).simplify()
    texte = texte + sympy.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{24}"
    display(Math(texte))
    Eq_4_4 = sympy.Eq( sum([b[i]*A[i,j]*A[j,k]*c[k] for k in range(s) for j in range(s) for i in range(s)]).simplify(), sympy.Rational(1,24) )
    return [Eq_4_1, Eq_4_2, Eq_4_3, Eq_4_4]



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

s = 2

# paramètres dans le cas générale
# b_0, b_1  = sympy.symbols('b_0 b_1', real=True)
# c_0, c_1 = sympy.symbols('c_0 c_1', real=True)
# a_00, a_01, a_10, a_11 = sympy.symbols('a_00 a_01 a_10 a_11', real=True)

# matrice dans le cas général
# b = [b_0, b_1]
# c = [c_0, c_1]
# A = [[a_00, a_01], [a_10, a_11]]

# matrice dans notre cas particulier
c = [sympy.Rational(1,3),1]
b = [sympy.Rational(3,4),sympy.Rational(1,4)]
A = [[sympy.Rational(1,3),0],[1,0]]

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

Matrice de Butcher

\(\displaystyle \left[\begin{matrix}\frac{1}{3} & \frac{1}{3} & 0\\1 & 1 & 0\\ & \frac{3}{4} & \frac{1}{4}\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}=\frac{1}{3}\text{ doit être égale à }\frac{1}{3}\)

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

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

\(\displaystyle \sum_{j=1}^s b_j c_j=\frac{1}{2}\text{ doit être égale à }\frac{1}{2}\)

On doit ajouter 2 conditions pour être d’ordre 3

\(\displaystyle \sum_{j=1}^s b_j c_j^2=\frac{1}{3}\text{ doit être égale à }\frac{1}{3}\)

\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j=\frac{1}{6}\text{ doit être égale à }\frac{1}{6}\)

On doit ajouter 4 conditions pour être d’ordre 4

\(\displaystyle \sum_{j=1}^s b_j c_j^3=\frac{5}{18}\text{ doit être égale à }\frac{1}{4}\)

\(\displaystyle \sum_{i,j=1}^s b_i c_i a_{ij} c_j=\frac{1}{9}\text{ doit être égale à }\frac{1}{8}\)

\(\displaystyle \sum_{i,j=1}^s b_i a_{ij} c_j^2=\frac{1}{18}\text{ doit être égale à }\frac{1}{12}\)

\(\displaystyle \sum_{i,j,k=1}^s b_i a_{ij} a_{jk} c_k=\frac{1}{18}\text{ doit être égale à }\frac{1}{24}\)

D’après ces résultats, le schéma est d’ordre \(3\).

Q7 [2 points] Étudier théoriquement la A-stabilité du schéma.

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_0 = -\beta\left(u_n+\frac{h}{3}K_0\right)\\ K_1 = -\beta\left(u_n+hK_0\right)\\ u_{n+1} = u_n + \frac{h}{4}\left(3K_0+K_1\right) & n=0,1,\dots N-1 \end{cases}\)

L’équation \(K_0 = -\beta\left(u_n+\frac{h}{3}K_0\right)\) donne \(K_0=\frac{-3\beta}{\beta h+3}u_n\) ainsi \(u_{n+1} = u_n + \frac{h}{4}\left(3K_0+K_1\right) = u_n + \frac{h}{4}\left(3K_0-\beta u_n-\beta h K_0\right) = \left(1-\frac{\beta h}{4}\right)u_n + \frac{h(3-\beta h)}{4}K_0\) et enfin \( u_{n+1} = \left(1-\frac{\beta h}{4}\right)u_n + \frac{h(3-\beta h)}{4}K_0 = \left(1-\frac{\beta h}{4}\right)u_n + \frac{h(3-\beta h)}{4}\frac{-3\beta}{\beta h+3}u_n = \left(1-\frac{\beta h}{4} - \frac{3\beta h(3-\beta h)}{4(3 + \beta h)}\right)u_n % u_n + \frac{h}{4}\left( \frac{-9\beta}{\beta h+3}+\beta\frac{2\beta h-3}{\beta h+3}\right)u_ n = \frac{(\beta h)^2-4\beta h+6}{2(\beta h+3)}u_n \)

Vérifions ce calcul:

Code 
# pour ne pas effacer l'affectation de "h", ici on vas l'appeler "dt" mais on affiche "h"
dt = sympy.Symbol('h')
u_np1 = sympy.Symbol('u_{n+1}')
beta = sympy.Symbol(r'\beta')
sympy.var('u_n,x')

K0 = sympy.solve(-x-beta*(u_n+dt/3*x),x)[0]
display(Math('K_0='+sympy.latex(K0)))

K1 = -beta*(u_n+dt*K0).factor()
display(Math('K_1='+sympy.latex(K1)))

u_np1 = (u_n+dt/4*(3*K0+K1)).factor()
display(Math('u_{n+1}='+sympy.latex(u_np1)))

\(\displaystyle K_0=- \frac{3 \beta u_{n}}{\beta h + 3}\)

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

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

On note \(x=\beta h\) et on étudie la fonction \(q\colon \mathbb{R}^+\to\mathbb{R}\) définie par \(q(x)=\dfrac{x^2-4x+6}{2(x+3)}=\dfrac{1}{2}\left(x-7+\dfrac{27}{x+3}\right)\). Le schéma est A-stable ssi \(|q(x)|<1\). D’après les calculs ci-dessous on conclut que le schéma est A-stable ssi \(h<\dfrac{6}{\beta}\).

Code 
sympy.var('x',nonnegative=True)

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

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("Le signe de q' est le signe du numérateur qui est une parabole. On cherche où elle s'annulle pour x>0")
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)+"\\approx"+sympy.latex(sympy.N(qmin))))
print("Conclusion: q(x)>0 pour tout x>0. Vérifions ce calcul:")
print("q(x)>0")
display(sympy.solve(q>0))
print("On resout maintenant q(x)=1:")
display(sympy.solve(q-1))
print('On trouve q(0)=q(6)=1')
print("Conclusion: q(x)<1 pour x<6. Vérifions ce calcul:")
print("q(x)<1")
display(sympy.solve(q<1))


sympy.plot(q,1,-1,xlim=[-1,15],ylim=[-1.5,1.5]);

\(\displaystyle q(x)=\frac{x^{2} - 4 x + 6}{2 x + 6}\)

\(\displaystyle q(0)=1\)

\(\displaystyle \lim_{x \to \infty}\left(\frac{x^{2} - 4 x + 6}{2 x + 6}\right)=\infty\)

On cherche les points stationnaires dans R^+
Le signe de q' est le signe du numérateur qui est une parabole. On cherche où elle s'annulle pour x>0
Conclusion: q(x)>0 pour tout x>0. Vérifions ce calcul:
q(x)>0
On resout maintenant q(x)=1:
On trouve q(0)=q(6)=1
Conclusion: q(x)<1 pour x<6. Vérifions ce calcul:
q(x)<1

\(\displaystyle q'(x)=\frac{x^{2} + 6 x - 18}{2 \left(x + 3\right)^{2}}\)

\(\displaystyle q'(x)=0 \iff x\in\left[ -3 + 3 \sqrt{3}\right]\)

\(\displaystyle x=-3 + 3 \sqrt{3}\text{ est un minimum et}\)

\(\displaystyle q(-3 + 3 \sqrt{3})=\frac{\sqrt{3} \left(- 12 \sqrt{3} + \left(-3 + 3 \sqrt{3}\right)^{2} + 18\right)}{18}\approx0.196152422706632\)

\(\displaystyle \text{True}\)

\(\displaystyle \left[ 0, \ 6\right]\)

\(\displaystyle 0 < x \wedge x < 6\)

Q8 [2 points] Implémenter le schéma et approcher la solution du problème de Cauchy suivant avec \(N+1\) points et \(N=8\). \( \begin{cases} y'(t)=-6y(t), &t\in[0;1]\\ y(0)=1 \end{cases} \)

Dans chaque point \(t_i\), il faut approcher \((K_0)_i\) en resolvant une équation. Si on utilise la fonction fsolve du module scipy.optimize, il faut initialiser fsolve avec une approximation de \((K_0)_i\). On choisira donc \(\varphi(t_i,u_i)\):

Code 
from scipy.optimize import fsolve

# version implicite
def RK(phi,tt):
    h = tt[1]-tt[0] 
    uu =np.zeros(len(tt))
    uu[0] = y0
    for n in range(len(tt)-1):
        k0 = fsolve( lambda x : -x+phi( tt[n]+h/3, uu[n]+h/3*x ) , phi(tt[n],uu[n]) )[0]
        k1 = phi( tt[n+1], uu[n]+h*k0 )
        uu[n+1] = uu[n] + h/4*(3*k0+k1 )
    return uu

# version explicite car phi = - beta h
# def RK(beta,tt):
#     h = tt[1]-tt[0]
#     x = beta*h
#     q = lambda x : (x**2-4*x+6)/(2*x+6)
#     uu = [y0*q(x)**n for n in range(len(tt))]
#     return uu




t0, y0, tfinal = 0, 1, 3
beta = 6
sol_exacte = lambda t : np.exp(-beta*t)
phi = lambda t,y : -beta*y

print('A-stable ssi h <',6/beta)
plt.figure(figsize=(10,5))
N = 8
tt = np.linspace(t0,tfinal,N+1)
h = tt[1]-tt[0]
yy = sol_exacte(tt)    
uu = RK(phi,tt)  # version avec fsolve
# uu = RK(beta,tt) # version explicite
plt.plot(tt,yy,'b-',label=("Exacte"))
plt.plot(tt,uu,'ro',label=("RK"))
plt.title(r' $N$={}, $h$={}'.format(N,h))
plt.xlabel('$t$')
plt.ylabel('$u$')
plt.legend() 
plt.grid(True)
A-stable ssi h < 1.0

1.3 Exercice : étude d’un schéma Predictor-Corrector multipas

On notera \(\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)\). Soit les méthodes multipas

\( \begin{aligned} u_{n+1}&=u_{n-3}+\frac{4h}{3}\left(2\varphi_{n}-\varphi_{n-1}+2\varphi_{n-2}\right)&\text{[P]} \\ u_{n+1}&=u_{n-1}+\frac{h}{3}\left(\varphi_{n+1}+4\varphi_{n}+\varphi_{n-1}\right)&\text{[C]} \end{aligned} \)

  1. Étudier consistance, ordre et zéro-stabilité de ces deux méthodes.
  2. Calculer empiriquement l’ordre de convergence des méthodes [P] et [C] ainsi que de la méthode Predictor-Corrector associée à ces deux méthodes sur le problème de Cauchy

\( \begin{cases} y'(t) = 1+\big(t-y(t)\big)^2, &\forall t \in I=[2,3],\\ y(2) = 1. \end{cases} \)

Code 
%reset -f

import numpy as np
import matplotlib.pyplot as plt
plt.rcParams.update({<span class="st">'font.size'</span>: <span class="dv">16</span>})

import sympy as sym
sym.init_printing()

from IPython.display import display, Math, Markdown



def xi(i, q, aa=None, bb=None, bm1=None):
    if aa is None:
        aa = sym.symbols(f'a_0:{</span>q<span class="sc">}')
    if bb is None:
        bb = sym.symbols(f'b_0:{</span>q<span class="sc">}')
    if bm1 is None:
        bm1 = sym.symbols('b_{-1}')
    p = len(aa)
    sa = sum((-j)**i * aa[j] for j in range(p))
    sb = bm1 + sum((-j)**(i-1) * bb[j] for j in range(1, p))
    # ??? si j=0 et i=0, il calcule (-j)**i =1 et on a bien a_0 dans la somme
    # mais si j=0 et i=1 il ne calcule pas b_0 et il faut l'ajouter à la main
    if i == 1:
        sb += bb[0]
    return (sa) + (i * sb)

def afficahge_xi(q, aa, bb, bm1):
    p = q-1
    omega = q + 2 if q%2==0 else q + 1 # i = 0, ..., omega
    display(Markdown(f"Pour une méthode à q={</span>q<span class="sc">} pas, on affiche $\\xi(i)$ pour i = 0, ..., {</span>q<span class="op">+</span><span class="dv">2</span><span class="sc">}"))
    for i in range(q+3):
        display(Markdown(f"$\\xi({</span>i<span class="sc">}) = {</span>sym<span class="sc">.</span>latex(xi(i, q, <span class="va">None</span>, <span class="va">None</span>, <span class="va">None</span>))<span class="sc">} = {</span>sym<span class="sc">.</span>latex(xi(i, q, aa, bb, bm1))<span class="sc">}$"))

Q9 [3 point] Étudier consistance, ordre et zéro-stabilité de la méthode P.

[P] est une méthode à \(q=4\) pas explicite :

  • \(p=3\)
  • \(b_{-1}=0\)
  • \(a_0=a_1=a_2=0\) et \(a_3=1\)
  • \(b_0=\frac{8}{3}\), \(b_1=\frac{-4}{3}\), \(b_2=\frac{8}{3}\) et \(b_3=0\)
  • La méthode est consistante d’ordre au plus \(q\) car explicite.
  • Elle est d’ordre \(\omega=4\) car \( \begin{cases} \xi(0)=\displaystyle\sum_{j=0}^p a_j= a_{0} + a_{1} + a_{2} + a_{3} = 1 \\ \xi(1)=\displaystyle\sum_{j=0}^p (-j)^{1}a_j+1\sum_{j=-1}^p (-j)^{0}b_j = - a_{1} - 2 a_{2} - 3 a_{3} + b_{0} + b_{1} + b_{2} + b_{3} + b_{-1}=1 \\ \xi(2)=\displaystyle\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j = - a_{1} - 2 a_{2} - 3 a_{3} + b_{0} + b_{1} + b_{2} + b_{3} + b_{-1}=1 \\ \xi(3)=\displaystyle\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j = - a_{1} - 8 a_{2} - 27 a_{3} + 3 \left(b_{1} + 4 b_{2} + 9 b_{3} + b_{-1}\right)=1 \\ \xi(4)=\displaystyle\sum_{j=0}^p (-j)^{4}a_j+4\sum_{j=-1}^p (-j)^{3}b_j = a_{1} + 16 a_{2} + 81 a_{3} - 4 \left(b_{1} + 8 b_{2} + 27 b_{3} - b_{-1}\right)=1 \end{cases} \) Rappel: \((-0)^0=1\) dans les formules précédentes.

Vérifions ces calculs ci-dessous:

Code 
# Paramètres
q = 3

a0 = 0; a1 = 0; a2 = 0; a3 = 1
b0 = sym.Rational(8,3); b1 = sym.Rational(-4,3); b2 = sym.Rational(8,3); b3 = 0
bm1 = 0

aa = [a0,a1,a2,a3]
bb = [b0,b1,b2,b3]
afficahge_xi(q, aa, bb, bm1)

Pour une méthode à q=3 pas, on affiche \(\xi(i)\) pour i = 0, …, 5

\(\xi(0) = a_{0} + a_{1} + a_{2} = 1\)

\(\xi(1) = - a_{1} - 2 a_{2} + b_{0} + b_{1} + b_{2} + b_{-1} = 1\)

\(\xi(2) = a_{1} + 4 a_{2} - 2 b_{1} - 4 b_{2} + 2 b_{-1} = 1\)

\(\xi(3) = - a_{1} - 8 a_{2} + 3 b_{1} + 12 b_{2} + 3 b_{-1} = 1\)

\(\xi(4) = a_{1} + 16 a_{2} - 4 b_{1} - 32 b_{2} + 4 b_{-1} = 1\)

\(\xi(5) = - a_{1} - 32 a_{2} + 5 b_{1} + 80 b_{2} + 5 b_{-1} = - \frac{109}{3}\)

  • Le premier polynôme caractéristique est \( \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^4-a_0r^3-a_1r^2-a_2r-a_3=r^4-1 \) dont les racines sont \( r_0=1,\ r_1=-1,\ r_2=i,\ r_3=-i \) La méthode est zéro-stable ssi \( \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ %|r_j|=1 \quad\implies \text{ $r_j$ a multiplicité 1, c'est-à-dire }\varrho'(r_j)\neq0, r_j\text{ a multiplicité}>1\quad\implies \quad |r_j|<1, \end{cases} \) ce qui est bien vérifié.

  • La méthode est convergente car consistante et zéro-stable.

Q10 [3 point] Étudier consistance, ordre et zéro-stabilité de la méthode C.

[C] est une méthode à \(q=2\) pas implicite :

  • \(p=1\)
  • \(b_{-1}=\frac{1}{3}\)
  • \(a_0=0\) et \(a_1=1\)
  • \(b_0=\frac{4}{3}\) et \(b_1=\frac{1}{3}\)
  • La méthode est consistante d’ordre au plus \(q+2=4\) car implicite avec un nombre pair de pas.
  • Elle est d’ordre \(\omega=4\) car

\( \begin{cases} \xi(0)=\displaystyle\sum_{j=0}^p a_j = a_{0} + a_{1} = 1, \\ \xi(1)=\displaystyle\sum_{j=0}^p (-j)^{1}a_j+1\sum_{j=-1}^p (-j)^{0}b_j = - a_{1} + b_{0} + b_{1} + b_{-1} =1, \\ \xi(2)=\displaystyle\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j = a_{1} - 2 b_{1} + 2 b_{-1} = 1, \\ \xi(3)=\displaystyle\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j = - a_{1} + 3 b_{1} + 3 b_{-1} = 1, \\ \xi(4)=\displaystyle\sum_{j=0}^p (-j)^{4}a_j+4\sum_{j=-1}^p (-j)^{3}b_j = a_{1} - 4 b_{1} + 4 b_{-1} = 1. \end{cases} \)

Rappel: \((-0)^0=1\) dans les formules précédentes.

Vérifions ces calculs ci-dessous:

Code 
# Paramètres
q = 2

a0 = 0; a1 = 1
b0 = sym.Rational(4,3); b1 = sym.Rational(1,3)
bm1 = sym.Rational(1,3)

aa = [a0,a1]
bb = [b0,b1]
afficahge_xi(q, aa, bb, bm1)

Pour une méthode à q=2 pas, on affiche \(\xi(i)\) pour i = 0, …, 4

\(\xi(0) = a_{0} + a_{1} = 1\)

\(\xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1} = 1\)

\(\xi(2) = a_{1} - 2 b_{1} + 2 b_{-1} = 1\)

\(\xi(3) = - a_{1} + 3 b_{1} + 3 b_{-1} = 1\)

\(\xi(4) = a_{1} - 4 b_{1} + 4 b_{-1} = 1\)

  • Le premier polynôme caractéristique est \( \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^2-a_0r-a_1=r^2-1 \) dont les racines sont \( r_0=1,\ r_1=-1 \) La méthode est zéro-stable ssi \( \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ r_j\text{ a multiplicité}>1\quad\implies \quad |r_j|<1, \end{cases} \) ce qui est bien vérifié.

  • La méthode est convergente car consistante et zéro-stable.

Q11 [2 point] Calculer empiriquement l’ordre de convergence des méthodes P et C ainsi que de la méthode Predictor-Corrector associée aux deux méthodes P et C sur le problème de Cauchy \( \begin{cases} y'(t) = 1+\big(t-y(t)\big)^2, &\forall t \in I=[2,3],\\ y(2) = 1. \end{cases} \)

Il est d’ordre \(4\) car le prédictor est d’ordre \(4\), le corrector d’ordre \(4\) et l’on a \( \omega_{[PC]} = \min \left\{ \omega_{[C]} , \omega_{[P]}+1 \right\} \)

On définit la solution exacte pour estimer les erreurs et on calcule la solution approchée pour différentes valeurs de \(N\) pour les trois schémas:

Code 
t0     = 2
tfinal = 3
y0     = 1

phi = lambda t,y : 1+(t-y)**2


t = sym.Symbol('t')
y = sym.Function('y')
edo = sym.Eq( sym.diff(y(t),t) , phi(t,y(t)) )
display(edo)
solpar = sym.dsolve( edo, y(t), ics={y(t0):y0}) 
display(solpar)

sol_exacte = sym.lambdify(t,solpar.rhs,'numpy')
# sol_exacte = lambda t : t+1/(1-t)

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


from scipy.optimize import fsolve
def multipas_C(phi, tt):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:2] = sol_exacte(tt[:2])
    for n in range(1,len(tt) - 1):
        k0 = phi(tt[n],uu[n])
        k1 = phi(tt[n-1],uu[n-1])
        temp = fsolve ( lambda x : -x + uu[n-1] + h/3*( phi(tt[n+1],x) + 4*k0 + k1 ), uu[n])
        uu[n+1] = temp[0]
    return uu

def multipas_PC(phi, tt):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:4] = sol_exacte(tt[:4])
    for n in range(3,len(tt) - 1):
        k0 = phi(tt[n],uu[n])
        k1 = phi(tt[n-1],uu[n-1])
        k2 = phi(tt[n-2],uu[n-2])
        u_tilde = uu[n-3]+ 4*h/3*( 2*k0 - k1 + 2*k2 ) 
        uu[n+1] = uu[n-1]+h/3* ( phi(tt[n+1],u_tilde) + 4*k0 + k1 ) 
    return uu


H = []
err_p = []
err_c = []
err_pc = []
N = 10
for k in range(7):
    N+=20
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu_p = multipas_P(phi, tt)
    uu_c = multipas_C(phi, tt)
    uu_pc = multipas_PC(phi, tt)
    H.append(h)
    err_p.append(max([abs(uu_p[i] - yy[i]) for i in range(len(yy))]))
    err_c.append(max([abs(uu_c[i] - yy[i]) for i in range(len(yy))]))
    err_pc.append(max([abs(uu_pc[i] - yy[i]) for i in range(len(yy))]))

plt.figure(figsize=(8,5))
plt.plot(np.log(H), np.log(err_p), 'r-o', label=f'Multipas  [P], ordre = {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err_p ), <span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:g}')
plt.plot(np.log(H), np.log(err_c), 'b-o', label=f'Multipas  [C], ordre = {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err_c ), <span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:g}')
plt.plot(np.log(H), np.log(err_pc),'c-o', label=f'Multipas [PC], ordre = {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err_pc), <span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:g}')
plt.xlabel('$\ln(h)$')
plt.ylabel('$\ln(e)$')
plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
plt.grid(True);

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

\(\displaystyle y{\left(t \right)} = \frac{t^{2} - t - 1}{t - 1}\)

On a bien

  • ordre du corrector \(\omega_{[C]}=4\)
  • ordre du predictor \(\omega_{[P]}=4\)

et donc

  • ordre du predictor-corrector \(\omega_{[PC]}=\min\{\omega_{[C]},\omega_{[P]}+1\}=4\)
Retour au sommet