In [1]:
from IPython.display import display, Latex
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read()) 
Out[1]:
In [2]:
import sys #only needed to determine Python version number
print('Python version ' + sys.version)
Python version 3.8.5 (default, Jan 27 2021, 15:41:15) 
[GCC 9.3.0]

M62_TD2 : Exercices schémas multipas et Runge-Kutta

Rappels schémas multipas

Considérons une méthode numérique à $q=p+1$ étapes (=pas) linéaire. Elle s'écrit sous la forme générale $$ u_{n+1} = \sum_{j=0}^p a_ju_{n-j} + h\sum_{j=0}^p b_j\varphi(t_{n-j},u_{n-j}) + hb_{-1}\varphi(t_{n+1},u_{n+1}), \qquad n=p,p+1,\dots,N-1 $$ où les $\{a_k\}$ et $\{b_k\}$ sont des coefficients donnés et $p\ge0$ un entier. En notant $\varphi_{n-j}\stackrel{\text{déf}}{=}\varphi(t_{n-j},u_{n-j})$ et en écrivant explicitement les séries, on a $$ u_{n+1} = a_0u_{n}+a_1u_{n-1}+\cdots+a_pu_{n-p} + h b_0\varphi_{n}+h b_1\varphi_{n-1}+\cdots+h b_{p}\varphi_{n-p} + hb_{-1}\varphi_{n+1}. $$
  1. La méthode est consistante ssi $$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=1, \\ \displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=1 \end{cases} $$

  2. La méthode est zéro-stable ssi $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1. \end{cases} $$ où $\varrho$ est le premier polynôme caractéristique défini par $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j} $$

  3. Consistance + zéro-stabilité = convergence

  4. Une méthode convergente est d'ordre $\omega\ge1$ si, de plus, $$ \sum_{j=0}^p (-j)^{i}a_j+i\sum_{j=-1}^p (-j)^{i-1}b_j=1 \quad i=2,\dots,\omega. $$

  5. 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}\xrightarrow[t\to+\infty]{}0.$ 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.

Exercice : consistance

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. La méthode multipas suivante est-elle consistante? $$ u_{n+1}=u_{n-1}+ \dfrac{h}{4}(3\varphi_{n}-\varphi_{n-1}) $$

Correction

C'est une méthode à 2 pas explicite :

  • $p=1$
  • $b_{-1}=0$
  • $a_0=0$ et $a_1=1$
  • $b_0=\frac{3}{4}$ et $b_1=-\frac{1}{4}$

La méthode n'est pas consistante car $$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=1, \\ \displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-1+\frac{3}{4}-\frac{1}{4}\neq1 \end{cases} $$

Exercice : convergence

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. La méthode multipas suivante est-elle convergente? $$ u_{n+1}=-4u_{n}+5u_{n-1}+ h(4\varphi_{n+1}+2\varphi_{n}) $$

Correction

C'est une méthode à 2 pas implicite :

  • $p=1$
  • $b_{-1}=4$
  • $a_0=-4$ et $a_1=5$
  • $b_0=2$ et $b_1=0$

La méthode est consistante car $$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=-4+5=1, \\ \displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-\big(0\times a_0+1\times a_1\big)+\big(b_{-1}+b_0+b_1\big)=-5+4+2+0=1. \end{cases} $$

Le premier polynôme caractéristique est $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{n-j}=r^2-\big(a_0 r^1+a_1r^0\big)=r^2-(-4r+5)=r^2+4r-5=(r-1)(r+5) $$ dont les racines sont $r_0=1$ et $r_1=-5$. La méthode n'est pas zéro-stable car $|r_1|>1$, donc la méthode ne converge pas.

Exercice : convergence

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. La méthode multipas suivante est-elle convergente? $$ u_{n+1}=-u_{n}+u_{n-1}+u_{n-2}+ 4h\varphi_{n} $$

Correction

C'est une méthode à 3 pas explicite :

  • $p=2$
  • $b_{-1}=0$
  • $a_0=-1$, $a_1=1$ et $a_2=1$
  • $b_0=4$, $b_1=0$ et $b_2=0$

La méthode est consistante car $$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=-1+1+1=1, \\ \displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-(0\times(-1)+1\times1+2\times1)+(0+4+0+0)=-3+4=1. \end{cases} $$

Le premier polynôme caractéristique est $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{n-j}=r^3-\big(-r^2+r+1\big)=r^3+r^2-r-1=(r-1)(r^2-1)=(r-1)^2(r+1) $$ dont les racines sont $r_0=1$ de multiplicité $2$ et $r_1=-1$. La méthode n'est pas zéro-stable car $\varrho'(r_0)\neq0$, donc la méthode ne converge pas.

In [3]:
import sympy as sym
sym.var('r')
rho=r**3+r**2-r-1
sym.factor(rho)
Out[3]:
$\displaystyle \left(r - 1\right) \left(r + 1\right)^{2}$

Exercice : schéma à deux pas d'ordre 4

La première barrière de Dahlquist affirme qu'un schéma implicite à $q=2$ étapes consistante et zéro-stable peut être d'ordre $\omega=q+2=4$.

Construire un schéma implicite à $q=2$ étapes consistante et zéro-stable d'ordre $\omega=q+2=4$.

Correction

Un schéma implicite à $q=2$ étapes s'écrit $$ u_{n+1}=a_0u_n+a_1u_{n-1}+h\big(b_{-1}\varphi_{n+1}+b_0\varphi_n+b_1\varphi_{n-1}\big) $$

La méthode est d'ordre 4 ssi $$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=1, \\ \displaystyle \sum_{j=0}^p (-j)a_j+\sum_{j=-1}^p b_j=1 \\ \displaystyle \sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j=1 \\ \displaystyle \sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j=1 \\ \displaystyle \sum_{j=0}^p (-j)^{4}a_j+4\sum_{j=-1}^p (-j)^{3}b_j=1 \end{cases} $$ ssi $$ \begin{cases} a_0+a_1=1, \\ {}-\big(0a_0+1a_1\big)+\big(b_{-1}+b_0+b_1\big)=1 \\ \big(0^2a_0+1^2a_1\big)-2\big(-1b_{-1}+0b_0+1b_1\big)=1 \\ {}-\big(0^3a_0+1^3a_1\big)+3\big((-1)^2b_{-1}+0^2b_0+1^2b_1\big)=1 \\ \big(0^4a_0+1^4a_1\big)-4\big((-1)^3b_{-1}+0^3b_0+1^3b_1\big)=1 \end{cases} $$ ssi $$ \begin{pmatrix} 1&1&0&0&0\\ 0&-1&1&1&1\\ 0&1&2&0&-2\\ 0&-1&3&0&3\\ 0&1&4&0&-4 \end{pmatrix} \begin{pmatrix} a_0\\a_1\\b_{-1}\\b_0\\b_1 \end{pmatrix} {}= \begin{pmatrix} 1\\1\\1\\1\\1 \end{pmatrix} $$

In [4]:
import sympy as sym
sym.init_printing()

sym.var('a_0,a_1,b_m1,b_0,b_1')

eq1 = sym.Eq(a_0 + a_1, 1) # eq1 = a_0+a_1-1
eq2 = sym.Eq(-a_1 + b_m1 + b_0 + b_1, 1)
eq3 = sym.Eq(a_1 + 2 * b_m1 - 2 * b_1, 1)
eq4 = sym.Eq(-a_1 + 3 * b_m1 + 3 * b_1, 1)
eq5 = sym.Eq(a_1 + 4 * b_m1 - 4 * b_1, 1)
sym.solve([eq1, eq2, eq3, eq4, eq5])
Out[4]:
$\displaystyle \left\{ a_{0} : 0, \ a_{1} : 1, \ b_{0} : \frac{4}{3}, \ b_{1} : \frac{1}{3}, \ b_{m1} : \frac{1}{3}\right\}$

On reconnait la méthode MS$_2$.

Le premier polynôme caractéristique est $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^2-a_0r-a_1=(r-1)(r+1) $$ dont les racines sont $r_0=1$ et $r_1=-1$. La méthode est zéro-stable car $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1, \end{cases} $$

Exercice : convergence

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Soit la méthode multipas $$ u_{n+1}=(1-\gamma)u_{n}+\gamma u_{n-1}+ \frac{h}{4}((\gamma+3)\varphi_{n+1}+(3\gamma+1)\varphi_{n-1}). $$

  1. Pour quelles valeurs de $\gamma$ est-elle consistante?
  2. Pour quelles valeurs de $\gamma$ est-elle zéro-stable?
  3. Pour quelles valeurs de $\gamma$ est-elle convergente?
  4. Pour quelles valeurs de $\gamma$ est-elle convergente d'ordre $2$? et d'ordre $3$?
  5. Pour $\gamma=\frac{1}{2}$ vérifier empiriquement la convergence sur le problème de Cauchy $$ \begin{cases} y'(t) = -4y(t)+t^2, &\forall t \in I=[0,1],\\ y(0) = 1 \end{cases} $$

Correction

C'est une méthode à 2 pas implicite :

  • $p=1$
  • $b_{-1}=\frac{\gamma+3}{4}$
  • $a_0=1-\gamma$ et $a_1=\gamma$
  • $b_0=0$ et $b_1=\frac{3\gamma+1}{4}$
  1. La méthode est consistante pour tout $\gamma$ car $$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=1-\gamma+\gamma=1 \text{ pour tout }\gamma, \\ \displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-\big(0\times(1-\gamma)+1\times \gamma\big)+\big(\frac{\gamma+3}{4}+0+\frac{3\gamma+1}{4}\big)=1 \text{ pour tout }\gamma \end{cases} $$

  2. Le premier polynôme caractéristique est $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^2-(1-\gamma)r-\gamma=(r-1)(r+\gamma) $$ dont les racines sont $r_0=1$ et $r_1=-\gamma$. La méthode est zéro-stable ssi $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1, \end{cases} $$ donc ssi $-1\le -\gamma\le1$ et $-\gamma\neq1$, i.e. ssi $-1< \gamma\le 1$.

  3. La méthode est convergente ssi elle est consistante et zéro-stable donc ssi $-1< \gamma\le 1$.

  4. La méthode est d'ordre au moins 2 pour tout $\gamma$ car $$ \sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j {}= \big(0^2a_0+1^2a_1\big)-2\big(-1b_{-1}+0b_0+1b_1\big) {}= \gamma-2\big(-\frac{\gamma+3}{4}+\frac{3\gamma+1}{4}\big) {}=1 $$ Pour que la méthode soit d'ordre 3 il faut que $$ 1 {}= \sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j {}= {}-\big(0^3a_0+1^3a_1\big)+3\big((-1)^2b_{-1}+0^2b_0+1^2b_1\big) {}= {}-\gamma+3\big(\frac{\gamma+3}{4}+\frac{3\gamma+1}{4}\big) {}= 2\gamma+3 $$ donc ssi $\gamma=-1$. Cependant cette valeur n'appartient pas au domain de zéro-stabilité.
    On conclut que la méthode est convergente d'ordre $2$ pour $-1< \gamma\le 1$ et non convergente sinon.

  5. On définit la solution exacte (calculée en utilisant le module sympy) pour estimer les erreurs:

In [5]:
%reset -f
%matplotlib inline

from matplotlib.pylab import *

# variables globales
t0     = 0
tfinal = 1
y0     = 1

phi = lambda t,y : -4*y+t**2

gamma=1/2
In [6]:
import sympy as sym
sym.init_printing()

t  = sym.Symbol('t')
y  = sym.Function('y')
edo= sym.Eq( sym.diff(y(t),t) , phi(t,y(t)) )
display(edo)
solgen = sym.dsolve(edo)
display(solgen)

consts = sym.solve( sym.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
display(consts)
solpar=solgen.subs(consts).simplify()
display(solpar)

sol_exacte = sym.lambdify(t,solpar.rhs,'numpy')
$\displaystyle \frac{d}{d t} y{\left(t \right)} = t^{2} - 4 y{\left(t \right)}$
$\displaystyle y{\left(t \right)} = \left(C_{1} + \frac{\left(8 t^{2} - 4 t + 1\right) e^{4 t}}{32}\right) e^{- 4 t}$
$\displaystyle \left\{ C_{1} : \frac{31}{32}\right\}$
$\displaystyle y{\left(t \right)} = \frac{t^{2}}{4} - \frac{t}{8} + \frac{1}{32} + \frac{31 e^{- 4 t}}{32}$

On calcule la solution approchée pour différentes valeurs de $N$:

In [7]:
from scipy.optimize import fsolve

def multipas(phi, tt,sol_exacte):
    h  = tt[1] - tt[0]
    uu = [sol_exacte(tt[0])]
    uu.append(sol_exacte(tt[1]))
    for i in range(1,len(tt) - 1):
        temp = fsolve(lambda x: -x + (1-gamma)*uu[i]+gamma*uu[i-1] + h/4 * ( (gamma+3)*phi(tt[i+1],x)+(3*gamma+1)*phi(tt[i-1],uu[i-1]) ), uu[i])
        uu.append(temp)
    return uu


H = []
err_mp = []
N = 10
for k in range(7):
    N+=20
    tt = linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = [sol_exacte(t) for t in tt]
    uu_mp = multipas(phi,tt,sol_exacte)
    H.append(h)
    err_mp.append(max([abs(uu_mp[i] - yy[i]) for i in range(len(yy))]))

print ('Multipas %1.2f' %(polyfit(log(H),log(err_mp), 1)[0]))

figure(figsize=(8,5))
plot(log(H), log(err_mp), 'r-o', label='Multipas')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
legend(bbox_to_anchor=(1.04, 1), loc='upper left')
grid(True)
show()
Multipas 1.97

Exercice : convergence

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Soit la méthode multipas $$ u_{n+1} {}=\left(\frac{1}{3}b^2+\frac{1}{2}b\right)u_{n} {}+\left(\frac{2}{3}b^2+\frac{1}{2}b-1\right)u_{n-1} {}+h\left(\left(-\frac{1}{6}b^2-\frac{1}{4}b+2\right)\varphi_{n} {} +\left(-\frac{1}{6}b^2-\frac{1}{4}b\right)\varphi_{n-1}\right). $$

  1. Pour quelles valeurs de $b$ est-elle consistante?
  2. Est-elle convergente pour ces valeurs?
  3. Quel ordre de convergente a-t-elle?
  4. Pour les valeurs de $b$ qui donnent la convergence, vérifier empiriquement l'ordre de convergence sur le problème de Cauchy $$ \begin{cases} y'(t) = -10\big(y(t)-1\big)^2, &\forall t \in I=[0,1],\\ y(0) = 2. \end{cases} $$

Correction

C'est une méthode à $q=2$ pas explicite :

  • $p=1$
  • $b_{-1}=0$
  • $a_0=\frac{1}{3}b^2+\frac{1}{2}b$ et $a_1=\frac{2}{3}b^2+\frac{1}{2}b-1$
  • $b_0=-\frac{1}{6}b^2-\frac{1}{4}b+2$ et $b_1=-\frac{1}{6}b^2-\frac{1}{4}b$
  1. La méthode est consistante ssi $$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=b^2+b-1=1, \\ \displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-\big(0a_0+1a_1\big)+\big(b_{-1}+b_0+b_1\big)=b^2+b-1=1 \end{cases} $$ donc ssi $b=-2$ ou $b=1$.

  2. Le premier polynôme caractéristique est $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^2-a_0r-a_1 {}= \begin{cases} r^2-\frac{1}{3}r-\frac{2}{3}&\text{si }b=-2, \\ r^2-\frac{5}{6}r-\frac{1}{6}&\text{si }b=1, \end{cases} $$ dont les racines sont $$ \begin{cases} r_0=1,\ r_1=-\frac{2}{3}&\text{si }b=-2, \\ r_0=1,\ r_1=-\frac{1}{6}&\text{si }b=1, \end{cases} $$ La méthode est zéro-stable ssi $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1, \end{cases} $$ donc elle est vérifiée pour les deux cas.

    La méthode est convergente ssi elle est consistante et zéro-stable donc ssi $b=-2$ ou $b=1$.

  3. La méthode est d'ordre au moins 2 car $$ \sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j {}= \begin{cases} \frac{2}{3}+2\frac{1}{6}=1&\text{si }b=-2, \\ \frac{1}{6}+2\frac{5}{12}=1&\text{si }b=1, \end{cases} $$ La méthode n'est pas d'ordre 3 car $$ \sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j {}= \begin{cases} {}-\frac{2}{3}+3\frac{-1}{6}\neq1&\text{si }b=-2, \\ {}-\frac{1}{6}+3\frac{-5}{12}\neq1&\text{si }b=1, \end{cases} $$

  4. On définit la solution exacte (calculée en utilisant le module sympy) pour estimer les erreurs:

In [8]:
%reset -f
%matplotlib inline

from matplotlib.pylab import *

# variables globales
t0     = 0
tfinal = 1
y0     = 2

phi = lambda t,y : -10*(y-1)**2
In [9]:
import sympy as sym
sym.init_printing()

t  = sym.Symbol('t')
y  = sym.Function('y')
edo= sym.Eq( sym.diff(y(t),t) , phi(t,y(t)) )
display(edo)
solgen = sym.dsolve(edo)
display(solgen)

consts = sym.solve( sym.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
display(consts)
solpar=solgen.subs(consts).simplify()
display(solpar.apart())

sol_exacte = sym.lambdify(t,solpar.rhs,'numpy')
$\displaystyle \frac{d}{d t} y{\left(t \right)} = - 10 \left(y{\left(t \right)} - 1\right)^{2}$
$\displaystyle y{\left(t \right)} = \frac{C_{1} - t - \frac{1}{10}}{C_{1} - t}$
$\displaystyle \left\{ C_{1} : - \frac{1}{10}\right\}$
$\displaystyle y{\left(t \right)} = 1 + \frac{1}{10 t + 1}$

On calcule la solution approchée pour différentes valeurs de $N$:

In [10]:
from scipy.optimize import fsolve

def multipas(b, phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = [sol_exacte(tt[0])]
    uu.append(sol_exacte(tt[1]))
    for i in range(1,len(tt) - 1):
        uu.append( (b**2/3+b/2)*uu[i]+(2*b**2/3+b/2-1)*uu[i-1]+h*(-b**2/6-b/4+2)*phi(tt[i],uu[i])+h*(-b**2/6-b/4)*phi(tt[i-1],uu[i-1]) )
    return uu


H = []
err_mp1 = []
err_mp2 = []
N = 10
for k in range(7):
    N+=20
    tt = linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = [sol_exacte(t) for t in tt]
    uu_mp1 = multipas( 1, phi, tt, sol_exacte)
    uu_mp2 = multipas(-2, phi, tt, sol_exacte)
    H.append(h)
    err_mp1.append(max([abs(uu_mp1[i] - yy[i]) for i in range(len(yy))]))
    err_mp2.append(max([abs(uu_mp2[i] - yy[i]) for i in range(len(yy))]))

print ('Multipas b= 1 ordre = %1.2f' %(polyfit(log(H),log(err_mp1), 1)[0]))
print ('Multipas b=-2 ordre = %1.2f' %(polyfit(log(H),log(err_mp2), 1)[0]))

figure(figsize=(8,5))
plot(log(H), log(err_mp1), 'r-o', label='Multipas b=1')
plot(log(H), log(err_mp2), 'b-d', label='Multipas b=-2')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
legend(bbox_to_anchor=(1.04, 1), loc='upper left')
grid(True)
show()
Multipas b= 1 ordre = 2.01
Multipas b=-2 ordre = 2.22

Exercice : convergence

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Soit la méthode multipas $$ u_{n+1} =\left(b^2-\frac{9}{2}b+\frac{11}{2}\right)u_{n} +\left(b^2-\frac{7}{2}b+\frac{3}{2}\right)u_{n-1} +h\left(\frac{b}{2}(b-2)\varphi_{n} +\frac{b}{2}\left(b-\frac{10}{3}\right)\varphi_{n-1}\right). $$

  • Pour quelles valeurs de $b$ est-elle consistante?
  • Est-elle convergente pour ces valeurs?
  • Quel ordre de convergente a-t-elle?
  • Pour les valeurs de $b$ qui donnent la convergence, vérifier empiriquement l'ordre de convergence sur le problème de Cauchy $$ \begin{cases} y'(t) = y(t)(1-y(t)), &\forall t \in I=[0,1],\\ y(0) = 2. \end{cases} $$

Correction

C'est une méthode à $q=2$ pas explicite :

  • $p=1$
  • $b_{-1}=0$
  • $a_0=b^2-\frac{9}{2}b+\frac{11}{2}$ et $a_1=b^2-\frac{7}{2}b+\frac{3}{2}$
  • $b_0=\frac{b}{2}(b-2)$ et $b_1=\frac{b}{2}\left(b-\frac{10}{3}\right)$
  • La méthode est consistante ssi $$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=2b^2-8b+7=1, \\ \displaystyle-\sum_{j=0}^p ja_j+\sum_{j=-1}^p b_j=-\big(0a_0+1a_1\big)+\big(b_{-1}+b_0+b_1\big)=\frac{5}{6}b-\frac{19}{6}=1 \end{cases} $$ donc ssi $b=3$. Vérifions ces calculs ci-dessous:
In [11]:
import sympy as sym
sym.init_printing()

sym.var('b')
p=1
bm1=0
a0=b**2-sym.Rational(9,2)*b+sym.Rational(11,2)
a1=b**2-sym.Rational(7,2)*b+sym.Rational(3,2)
b0=sym.Rational(1,2)*b*(b-2)
b1=sym.Rational(1,2)*b*(b-sym.Rational(10,3))

cond1=a0+a1
cond2=-(0*a0+1*a1)+(bm1+b0+b1)
display(cond1)
display(cond2.factor())
display(sym.solve([sym.Eq(cond1,1),sym.Eq(cond2,1)],b))
$\displaystyle 2 b^{2} - 8 b + 7$
$\displaystyle \frac{5 b - 9}{6}$
$\displaystyle \left[ \left( 3\right)\right]$

Dans la suite on considerera $b=3$.

In [12]:
from IPython.display import display, Math
display(Math('a_0='+sym.latex(a0.subs(b,3))))
display(Math('a_1='+sym.latex(a1.subs(b,3))))
display(Math('b_0='+sym.latex(b0.subs(b,3))))
display(Math('b_1='+sym.latex(b1.subs(b,3))))
$\displaystyle a_0=1$
$\displaystyle a_1=0$
$\displaystyle b_0=\frac{3}{2}$
$\displaystyle b_1=- \frac{1}{2}$
  • Le premier polynôme caractéristique est $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^2-a_0r-a_1 \stackrel{b=3}{=} r^2-r $$ dont les racines sont $$ r_0=1,\ r_1=0 $$ La méthode est zéro-stable ssi $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1, \end{cases} $$ ce qui est bien vérifié.

    La méthode est convergente ssi elle est consistante et zéro-stable donc ssi $b=3$.

  • La méthode est d'ordre au moins 2 car $$ \sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j {}= \big(0^2a_0+1^2a_1\big)+2\big(-(-1)b_{-1}+0b_0+1b_1\big) \stackrel{b=3}{=} 1 $$ La méthode n'est pas d'ordre 3 car $$ \sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j {}= {}-a_1+3\big(b_{-1}+b_1\big) \stackrel{b=3}{\neq}1 $$
  • On définit la solution exacte (calculée en utilisant le module sympy) pour estimer les erreurs:
In [13]:
%reset -f
%matplotlib inline

from matplotlib.pylab import *

# variables globales
t0     = 0
tfinal = 1
y0     = 2

phi = lambda t,y : y*(1-y)
In [14]:
import sympy as sym
sym.init_printing()

t  = sym.Symbol('t')
y  = sym.Function('y')
edo= sym.Eq( sym.diff(y(t),t) , phi(t,y(t)) )
display(edo)
solgen = sym.dsolve(edo)
display(solgen)

consts = sym.solve( sym.Eq( y0, solgen.rhs.subs(t,t0)) , dict=True)[0]
display(consts)
solpar=solgen.subs(consts).simplify()
display(solpar.apart())

sol_exacte = sym.lambdify(t,solpar.rhs,'numpy')
$\displaystyle \frac{d}{d t} y{\left(t \right)} = \left(1 - y{\left(t \right)}\right) y{\left(t \right)}$
$\displaystyle y{\left(t \right)} = \frac{1}{C_{1} e^{- t} + 1}$
$\displaystyle \left\{ C_{1} : - \frac{1}{2}\right\}$
$\displaystyle y{\left(t \right)} = 1 + \frac{1}{2 e^{t} - 1}$

On calcule la solution approchée pour différentes valeurs de $N$:

In [15]:
from scipy.optimize import fsolve

def multipas(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = [sol_exacte(tt[0])]
    uu.append(sol_exacte(tt[1]))
    for i in range(1,len(tt) - 1):
        uu.append( uu[i]+h*3/2*phi(tt[i],uu[i])-h/2*phi(tt[i-1],uu[i-1]) )
    return uu


H = []
err_mp = []
N = 10
for k in range(7):
    N+=20
    tt = linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = [sol_exacte(t) for t in tt]
    uu_mp = multipas(phi, tt, sol_exacte)
    H.append(h)
    err_mp.append(max([abs(uu_mp[i] - yy[i]) for i in range(len(yy))]))

print ('Multipas b= 3 ordre=%1.2f' %(polyfit(log(H),log(err_mp), 1)[0]))

figure(figsize=(8,5))
plot(log(H), log(err_mp), 'r-o', label='Multipas b=3')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
legend(bbox_to_anchor=(1.04, 1), loc='upper left')
grid(True)
show()
Multipas b= 3 ordre=1.98

Exercice : convergence predictor-corrector

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Soit les méthodes multipas \begin{align*} u_{n+1}&=u_{n}+\frac{h}{12}\left(23\varphi_{n}-16\varphi_{n-1}+5\varphi_{n-2}\right)&\text{[P]} \\ u_{n+1}&=u_{n}+\frac{h}{24}\left(9\varphi_{n+1}+19\varphi_{n}-5\varphi_{n-1}+\varphi_{n-2}\right)&\text{[C]} \end{align*}

  • Étudier consistance, ordre et zéro-stabilité de la méthode P
  • Étudier consistance, ordre et zéro-stabilité de la méthode C
  • Soit la méthode PC suivante: $$ \begin{cases} \tilde u=u_{n}+\frac{h}{12}\left(23\varphi_{n}-16\varphi_{n-1}+5\varphi_{n-2}\right) \\ u_{n+1}=u_{n}+\frac{h}{24}\left(9\varphi(t_{n+1},\tilde u)+19\varphi_{n}-5\varphi_{n-1}+\varphi_{n-2}\right) \end{cases} $$ Calculer empiriquement l'ordre de convergence de la méthode PC sur le problème de Cauchy $$ \begin{cases} y'(t) = \frac{1}{1+t^2}-2y^2(t), &\forall t \in I=[0,1],\\ y(0) = 0, \end{cases} $$ dont la solution exacte est $y(t)=\frac{t}{1+t^2}$.

Correction schéma P

P est une méthode à $q=3$ pas explicite :

  • $p=2$
  • $b_{-1}=0$
  • $a_0=1$ et $a_1=a_2=0$
  • $b_0=\frac{23}{12}$, $b_1=\frac{-16}{12}$ et $b_2=\frac{5}{12}$
  • La méthode est consistante d'ordre $\omega=3$ car $$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=1, \\ \displaystyle\sum_{j=0}^p -ja_j+\sum_{j=-1}^p b_j=-a_1-2a_2+\big(b_{-1}+b_0+b_1+b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j=a_1+4a_2-2\big(-b_{-1}+b_1+2b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j=-a_1-8a_2+3\big(b_{-1}+b_1+4b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{4}a_j+4\sum_{j=-1}^p (-j)^{3}b_j=a_1+16a_2-4\big(-b_{-1}+b_1+8b_2\big)\neq1 \end{cases} $$ Vérifions ces calculs ci-dessous:
In [16]:
import sympy as sym
sym.init_printing()
from IPython.display import display, Math
 
p=1
bm1=0
a0=1
a1=0
a2=0
b0=sym.Rational(23,12)
b1=sym.Rational(-16,12)
b2=sym.Rational(5,12)


aa=[a0,a1,a2]
bb=[b0,b1,b2]

cond=[]
cond.append(sum(aa)) 
cond.append(sum([-j*aa[j]+bb[j] for j in range(len(aa))])+bm1)
if cond==[1,1]:
    display(Math("\omega\ge1"))
    for n in range(2,9):
        cond.append(sum( [ (-j)**n*aa[j]+n*(-j)**(n-1)*bb[j] for j in range(len(aa)) ])+n*bm1)
        if cond[n]==1:
            display(Math("\omega\ge"+str(n)))
        else:
            display(Math("\omega\le"+str(n)))
            break
$\displaystyle \omega\ge1$
$\displaystyle \omega\ge2$
$\displaystyle \omega\ge3$
$\displaystyle \omega\le4$
  • Le premier polynôme caractéristique est $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^3-a_0r^2-a_1r-a_2 {}= r^3-r^2=r^2(r-1) $$ dont les racines sont $$ r_0=1,\ r_1=r_2=0 $$ La méthode est zéro-stable ssi $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1, \end{cases} $$ ce qui est bien vérifié.

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

Correction schéma C

C est une méthode à $q=3$ pas implicite :

  • $p=2$
  • $b_{-1}=\frac{9}{24}$
  • $a_0=1$ et $a_1=a_2=0$
  • $b_0=\frac{19}{24}$, $b_1=\frac{-5}{24}$ et $b_2=\frac{1}{24}$
  • La méthode est consistante d'ordre $\omega=4$ car $$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=1, \\ \displaystyle\sum_{j=0}^p -ja_j+\sum_{j=-1}^p b_j=-a_1-2a_2+\big(b_{-1}+b_0+b_1+b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{2}a_j+2\sum_{j=-1}^p (-j)^{1}b_j=a_1+4a_2-2\big(-b_{-1}+b_1+2b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{3}a_j+3\sum_{j=-1}^p (-j)^{2}b_j=-a_1-8a_2+3\big(b_{-1}+b_1+4b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{4}a_j+4\sum_{j=-1}^p (-j)^{3}b_j=a_1+16a_2-4\big(-b_{-1}+b_1+8b_2\big)=1 \\ \displaystyle\sum_{j=0}^p (-j)^{5}a_j+5\sum_{j=-1}^p (-j)^{4}b_j=-a_1+32a_2+5\big(b_{-1}+b_1+16b_2\big)\neq1 \end{cases} $$ Vérifions ces calculs ci-dessous:
In [17]:
import sympy as sym
sym.init_printing()
from IPython.display import display, Math
 
p=1
bm1=sym.Rational(9,24)
a0=1
a1=0
a2=0
b0=sym.Rational(19,24)
b1=sym.Rational(-5,24)
b2=sym.Rational(1,24)


aa=[a0,a1,a2]
bb=[b0,b1,b2]

cond=[]
cond.append(sum(aa)) 
cond.append(sum([-j*aa[j]+bb[j] for j in range(len(aa))])+bm1)
if cond==[1,1]:
    display(Math("\omega\ge1"))
    for n in range(2,9):
        cond.append(sum( [ (-j)**n*aa[j]+n*(-j)**(n-1)*bb[j] for j in range(len(aa)) ])+n*bm1)
        if cond[n]==1:
            display(Math("\omega\ge"+str(n)))
        else:
            display(Math("\omega\le"+str(n)))
            break
$\displaystyle \omega\ge1$
$\displaystyle \omega\ge2$
$\displaystyle \omega\ge3$
$\displaystyle \omega\ge4$
$\displaystyle \omega\le5$
  • Le premier polynôme caractéristique est $$ \varrho(r)=r^{p+1}-\sum_{j=0}^p a_jr^{p-j}=r^3-a_0r^2-a_1r-a_2 {}= r^3-r^2=r^2(r-1) $$ dont les racines sont $$ r_0=1,\ r_1=r_2=0 $$ La méthode est zéro-stable ssi $$ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ \varrho'(r_j)\neq0 \text{ si } |r_j|=1, \end{cases} $$ ce qui est bien vérifié.

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

Correction schéma PC

Il est d'ordre $4$ car le prédictor est d'ordre $3$ et le corrector d'ordre $4$.

On définit la solution exacte pour estimer les erreurs:

In [18]:
%reset -f
%matplotlib inline

from matplotlib.pylab import *

# variables globales
t0     = 0
tfinal = 1
y0     = 0

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

sol_exacte = lambda t : t/(1+t**2)

On calcule la solution approchée pour différentes valeurs de $N$ pour les trois schémas:

In [19]:
from scipy.optimize import fsolve


def multipas_P(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = [sol_exacte(tt[0])]
    uu.append(sol_exacte(tt[1]))
    uu.append(sol_exacte(tt[2]))
    for i in range(2,len(tt) - 1):
        uu.append( uu[i]+h/12* ( 23*phi(tt[i],uu[i])-16*phi(tt[i-1],uu[i-1])+5*phi(tt[i-2],uu[i-2]) )  )
    return uu

def multipas_C(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = [sol_exacte(tt[0])]
    uu.append(sol_exacte(tt[1]))
    uu.append(sol_exacte(tt[2]))
    for i in range(2,len(tt) - 1):
        temp = fsolve ( lambda x : -x+uu[i]+h/24* ( 9*phi(tt[i+1],x) + 19*phi(tt[i],uu[i])-5*phi(tt[i-1],uu[i-1])+phi(tt[i-2],uu[i-2]) ), uu[i])
        uu.append(temp)
    return uu

def multipas_PC(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = [sol_exacte(tt[0])]
    uu.append(sol_exacte(tt[1]))
    uu.append(sol_exacte(tt[2]))
    for i in range(2,len(tt) - 1):
        u_tilde = uu[i]+h/12* ( 23*phi(tt[i],uu[i])-16*phi(tt[i-1],uu[i-1])+5*phi(tt[i-2],uu[i-2]) )
        uu.append(  uu[i]+h/24* ( 9*phi(tt[i+1],u_tilde) + 19*phi(tt[i],uu[i])-5*phi(tt[i-1],uu[i-1])+phi(tt[i-2],uu[i-2]) )  )
    return uu


H = []
err_p = []
err_c = []
err_pc = []
N = 10
for k in range(7):
    N+=20
    tt = linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = [sol_exacte(t) for t in tt]
    uu_p = multipas_P(phi, tt, sol_exacte)
    uu_c = multipas_C(phi, tt, sol_exacte)
    uu_pc = multipas_PC(phi, tt, sol_exacte)
    H.append(h)
    err_p.append(max([abs(uu_p[i] - yy[i]) for i in range(len(yy))]))
    err_c.append(max([abs(uu_c[i] - yy[i]) for i in range(len(yy))]))
    err_pc.append(max([abs(uu_pc[i] - yy[i]) for i in range(len(yy))]))
    
print ('Multipas P ordre  = %1.2f' %(polyfit(log(H),log(err_p), 1)[0]))
print ('Multipas C ordre  = %1.2f' %(polyfit(log(H),log(err_c), 1)[0]))
print ('Multipas PC ordre = %1.2f' %(polyfit(log(H),log(err_pc), 1)[0]))

figure(figsize=(8,5))
plot(log(H), log(err_p), 'r-o', label='Multipas P')
plot(log(H), log(err_c), 'b-o', label='Multipas C')
plot(log(H), log(err_pc), 'c-o', label='Multipas PC')
xlabel('$\ln(h)$')
ylabel('$\ln(e)$')
legend(bbox_to_anchor=(1.04, 1), loc='upper left')
grid(True);
Multipas P ordre  = 2.99
Multipas C ordre  = 3.86
Multipas PC ordre = 4.00

Rappels schémas Runge-Kutta

Une méthode de Runge-Kutta à $s\ge1$ étages s'écrit: $$ \begin{array}{c|c} \mathbf{c} & \mathbb{A}\\ \hline &\mathbf{b}^T \end{array} \quad\implies\quad \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}$$
  1. Soit $\omega$ l'ordre de la méthode.

    • Si $\begin{cases} \displaystyle\sum_{j=1}^s b_{i}=1& \\ \displaystyle c_i=\sum_{j=1}^s a_{ij}&i=1,\dots,s \end{cases}$ alors la méthode est consistante (i.e. $\omega\ge1$)

    • Si $\displaystyle\sum_{j=1}^s b_j c_j=\frac{1}{2}$ alors $\omega\ge2$

    • Si $\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}$ alors $\omega\ge3$
    • Si $\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}$ alors $\omega\ge4$
  2. La méthode est zéro-stable

  3. Consistance + zéro-stabilité = convergence

  4. 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}\xrightarrow[t\to+\infty]{}0.$ 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.

    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}$$

Exercice

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 indiquée, il est indiqué le nombre d'étages $s$, son ordre de convergence $\omega$ et si le schéma est explicite, semi-implicite ou implicite.

Exercice: choisir une méthode et vérifier si l'ordre de convergence est correct.

Matrice de Butcher s $\omega$ Exp/Imp/Semi-Implicite
Gauss-1 $\begin{array}{c cc} \frac{1}{2} & \frac{1}{2}\ \hline & 1 \end{array}$ 1 2 I
EE $\begin{array}{c cc} 1 & 1\ \hline & 1 \end{array}$ 1 1 E
EI $\begin{array}{c cc} 0 & 0\ \hline & 1 \end{array}$ 1 1 I
Matrice de Butcher s $\omega$ Exp/Imp/Semi-Implicite
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 I
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 SI
Hammer & 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 ? SI
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 I
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 I
Heun $\begin{array}{c cc} 0 & 0 & 0\ 1 & 1 & 0\ \hline & \frac{1}{2} & \frac{1}{2} \end{array}$ 2 2 E
EM $\begin{array}{c cc} 0 & 0 & 0\ \frac{1}{2} & \frac{1}{2} & 0\ \hline & 0 & 1 \end{array}$ 2 2 E
$\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 E
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 E
$\begin{array}{c cc} 0 & 0 & 0\ \frac{3}{4} & \frac{3}{4} & 0\ \hline & 0 & 1 \end{array}$ 2 1 E
$\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 E
Matrice de Butcher s $\omega$ Exp/Imp/Semi-Implicite
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
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
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
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
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
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
$\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{4}{3} & -\frac{1}{6} \end{array}$ 3 3 E
$\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
$\begin{array}{c ccc} 0 & 0 & 0 & 0 \ 1 & 1 & 0&0\ \frac{1}{2} & \frac{1}{4} &\frac{1}{4}&0\ \hline & \frac{1}{6} & \frac{1}{6} & \frac{2}{3} \end{array}$ 3 3 E
$\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
$\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
$\begin{array}{c ccc} 0 & 0 & 0 & 0 \ \frac{2}{3} & \frac{2}{3} & 0&0\ 0 & -1 &1&0\ \hline & 0 & \frac{3}{4} & \frac{1}{4} \end{array}$ 3 3 E
$\begin{array}{c ccc} 0 & 0 & 0 & 0 \ \frac{1}{3} & \frac{1}{3} & 0&0\ 1 & -1 &2&0\ \hline & 0 & \frac{3}{4} & \frac{1}{4} \end{array}$ 3 3 E
$\begin{array}{c ccc} 0 & 0 & 0 & 0 \ \frac{1}{3} & \frac{1}{3} & 0&0\ \frac{2}{3} & \frac{1}{3} &\frac{1}{3}&0\ \hline & \frac{1}{3} & \frac{1}{3} & \frac{1}{3} \end{array}$ 3 1 E
Matrice de Butcher s $\omega$ Exp/Imp/Semi-Implicite
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 $\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
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
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
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
$\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
$\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
$\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
$\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
$\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
$\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
Matrice de Butcher s $\omega$ Exp/Imp/Semi-Implicite
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
Matrice de Butcher s $\omega$ Exp/Imp/Semi-Implicite
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
Matrice de Butcher s $\omega$ Exp/Imp/Semi-Implicite
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 |