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

Auteur·rice

Gloria Faccanoni

Date de publication

29 mars 2026

Code

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

import sympy as sp

from IPython.display import display, Math, Markdown

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

Pour \(\omega \geq 1\) il faut \( \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} \)
Pour \(\omega \geq 2\) il faut les conditions pour \(\omega \geq 1\) et \( \sum_{j=1}^s b_j c_j = \frac{1}{2} \quad\text{ i.e. }\quad \mathbf{b}^T\mathbf{c} = \frac{1}{2}. \)
Pour \(\omega \geq 3\) il faut 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} \quad\text{ i.e. }\quad \mathbf{b}^T\mathbb{A}\mathbf{c} = \frac{1}{6} \end{cases} \)
Pour \(\omega \geq 4\) il faut 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} \quad\text{ i.e. }\quad \mathbf{b}^T\mathbb{A}^2\mathbf{c} = \frac{1}{24} \end{cases} \)
Pour \(\omega \geq 5\) il faut les conditions pour \(\omega \geq 4\) et
\( \begin{cases} \displaystyle\sum_{j=1}^s b_j c_j^4 = \frac{1}{5} \\ \displaystyle\sum_{i,j=1}^s b_i c_i^2 a_{ij} c_j = \frac{1}{10} \\ \displaystyle\sum_{i=1}^s b_i \!\left(\sum_{j=1}^s a_{ij} c_j\right)^{2} = \frac{1}{20} \\ \displaystyle\sum_{i,j=1}^s b_i a_{ij} c_j^3 = \frac{1}{20} \\ \displaystyle\sum_{i,j=1}^s b_i c_i a_{ij} c_j^2 = \frac{1}{15} \\ \displaystyle\sum_{i,j,k=1}^s b_i c_i a_{ij}a_{jk} c_k = \frac{1}{30} \\ \displaystyle\sum_{i,j,k=1}^s b_i a_{ij} c_j a_{jk} c_k = \frac{1}{40} \\ \displaystyle\sum_{i,j,k=1}^s b_i a_{ij}a_{jk} c_k^2 = \frac{1}{60} \\ \displaystyle\sum_{i,j,k,\ell=1}^s b_i a_{ij}a_{jk}a_{k\ell} c_\ell = \frac{1}{120} \quad\text{ i.e. }\quad \mathbf{b}^T\mathbb{A}^3\mathbf{c} = \frac{1}{120}\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.

Étant une méthode à un pas, on peut réécrire \(u_{n+1}\) en fonction de \(u_n\) : \( u_{n+1} = q(h\beta)u_n. \) Le schéma est A-stable ssi \(|q(h\beta)|<1\).

Remarque : en introduisant les notations vectorielles \(\mathbf{K} = (K_1, \dots, K_s)^T\) et \(\mathbf{1} = (1, \dots, 1)^T\), le coefficient d’amplification se réécrit de manière compacte : \( q(\beta h) = \left( 1 - \beta h \mathbf{b}^T (\mathbb{I} + \beta h \mathbb{A})^{-1} \mathbf{1} \right) = \frac{\det(\mathbb{I}_s+\beta h\mathbb{A}-\beta h\mathbf{1}_s\mathbf{b}^T)}{\det(\mathbb{I}_s+\beta h\mathbb{A})}. \)

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

# =============================================================================
# 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
#
# =============================================================================


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

    j = sp.symbols('j')
    
    if type is not None:
        if type == "implicite":       ordre_max = 2*s
        elif type=="semi-implicite":  ordre_max = 2*s
        elif type=="explicite":       ordre_max = s
    else:
        ordre_max = 5 # par défaut, on affiche les conditions jusqu'à l'ordre 5

    if c is None: c = sp.symbols(f'c_0:{s}')    
    if b is None: b = sp.symbols(f'b_0:{s}')
    if A is None:
        if type == "implicite":       A = sp.MatrixSymbol('a', s, s)
        elif type=="semi-implicite":  A = sp.Matrix(s, s, lambda i,j: sp.Symbol('a{}{}'.format(i, j)) if i >= j else 0)
        elif type=="explicite":       A = sp.Matrix(s, s, lambda i,j: sp.Symbol('a{}{}'.format(i, j)) if i > j else 0)
    else: 
        A = sp.Matrix(A)
    
    display(Markdown("**Matrice de Butcher**"))
    But = matrice_Butcher(s, A, b, c)
    
    Eqs = []

    display(Markdown(f"**On a {s+1} 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)
    
    if ordre_max < 2:
        display(Markdown("**On ne peut pas avoir d'ordre 2 avec ce schéma**"))
    else:
        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)

    if ordre_max < 3:
        display(Markdown("**On ne peut pas avoir d'ordre 3 avec ce schéma**"))
    else:
        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)  

    if ordre_max < 4: 
        display(Markdown("**On ne peut pas avoir d'ordre 4 avec ce schéma**"))
    else:  
        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)

    if ordre_max < 5:
        display(Markdown("**On ne peut pas avoir d'ordre 5 avec ce schéma**"))
    else:
        display(Markdown("**On doit ajouter 6 conditions pour être d'ordre 5**"))
        Eq_5 = ordre_5(s, A, b, c, j)
        Eqs.append(Eq_5)
        display(*Eq_5)

    return Eqs


def matrice_Butcher(s, A, b, c):
    But = sp.Matrix(A)
    But = But.col_insert(0, sp.Matrix(c))   
    last = [sp.Symbol(" ")]
    last.extend(b)
    But = But.row_insert(s,sp.Matrix(last).T)
    display(Math(sp.latex(sp.Matrix(But))))
    return But


def ordre_1(s, A, b, c, j):
    texte = r"\sum_{j=1}^s b_j =" + f"{sum(b).simplify()}"
    texte += r"\text{ doit être égale à }1"
    # display(Math(texte))
    Eq_1 = [sp.Eq( sum(b).simplify(), 1 )]
    for i in range(s):
        somma = sp.summation(A[i,j],(j,0,s-1)).simplify()
        texte = r'\sum_{j=1}^s a_{'+str(i)+r'j}=' + sp.latex( somma )
        texte += r"\text{ doit être égale à }"+sp.latex(c[i])
        # display(Math( texte ))
        Eq_1.append(sp.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=' +sp.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 = sp.Eq( sum([b[i]*c[i] for i in range(s)]).simplify(), sp.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 += sp.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 = sp.Eq( sum([b[i]*c[i]**2 for i in range(s)]).simplify(), sp.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 + sp.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{6}"
    # display(Math(texte))
    Eq_3_2 = sp.Eq( sum([b[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify(), sp.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 += sp.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 = sp.Eq( sum([b[i]*c[i]**3 for i in range(s)]).simplify(), sp.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 + sp.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{8}"
    # display(Math(texte))
    Eq_4_2 = sp.Eq( sum([b[i]*c[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify(), sp.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 + sp.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{12}"
    # display(Math(texte))
    Eq_4_3 = sp.Eq( sum([b[i]*A[i,j]*c[j]**2 for j in range(s) for i in range(s)]).simplify(), sp.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 + sp.latex(somma)
    texte += r"\text{ doit être égale à }\frac{1}{24}"
    # display(Math(texte))
    Eq_4_4 = sp.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(), sp.Rational(1,24) )
    return [Eq_4_1, Eq_4_2, Eq_4_3, Eq_4_4]



def ordre_5(s, A, b, c, j):
    Eq_5_1 = sp.Eq(sum([b[i]*c[i]**4 for i in range(s)]).simplify(), sp.Rational(1, 5))

    somma = sum([b[i]*c[i]**2*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify()
    Eq_5_2 = sp.Eq(somma, sp.Rational(1, 10))

    somma = sum([b[i]*(sum([A[i,j]*c[j] for j in range(s)]))**2 for i in range(s)]).simplify()
    Eq_5_3 = sp.Eq(somma, sp.Rational(1, 20))

    somma = sum([b[i]*A[i,j]*c[j]**3 for j in range(s) for i in range(s)]).simplify()
    Eq_5_4 = sp.Eq(somma, sp.Rational(1, 20))

    somma = sum([b[i]*c[i]*A[i,j]*c[j]**2 for j in range(s) for i in range(s)]).simplify()
    Eq_5_5 = sp.Eq(somma, sp.Rational(1, 15))

    somma = sum([b[i]*c[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()
    Eq_5_6 = sp.Eq(somma, sp.Rational(1, 30))

    somma = sum([b[i]*A[i,j]*c[j]*A[j,k]*c[k]
                 for k in range(s) for j in range(s) for i in range(s)]).simplify()
    Eq_5_7 = sp.Eq(somma, sp.Rational(1, 40))

    somma = sum([b[i]*A[i,j]*A[j,k]*c[k]**2
                 for k in range(s) for j in range(s) for i in range(s)]).simplify()
    Eq_5_8 = sp.Eq(somma, sp.Rational(1, 60))

    somma = sum([b[i]*A[i,j]*A[j,k]*A[k,l]*c[l]
                 for l in range(s) for k in range(s)
                 for j in range(s) for i in range(s)]).simplify()
    Eq_5_9 = sp.Eq(somma, sp.Rational(1, 120))

    return [Eq_5_1, Eq_5_2, Eq_5_3, Eq_5_4, Eq_5_5, Eq_5_6, Eq_5_7, Eq_5_8, Eq_5_9]

Code

# TEST 1
# =============================================================================
# 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
# Affichage des conditions de consistance : attention aux indices qui vont de 0 à s-1
# =============================================================================
s = 5
type = "explicite" # implicite, semi-implicite, explicite
display(Markdown(f"<mark>**Conditions pour un schéma {type} avec s = {s} étages**</mark>"))
Eqs = ordre_RK(s=s, type=type)

Conditions pour un schéma explicite avec s = 5 étages

Matrice de Butcher

\(\displaystyle \left[\begin{matrix}c_{0} & 0 & 0 & 0 & 0 & 0\\c_{1} & a_{10} & 0 & 0 & 0 & 0\\c_{2} & a_{20} & a_{21} & 0 & 0 & 0\\c_{3} & a_{30} & a_{31} & a_{32} & 0 & 0\\c_{4} & a_{40} & a_{41} & a_{42} & a_{43} & 0\\ & b_{0} & b_{1} & b_{2} & b_{3} & b_{4}\end{matrix}\right]\)

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

\(\displaystyle b_{0} + b_{1} + b_{2} + b_{3} + b_{4} = 1\)

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

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

\(\displaystyle a_{20} + a_{21} = c_{2}\)

\(\displaystyle a_{30} + a_{31} + a_{32} = c_{3}\)

\(\displaystyle a_{40} + a_{41} + a_{42} + a_{43} = c_{4}\)

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

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

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

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

\(\displaystyle a_{10} b_{1} c_{0} + a_{20} b_{2} c_{0} + a_{21} b_{2} c_{1} + a_{30} b_{3} c_{0} + a_{31} b_{3} c_{1} + a_{32} b_{3} c_{2} + a_{40} b_{4} c_{0} + a_{41} b_{4} c_{1} + a_{42} b_{4} c_{2} + a_{43} b_{4} c_{3} = \frac{1}{6}\)

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

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

\(\displaystyle a_{10} b_{1} c_{0} c_{1} + a_{20} b_{2} c_{0} c_{2} + a_{21} b_{2} c_{1} c_{2} + a_{30} b_{3} c_{0} c_{3} + a_{31} b_{3} c_{1} c_{3} + a_{32} b_{3} c_{2} c_{3} + a_{40} b_{4} c_{0} c_{4} + a_{41} b_{4} c_{1} c_{4} + a_{42} b_{4} c_{2} c_{4} + a_{43} b_{4} c_{3} c_{4} = \frac{1}{8}\)

\(\displaystyle a_{10} b_{1} c_{0}^{2} + a_{20} b_{2} c_{0}^{2} + a_{21} b_{2} c_{1}^{2} + a_{30} b_{3} c_{0}^{2} + a_{31} b_{3} c_{1}^{2} + a_{32} b_{3} c_{2}^{2} + a_{40} b_{4} c_{0}^{2} + a_{41} b_{4} c_{1}^{2} + a_{42} b_{4} c_{2}^{2} + a_{43} b_{4} c_{3}^{2} = \frac{1}{12}\)

\(\displaystyle a_{10} a_{21} b_{2} c_{0} + a_{10} a_{31} b_{3} c_{0} + a_{10} a_{41} b_{4} c_{0} + a_{20} a_{32} b_{3} c_{0} + a_{20} a_{42} b_{4} c_{0} + a_{21} a_{32} b_{3} c_{1} + a_{21} a_{42} b_{4} c_{1} + a_{30} a_{43} b_{4} c_{0} + a_{31} a_{43} b_{4} c_{1} + a_{32} a_{43} b_{4} c_{2} = \frac{1}{24}\)

On doit ajouter 6 conditions pour être d’ordre 5

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

\(\displaystyle a_{10} b_{1} c_{0} c_{1}^{2} + a_{20} b_{2} c_{0} c_{2}^{2} + a_{21} b_{2} c_{1} c_{2}^{2} + a_{30} b_{3} c_{0} c_{3}^{2} + a_{31} b_{3} c_{1} c_{3}^{2} + a_{32} b_{3} c_{2} c_{3}^{2} + a_{40} b_{4} c_{0} c_{4}^{2} + a_{41} b_{4} c_{1} c_{4}^{2} + a_{42} b_{4} c_{2} c_{4}^{2} + a_{43} b_{4} c_{3} c_{4}^{2} = \frac{1}{10}\)

\(\displaystyle a_{10}^{2} b_{1} c_{0}^{2} + b_{2} \left(a_{20} c_{0} + a_{21} c_{1}\right)^{2} + b_{3} \left(a_{30} c_{0} + a_{31} c_{1} + a_{32} c_{2}\right)^{2} + b_{4} \left(a_{40} c_{0} + a_{41} c_{1} + a_{42} c_{2} + a_{43} c_{3}\right)^{2} = \frac{1}{20}\)

\(\displaystyle a_{10} b_{1} c_{0}^{3} + a_{20} b_{2} c_{0}^{3} + a_{21} b_{2} c_{1}^{3} + a_{30} b_{3} c_{0}^{3} + a_{31} b_{3} c_{1}^{3} + a_{32} b_{3} c_{2}^{3} + a_{40} b_{4} c_{0}^{3} + a_{41} b_{4} c_{1}^{3} + a_{42} b_{4} c_{2}^{3} + a_{43} b_{4} c_{3}^{3} = \frac{1}{20}\)

\(\displaystyle a_{10} b_{1} c_{0}^{2} c_{1} + a_{20} b_{2} c_{0}^{2} c_{2} + a_{21} b_{2} c_{1}^{2} c_{2} + a_{30} b_{3} c_{0}^{2} c_{3} + a_{31} b_{3} c_{1}^{2} c_{3} + a_{32} b_{3} c_{2}^{2} c_{3} + a_{40} b_{4} c_{0}^{2} c_{4} + a_{41} b_{4} c_{1}^{2} c_{4} + a_{42} b_{4} c_{2}^{2} c_{4} + a_{43} b_{4} c_{3}^{2} c_{4} = \frac{1}{15}\)

\(\displaystyle a_{10} a_{21} b_{2} c_{0} c_{2} + a_{10} a_{31} b_{3} c_{0} c_{3} + a_{10} a_{41} b_{4} c_{0} c_{4} + a_{20} a_{32} b_{3} c_{0} c_{3} + a_{20} a_{42} b_{4} c_{0} c_{4} + a_{21} a_{32} b_{3} c_{1} c_{3} + a_{21} a_{42} b_{4} c_{1} c_{4} + a_{30} a_{43} b_{4} c_{0} c_{4} + a_{31} a_{43} b_{4} c_{1} c_{4} + a_{32} a_{43} b_{4} c_{2} c_{4} = \frac{1}{30}\)

\(\displaystyle a_{10} a_{21} b_{2} c_{0} c_{1} + a_{10} a_{31} b_{3} c_{0} c_{1} + a_{10} a_{41} b_{4} c_{0} c_{1} + a_{20} a_{32} b_{3} c_{0} c_{2} + a_{20} a_{42} b_{4} c_{0} c_{2} + a_{21} a_{32} b_{3} c_{1} c_{2} + a_{21} a_{42} b_{4} c_{1} c_{2} + a_{30} a_{43} b_{4} c_{0} c_{3} + a_{31} a_{43} b_{4} c_{1} c_{3} + a_{32} a_{43} b_{4} c_{2} c_{3} = \frac{1}{40}\)

\(\displaystyle a_{10} a_{21} b_{2} c_{0}^{2} + a_{10} a_{31} b_{3} c_{0}^{2} + a_{10} a_{41} b_{4} c_{0}^{2} + a_{20} a_{32} b_{3} c_{0}^{2} + a_{20} a_{42} b_{4} c_{0}^{2} + a_{21} a_{32} b_{3} c_{1}^{2} + a_{21} a_{42} b_{4} c_{1}^{2} + a_{30} a_{43} b_{4} c_{0}^{2} + a_{31} a_{43} b_{4} c_{1}^{2} + a_{32} a_{43} b_{4} c_{2}^{2} = \frac{1}{60}\)

\(\displaystyle a_{10} a_{21} a_{32} b_{3} c_{0} + a_{10} a_{21} a_{42} b_{4} c_{0} + a_{10} a_{31} a_{43} b_{4} c_{0} + a_{20} a_{32} a_{43} b_{4} c_{0} + a_{21} a_{32} a_{43} b_{4} c_{1} = \frac{1}{120}\)

Code

# TEST 2 :  schéma de Cash-Karp, un classique explicite à 6 étages d'ordre 5 
# =============================================================================

s = 6
type = "explicite"

# Tableau de Butcher de Cash-Karp (1990) — explicite, ordre 5
A = sp.Matrix([
    [0, 0, 0, 0, 0, 0],
    [sp.Rational(1,5),       0, 0, 0, 0, 0],
    [sp.Rational(3,40),      sp.Rational(9,40),       0, 0, 0, 0],
    [sp.Rational(3,10),     -sp.Rational(9,10),       sp.Rational(6,5),       0, 0, 0],
    [-sp.Rational(11,54),    sp.Rational(5,2),        -sp.Rational(70,27),    sp.Rational(35,27),      0, 0],
    [sp.Rational(1631,55296), sp.Rational(175,512),   sp.Rational(575,13824), sp.Rational(44275,110592), sp.Rational(253,4096), 0]
])

# c_i = somme des lignes de A (condition de consistance)
c = [sum(A[i, j] for j in range(s)) for i in range(s)]
# c = [0, 1/5, 3/10, 3/5, 1, 7/8]

b = [
    sp.Rational(37, 378),
    0,
    sp.Rational(250, 621),
    sp.Rational(125, 594),
    0,
    sp.Rational(512, 1771)
]

display(Markdown(f"### Conditions pour un schéma {type} avec s = {s} étages"))
Eqs = ordre_RK(s=s, type=type, A=A, b=b, c=c)

Conditions pour un schéma explicite avec s = 6 étages

Matrice de Butcher

\(\displaystyle \left[\begin{matrix}0 & 0 & 0 & 0 & 0 & 0 & 0\\\frac{1}{5} & \frac{1}{5} & 0 & 0 & 0 & 0 & 0\\\frac{3}{10} & \frac{3}{40} & \frac{9}{40} & 0 & 0 & 0 & 0\\\frac{3}{5} & \frac{3}{10} & - \frac{9}{10} & \frac{6}{5} & 0 & 0 & 0\\1 & - \frac{11}{54} & \frac{5}{2} & - \frac{70}{27} & \frac{35}{27} & 0 & 0\\\frac{7}{8} & \frac{1631}{55296} & \frac{175}{512} & \frac{575}{13824} & \frac{44275}{110592} & \frac{253}{4096} & 0\\ & \frac{37}{378} & 0 & \frac{250}{621} & \frac{125}{594} & 0 & \frac{512}{1771}\end{matrix}\right]\)

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

\(\displaystyle \text{True}\)

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

\(\displaystyle \text{True}\)

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

\(\displaystyle \text{True}\)

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

\(\displaystyle \text{True}\)

On doit ajouter 6 conditions pour être d’ordre 5

\(\displaystyle \text{True}\)

Le calcul des conditions de consistance peut aussi s’averer fastidieux pour les méthodes à plusieurs étages. On peut aussi utiliser module sympy pour les obtenir :

Code

# =============================================================================
# Calcul de la fonction de stabilité q(x) pour un schéma de Runge-Kutta
# =============================================================================

def calculer_q_rk_general(s, A=None, b=None):
    x = sp.symbols('x', positive=True)
    # beta, h = sp.symbols('beta h')
    # x = beta * h
    
    if A is None: A = sp.Matrix(s, s, lambda i, j: sp.symbols(f'a_{i+1}{j+1}'))
    if b is None: b = sp.Matrix(s, 1, lambda i, j: sp.symbols(f'b_{i+1}'))
    
    ones = sp.Matrix.ones(s, 1)
    I = sp.eye(s)
    q = sp.det(I + x * A - x * ones * b.T) / sp.det(I + x * A)
    return sp.simplify(q)

Code

# TEST 
# =============================================================================

A = sp.Matrix([
    [0, 0, 0],
    [sp.Rational(1, 2), 0, 0],
    [0, sp.Rational(1, 2), 0]
])
b = sp.Matrix([sp.Rational(1, 6), sp.Rational(2, 3), sp.Rational(1, 6)])    
s = len(b)
q = calculer_q_rk_general(s=s, A=A, b=b)
display(Markdown(f"\(q(x) = {sp.latex(q)}\)"))

\(q(x) = - \frac{x^{3}}{24} + \frac{5 x^{2}}{12} - x + 1\)

Exercices

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

On considère la méthode RK donnée par la matrice 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. Déterminer sous quelle condition sur \(h\) la convergence est monotone.
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]\).

\( \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.

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

s = 2

b = [ sp.Rational(1,4), sp.Rational(3,4) ]
c = [ 0, sp.Rational(2,3) ]
A = [ [ 0, 0] , [ sp.Rational(2,3),  0] ]

Eqs = ordre_RK(s=s, A=A, b=b, c=c, type="explicite")

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 \text{True}\)

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

\(\displaystyle \text{True}\)

On ne peut pas avoir d’ordre 3 avec ce schéma

On ne peut pas avoir d’ordre 4 avec ce schéma

On ne peut pas avoir d’ordre 5 avec ce schéma

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

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 = sp.symbols(r'u_n')
u_np1 = sp.symbols(r'u_{n+1}')

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

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

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

Eqs = [ sp.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 = sp.solve( Eqs, kk ) # c'est un dictionnaire
KK = [sol[kk[i]] for i in range(s)]
schema = sp.Eq( u_np1 , (u_n + dt*sum([KK[j]*b[j] for j in range(s)])).factor(u_n) )
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} = u_{n} \left(\frac{\beta^{2} h^{2}}{2} - \beta h + 1\right)\)

Code

# sinon directement avec la fonction calculer_q_rk_general
b = sp.Matrix(b)
A = sp.Matrix(A)  
q = calculer_q_rk_general(s=s,A=A,b=b)
display(Markdown(f"\(q(x) = {sp.latex(q)}\)"))

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

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(sp.latex(sp.Eq(sp.Symbol('q(x)'),q.simplify()))))
display( sp.solve(abs(q)<1)  )

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

\(\displaystyle x < 2\)

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

q_np = sp.lambdify(x, q, 'numpy')

x_max = 2.5
xx = np.linspace(0, x_max, 500)
plt.figure(figsize=(10, 4))
plt.plot(xx, q_np(xx), label='q(x)', color='blue')

# Tracé des limites de stabilité (|q(x)| = 1)
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(xx, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')

# Configuration des axes
plt.xlabel(r'$x = \beta h$')
plt.ylabel(r'$q(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(0, x_max)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');


plt.tight_layout()
plt.show()

Code

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

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

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

sol_exacte = sp.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)

for N in range(10,100,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) )
    if err[-1]<1.e-10:
        break
    ax1.plot(tt,uu,label=f'Approchée avec N={N}')

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 {np.polyfit(np.log(H),np.log(err),1)[0]:1.2f}')
ax2.grid(True)

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

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

On considère la méthode RK donnée par la matri suivante : \( \def\arraystretch{2} \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]\).

\( \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.

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

s = 2
b = [ sp.Rational(1,4), sp.Rational(3,4) ]
c = [ 0, sp.Rational(2,3) ]
A = [ [ sp.Rational(1,4), -sp.Rational(1,4)] , [ sp.Rational(1,4),  sp.Rational(5,12) ] ]

Eqs = ordre_RK(s=s, type="implicite", A=A, b=b, c=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 \text{True}\)

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

\(\displaystyle \text{True}\)

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

\(\displaystyle \text{True}\)

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

\(\displaystyle \text{False}\)

On ne peut pas avoir d’ordre 5 avec ce schéma

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

Code

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

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

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

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

# Eqs = [ sp.Eq( k_0, phi(u_n+dt/4*(k_0-k_1))   ) ,
#         sp.Eq( k_1, phi(u_n+dt/4*(3*k_0+k_1)) ) ]

Eqs = [ sp.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 = sp.solve( Eqs, kk ) # c'est un dictionnaire
for key,val in sol.items():
    display(sp.Eq(key, val))

KK = [sol[kk[i]] for i in range(s)]
schema = sp.Eq( u_np1 , (u_n + dt*sum([KK[j]*b[j] for j in range(s)])).factor() )
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 k_{0} = \frac{- 4 \beta^{2} h u_{n} - 6 \beta u_{n}}{\beta^{2} h^{2} + 4 \beta h + 6}\)

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

\(\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(Markdown(f"\(q(x) = {sp.latex(q)}\)"))

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

Code

# sinon directement avec la fonction calculer_q_rk_general
b = sp.Matrix(b)
A = sp.Matrix(A)  
q = calculer_q_rk_general(s=s,A=A,b=b)
display(Markdown(f"\(q(x) = {sp.latex(q)}\)"))

\(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 \(-1<q(x)<1\) pour tout \(x>0\) :

Code

display(Markdown(f"Est-ce que $|q(x)|<1$ pour tout $x>0$ ? {sp.solve(abs(q)<1)}"))

Est-ce que \(|q(x)|<1\) pour tout \(x>0\) ? True

Code

q_np = sp.lambdify(x, q, 'numpy')

x_max = 15
xx = np.linspace(0, x_max, 500)
plt.figure(figsize=(10, 4))
plt.plot(xx, q_np(xx), label='q(x)', color='blue')

# Tracé des limites de stabilité (|q(x)| = 1)
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(xx, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')

# Configuration des axes
plt.xlabel(r'$x = \beta h$')
plt.ylabel(r'$q(x)$')
plt.ylim(-1.2, 1.2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(0, x_max)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');


plt.tight_layout()
plt.show()

Code

# Justifions ce calcul de sympy par une rapide étude de q(x)
# ==============================================================
# q(0)
display(Markdown(f"$q(0) = {q.subs(x,0)}$"))
# q(x) -> 0^- si x -> +\infty
display(Markdown( r"$\lim_{x\to +\infty} q(x) = " + sp.latex(sp.limit(q,x,sp.oo)) + "$"))
# q'(x)
q_prime = sp.diff(q,x).factor()
display(Markdown(f"$q'(x) = {sp.latex(q_prime)}$"))
# Elle admer un seul minimum en
# q est décroissante sur [0,x_min] et croissante sur [x_min,+\infty]
x_min = sp.solve(q_prime,x)
display(Markdown(f"Le minimum de q(x) est atteint en $x = {sp.latex(x_min[0])}$"))
display(Markdown(f"$q'(x)>0$ ssi ${sp.latex(sp.solve(q_prime>0))}$"))
display(Markdown(f"$q'(x)<0$ ssi ${sp.latex(sp.solve(q_prime<0))}$"))
display(Markdown(f"Valeur du minimum : $q({sp.latex(x_min[0])}) = {sp.latex(q.subs(x,x_min[0]))}$ qui est >-1"))

\(q(0) = 1\)

\(\lim_{x\to +\infty} q(x) = 0\)

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

Le minimum de q(x) est atteint en \(x = 3 + 3 \sqrt{3}\)

\(q'(x)>0\) ssi \(3 + 3 \sqrt{3} < x\)

\(q'(x)<0\) ssi \(x < 3 + 3 \sqrt{3}\)

Valeur du minimum : \(q(3 + 3 \sqrt{3}) = - \frac{6 \sqrt{3}}{18 + 12 \sqrt{3} + \left(3 + 3 \sqrt{3}\right)^{2}}\) qui est >-1

Code

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

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

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

sol_exacte = sp.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)

for N in range(10,100,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={N}')
    if err[-1]<1.e-10:
        break

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 {np.polyfit(np.log(H),np.log(err),1)[0]:1.2f}')
ax2.grid(True)

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

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

On considère la méthode RK donnée par la matri suivante : \( \def\arraystretch{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} \)

É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]\).

\( \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.

Code

s = 2

b = [ sp.Rational(3,4), sp.Rational(1,4) ]
c = [ sp.Rational(1,3) , 1 ]
A = [ [ sp.Rational(5,12), -sp.Rational(1,12)] , [ sp.Rational(3,4),  sp.Rational(1,4) ] ]

Eqs = ordre_RK(s=s, type="implicite", A=A, b=b, c=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 \text{True}\)

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

\(\displaystyle \text{True}\)

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

\(\displaystyle \text{True}\)

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

\(\displaystyle \text{False}\)

On ne peut pas avoir d’ordre 5 avec ce schéma

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

Code

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

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

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

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

Eqs = [ sp.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 = sp.solve( Eqs, kk ) # c'est un dictionnaire
KK = [sol[kk[i]] for i in range(s)]
schema = sp.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(sp.latex(sp.Eq(sp.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(sp.latex(sp.Eq(sp.Symbol('q(x)'),q.simplify()))))

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

Code

# sinon directement avec la fonction calculer_q_rk_general
b = sp.Matrix(b)
A = sp.Matrix(A)  
q = calculer_q_rk_general(s=s,A=A,b=b)
display(Markdown(f"\(q(x) = {sp.latex(q)}\)"))

\(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

sp.solve(abs(q)<1)

\(\displaystyle \text{True}\)

Code

q_np = sp.lambdify(x, q, 'numpy')

x_max = 10
xx = np.linspace(0, x_max, 500)
plt.figure(figsize=(10, 4))
plt.plot(xx, q_np(xx), label='q(x)', color='blue')

# Tracé des limites de stabilité (|q(x)| = 1)
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(xx, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')

# Configuration des axes
plt.xlabel(r'$x = \beta h$')
plt.ylabel(r'$q(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(0, x_max)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');


plt.tight_layout()
plt.show()

Code

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

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

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

sol_exacte = sp.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={N}')

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 {np.polyfit(np.log(H),np.log(err),1)[0]:1.2f}')
ax2.grid(True)

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

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

On considère la méthode RK donnée par la matrice suivante : \( \def\arraystretch{2} \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]\).

\( \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.

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

s = 3
type = "explicite"
b = [ sp.Rational(1,3), sp.Rational(1,3), sp.Rational(1,3)]
c = [ 0, sp.Rational(1,2), 1 ]
A = [[ 0, 0, 0], 
    [ sp.Rational(1,2), 0, 0],
    [ 0, 1, 0]]

Eqs = ordre_RK(s=s, type=type, A=A, b=b, c=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 \text{True}\)

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

\(\displaystyle \text{True}\)

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

\(\displaystyle \text{False}\)

\(\displaystyle \text{True}\)

On ne peut pas avoir d’ordre 4 avec ce schéma

On ne peut pas avoir d’ordre 5 avec ce schéma

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

Code

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

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

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

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

Eqs = [ sp.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 = sp.solve( Eqs, kk ) # c'est un dictionnaire
KK = [sol[kk[i]] for i in range(s)]
schema = sp.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(sp.latex(sp.Eq(sp.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(sp.latex(sp.Eq(sp.Symbol('q(x)'),q.simplify()))))

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

Code

# sinon directement avec la fonction calculer_q_rk_general
b = sp.Matrix(b)
A = sp.Matrix(A)  
q = calculer_q_rk_general(s=s,A=A,b=b)
display(Markdown(f"\(q(x) = {sp.latex(q)}\)"))

\(q(x) = - \frac{x^{3}}{6} + \frac{x^{2}}{2} - x + 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

sol = sp.solveset(q+1, x, domain=sp.S.Reals)  
display(Markdown(f"$\\beta h<{sol.args[0].evalf()}$"))

\(\beta h<2.51274532661833\)

Code

q_np = sp.lambdify(x, q, 'numpy')

x_max = 3
xx = np.linspace(0, x_max, 500)
plt.figure(figsize=(10, 4))
plt.plot(xx, q_np(xx), label='q(x)', color='blue')

# Tracé des limites de stabilité (|q(x)| = 1)
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(xx, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')

# Configuration des axes
plt.xlabel(r'$x = \beta h$')
plt.ylabel(r'$q(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(0, x_max)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');


plt.tight_layout()
plt.show()

Code

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

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

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

sol_exacte = sp.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={N}')

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 {np.polyfit(np.log(H),np.log(err),1)[0]:1.2f}')
ax2.grid(True)

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

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

On considère la méthode RK donnée par la matrice suivante : \( \def\arraystretch{2} \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]\).

\( \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.

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

s = 3
type = "implicite"

b = [ sp.Rational(1,6), sp.Rational(4,6), sp.Rational(1,6)]
c = [ 0, sp.Rational(1,2), 1 ]
A = [[ 0, 0, 0], 
    [ sp.Rational(5,24), sp.Rational(1,3), -sp.Rational(1,24)], 
    [ 0, 1, 0]]

Eqs = ordre_RK(s=s, type=type, A=A, b=b, c=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 \text{True}\)

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

\(\displaystyle \text{True}\)

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

\(\displaystyle \text{True}\)

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

\(\displaystyle \text{True}\)

\(\displaystyle \text{False}\)

On doit ajouter 6 conditions pour être d’ordre 5

\(\displaystyle \text{False}\)

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

Code

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

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

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

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

Eqs = [ sp.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 = sp.solve( Eqs, kk ) # c'est un dictionnaire
KK = [sol[kk[i]] for i in range(s)]
schema = sp.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(sp.latex(sp.Eq(sp.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(sp.latex(sp.Eq(sp.Symbol('q(x)'),q.simplify()))))

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

Code

# sinon directement avec la fonction calculer_q_rk_general
b = sp.Matrix(b)
A = sp.Matrix(A)  
q = calculer_q_rk_general(s=s,A=A,b=b)
display(Markdown(f"\(q(x) = {sp.latex(q)}\)"))

\(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

sol = sp.solveset(q+1, x, domain=sp.S.Reals)  
display(Markdown(f"$\\beta h<{sol.args[0]}$"))

\(\beta h<6\)

Code

q_np = sp.lambdify(x, q, 'numpy')

x_max = 10
xx = np.linspace(0, x_max, 500)
plt.figure(figsize=(10, 4))
plt.plot(xx, q_np(xx), label='q(x)', color='blue')

# Tracé des limites de stabilité (|q(x)| = 1)
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(xx, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')

# Configuration des axes
plt.xlabel(r'$x = \beta h$')
plt.ylabel(r'$q(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(0, x_max)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');


plt.tight_layout()
plt.show()

Code

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

phi = lambda t,y: -y

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

sol_exacte = sp.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={N}')

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 {np.polyfit(np.log(H),np.log(err),1)[0]:1.2f}')
ax2.grid(True)

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

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

On considère la méthode RK donnée par la matrice suivante : \( \def\arraystretch{2} \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]\).

\( \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.

Code

s = 3

type = "explicite"

b_0, b_1 = sp.symbols('b_0 b_1', real=True)
b = [ b_0, b_1, 0 ]
c = [ 0, sp.Rational(1,3), sp.Rational(2,3) ]
A = [[ 0, 0, 0], 
    [ sp.Rational(1,3), 0, 0],
    [ 0, sp.Rational(2,3), 0]]

Eqs = ordre_RK(s=s, type=type, A=A, b=b, c=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 b_{0} + b_{1} = 1\)

\(\displaystyle \text{True}\)

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

\(\displaystyle \frac{b_{1}}{3} = \frac{1}{2}\)

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

\(\displaystyle \frac{b_{1}}{9} = \frac{1}{3}\)

\(\displaystyle \text{False}\)

On ne peut pas avoir d’ordre 4 avec ce schéma

On ne peut pas avoir d’ordre 5 avec ce schéma

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

Code

b_1 = sp.Rational(3,2)
b_0 = 1-b_1
b = [ b_0, b_1, 0 ]
Eqs = ordre_RK(s=s, type=type, A=A, b=b, c=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 \text{True}\)

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

\(\displaystyle \text{True}\)

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

\(\displaystyle \text{False}\)

On ne peut pas avoir d’ordre 4 avec ce schéma

On ne peut pas avoir d’ordre 5 avec ce schéma

Code

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

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

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

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

Eqs = [ sp.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 = sp.solve( Eqs, kk ) # c'est un dictionnaire
KK = [sol[kk[i]] for i in range(s)]
schema = sp.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(sp.latex(sp.Eq(sp.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(sp.latex(sp.Eq(sp.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

sp.solve(abs(q)<1)

\(\displaystyle x < 2\)

Code

q_np = sp.lambdify(x, q, 'numpy')

x_max = 3
xx = np.linspace(0, x_max, 500)
plt.figure(figsize=(10, 4))
plt.plot(xx, q_np(xx), label='q(x)', color='blue')

# Tracé des limites de stabilité (|q(x)| = 1)
plt.axhline(y=1, color='black', linestyle='-', alpha=0.5)
plt.axhline(y=-1, color='black', linestyle='-', alpha=0.5)
plt.fill_between(xx, -1, 1, color='gray', alpha=0.1, label='Zone de A-Stabilité')

# Configuration des axes
plt.xlabel(r'$x = \beta h$')
plt.ylabel(r'$q(x)$')
plt.ylim(-2, 2) # On zoome sur la zone intéressante [-1, 1]
plt.xlim(0, x_max)
plt.grid(True, which='both', linestyle='--', linewidth=0.5)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');


plt.tight_layout()
plt.show()

Code

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

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

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

sol_exacte = sp.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={N}')

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 {np.polyfit(np.log(H),np.log(err),1)[0]: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}}}\)

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.

Cf. Numerical Methods for Ordinary Differential Equations, J. C. Butcher, John Wiley & Sons, 2004

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.

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

\( \def\arraystretch{2} \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} \)

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

\( \def\arraystretch{2} \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{Lobatto IIIC-2} & \begin{array}{c|cc} 0 & \frac{1}{2} & -\frac{1}{2}\\ 1 & \frac{1}{2} & \frac{1}{2}\\ \hline & \frac{1}{2} & \frac{1}{2} \end{array} & 2 & 2 & \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-I-2} & \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 & 3 & \text{SI} \\ \hline \text{Randau-Ia-2} & \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} & 2 & 3 & \text{I} \\ \hline \text{Randau-II-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{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 (=Lobatto IIIa)} & \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{SI} \\ \hline \text{Heun (=Lobatto III)} & \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} \)

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

\( \def\arraystretch{2} \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-I-3} & \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \frac{6-\sqrt{6}}{10} & \frac{9+\sqrt{6}}{75} & \frac{24+\sqrt{6}}{120}&\frac{168-73\sqrt{6}}{600}\\ \frac{6+\sqrt{6}}{10} & \frac{9-\sqrt{6}}{75} &\frac{168+73\sqrt{6}}{600}&\frac{24-\sqrt{6}}{120}\\ \hline & \frac{1}{9} &\frac{16+\sqrt{6}}{36}& \frac{16-\sqrt{6}}{36} \end{array} & 3 & 5 & I \\ \hline \text{Radau-Ia-3} & \begin{array}{c|ccc} 0 & \frac{1}{9} & \frac{-1-\sqrt{6}}{18} & \frac{-1+\sqrt{6}}{18} \\ \frac{6-\sqrt{6}}{10} & \frac{1}{9} & \frac{88+7\sqrt{6}}{360}&\frac{88-43\sqrt{6}}{360}\\ \frac{6+\sqrt{6}}{10} & \frac{1}{9} &\frac{88+43\sqrt{6}}{360}&\frac{88-7\sqrt{6}}{360}\\ \hline & \frac{1}{9} &\frac{16+\sqrt{6}}{36}& \frac{16-\sqrt{6}}{36} \end{array} & 3 & 5 & I \\ \hline \text{Radau-II-3} & \begin{array}{c|ccc} \frac{4-\sqrt{6}}{10} & \frac{24-\sqrt{6}}{120} & \frac{24-11\sqrt{6}}{120} & 0 \\ \frac{4+\sqrt{6}}{10} & \frac{24+11\sqrt{6}}{120} & \frac{24+\sqrt{6}}{120}&0\\ 1 & \frac{6-\sqrt{6}}{12} &\frac{6+6\sqrt{6}}{12}&0\\ \hline & \frac{16-\sqrt{6}}{36} &\frac{16+\sqrt{6}}{36}& \frac{1}{9} \end{array} & 3 & 5 & 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 III} & \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 & 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 IIIb} & \begin{array}{c|ccc} 0 & \frac{1}{6} & -\frac{1}{6} & 0 \\ \frac{1}{2} & \frac{1}{6} & \frac{1}{3}&0\\ 1 & \frac{1}{6} &\frac{5}{6}&0\\ \hline & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \end{array} & 3 & 4 & I \\ \hline \text{Lobatto IIIc} & \begin{array}{c|ccc} 0 & \frac{1}{6} & \frac{-1}{3} & \frac{1}{6} \\ \frac{1}{2} & \frac{1}{6} & \frac{5}{12}&\frac{-1}{12}\\ 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{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 (RK32)} & \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 \text{RK31} & \begin{array}{c|ccc} 0 & 0 & 0 & 0 \\ \frac{2}{3} & \frac{2}{3} & 0&0\\ \frac{2}{3} & \frac{1}{3} &\frac{1}{3}&0\\ \hline & \frac{1}{4} & 0 & \frac{3}{4} \end{array} & 3 & 3 & E \\ \hline \end{array} \)

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

\( \def\arraystretch{2} \begin{array}{|l|l|l|l|l|} \hline & \text{Matrice de Butcher} & s & \omega & \text{Exp/Imp/Semi-Implicite} \\ \hline \text{Lobatto III} & \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 \text{Lobatto IIIa} & \begin{array}{c|cccc} 0 & 0 & 0 & 0 & 0\\ \frac{5-\sqrt{5}}{10} & \frac{11+\sqrt{5}}{120} & \frac{25-\sqrt{5}}{120} & \frac{25-13\sqrt{5}}{120} & \frac{-1+\sqrt{5}}{120}\\ \frac{5+\sqrt{5}}{10} & \frac{11-\sqrt{5}}{120} &\frac{25+13\sqrt{5}}{120} &\frac{25+\sqrt{5}}{120}&\frac{-1-\sqrt{5}}{120}\\ 1 & \frac{1}{12} & \frac{5}{12} & \frac{5}{12}& \frac{1}{12} \\ \hline & \frac{1}{12} & \frac{5}{12} & \frac{5}{12} & \frac{1}{12} \end{array} & 4 & 6 & I \\ \hline \text{Lobatto IIIb} & \begin{array}{c|cccc} 0 & \frac{1}{12} & \frac{-1-\sqrt{5}}{24} & \frac{-1+\sqrt{5}}{24} & 0\\ \frac{5-\sqrt{5}}{10} & \frac{1}{12} & \frac{25+\sqrt{5}}{120} & \frac{25-13\sqrt{5}}{120} & 0\\ \frac{5+\sqrt{5}}{10} & \frac{1}{12} &\frac{25+13\sqrt{5}}{120} &\frac{25-\sqrt{5}}{120}&0\\ 1 & \frac{1}{12} & \frac{11-\sqrt{5}}{120} & \frac{11+\sqrt{5}}{120} & 0 \\ \hline & \frac{1}{12} & \frac{5}{12} & \frac{5}{12} & \frac{1}{12} \end{array} & 4 & 6 & I \\ \hline \text{Lobatto IIIc} & \begin{array}{c|cccc} 0 & \frac{1}{12} & \frac{-\sqrt{5}}{12} & \frac{\sqrt{5}}{12} & -\frac{1}{12}\\ \frac{5-\sqrt{5}}{10} & \frac{1}{12} & \frac{1}{4} & \frac{10-7\sqrt{5}}{60} & \frac{\sqrt{5}}{60}\\ \frac{5+\sqrt{5}}{10} & \frac{1}{12} &\frac{10+7\sqrt{5}}{60} &\frac{1}{4}&-\frac{\sqrt{5}}{60}\\ 1 & \frac{1}{12} & \frac{5}{12} & \frac{5}{12} & \frac{1}{12} \\ \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} \)

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

\( \def\arraystretch{2} \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} \)

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

\( \def\arraystretch{2} \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{RK5} & \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 &0&\frac{1}{2}&0&0&0\\ \frac{3}{4} & \frac{3}{16} &-\frac{3}{8}&\frac{3}{8}&\frac{9}{16}&0&0\\ 1 & \frac{-3}{7} &\frac{8}{7}&\frac{6}{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 & \begin{array}{c|cccccc} 0 & 0 & 0 & 0 & 0 &0&0 \\ \frac{1}{2} & \frac{1}{2} & 0&0&0&0&0\\ 1 & \frac{-9}{4} &\frac{13}{4}&0&0&0&0\\ \frac{1}{4} & \frac{9}{64} & \frac{5}{32} &\frac{-3}{64}&0&0&0\\ \frac{7}{10} & \frac{63}{625} & \frac{259}{2500}&\frac{231}{2500}&\frac{252}{625}&0&0\\ 1 & \frac{-27}{50} &\frac{-139}{50}&\frac{-21}{50}&\frac{56}{25}&\frac{5}{2}&0\\ \hline & \frac{1}{14} & 0& 0 & \frac{32}{81} & \frac{250}{567} & \frac{5}{54} \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} \)

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

\( \def\arraystretch{2} \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} \)