TD3 - Exercices sur les schémas de Runge-Kutta

Auteur·rice

Gloria FACCANONI

Date de publication

17 mars 2025

1 Aide-mémoire

Définitions

Soit \(s\ge1\) un entier. On se donne les \(2s+s^2\) coefficients suivants :

\(\mathbf{c}=(c_1,c_2,\dots,c_s)\),
\(\mathbf{b}=(b_1,b_2,\dots,b_s)\).
\(\mathbb{A}=(a_{ij})_{1\le i,j\le s}\),

qu’on regroupe dans un tableau de Butcher :

\(\begin{array}{c|c} \mathbf{c} & \mathbb{A}\\ \hline &\mathbf{b}^T \end{array}.\)

La méthode de Runge-Kutta à \(s\) étages associée à ce tableau est définie par :

\(\begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_{n}+h \displaystyle\sum_{i=1}^sb_iK_i& n=0,1,\dots N-1\\ K_i=\displaystyle\varphi\left( t_n+hc_i,u_n+h\sum_{j=1}^{s}a_{ij}K_j \right) & i=1,\dots s. \end{cases}\)

Méthode implicite : s’il existe un terme \(a_{ij}\neq0\) pour \(j> i\), alors la méthode est dite implicite. Pour calculer les \(K_i\), il faut résoudre un système d’équations non linéaires.
Méthode semi-implicite : si \(a_{ij}=0\) pour \(j>i\) (la matrice \(\mathbb{A}\) est triangulaire inférieure mais au moins un terme de la diagonale est non nul), alors la méthode est dite semi-implicite. Pour calculer les \(K_i\), il suffit de résoudre l’une après l’autre les équations non linéaires (et non pas un système).
Méthode explicite : si \(a_{ij}=0\) pour \(j\ge i\) (la matrice \(\mathbb{A}\) est triangulaire inférieure stricte), alors la méthode est dite explicite. Le calcul des \(K_i\) est alors direct.

Théorème de Lax-Richtmyer

Soit \(\omega\ge1\). Le schéma est convergent (d’ordre \(\omega\)) ssi il est zéro-stable et consistant (d’ordre \(\omega\ge1\)).

Zéro-stabilité

Les méthodes de Runge-Kutta sont zéro-stables.

Consistance (et ordre de consistance \(\omega\))

Une méthode de Runge-Kutta est consistante d’ordre \(\omega\ge1\) si elle satisfait les conditions suivantes :

Ordre minimum	Conditions à satisfaire
\(\omega \geq 1\)	\(\begin{cases} \displaystyle 1 = \sum_{j=1}^s b_j \\ \displaystyle c_i = \sum_{j=1}^s a_{ij}, \quad i=1,\dots,s \end{cases}\)
\(\omega \geq 2\)	Les conditions pour \(\omega \geq 1\) et \(\displaystyle\sum_{j=1}^s b_j c_j = \frac{1}{2}\)
\(\omega \geq 3\)	Les conditions pour \(\omega \geq 2\) et \(\begin{cases} \displaystyle\sum_{j=1}^s b_j c_j^2 = \frac{1}{3} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i a_{ij} c_j = \frac{1}{6} \end{cases}\)
\(\omega \geq 4\)	Les conditions pour \(\omega \geq 3\) et \(\begin{cases} \displaystyle\sum_{j=1}^s b_j c_j^3 = \frac{1}{4} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i c_i a_{ij} c_j = \frac{1}{8} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i a_{ij} c_j^2 = \frac{1}{12} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s\sum_{k=1}^s b_i a_{ij}a_{jk} c_k = \frac{1}{24} \end{cases}\)

Barrière de Butcher

L’ordre \(\omega\ge1\) d’une méthode de Runge-Kutta à \(s\ge1\) étapes convergente satisfait les inégalités suivantes :

\( \omega \leq \begin{cases} 2s, & \text{si la méthode est implicite}, \\ s, & \text{si la méthode est explicite et } s < 5, \\ s - 1, & \text{si la méthode est explicite et } s \geq 5. \end{cases} \)

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)=1 \end{cases}\)

Sa solution est \(y(t)=e^{-\beta t}\) et l’on a \(y(t)\xrightarrow[t\to+\infty]{}0.\)

Une méthode de Runge-Kutta à \(s\ge1\) étages pour \(y'(t)=-\beta y(t)\) s’écrit: \( \begin{cases} u_0=y_0, \\ u_{n+1}=u_{n}+h \displaystyle\sum_{i=1}^sb_iK_i& n=0,1,\dots N-1\\ K_i=\displaystyle -\beta \left( u_n+h\sum_{j=1}^{s}a_{ij}K_j \right) & i=1,\dots s. \end{cases} \)

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

Le calcul des conditions de consistance est un peu long et fastidieux dès que \(s\) est grand ou \(\omega\) est élevé. On peut utiliser un logiciel de calcul formel pour les obtenir, par exemple le module sympy.

On doit simplement indiquer le nombre d’étages \(s\) et le type de méthode (implicite, explicite, semi-implicite) pour obtenir les conditions de consistance. Indiquer le type de méthode permet de simplifier les calculs car de nombreux coefficients de la matrice sont nuls.

Code

# =============================================================================
# Calcul des conditions de consistance pour un schéma RK
# On doit indiquer le nombre d'étages s et le type de schéma implicite, semi-implicite ou explicite
# =============================================================================
s = 2
type = "semi-implicite" # implicite, semi-implicite, explicite






# =============================================================================
# Fonctions pour calculer les conditions de consistance pour un schéma RK.
#
# La fonction ordre_RK calcule les conditions pour les ordres 1, 2, 3 et 4 (même si théoriquement cela n'est pas possible) 
# - affiche la matrice de Butcher
# - affiche les conditions de consistance obtenues pour chaque ordre par la fonction ordre_i
# - renvoie une liste d'équations pour une éventuelle résolution
#
# Chaque fonction ordre_i calcule les conditions pour être d'ordre i
# - peut afficher les conditions de consistance [actuellement en commentaire]
# - renvoie la liste d'équations
#
# =============================================================================

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




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

    j = sympy.symbols('j')
    
    if A is None:
        if type == "implicite":
            A = sympy.MatrixSymbol('a', s, s)
        elif type=="semi-implicite":
            A = sympy.Matrix(s, s, lambda i,j: sympy.Symbol('a{}{}'.format(i, j)) if i >= j else 0)
        elif type=="explicite":
            A = sympy.Matrix(s, s, lambda i,j: sympy.Symbol('a{}{}'.format(i, j)) if i > j else 0)
            
    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(*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(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(*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)
    display(*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 = r"\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 = r'\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 = r'\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 = r'\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]

    
# ===================================================================================    
# Affichage des conditions de consistance : attention aux indices qui vont de 0 à s-1
# ===================================================================================
display(Markdown(f"### Conditions pour un schéma {</span><span class="bu">type</span><span class="sc">} avec s = {</span>s<span class="sc">} étages "))
Eqs = ordre_RK(s, type)

1.0.1 Conditions pour un schéma semi-implicite avec s = 2 étages

Matrice de Butcher

\(\displaystyle \left[\begin{matrix}c_{0} & a_{00} & 0\\c_{1} & a_{10} & a_{11}\\ & b_{0} & b_{1}\end{matrix}\right]\)

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

\(\displaystyle b_{0} + b_{1} = 1\)

\(\displaystyle a_{00} = c_{0}\)

\(\displaystyle a_{10} + a_{11} = c_{1}\)

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

\(\displaystyle b_{0} c_{0} + b_{1} c_{1} = \frac{1}{2}\)

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

\(\displaystyle b_{0} c_{0}^{2} + b_{1} c_{1}^{2} = \frac{1}{3}\)

\(\displaystyle a_{00} b_{0} c_{0} + a_{10} b_{1} c_{0} + a_{11} b_{1} c_{1} = \frac{1}{6}\)

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

\(\displaystyle b_{0} c_{0}^{3} + b_{1} c_{1}^{3} = \frac{1}{4}\)

\(\displaystyle a_{00} b_{0} c_{0}^{2} + a_{10} b_{1} c_{0} c_{1} + a_{11} b_{1} c_{1}^{2} = \frac{1}{8}\)

\(\displaystyle a_{00} b_{0} c_{0}^{2} + a_{10} b_{1} c_{0}^{2} + a_{11} b_{1} c_{1}^{2} = \frac{1}{12}\)

\(\displaystyle a_{00}^{2} b_{0} c_{0} + a_{00} a_{10} b_{1} c_{0} + a_{10} a_{11} b_{1} c_{0} + a_{11}^{2} b_{1} c_{1} = \frac{1}{24}\)

2 Exercices

2.0.1 🎨 Exercice : schéma de Ralston (\(s=2\), explicite)

On considère la méthode RK donnée par la matri suivante : \( \begin{array}{c|cc} 0 & 0 & 0 \\ \frac{2}{3} & \frac{2}{3} & 0 \\ \hline & \frac{1}{4} & \frac{3}{4} \end{array} \)

Écrire le schéma correspondant.
Déterminer l’ordre de convergence.
Déterminer sous quelle condition sur \(h\) ce schéma est A-stable.
Tester empiriquement l’ordre de convergence de ce schéma sur le problème de Cauchy \(y'(t) = t(1-t^2)-y(t)\) et \(y(0) = 1\) pour \(t\in[0,2]\).

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

Correction point 1 : schéma

\( \begin{cases} u_0 = y_0 \\ k_0 = \varphi\Big(t_n,u_n\Big) \\ k_1 = \varphi\Big(t_n+\frac{2}{3}h,u_n+\frac{2}{3}hk_0\Big) \\ u_{n+1}= u_n + h\big(\frac{1}{4}k_0+\frac{3}{4}k_1\big) \end{cases} \)

Il s’agit d’un schéma à \(s=2\) étages explicite.

Correction point 2 : ordre de convergence

D’après la barrière de Butcher, un schéma RK à \(s=2\) étages explicite ne peut pas être d’ordre supérieur à \(s=2\).

On peut bien-sûr écrire les formules à la main, mais on peut aussi utiliser sympy comme suit :

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 = r"\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 = r'\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 = r'\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 = r'\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
b = [ sympy.Rational(1,4), sympy.Rational(3,4) ]
c = [ 0, sympy.Rational(2,3) ]
A = [ [ 0, 0] , [ sympy.Rational(2,3),  0] ]

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

Matrice de Butcher

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

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

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

En conclusion, le schéma est d’ordre 2.

Correction point 3 : A-stabilité

On applique le schéma au problème de Cauchy \(y'(t) = -\beta y(t)\) et \(y(0) = 1\) avec \(\beta>0\) dont la solution est \(y(t)=e^{-\beta t}\).

On obtient alors la relation de récurrence suivante :

Code

u_n = sympy.symbols(r'u_n')
u_np1 = sympy.symbols(r'u_{n+1}')

dt = sympy.symbols(r'h', positive=True)
x = sympy.symbols(r'x', positive=True)

k_0 = sympy.symbols(r'k_0')
k_1 = sympy.symbols(r'k_1')
kk = [k_0, k_1]

beta = sympy.symbols(r'\beta', positive=True)
phi = lambda y: -beta*y

Eqs = [ sympy.Eq( kk[i] ,  phi(u_n + dt*sum([kk[j]*A[i][j] for j in range(s)])) ) for i in range(s) ]
for eq in Eqs:
    display(eq)
    
sol = sympy.solve( Eqs, kk ) # c'est un dictionnaire
KK = [sol[kk[i]] for i in range(s)]
schema = sympy.Eq( u_np1 , u_n + dt*sum([KK[j]*b[j] for j in range(s)]) )
display(schema)

\(\displaystyle k_{0} = - \beta u_{n}\)

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

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

Code

display(Math(sympy.latex(sympy.Eq(sympy.Symbol('u_{n+1}'),schema.rhs.simplify()))))

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

On peut expliciter \(u_n\) en fonction de \(u_0\) et du produit \(\beta h\) (on notera \(x\) le produit \(\beta h\)) :

Code

q = (schema.rhs/u_n).simplify().subs(beta*dt,x)
display(Math(sympy.latex(sympy.Eq(sympy.Symbol('q(x)'),q.simplify()))))

\(\displaystyle q(x) = \frac{x \left(x - 2\right)}{2} + 1\)

La suite vérifie \(u_n \xrightarrow[n\to \infty]{} 0\) ssi \(\left|q(x)\right|<1\). En étudiant brièvement la fonction \(q(x)\), on trouve que \(q(x)<1\) ssi \(x<6\) donné ci-dessous:.

Code

sympy.plot(q, (x,0,2.5), ylim=(-1.5,1.5),
            ylabel='q(x)', xlabel='x', title='q(x) en fonction de x')

sympy.solve(abs(q)<1)

\(\displaystyle x < 2\)

Correction point 4 : test de l’ordre

Code

##################################################
# EDO
##################################################
t0     = 0
tfinal = 2
y0     = 1

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

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

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

##################################################
# schéma 
##################################################

# Paramètres du schéma recuperés depuis la matrice de Butcher obtenue précédemment
cc = np.array(c)
bb = np.array(b)
AA = np.array(A)

s = len(cc)
# Schéma de Runge-Kutta implicite à s étages
def RK(phi, tt, y0):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0 
    for n in range(len(tt) - 1):
        k_0 = phi(tt[n], uu[n])
        k_1 = phi(tt[n]+cc[1]*h, uu[n]+h*AA[1,0]*k_0)
        KK = [k_0, k_1]
        uu[n+1] = uu[n] + h*sum([KK[j]*b[j] for j in range(s)])
    return uu


##################################################
# ordre
##################################################
H = []
err = []


plt.figure(figsize=(16,5))

ax1 = plt.subplot(1,2,1)
plt.xlabel('$t$')
plt.ylabel('$y$')
ax1.grid(True)

N = 10
for k in range(6):
    N += 10
    tt = np.linspace(t0, tfinal, N + 1)
    H.append( tt[1] - tt[0] )
    yy = sol_exacte(tt)
    uu = RK(phi, tt, y0)
    err.append( np.linalg.norm(uu-yy,np.inf) )
    ax1.plot(tt,uu,label=f'Approchée avec N={</span>N<span class="sc">}')

ax1.plot(tt,yy,label='Exacte') # sur la grille la plus fine
ax1.legend()

ax2 = plt.subplot(1,2,2)
ax2.plot( np.log(H), np.log(err), 'r-o')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.title(f'Le schéma a ordre {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err),<span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:1.2f}')
ax2.grid(True)

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

2.0.2 🎨 Exercice : schéma de Radau IA (\(s=2\), implicite)

On considère la méthode RK donnée par la matri suivante : \( \begin{array}{c|cc} 0 & \frac{1}{4} & -\frac{1}{4} \\ \frac{2}{3} & \frac{1}{4} & \frac{5}{12} \\ \hline & \frac{1}{4} & \frac{3}{4} \end{array} \)

Écrire le schéma correspondant.
Déterminer l’ordre de convergence.
Déterminer sous quelle condition sur \(h\) ce schéma est A-stable.
Tester empiriquement l’ordre de convergence de ce schéma sur le problème de Cauchy \(y'(t) = t(1-t^2)-y(t)\) et \(y(0) = 1\) pour \(t\in[0,2]\).

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

Correction point 1 : schéma

\( \begin{cases} u_0 = y_0 \\ k_0 = \varphi\Big(t_n,u_n+\frac{h}{4}(k_0-k_1)\Big) \\ k_1 = \varphi\Big(t_n+\frac{2}{3}h,u_n+\frac{h}{12}(3k_0+ 5k_1)\Big) \\ u_{n+1}= u_n + \frac{h}{4}\big(k_0+3k_1\big) \end{cases} \)

Il s’agit d’un schéma à \(s=2\) étages implicite.

Correction point 2 : ordre de convergence

D’après la barrière de Butcher, un schéma RK à \(s=2\) étages implicite ne peut pas être d’ordre supérieur à \(s=4\).

On peut bien-sûr écrire les formules à la main, mais on peut aussi utiliser sympy comme suit :

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 = r"\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 = r'\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 = r'\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 = r'\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
b = [ sympy.Rational(1,4), sympy.Rational(3,4) ]
c = [ 0, sympy.Rational(2,3) ]
A = [ [ sympy.Rational(1,4), -sympy.Rational(1,4)] , [ sympy.Rational(1,4),  sympy.Rational(5,12) ] ]

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

Matrice de Butcher

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

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

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{2}{9}\text{ doit être égale à }\frac{1}{4}\)

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

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

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

En conclusion, le schéma est d’ordre 3.

Correction point 3 : A-stabilité

On applique le schéma au problème de Cauchy \(y'(t) = -\beta y(t)\) et \(y(0) = 1\) avec \(\beta>0\) dont la solution est \(y(t)=e^{-\beta t}\).

On obtient alors la relation de récurrence suivante :

Code

u_n = sympy.symbols(r'u_n')
u_np1 = sympy.symbols(r'u_{n+1}')

dt = sympy.symbols(r'h', positive=True)
x = sympy.symbols(r'x', positive=True)

k_0 = sympy.symbols(r'k_0')
k_1 = sympy.symbols(r'k_1')
kk = [k_0, k_1]

beta = sympy.symbols(r'\beta', positive=True)
phi = lambda y: -beta*y

# Eqs = [ sympy.Eq( k_0, phi(u_n+dt/4*(k_0-k_1))   ) ,
#         sympy.Eq( k_1, phi(u_n+dt/4*(3*k_0+k_1)) ) ]
Eqs = [ sympy.Eq( kk[i] ,  phi(u_n + dt*sum([kk[j]*A[i][j] for j in range(s)])) ) for i in range(s) ]
for eq in Eqs:
    display(eq)
    
sol = sympy.solve( Eqs, kk ) # c'est un dictionnaire
display(sol)
KK = [sol[kk[i]] for i in range(s)]
schema = sympy.Eq( u_np1 , u_n + dt*sum([KK[j]*b[j] for j in range(s)]) )
display(schema)

\(\displaystyle k_{0} = - \beta \left(h \left(\frac{k_{0}}{4} - \frac{k_{1}}{4}\right) + u_{n}\right)\)

\(\displaystyle k_{1} = - \beta \left(h \left(\frac{k_{0}}{4} + \frac{5 k_{1}}{12}\right) + u_{n}\right)\)

\(\displaystyle \left\{ k_{0} : \frac{- 4 \beta^{2} h u_{n} - 6 \beta u_{n}}{\beta^{2} h^{2} + 4 \beta h + 6}, \ k_{1} : - \frac{6 \beta u_{n}}{\beta^{2} h^{2} + 4 \beta h + 6}\right\}\)

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

Code

display(Math(sympy.latex(sympy.Eq(sympy.Symbol('u_{n+1}'),schema.rhs.simplify()))))

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

On peut expliciter \(u_n\) en fonction de \(u_0\) et du produit \(\beta h\) (on notera \(x\) le produit \(\beta h\)) :

Code

q = (schema.rhs/u_n).simplify().subs(beta*dt,x)
display(Math(sympy.latex(sympy.Eq(sympy.Symbol('q(x)'),q.simplify()))))

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

La suite vérifie \(u_n \xrightarrow[n\to \infty]{} 0\) ssi \(\left|q(x)\right|<1\). En étudiant brièvement la fonction \(q(x)\), on trouve que \(q(x)<1\) pour tout \(x>0\) :

Code

sympy.plot(q, (x,0,2.5), ylim=(-1.5,1.5),
            ylabel='q(x)', xlabel='x', title='q(x) en fonction de x')

sympy.solve(abs(q)<1)

\(\displaystyle \text{True}\)

Correction point 4 : test de l’ordre

Code

##################################################
# EDO
##################################################
t0     = 0
tfinal = 2
y0     = 1

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

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

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

##################################################
# schéma 
##################################################

# Paramètres du schéma recuperés depuis la matrice de Butcher obtenue précédemment
cc = np.array(c)
bb = np.array(b)
AA = np.array(A)

s = len(cc)
# Schéma de Runge-Kutta implicite à s étages
def RK(phi, tt, y0):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0 
    for n in range(len(tt) - 1):
        k_init = phi(tt[n], uu[n])
        eqs = lambda KK: [ -KK[0] + phi(tt[n]        , uu[n] + h/4*( KK[0] - KK[1]) ) ,
                           -KK[1] + phi(tt[n] + 2/3*h, uu[n] + h/12*( 3*KK[0] + 5*KK[1]) ) 
                           ]
        KK = fsolve( eqs, [k_init,k_init] )
        uu[n+1] = uu[n] + h/4*(KK[0]+ 3*KK[1])
        # eqs = lambda KK: [ -KK[i] + phi(tt[n] + cc[i]*h, uu[n] + h*sum([AA[i,j]*KK[j] for j in range(s)])) for i in range(s)]
        # KK = fsolve( eqs , k_init*np.ones(s) )
        # uu[n+1] = uu[n] + h*sum([KK[j]*b[j] for j in range(s)])
    return uu


##################################################
# ordre
##################################################
H = []
err = []


plt.figure(figsize=(16,5))

ax1 = plt.subplot(1,2,1)
plt.xlabel('$t$')
plt.ylabel('$y$')
ax1.grid(True)

N = 10
for k in range(6):
    N += 10
    tt = np.linspace(t0, tfinal, N + 1)
    yy = sol_exacte(tt)
    uu = RK(phi, tt, y0)
    H.append( tt[1] - tt[0] )
    err.append( np.linalg.norm(uu-yy,np.inf) )
    ax1.plot(tt,uu,label=f'Approchée avec N={</span>N<span class="sc">}')

ax1.plot(tt,yy,label='Exacte') # sur la grille la plus fine
ax1.legend()

ax2 = plt.subplot(1,2,2)
ax2.plot( np.log(H), np.log(err), 'r-o')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.title(f'Le schéma a ordre {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err),<span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:1.2f}')
ax2.grid(True)

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

2.0.3 🎨 Exercice : schéma de Radau IIA (\(s=2\), implicite)

On considère la méthode RK donnée par la matri suivante : \( \begin{array}{c|cc} \frac{1}{3} & \frac{5}{12} & -\frac{1}{12} \\ 1 & \frac{3}{4} & \frac{1}{4} \\ \hline & \frac{3}{4} & \frac{1}{4} \end{array} \)

Écrire le schéma correspondant.
Déterminer l’ordre de convergence.
Déterminer sous quelle condition sur \(h\) ce schéma est A-stable.
Tester empiriquement l’ordre de convergence de ce schéma sur le problème de Cauchy \(y'(t) = t(1-t^2)-y(t)\) et \(y(0) = 1\) pour \(t\in[0,2]\).

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

Correction point 1 : schéma

\( \begin{cases} u_0 = y_0 \\ k_0 = \varphi\Big(t_n+\frac{1}{3}h ,u_n+\frac{h}{12}(5k_0-k_1)\Big) \\ k_1 = \varphi\Big(t_{n+1}, u_n+\frac{h}{4}(3k_0+k_1)\Big)\\ u_{n+1}= u_n + \frac{h}{4}\big(3k_0+k_1\big) \end{cases} \)

Il s’agit d’un schéma à \(s=2\) étages implicite.

Correction point 2 : ordre de convergence

D’après la barrière de Butcher, un schéma RK à \(s=2\) étages implicite ne peut pas être d’ordre supérieur à \(s=4\).

On peut bien-sûr écrire les formules à la main, mais on peut aussi utiliser sympy comme suit :

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 = r"\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 = r'\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 = r'\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 = r'\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
b = [ sympy.Rational(3,4), sympy.Rational(1,4) ]
c = [ sympy.Rational(1,3) , 1 ]
A = [ [ sympy.Rational(5,12), -sympy.Rational(1,12)] , [ sympy.Rational(3,4),  sympy.Rational(1,4) ] ]

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

Matrice de Butcher

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

En conclusion, le schéma est d’ordre 3.

Correction point 3 : A-stabilité

On applique le schéma au problème de Cauchy \(y'(t) = -\beta y(t)\) et \(y(0) = 1\) avec \(\beta>0\) dont la solution est \(y(t)=e^{-\beta t}\).

On obtient alors la relation de récurrence suivante :

Code

u_n = sympy.symbols(r'u_n')
u_np1 = sympy.symbols(r'u_{n+1}')

dt = sympy.symbols(r'h', positive=True)
x = sympy.symbols(r'x', positive=True)

k_0 = sympy.symbols(r'k_0')
k_1 = sympy.symbols(r'k_1')
kk = [k_0, k_1]

beta = sympy.symbols(r'\beta', positive=True)
phi = lambda y: -beta*y

Eqs = [ sympy.Eq( kk[i] ,  phi(u_n + dt*sum([kk[j]*A[i][j] for j in range(s)])) ) for i in range(s) ]
for eq in Eqs:
    display(eq)
    
sol = sympy.solve( Eqs, kk ) # c'est un dictionnaire
KK = [sol[kk[i]] for i in range(s)]
schema = sympy.Eq( u_np1 , u_n + dt*sum([KK[j]*b[j] for j in range(s)]) )
display(schema)

\(\displaystyle k_{0} = - \beta \left(h \left(\frac{5 k_{0}}{12} - \frac{k_{1}}{12}\right) + u_{n}\right)\)

\(\displaystyle k_{1} = - \beta \left(h \left(\frac{3 k_{0}}{4} + \frac{k_{1}}{4}\right) + u_{n}\right)\)

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

Code

display(Math(sympy.latex(sympy.Eq(sympy.Symbol('u_{n+1}'),schema.rhs.simplify()))))

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

On peut expliciter \(u_n\) en fonction de \(u_0\) et du produit \(\beta h\) (on notera \(x\) le produit \(\beta h\)) :

Code

q = (schema.rhs/u_n).simplify().subs(beta*dt,x)
display(Math(sympy.latex(sympy.Eq(sympy.Symbol('q(x)'),q.simplify()))))

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

La suite vérifie \(u_n \xrightarrow[n\to \infty]{} 0\) ssi \(\left|q(x)\right|<1\). En étudiant brièvement la fonction \(q(x)\), on trouve que \(q(x)<1\) pour tout \(x>0\) :

Code

sympy.plot(q, (x,0,2.5), ylim=(-1.5,1.5),
            ylabel='q(x)', xlabel='x', title='q(x) en fonction de x')

sympy.solve(abs(q)<1)

\(\displaystyle \text{True}\)

Correction point 4 : test de l’ordre

Code

##################################################
# EDO
##################################################
t0     = 0
tfinal = 2
y0     = 1

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

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

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

##################################################
# schéma 
##################################################

# Paramètres du schéma recuperés depuis la matrice de Butcher obtenue précédemment
cc = np.array(c)
bb = np.array(b)
AA = np.array(A)

s = len(cc)
# Schéma de Runge-Kutta implicite à s étages
def RK(phi, tt, y0):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0 
    for n in range(len(tt) - 1):
        eqs = lambda KK: [ -KK[i] + phi(tt[n] + cc[i]*h, uu[n] + h*sum([AA[i,j]*KK[j] for j in range(s)])) for i in range(s)]
        k_init = phi(tt[n], uu[n])
        KK = fsolve( eqs , k_init*np.ones(s) )
        uu[n+1] = uu[n] + h*sum([KK[j]*b[j] for j in range(s)])
    return uu


##################################################
# ordre
##################################################
H = []
err = []


plt.figure(figsize=(16,5))

ax1 = plt.subplot(1,2,1)
plt.xlabel('$t$')
plt.ylabel('$y$')
ax1.grid(True)

N = 10
for k in range(6):
    N += 10
    tt = np.linspace(t0, tfinal, N + 1)
    H.append( tt[1] - tt[0] )
    yy = sol_exacte(tt)
    uu = RK(phi, tt, y0)
    err.append( np.linalg.norm(uu-yy,np.inf) )
    ax1.plot(tt,uu,label=f'Approchée avec N={</span>N<span class="sc">}')

ax1.plot(tt,yy,label='Exacte') # sur la grille la plus fine
ax1.legend()

ax2 = plt.subplot(1,2,2)
ax2.plot( np.log(H), np.log(err), 'r-o')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.title(f'Le schéma a ordre {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err),<span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:1.2f}')
ax2.grid(True)

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

2.0.4 🎨 Exercice (\(s=3\), explicite)

On considère la méthode RK donnée par la matrice suivante : \( \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \frac{1}{2} & \frac{1}{2} & 0 & 0 \\ 1 & 0 & 1 & 0 \\ \hline & \frac{1}{3} & \frac{1}{3} & \frac{1}{3} \end{array} \)

Écrire le schéma correspondant.
Déterminer l’ordre de convergence.
Déterminer sous quelle condition sur \(h\) ce schéma est A-stable.
Tester empiriquement l’ordre de convergence de ce schéma sur le problème de Cauchy \(y'(t) = t(1-t^2)-y(t)\) et \(y(0) = 1\) pour \(t\in[0,2]\).

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

Correction point 1 : schéma

\( \begin{cases} u_0 = y_0 \\ k_0 = \varphi\Big(t_n,u_n\Big) \\ k_1 = \varphi\Big(t_n+\frac{h}{2},u_n+\frac{h}{2}k_0\Big) \\ k_2 = \varphi\Big(t_n+h,u_n+hk_1\Big) \\ u_{n+1}= u_n + h\big(\frac{1}{3}k_0+\frac{1}{3}k_1+\frac{1}{3}k_2\big) \end{cases} \)

Il s’agit d’un schéma à \(s=3\) étages explicite.

Correction point 2 : ordre de convergence

D’après la barrière de Butcher, un schéma RK à \(s=3\) étages explicite ne peut pas être d’ordre supérieur à \(s=3\).

On peut bien-sûr écrire les formules à la main, mais on peut aussi utiliser un solveur pour résoudre le système d’équations non-linéaires comme suit :

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 = r"\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 = r'\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 = r'\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 = r'\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 = 3

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

# matrice dans le cas général
# b = [b_0, b_1, b_2]
# c = [c_0, c_1, c_2]
# A = [[a_00, a_01, a_02], [a_10, a_11, a_12], [a_20, a_21, a_22]]

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

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

Matrice de Butcher

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\\frac{1}{2} & \frac{1}{2} & 0 & 0\\1 & 0 & 1 & 0\\ & \frac{1}{3} & \frac{1}{3} & \frac{1}{3}\end{matrix}\right]\)

On a 4 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{1}{2}\text{ doit être égale à }\frac{1}{2}\)

\(\displaystyle \sum_{j=1}^s a_{2j}=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{5}{12}\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}\)

En conclusion, le schéma est d’ordre 2.

Correction point 3 : A-stabilité

On applique le schéma au problème de Cauchy \(y'(t) = -\beta y(t)\) et \(y(0) = 1\) avec \(\beta>0\) dont la solution est \(y(t)=e^{-\beta t}\).

On obtient alors la relation de récurrence suivante :

Code

u_n = sympy.symbols(r'u_n')
u_np1 = sympy.symbols(r'u_{n+1}')

dt = sympy.symbols(r'h', positive=True)
x = sympy.symbols(r'x', positive=True)

k_0 = sympy.symbols(r'k_0')
k_1 = sympy.symbols(r'k_1')
k_2 = sympy.symbols(r'k_2')
kk = [k_0, k_1, k_2]

beta = sympy.symbols(r'\beta', positive=True)
phi = lambda y: -beta*y

Eqs = [ sympy.Eq( kk[i] ,  phi(u_n + dt*sum([kk[j]*A[i][j] for j in range(s)])) ) for i in range(s) ]
for eq in Eqs:
    display(eq)
    
sol = sympy.solve( Eqs, kk ) # c'est un dictionnaire
KK = [sol[kk[i]] for i in range(s)]
schema = sympy.Eq( u_np1 , u_n + dt*sum([KK[j]*b[j] for j in range(s)]) )
display(schema)

\(\displaystyle k_{0} = - \beta u_{n}\)

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

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

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

Code

display(Math(sympy.latex(sympy.Eq(sympy.Symbol('u_{n+1}'),schema.rhs.simplify()))))

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

On peut expliciter \(u_n\) en fonction de \(u_0\) et du produit \(\beta h\) (on notera \(x\) le produit \(\beta h\)) :

Code

q = (schema.rhs/u_n).simplify().subs(beta*dt,x)
display(Math(sympy.latex(sympy.Eq(sympy.Symbol('q(x)'),q.simplify()))))

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

La suite vérifie \(u_n \xrightarrow[n\to \infty]{} 0\) ssi \(\left|q(x)\right|<1\). En étudiant brièvement la fonction \(q(x)\), on trouve que \(q(x)<1\) ssi \(x<6\) donné ci-dessous:.

Code

sympy.plot(q, (x,0,4), ylim=(-1.5,1.5),
            ylabel='q(x)', xlabel='x', title='q(x) en fonction de x')

sol = sympy.solveset(q+1, x, domain=sympy.S.Reals)  
display(Markdown(f"$\\beta h<{</span>sol<span class="sc">.</span>args[<span class="dv">0</span>]<span class="sc">.</span>evalf()<span class="sc">}$"))

\(\beta h<2.51274532661833\)

Correction point 4 : test de l’ordre

Code

##################################################
# EDO
##################################################
t0     = 0
tfinal = 2
y0     = 1

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

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

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

##################################################
# schéma 
##################################################

# Paramètres du schéma recuperés depuis la matrice de Butcher obtenue précédemment
cc = np.array(c)
bb = np.array(b)
AA = np.array(A)

s = len(cc)
# Schéma de Runge-Kutta implicite à s étages
def RK(phi, tt, y0):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0 
    for n in range(len(tt) - 1):
        k_0 = phi(tt[n], uu[n])
        k_1 = phi(tt[n]+cc[1]*h, uu[n]+h*AA[1,0]*k_0)
        k_2 = phi(tt[n]+cc[2]*h, uu[n]+h*(AA[2,0]*k_0+AA[2,1]*k_1))
        KK = [k_0, k_1, k_2]
        uu[n+1] = uu[n] + h*sum([KK[j]*b[j] for j in range(s)])
    return uu


##################################################
# ordre
##################################################
H = []
err = []


plt.figure(figsize=(16,5))

ax1 = plt.subplot(1,2,1)
plt.xlabel('$t$')
plt.ylabel('$y$')
ax1.grid(True)

N = 10
for k in range(6):
    N += 10
    tt = np.linspace(t0, tfinal, N + 1)
    H.append( tt[1] - tt[0] )
    yy = sol_exacte(tt)
    uu = RK(phi, tt, y0)
    err.append( np.linalg.norm(uu-yy,np.inf) )
    ax1.plot(tt,uu,label=f'Approchée avec N={</span>N<span class="sc">}')

ax1.plot(tt,yy,label='Exacte') # sur la grille la plus fine
ax1.legend()

ax2 = plt.subplot(1,2,2)
ax2.plot( np.log(H), np.log(err), 'r-o')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.title(f'Le schéma a ordre {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err),<span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:1.2f}')
ax2.grid(True)

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

2.0.5 🎨 Exercice (\(s=3\), implicite)

On considère la méthode RK donnée par la matrice suivante : \( \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \frac{1}{2} & \frac{5}{24} & \frac{1}{3} & -\frac{1}{24} \\ 1 & 0 & 1 & 0 \\ \hline & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \end{array} \)

Écrire le schéma correspondant.
Déterminer l’ordre de convergence.
Déterminer sous quelle condition sur \(h\) ce schéma est A-stable.
Tester empiriquement l’ordre de convergence de ce schéma sur le problème de Cauchy \(y'(t) = -y(t)\) et \(y(0) = 1\) pour \(t\in[0,10]\).

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

Correction point 1 : schéma

\( \begin{cases} u_0 = y_0 \\ k_0 = \varphi\Big(t_n,u_n\Big) \\ k_1 = \varphi\Big(t_n+\frac{h}{2},u_n+h\big(\frac{5}{24}k_0+\frac{1}{3}k_1-\frac{1}{24}k_2\big)\Big) \\ k_2 = \varphi\Big(t_n+h,u_n+hk_1\Big) \\ u_{n+1}= u_n + h\left(\frac{1}{6}k_0+\frac{2}{3}k_1+\frac{1}{6}k_2\right) \end{cases} \)

Il s’agit d’un schéma à \(s=3\) étages implicite.

Correction point 2 : ordre de convergence

D’après la barrière de Butcher, un schéma RK à \(s=3\) étages implicite ne peut pas être d’ordre supérieur à \(s=6\).

On peut bien-sûr écrire les formules à la main, mais on peut aussi utiliser un solveur pour résoudre le système d’équations non-linéaires comme suit :

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 = r"\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 = r'\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 = r'\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 = r'\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 = 3

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

# matrice dans le cas général
# b = [b_0, b_1, b_2]
# c = [c_0, c_1, c_2]
# A = [[a_00, a_01, a_02], [a_10, a_11, a_12], [a_20, a_21, a_22]]

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

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

Matrice de Butcher

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\\frac{1}{2} & \frac{5}{24} & \frac{1}{3} & - \frac{1}{24}\\1 & 0 & 1 & 0\\ & \frac{1}{6} & \frac{2}{3} & \frac{1}{6}\end{matrix}\right]\)

On a 4 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{1}{2}\text{ doit être égale à }\frac{1}{2}\)

\(\displaystyle \sum_{j=1}^s a_{2j}=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{1}{4}\text{ doit être égale à }\frac{1}{4}\)

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

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

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

En conclusion, le schéma est d’ordre 3.

Correction point 3 : A-stabilité

On applique le schéma au problème de Cauchy \(y'(t) = -\beta y(t)\) et \(y(0) = 1\) avec \(\beta>0\) dont la solution est \(y(t)=e^{-\beta t}\).

On obtient alors la relation de récurrence suivante :

Code

u_n = sympy.symbols(r'u_n')
u_np1 = sympy.symbols(r'u_{n+1}')

dt = sympy.symbols(r'h', positive=True)
x = sympy.symbols(r'x', positive=True)

k_0 = sympy.symbols(r'k_0')
k_1 = sympy.symbols(r'k_1')
k_2 = sympy.symbols(r'k_2')
kk = [k_0, k_1, k_2]

beta = sympy.symbols(r'\beta', positive=True)
phi = lambda y: -beta*y

Eqs = [ sympy.Eq( kk[i] ,  phi(u_n + dt*sum([kk[j]*A[i][j] for j in range(s)])) ) for i in range(s) ]
for eq in Eqs:
    display(eq)
    
sol = sympy.solve( Eqs, kk ) # c'est un dictionnaire
KK = [sol[kk[i]] for i in range(s)]
schema = sympy.Eq( u_np1 , u_n + dt*sum([KK[j]*b[j] for j in range(s)]) )
display(schema)

\(\displaystyle k_{0} = - \beta u_{n}\)

\(\displaystyle k_{1} = - \beta \left(h \left(\frac{5 k_{0}}{24} + \frac{k_{1}}{3} - \frac{k_{2}}{24}\right) + u_{n}\right)\)

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

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

Code

display(Math(sympy.latex(sympy.Eq(sympy.Symbol('u_{n+1}'),schema.rhs.simplify()))))

\(\displaystyle u_{n+1} = \frac{u_{n} \left(- \beta^{3} h^{3} + 5 \beta^{2} h^{2} - 16 \beta h + 24\right)}{\beta^{2} h^{2} + 8 \beta h + 24}\)

On peut expliciter \(u_n\) en fonction de \(u_0\) et du produit \(\beta h\) (on notera \(x\) le produit \(\beta h\)) :

Code

q = (schema.rhs/u_n).simplify().subs(beta*dt,x)
display(Math(sympy.latex(sympy.Eq(sympy.Symbol('q(x)'),q.simplify()))))

\(\displaystyle q(x) = \frac{- x^{3} + 5 x^{2} - 16 x + 24}{x^{2} + 8 x + 24}\)

La suite vérifie \(u_n \xrightarrow[n\to \infty]{} 0\) ssi \(\left|q(x)\right|<1\). En étudiant brièvement la fonction \(q(x)\), on trouve que \(q(x)<1\) ssi \(x<6\) donné ci-dessous:.

Code

sympy.plot(q, (x,0,10), ylim=(-1.5,1.5),
            ylabel='q(x)', xlabel='x', title='q(x) en fonction de x')

sol = sympy.solveset(q+1, x, domain=sympy.S.Reals)  
display(Markdown(f"$\\beta h<{</span>sol<span class="sc">.</span>args[<span class="dv">0</span>]<span class="sc">}$"))

\(\beta h<6\)

Correction point 4 : test de l’ordre

Code

##################################################
# EDO
##################################################
t0     = 0
tfinal = 10
y0     = 1

phi = lambda t,y: -y

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

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

##################################################
# schéma 
##################################################

# Paramètres du schéma recuperés depuis la matrice de Butcher obtenue précédemment
cc = np.array(c)
bb = np.array(b)
AA = np.array(A)

s = len(cc)
# Schéma de Runge-Kutta implicite à s étages
def RK(phi, tt, y0):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0 
    for n in range(len(tt) - 1):
        k_init = phi(tt[n],uu[n])
        KK = fsolve( lambda kk: [ -kk[i] +  phi( tt[n]+cc[i]*h ,  uu[n] + h*sum([kk[j]*AA[i,j] for j in range(s)]) ) for i in range(s) ], k_init*np.ones((s,1)) )
        uu[n+1] = uu[n] + h*sum([KK[j]*b[j] for j in range(s)])
    return uu


##################################################
# ordre
##################################################
H = []
err = []


plt.figure(figsize=(16,5))

ax1 = plt.subplot(1,2,1)
plt.xlabel('$t$')
plt.ylabel('$y$')
ax1.grid(True)

N = 10
for k in range(6):
    N += 10
    tt = np.linspace(t0, tfinal, N + 1)
    H.append( tt[1] - tt[0] )
    yy = sol_exacte(tt)
    uu = RK(phi, tt, y0)
    err.append( np.linalg.norm(uu-yy,np.inf) )
    ax1.plot(tt,uu,label=f'Approchée avec N={</span>N<span class="sc">}')

ax1.plot(tt,yy,label='Exacte') # sur la grille la plus fine
ax1.legend()

ax2 = plt.subplot(1,2,2)
ax2.plot( np.log(H), np.log(err), 'r-o')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.title(f'Le schéma a ordre {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err),<span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:1.2f}')
ax2.grid(True)

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

2.0.6 🎨 Exercice (\(s=3\), explicite avec paramètres)

On considère la méthode RK donnée par la matrice suivante : \( \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \frac{1}{3} & \frac{1}{3} & 0 & 0 \\ \frac{2}{3} & 0 & \frac{2}{3} & 0 \\ \hline & b_0 & b_1 & 0 \end{array} \)

Écrire le schéma correspondant.
Déterminer \(b_0\) et \(b_1\) pour que l’ordre de convergence soit maximal. Dans la suite de l’exercice on prendra ces valeurs.
Déterminer sous quelle condition sur \(h\) ce schéma est A-stable.
Tester empiriquement l’ordre de convergence de ce schéma sur le problème de Cauchy \(y'(t) = y(t)(1-y(t))t\) et \(y(0) = 2\) pour \(t\in[0,2]\).

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

Correction point 1 : schéma

\( \begin{cases} u_0 = y_0 \\ k_0 = \varphi\Big(t_n,u_n\Big) \\ k_1 = \varphi\Big(t_n+\frac{h}{3},u_n+\frac{h}{3}k_0\Big) \\ k_2 = \varphi\Big(t_n+\frac{2}{3}h,u_n+\frac{2}{3}hk_1\Big) \\ u_{n+1}= u_n + h\big(b_0k_0+b_1k_1\big) \end{cases} \)

Il s’agit d’un schéma à \(s=3\) étages explicite.

Correction point 2 : ordre de convergence

D’après la barrière de Butcher, un schéma RK à \(s=3\) étages explicite ne peut pas être d’ordre supérieur à \(s=3\).

On peut bien-sûr écrire les formules à la main, mais on peut aussi utiliser sympy comme suit :

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 = r"\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 = r'\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 = r'\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 = r'\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 = 3

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

# matrice dans le cas général
# b = [b_0, b_1, b_2]
# c = [c_0, c_1, c_2]
# A = [[a_00, a_01, a_02], [a_10, a_11, a_12], [a_20, a_21, a_22]]

# matrice dans notre cas particulier
b_0, b_1 = sympy.symbols('b_0 b_1', real=True)
b = [ b_0, b_1, 0 ]
c = [ 0, sympy.Rational(1,3), sympy.Rational(2,3) ]
A = [[ 0, 0, 0], 
    [ sympy.Rational(1,3), 0, 0],
    [ 0, sympy.Rational(2,3), 0]]

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

Matrice de Butcher

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\\frac{1}{3} & \frac{1}{3} & 0 & 0\\\frac{2}{3} & 0 & \frac{2}{3} & 0\\ & b_{0} & b_{1} & 0\end{matrix}\right]\)

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

\(\displaystyle \sum_{j=1}^s b_j =b_0 + b_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{1}{3}\text{ doit être égale à }\frac{1}{3}\)

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

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

\(\displaystyle \sum_{j=1}^s b_j c_j=\frac{b_{1}}{3}\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{b_{1}}{9}\text{ doit être égale à }\frac{1}{3}\)

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

En conclusion, si \(b_1=\frac{3}{2}\) et \(b_0=-\frac{1}{2}\) le schéma est d’ordre 2 :

Code

b_1 = sympy.Rational(3,2)
b_0 = 1-b_1
b = [ b_0, b_1, 0 ]
Eqs = ordre_RK(s, A, b, c)

Matrice de Butcher

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0\\\frac{1}{3} & \frac{1}{3} & 0 & 0\\\frac{2}{3} & 0 & \frac{2}{3} & 0\\ & - \frac{1}{2} & \frac{3}{2} & 0\end{matrix}\right]\)

On a 4 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{1}{3}\text{ doit être égale à }\frac{1}{3}\)

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

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}{6}\text{ doit être égale à }\frac{1}{3}\)

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

Correction point 3 : A-stabilité

On applique le schéma au problème de Cauchy \(y'(t) = -\beta y(t)\) et \(y(0) = 1\) avec \(\beta>0\) dont la solution est \(y(t)=e^{-\beta t}\).

On obtient alors la relation de récurrence suivante :

Code

u_n = sympy.symbols(r'u_n')
u_np1 = sympy.symbols(r'u_{n+1}')

dt = sympy.symbols(r'h', positive=True)
x = sympy.symbols(r'x', positive=True)

k_0 = sympy.symbols(r'k_0')
k_1 = sympy.symbols(r'k_1')
k_2 = sympy.symbols(r'k_2')
kk = [k_0, k_1, k_2]

beta = sympy.symbols(r'\beta', positive=True)
phi = lambda y: -beta*y

Eqs = [ sympy.Eq( kk[i] ,  phi(u_n + dt*sum([kk[j]*A[i][j] for j in range(s)])) ) for i in range(s) ]
for eq in Eqs:
    display(eq)
    
sol = sympy.solve( Eqs, kk ) # c'est un dictionnaire
KK = [sol[kk[i]] for i in range(s)]
schema = sympy.Eq( u_np1 , u_n + dt*sum([KK[j]*b[j] for j in range(s)]) )
display(schema)

\(\displaystyle k_{0} = - \beta u_{n}\)

\(\displaystyle k_{1} = - \beta \left(\frac{h k_{0}}{3} + u_{n}\right)\)

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

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

Code

display(Math(sympy.latex(sympy.Eq(sympy.Symbol('u_{n+1}'),schema.rhs.simplify()))))

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

On peut expliciter \(u_n\) en fonction de \(u_0\) et du produit \(\beta h\) (on notera \(x\) le produit \(\beta h\)) :

Code

q = (schema.rhs/u_n).simplify().subs(beta*dt,x)
display(Math(sympy.latex(sympy.Eq(sympy.Symbol('q(x)'),q.simplify()))))

\(\displaystyle q(x) = \frac{x \left(x - 2\right)}{2} + 1\)

La suite vérifie \(u_n \xrightarrow[n\to \infty]{} 0\) ssi \(\left|q(x)\right|<1\). En étudiant brièvement la fonction \(q(x)\), on trouve que \(q(x)<1\) ssi \(x<6\) donné ci-dessous:.

Code

sympy.plot(q, (x,0,2.5), ylim=(-1.5,1.5),
            ylabel='q(x)', xlabel='x', title='q(x) en fonction de x')

sympy.solve(abs(q)<1)

\(\displaystyle x < 2\)

Correction point 4 : test de l’ordre

Code

##################################################
# EDO
##################################################
t0     = 0
tfinal = 2
y0     = 2

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

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

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

##################################################
# schéma 
##################################################

# Paramètres du schéma recuperés depuis la matrice de Butcher obtenue précédemment
cc = np.array(c)
bb = np.array(b)
AA = np.array(A)

s = len(cc)
# Schéma de Runge-Kutta implicite à s étages
def RK(phi, tt, y0):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[0] = y0 
    for n in range(len(tt) - 1):
        k_0 = phi(tt[n], uu[n])
        k_1 = phi(tt[n]+cc[1]*h, uu[n]+h*AA[1,0]*k_0)
        k_2 = phi(tt[n]+cc[2]*h, uu[n]+h*(AA[2,0]*k_0+AA[2,1]*k_1))
        KK = [k_0, k_1, k_2]
        uu[n+1] = uu[n] + h*sum([KK[j]*b[j] for j in range(s)])
    return uu


##################################################
# ordre
##################################################
H = []
err = []


plt.figure(figsize=(16,5))

ax1 = plt.subplot(1,2,1)
plt.xlabel('$t$')
plt.ylabel('$y$')
ax1.grid(True)

N = 10
for k in range(6):
    N += 10
    tt = np.linspace(t0, tfinal, N + 1)
    H.append( tt[1] - tt[0] )
    yy = sol_exacte(tt)
    uu = RK(phi, tt, y0)
    err.append( np.linalg.norm(uu-yy,np.inf) )
    ax1.plot(tt,uu,label=f'Approchée avec N={</span>N<span class="sc">}')

ax1.plot(tt,yy,label='Exacte') # sur la grille la plus fine
ax1.legend()

ax2 = plt.subplot(1,2,2)
ax2.plot( np.log(H), np.log(err), 'r-o')
plt.xlabel(r'$\ln(h)$')
plt.ylabel(r'$\ln(e)$')
plt.title(f'Le schéma a ordre {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err),<span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:1.2f}')
ax2.grid(True)

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

3 Catalogue de quelques schémas 🎨 Runge-Kutta rencontrés dans la littérature

Dans la littérature on trouve beaucoup de schémas de Runge-Kutta. Ci-dessous une petite liste de ces méthodes. Pour chaque méthode de Runge-Kutta donné à travers sa matrice de Butcher, il est indiqué le nombre d’étages \(s\), son ordre de convergence \(\omega\) comme reporté dans la littérature et si le schéma est explicite, semi-implicite ou implicite.

Exercice : étudier l’ordre théorique de convergence (au plus jusqu’à 4) et vérifier empiriquement si l’ordre annoncé coïncide avec l’ordre observé sur un problème de Cauchy dont on connaît la solution exacte.

3.1 Schémas à \(s=1\) étage

\( \begin{array}{|c|l|c|c|c|} \hline & \text{Matrice de Butcher} & s & \omega & \text{Exp/Imp/Semi-Implicite} \\ \hline \text{Gauss-1} & \begin{array}{c|cc} \frac{1}{2} & \frac{1}{2} \\ \hline & 1 \end{array} & 1 & 2 & I \\ \hline \text{EE} & \begin{array}{c|cc} 1 & 1 \\ \hline & 1 \end{array} & 1 & 1 & E \\ \hline \text{EI} & \begin{array}{c|cc} 0 & 0 \\ \hline & 1 \end{array} & 1 & 1 & I \\ \hline \end{array} \)

3.2 Schémas à \(s=2\) étages

\( \begin{array}{|l|l|l|l|l|} \hline & \text{Matrice de Butcher} & s & \omega & \text{Exp/Imp/Semi-Implicite} \\ \hline \text{SDIRK order 3} & \begin{array}{c|cc} \gamma & \gamma & 0\\ 1-\gamma & 1-2\gamma & \gamma\\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} \text{avec } \gamma=\frac{3\pm\sqrt{3}}{6} & 2 & 3 & \text{SI} \\ \hline \text{Gauss-2} & \begin{array}{c|cc} \frac{1}{2}-\frac{\sqrt{3}}{6} & \frac{1}{4} & \frac{1}{4}-\frac{\sqrt{3}}{6}\\ \frac{1}{2}+\frac{\sqrt{3}}{6} & \frac{1}{4}+\frac{\sqrt{3}}{6} & \frac{1}{4}\\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} & 2 & 4 & \text{I} \\ \hline \text{DIRK-2} & \begin{array}{c|cc} \frac{1}{3} & \frac{1}{3} & 0\\ 1 & 1 & 0\\ \hline & \frac{3}{4} & \frac{1}{4} \end{array} & 2 & 3 & \text{SI} \\ \hline \text{Hammer and Hollingsworth} & \begin{array}{c|cc} 0 & 0 & 0\\ \frac{2}{3} & \frac{1}{3} & \frac{1}{3}\\ \hline & \frac{1}{4} & \frac{3}{4} \end{array} & 2 & ? & \text{SI} \\ \hline \text{Randau-IIa-2} & \begin{array}{c|cc} \frac{1}{3} & \frac{5}{12} & -\frac{1}{12}\\ 1 & \frac{3}{4} & \frac{1}{4}\\ \hline & \frac{3}{4} & \frac{1}{4} \end{array} & 2 & 3 & \text{I} \\ \hline \text{CN} & \begin{array}{c|cc} 0 & 0 & 0\\ 1 & \frac{1}{2} & \frac{1}{2}\\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} & 2 & 2 & \text{I} \\ \hline \text{Heun} & \begin{array}{c|cc} 0 & 0 & 0\\ 1 & 1 & 0\\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} & 2 & 2 & \text{E} \\ \hline \text{EM} & \begin{array}{c|cc} 0 & 0 & 0\\ \frac{1}{2} & \frac{1}{2} & 0\\ \hline & 0 & 1 \end{array} & 2 & 2 & \text{E} \\ \hline & \begin{array}{c|cc} 0 & 0 & 0\\ \frac{2}{3} & \frac{2}{3} & 0\\ \hline & \frac{1}{4} & \frac{3}{4}\end{array} & 2 & 2 & \text{E} \\ \hline \text{Ralston} & \begin{array}{c|cc} 0 & 0 & 0\\ \frac{3}{4} & \frac{3}{4} & 0\\ \hline & \frac{1}{3} & \frac{2}{3} \end{array} & 2 & 2 & \text{E} \\ \hline & \begin{array}{c|cc} 0 & 0 & 0\\ \frac{3}{4} & \frac{3}{4} & 0\\ \hline & 0 & 1 \end{array} & 2 & 1 & \text{E} \\ \hline & \begin{array}{c|cc} 0 & 0 & 0\\ \frac{2}{3} & \frac{2}{3} & 0\\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} & 2 & 1 & \text{E} \\ \hline \end{array} \)

3.3 Schémas à \(s=3\) étages

\( \begin{array}{|l|l|l|l|l|} \hline & \text{Matrice de Butcher} & s & \omega & \text{Exp/Imp/Semi-Implicite} \\ \hline \text{Gauss-3} & \begin{array}{c|ccc} \frac{1}{2}-\frac{\sqrt{15}}{10} & \frac{5}{36} & \frac{2}{9}-\frac{\sqrt{15}}{15} & \frac{5}{36}-\frac{\sqrt{15}}{30} \\ \frac{1}{2} & \frac{5}{36}+\frac{\sqrt{15}}{24} & \frac{2}{9}&\frac{5}{36}-\frac{\sqrt{15}}{24}\\ \frac{1}{2}+\frac{\sqrt{15}}{10} & \frac{5}{36}+\frac{\sqrt{15}}{30} &\frac{2}{9}+\frac{\sqrt{15}}{15}&\frac{5}{36}\\ \hline & \frac{5}{18} &\frac{4}{9}& \frac{5}{18} \end{array} & 3 & 6 & I \\ \hline \text{Radau-IIa-3} & \begin{array}{c|ccc} \frac{4-\sqrt{6}}{10} & \frac{88-7\sqrt{6}}{360} & \frac{296-169\sqrt{6}}{1800} & \frac{-2+3\sqrt{6}}{225} \\ \frac{4+\sqrt{6}}{10} & \frac{296+169\sqrt{6}}{1800} & \frac{88+7\sqrt{6}}{360}&\frac{-2-3\sqrt{6}}{225}\\ 1 & \frac{16-\sqrt{6}}{36} &\frac{16+\sqrt{6}}{36}&\frac{1}{9}\\ \hline & \frac{16-\sqrt{6}}{36} &\frac{16+\sqrt{6}}{36}& \frac{1}{9} \end{array} & 3 & 5 & I \\ \hline \text{Lobatto IIIa} & \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \frac{1}{2} & \frac{5}{24} & \frac{1}{3}&-\frac{1}{24}\\ 1 & \frac{1}{6} &\frac{2}{3}&\frac{1}{6}\\ \hline & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \end{array} & 3 & 4 & I \\ \hline \text{Lobatto} & \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \frac{1}{2} & \frac{1}{4} & \frac{1}{4}&0\\ 1 & 0 &1&0\\ \hline & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \end{array} & 3 & 4 & SI \\ \hline \text{Heun3} & \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \frac{1}{3} & \frac{1}{3} & 0&0\\ \frac{2}{3} & 0 &\frac{2}{3}&0\\ \hline & \frac{1}{4} & 0 & \frac{3}{4} \end{array} & 3 & 3 & E \\ \hline \text{Kutta} & \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \frac{1}{2} & \frac{1}{2} & 0&0\\ 1 & -1 &2&0\\ \hline & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \end{array} & 3 & 3 & E \\ \hline \text{Nystron} & \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \frac{2}{3} & \frac{2}{3} & 0&0\\ \frac{2}{3} & 0 &\frac{2}{3}&0\\ \hline & \frac{1}{4} & \frac{3}{8} & \frac{3}{8} \end{array} & 3 & 3 & E \\ \hline \end{array} \)

3.4 Schémas à \(s=4\) étages

\( \begin{array}{|l|l|l|l|l|} \hline & \text{Matrice de Butcher} & s & \omega & \text{Exp/Imp/Semi-Implicite} \\ \hline \text{Lobatto} & \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\ \frac{5-\sqrt{5}}{10} & \frac{5+\sqrt{5}}{60} & \frac{1}{6} & \frac{15-7\sqrt{5}}{60} & 0\\ \frac{5+\sqrt{5}}{10} & \frac{5-\sqrt{5}}{60} &\frac{15+7\sqrt{5}}{60} &\frac{1}{6}&0\\ 1 & \frac{1}{6} & \frac{5-\sqrt{5}}{12} & \frac{5+\sqrt{5}}{12}& 0 \\ \hline & \frac{1}{12} & \frac{5}{12} & \frac{5}{12} & \frac{1}{12} \end{array} & 4 & 6 & I \\ \hline & \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\ \frac{1}{4} & \frac{1}{8} & \frac{1}{8} & 0 & 0\\ \frac{7}{10} & -\frac{1}{100} &\frac{14}{25} &\frac{3}{20}&0\\ 1 & \frac{2}{7} & 0 & \frac{5}{7} & 0\\ \hline & \frac{1}{14} & \frac{32}{81} & \frac{250}{567} & \frac{5}{54} \end{array} & 4 & 5 & SI \\ \hline \text{RK4-1} & \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\ \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0\\ \frac{1}{2} & 0 &\frac{1}{2} &0&0\\ 1 & 0 & 0 & 1 & 0\\ \hline & \frac{1}{6} & \frac{1}{3} & \frac{1}{3} & \frac{1}{6} \end{array} & 4 & 4 & E \\ \hline \text{RK4-2} & \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\ \frac{1}{4} & \frac{1}{4} & 0 & 0 & 0\\ \frac{1}{2} & 0 &\frac{1}{2} &0&0\\ 1 & 1 & -2 & 2 & 0\\ \hline & \frac{1}{6} & 0 & \frac{2}{3} & \frac{1}{6} \end{array} & 4 & 4 & E \\ \hline \text{Règle 3/8} & \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\ \frac{1}{3} & \frac{1}{3} & 0 & 0 & 0\\ \frac{2}{3} &-\frac{1}{3} & 1 &0&0\\ 1 & 1 & -1 & 1 & 0\\ \hline & \frac{1}{8} & \frac{3}{8} & \frac{3}{8} & \frac{1}{8} \end{array} & 4 & 4 & E \\ \hline & \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\ \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0\\ 0 &-1 & 1 &0&0\\ 1 & -1 & \frac{3}{2} & \frac{1}{2} & 0\\ \hline & \frac{1}{12} & \frac{2}{3} & \frac{1}{12} & \frac{1}{6} \end{array} & 4 & 4 & E \\ \hline & \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\ 1 & 1 & 0 & 0 & 0\\ \frac{1}{2} &\frac{3}{8} & \frac{1}{8} &0&0\\ 1 & -2 & -1 & 4 & 0\\ \hline & \frac{1}{6} & \frac{1}{12} & \frac{2}{3} & \frac{1}{12} \end{array} & 4 & 4 & E \\ \hline & \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\ \frac{2}{3} & \frac{2}{3} & 0 & 0 & 0\\ \frac{1}{3} &\frac{1}{12} & \frac{1}{4} &0&0\\ 1 &-\frac{5}{4} & \frac{1}{4} & 2 & 0\\ \hline & \frac{1}{8} & \frac{3}{8} & \frac{3}{8} & \frac{1}{8} \end{array} & 4 & 4 & E \\ \hline & \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\ \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0\\ 1 &0 & 1 &0&0\\ 1 &0 & 0 & 1 & 0\\ \hline & \frac{1}{6} & \frac{2}{3} & 0 & \frac{1}{6} \end{array} & 4 & 3 & E \\ \hline & \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\ \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0\\ \frac{1}{3} &\frac{1}{3} & 0 &0&0\\ \frac{2}{3} &\frac{1}{3} & 0 & \frac{1}{3} & 0\\ \hline & 0 & -2 & \frac{3}{2} & \frac{3}{2} \end{array} & 4 & 3 & E \\ \hline & \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\ \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0\\ \frac{1}{3} &\frac{1}{3} & 0 &0&0\\ \frac{2}{3} &\frac{1}{3} & 0 & \frac{1}{3} & 0\\ \hline & 0 & -1 & 1 & 1 \end{array} & 4 & 2 & E \\ \hline \end{array} \)

3.5 Schémas à \(s=5\) étages

\( \begin{array}{|l|l|l|l|l|} \hline & \text{Matrice de Butcher} & s & \omega & \text{Exp/Imp/Semi-Implicite} \\ \hline \text{Merson} & \begin{array}{c|ccccc} 0 & 0 & 0 & 0 & 0 &0 \\ \frac{1}{3} & \frac{1}{3} & 0&0&0&0\\ \frac{1}{3} & \frac{1}{6} &\frac{1}{6}&0&0&0\\ \frac{1}{2} & \frac{1}{8} &0&\frac{3}{8}&0&0\\ 1 & \frac{1}{2} &0&-\frac{3}{2}&2&0\\ \hline & \frac{1}{6} & 0 & 0 & \frac{2}{3} & \frac{1}{6} \end{array} & 5 & 5 & E \\ \hline \end{array} \)

3.6 Schémas à \(s=6\) étages

\( \begin{array}{|l|l|l|l|l|} \hline & \text{Matrice de Butcher} & s & \omega & \text{Exp/Imp/Semi-Implicite} \\ \hline \text{Fehlberg} | & \begin{array}{c|cccccc} 0 & 0 & 0 & 0 & 0 &0&0 \\ \frac{1}{4} & \frac{1}{4} & 0&0&0&0&0\\ \frac{3}{8} & \frac{3}{32} &\frac{9}{32}&0&0&0&0\\ \frac{12}{13} & \frac{1932}{2197} &-\frac{7200}{2197}&\frac{7269}{2197}&0&0&0\\ 1 & \frac{439}{216} &-8&\frac{3680}{513}&-\frac{845}{4104}&0&0\\ \frac{1}{2} & -\frac{8}{27} &2&-\frac{3544}{2565}&\frac{1859}{4104}&-\frac{11}{40}&0\\ \hline & \frac{25}{216} & 0&\frac{1408}{2565} & \frac{2197}{4104} & -\frac{1}{5} & 0 \end{array} & 6 & 4 & E \\ \hline \text{Butcher} | & \begin{array}{c|cccccc} 0 & 0 & 0 & 0 & 0 &0&0 \\ \frac{1}{4} & \frac{1}{4} & 0&0&0&0&0\\ \frac{1}{4} & \frac{1}{8} &\frac{1}{8}&0&0&0&0\\ \frac{1}{2} & 0 &-\frac{1}{2}&1&0&0&0\\ \frac{3}{4} & \frac{3}{16} &0&0&\frac{9}{16}&0&0\\ 1 & \frac{-3}{7} &\frac{2}{7}&\frac{12}{7}&\frac{-12}{7}&\frac{8}{7}&0\\ \hline & \frac{7}{90} & 0&\frac{32}{90} & \frac{12}{90} & \frac{32}{90} & \frac{7}{90} \end{array} & 6 & 5 & E \\ \hline \text{Nystrøm} | & \begin{array}{c|cccccc} 0 & 0 & 0 & 0 & 0 &0&0 \\ \frac{1}{3} & \frac{1}{3} & 0&0&0&0&0\\ \frac{2}{5} & \frac{4}{25} &\frac{6}{25}&0&0&0&0\\ 1 & \frac{1}{4} &-3&\frac{15}{4}&0&0&0\\ \frac{2}{3} & \frac{2}{27} &\frac{10}{9}&-\frac{50}{81}&\frac{8}{81}&0&0\\ \frac{4}{5} & \frac{2}{25} &\frac{12}{25}&\frac{2}{15}&\frac{8}{75}&0&0\\ \hline & \frac{23}{192} & 0&\frac{125}{192} & 0 & -\frac{27}{64} & \frac{125}{192} \end{array} & 6 & 5 & E \\ \hline \end{array} \)

3.7 Schémas à \(s=7\) étages

\( \begin{array}{|l|l|l|l|l|} \hline & \text{Matrice de Butcher} & s & \omega & \text{Exp/Imp/Semi-Implicite} \\ \hline \text{Butcher} | & \begin{array}{c|ccccccc} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{2} & \frac{1}{2} & 0&0&0&0&0&0\\ \frac{2}{3} & \frac{2}{9} &\frac{4}{9}&0&0&0&0&0\\ \frac{1}{3} & \frac{7}{36} &\frac{8}{36}&-\frac{3}{36}&0&0&0&0\\ \frac{5}{6} & -\frac{35}{144} &-\frac{220}{144}&\frac{105}{144}&\frac{270}{144}&0&0&0\\ \frac{1}{6} & -\frac{1}{360} &-\frac{110}{360}&-\frac{45}{360}&\frac{180}{360}&\frac{36}{360}&0&0\\ 1 & -\frac{41}{160} &\frac{22}{13}&\frac{43}{156}&-\frac{118}{39}&\frac{32}{195}&\frac{80}{39}&0\\ \hline & \frac{13}{200} & 0 &\frac{11}{40} & \frac{11}{40} & \frac{4}{25} & \frac{4}{25} & \frac{13}{200} \end{array} & 7 & 6 & E \\ \hline \end{array} \)

Retour au sommet