TD2 - Exercices schémas multipas et Runge-Kutta

Auteur·rice

Gloria FACCANONI

Date de publication

15 avril 2024

1 Schémas 👣 multipas 👣 linéaires

Aide-mémoire

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

Pour $i=0,1,\dots$ introduisons la quantité $ \xi(i) = \sum_{j=0}^p (-j)^{i}a_j+i\sum_{j=-1}^p (-j)^{i-1}b_j, $ ayant posé $(-j)^0=1$ même pour $j=0$.

La méthode est consistante ssi $\xi(0)=\xi(1)=1$, autrement dit 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} $
La méthode est zéro-stable ssi $ \begin{cases} |r_j|\le1 \quad\text{pour tout }j=0,\dots,p \\ r_j\text{ a multiplicité}>1\quad\implies \quad |r_j|<1, \end{cases} $ 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}. $
Consistance + zéro-stabilité = convergence
Une méthode convergente est d’ordre $\omega$ ssi $\xi(0)=\xi(1)=\dots=\xi(\omega)=1$, autrement dit ssi elle est consistante et, 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. $ Première barrière de Dahlquist : l’ordre $\omega$ d’une méthode multipas à $q$ étapes convergente satisfait les inégalités suivantes:

si $b_{-1}=0$ (la méthode est explicite) alors $\omega\le q$
si $b_{-1}\neq0$ (la méthode est implicite) on doit considérér deux cas
- si $q$ est impair alors $\omega\le q+1$
- si $q$ est pair alors $\omega\le q+2$

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.

Code

from sympy import symbols, Symbol, latex, factor, solve
from IPython.display import display, Math, Markdown, Latex

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


# Paramètres
p = 1     # j = 0...p # coeffs du schema
q = p + 1
omega = q + 2 if q%2==0 else q + 1 # i = 0, ..., omega

display(Markdown(f"Pour une méthode à q={</span>p <span class="op">+</span> <span class="dv">1</span><span class="sc">} pas, on affiche $\\xi(i)$ pour i = 0, ..., {</span>omega<span class="sc">}"))
for i in range(omega + 1):
    display(Markdown(f"$\\xi({</span>i<span class="sc">}) = {</span>latex(xi(i, q))<span class="sc">}$"))


# systeme = [xi(i, q) for i in range(omega + 2)]
# solve(systeme)

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

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

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

$\xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right)$

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

$\xi(4) = a_{1} - 4 \left(b_{1} - b_{-1}\right)$

1.1 👣 Exercice : consistance d’un schéma à deux pas explicite

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

1.1.1 Correction

C’est une méthode à q=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} $

Voyons sur un exemple ce qui se passe lorsqu’on utilise un schéma non consistant. On compare la solution exacte et la solution approchée du problème de Cauchy $y'(t)=2(y(t)-\cos(t))-\sin(t)$, $y(0)=1$ sur l’intervalle $[0,2\pi]$ dont la solution exacte est $y(t) = cos(t)$. Comme visible sur la figure suivante, la solution numérique ne correspond pas à la solution exacte et on peut vérifier qu’un comportement similaire se produit également si on raffine la discrétisation. Il s’agit d’une situation typique de manque de consistance : l’application d’une méthode non consistante produit l’allure d’une autre fonction plutôt que celle reproduisant la solution exacte.

Code

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

t0, y0 = 0,1 

# Solution exacte
t = sp.symbols('t')
y = sp.Function('y')
edo = sp.Eq(sp.diff(y(t),t) , 2*(y(t)-sp.cos(t))-sp.sin(t) )
# display(edo)
exacte_sp = sp.dsolve( edo , y(t),  ics={y(t0):y0})
# display(exacte_sp)


# Solution approchée
def schema(tt,phi,sol_exacte):
    uu = np.zeros_like(tt)
    q = 2
    uu[:q] = sol_exacte(tt[:q])
    h = tt[1]-tt[0]
    for n in range(q-1,len(tt)-1):
        k0 = phi(tt[n],uu[n])
        k1 = phi(tt[n-1],uu[n-1])
        uu[n+1] = uu[n-1] + h/4*(3*k0-k1)
    return uu


# Comparaison
sol_exacte = sp.lambdify(t, exacte_sp.rhs, 'numpy')
tfin = 2*np.pi
phi =  lambda t,y: 2*(y-np.cos(t))-np.sin(t)
N = 100
tt = np.linspace(t0,tfin,1+N)
yy = sol_exacte(tt)
uu = schema(tt,phi,sol_exacte)
plt.plot(tt,yy,label='exacte')
plt.plot(tt,uu,label='schema')
plt.legend();

1.2 👣 Exercice : convergence d’un schéma à 2 pas implicite

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

1.2.1 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 \xi(0)=\sum_{j=0}^p a_j=-4+5=1, \\ \displaystyle \xi(1)=-\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 $ \begin{aligned} \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) \end{aligned} $ 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.

Voyons sur un exemple ce qui se passe lorsqu’on utilise un schéma consistant mais non zéro-stable. On compare la solution exacte et la solution approchée du problème de Cauchy $y'(t)=2(y(t)-\cos(t))-\sin(t)$, $y(0)=1$ sur l’intervalle $[0,2\pi]$ dont la solution exacte est $y(t) = cos(t)$.

La cosistance est une propriété de précision locale, tandis que la zéro-stabilité et la convergence sont des propriétés de précision globale. La zéro-stabilité assure que la solution numérique ne diverge pas lorsque $h$ tend vers $0$.

Code

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

t0, y0 = 0,1 

# Solution exacte
t = sp.symbols('t')
y = sp.Function('y')
edo = sp.Eq(sp.diff(y(t),t) , 2*(y(t)-sp.cos(t))-sp.sin(t) )
# display(edo)
exacte_sp = sp.dsolve( edo , y(t),  ics={y(t0):y0})
# display(exacte_sp)

# Solution approchée
def schema(tt,phi,sol_exacte):
    uu = np.zeros_like(tt)
    q = 2
    uu[:q] = sol_exacte(tt[:q])
    h = tt[1]-tt[0]
    for n in range(q-1,len(tt)-1):
        eq = lambda x : -x -4*uu[n]+5*uu[n-1] + h*(2*phi(tt[n],uu[n])+phi(tt[n+1],x))
        uu[n+1] =  fsolve(eq,uu[n])[0]
    return uu


# Comparaison
sol_exacte = sp.lambdify(t, exacte_sp.rhs, 'numpy')
tfin = 2*np.pi
phi =  lambda t,y: 2*(y-np.cos(t))-np.sin(t)
N = 20
tt = np.linspace(t0,tfin,1+N)
yy = sol_exacte(tt)
uu = schema(tt,phi,sol_exacte)
plt.plot(tt,yy,label='exacte')
plt.plot(tt,uu,label='schema')
plt.ylim(-3,3)
plt.legend();

1.3 👣 Exercice : convergence d’un schéma à 3 pas explicite

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

1.3.1 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.

Code

import sympy as sym
sym.var('r')
rho = r**3+r**2-r-1
sym.factor(rho)

$\displaystyle \left(r - 1\right) \left(r + 1\right)^{2}$

1.4 👣 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$.

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

Code

import sympy as sym
sym.init_printing()

a_0 = sym.symbols('a_0')
a_1 = sym.symbols('a_1')
b_m1 = sym.symbols('b_{-1}')
b_0 = sym.symbols('b_0')
b_1 = sym.symbols('b_1')

# p=1 => q=2 et omega=4
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) # eq2 = -a_1+b_{-1}+b_0+b_1-1
eq3 = sym.Eq(a_1 + 2 * b_m1 - 2 * b_1, 1) # eq3 = a_1+2*b_{-1}-2*b_1-1
eq4 = sym.Eq(-a_1 + 3 * b_m1 + 3 * b_1, 1) # eq4 = -a_1+3*b_{-1}+3*b_1-1
eq5 = sym.Eq(a_1 + 4 * b_m1 - 4 * b_1, 1) # eq5 = a_1+4*b_{-1}-4*b_1-1
sol = sym.solve([eq1, eq2, eq3, eq4, eq5]) # résolution du système d'équations
display(sol)

# Résultat
u_np1 = sym.symbols('u_{n+1}')
u_n = sym.symbols('u_n')
u_nm1 = sym.symbols('u_{n-1}')
phi_np1 = sym.symbols(r'\varphi_{n+1}')
phi_n = sym.symbols(r'\varphi_n')
phi_nm1 = sym.symbols(r'\varphi_{n-1}')
h = sym.symbols('h')
eq = sym.Eq(u_np1, (a_0 * u_n + a_1 * u_nm1 + h * (b_m1 * phi_np1 + b_0 * phi_n + b_1 * phi_nm1)).subs(sol))
display(eq)

$\displaystyle \left\{ a_{0} : 0, \ a_{1} : 1, \ b_{0} : \frac{4}{3}, \ b_{1} : \frac{1}{3}, \ b_{-1} : \frac{1}{3}\right\}$

$\displaystyle u_{n+1} = h \left(\frac{4 \varphi_{n}}{3} + \frac{\varphi_{n+1}}{3} + \frac{\varphi_{n-1}}{3}\right) + u_{n-1}$

On reconnait la méthode $\text{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} $

1.5 👣 Exercice : convergence d’un schéma à 2 pas implicite

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

Pour quelles valeurs de $\gamma$ est-elle consistante?
Pour quelles valeurs de $\gamma$ est-elle zéro-stable?
Pour quelles valeurs de $\gamma$ est-elle convergente?
Pour quelles valeurs de $\gamma$ est-elle convergente d’ordre $2$? et d’ordre $3$?
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} $

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

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

Zéro-stabilité 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$.

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

Ordre de convergence 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 domaine de zéro-stabilité.
On conclut que la méthode est convergente d’ordre $2$ pour $-1< \gamma\le 1$ et non convergente sinon.

Vérification empirique On définit la solution exacte (calculée en utilisant le module sympy) pour estimer les erreurs, puis on définit le schéma et enfin on calcule les erreurs pour différentes valeurs de $h$.

Code

%reset -f
%matplotlib inline

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

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

Code

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)

# cond_init = sym.Eq( y0, solgen.rhs.subs(t,t0))
# display(cond_init)

# consts = sym.solve( cond_init , dict=True)[0]
# display(consts)

# solpar = solgen.subs(consts).simplify()
# display(solpar)

solpar = sym.dsolve(edo, ics={y(t0):y0})
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)} = \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$:

Code

import numpy as np
from scipy.optimize import fsolve

def multipas(phi, tt, yy, gamma):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    # initialisation
    q = 2
    uu[:q] = yy[:q]
    # recurrence
    for n in range(q-1,len(tt) - 1):
        temp = fsolve(lambda x: -x + (1-gamma)*uu[n]+gamma*uu[n-1] + h/4 * ( (gamma+3)*phi(tt[n+1],x)+(3*gamma+1)*phi(tt[n-1],uu[n-1]) ), uu[n])
        uu[n+1] = temp[0]
    return uu

Code

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

H = []
err_mp = []
N = 10
gamma = 1/2

for k in range(7):
    N += 20
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu_mp = multipas(phi, tt, yy, gamma)
    H.append(h)
    err_mp.append(  np.linalg.norm(uu_mp-yy, np.inf) )
         
         

plt.figure(figsize=(8,5))
plt.plot(np.log(H), np.log(err_mp), 'r-o')
plt.xlabel('$\ln(h)$')
plt.ylabel('$\ln(e)$')
plt.title(f'Ordre de convergence = pente de la droite = {</span>(np.polyfit(np.log(H),np.log(err_mp), <span class="dv">1</span>)[<span class="dv">0</span>])<span class="sc">}')
plt.grid(True)

1.6 👣 Exercice : convergence d’un schéma à 2 pas explicite

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}\gamma^2+\frac{1}{2}\gamma\right)u_{n}+ \left(\frac{2}{3}\gamma^2+\frac{1}{2}\gamma-1\right)u_{n-1}+ h\left[\left(-\frac{1}{6}\gamma^2-\frac{1}{4}\gamma+2\right)\varphi_{n}+ \left(-\frac{1}{6}\gamma^2-\frac{1}{4}\gamma\right)\varphi_{n-1}\right]. $ 1. Pour quelles valeurs de $\gamma$ est-elle consistant e? 1. Est-elle convergente pour ces valeurs? 1. Quel ordre de convergente a-t-elle? 1. Pour les valeurs de $\gamma$ 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} $

1.6.1 Correction

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

$p=1$
$b_{-1}=0$
$a_0=\frac{1}{3}\gamma^2+\frac{1}{2}\gamma$ et $a_1=\frac{2}{3}\gamma^2+\frac{1}{2}\gamma-1$
$b_0=-\frac{1}{6}\gamma^2-\frac{1}{4}\gamma+2$ et $b_1=-\frac{1}{6}\gamma^2-\frac{1}{4}\gamma$

Consistance
La méthode est consistante ssi
$ \begin{cases} \displaystyle\sum_{j=0}^p a_j=\gamma^2+\gamma-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)=\gamma^2+\gamma-1=1 \end{cases} $ donc ssi $\gamma=-2$ ou $\gamma=1$.

Zéro-stabilité
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 }\gamma=-2, \\ r^2-\frac{5}{6}r-\frac{1}{6}&\text{si }\gamma=1, \end{cases} $
dont les racines sont
$ \begin{cases} r_0=1,\ r_1=-\frac{2}{3}&\text{si }\gamma=-2, \\ r_0=1,\ r_1=-\frac{1}{6}&\text{si }\gamma=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} $
Cette condition est vérifiée pour les deux cas.

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

Ordre
La méthode est d’ordre au moins 2 pour les deux choix de $\gamma$ 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 }\gamma=-2, \\ \frac{1}{6}+2\frac{5}{12}=1&\text{si }\gamma=1, \end{cases} $
La méthode n’est pas d’ordre 3 pour aucun des deux choix de $\gamma$ 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 }\gamma=-2, \\ {}-\frac{1}{6}+3\frac{-5}{12}\neq1&\text{si }\gamma=1, \end{cases} $

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

Code

%reset -f
%matplotlib inline


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

from scipy.optimize import fsolve

import sympy as sym
sym.init_printing()

from IPython.display import display, Math

Code

t0     = 0
tfinal = 1
y0     = 2

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

Code

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

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{- t - \frac{1}{5}}{- t - \frac{1}{10}}$

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

Code

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


H = []
err_mp1 = []
err_mp2 = []
N = 10
for k in range(7):
    N+=20
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu_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[n] - yy[n]) for n in range(len(yy))]))
    err_mp2.append(max([abs(uu_mp2[n] - yy[n]) for n in range(len(yy))]))

ordre_gamma_1 = np.polyfit(np.log(H),np.log(err_mp1), 1)[0]
ordre_gamma_m2 = np.polyfit(np.log(H),np.log(err_mp2), 1)[0]

plt.figure(figsize=(8,5))
plt.plot(np.log(H), np.log(err_mp1), 'r-o', label=f'Multipas b=1 : ordre={</span>ordre_gamma_1<span class="sc">:.2f}')
plt.plot(np.log(H), np.log(err_mp2), 'b-d', label=f'Multipas b=-2: ordre={</span>ordre_gamma_m2<span class="sc">:.2f}')
plt.xlabel('$\ln(h)$')
plt.ylabel('$\ln(e)$')
plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
plt.grid(True)

1.7 👣 Exercice : convergence d’un schéma à 2 pas explicite

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

Pour quelles valeurs de $\gamma$ est-elle consistante?
Est-elle convergente pour ces valeurs?
Quel ordre de convergente a-t-elle?
Pour les valeurs de $\gamma$ 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} $

1.7.1 Correction

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

$p=1$
$b_{-1}=0$
$a_0=\gamma^2-\frac{9}{2}\gamma+\frac{11}{2}$ et $a_1=\gamma^2-\frac{7}{2}\gamma+\frac{3}{2}$
$b_0=\frac{\gamma}{2}(\gamma-2)$ et $b_1=\frac{\gamma}{2}\left(\gamma-\frac{10}{3}\right)$
La méthode est consistante ssi $\xi(0)=\xi(1)=1$: $ \begin{cases} \xi(0)=\displaystyle\sum_{j=0}^p a_j=2\gamma^2-8\gamma+7=1, \\ \xi(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}\gamma-\frac{19}{6}=1 \end{cases} $ donc ssi $\gamma=3$. Vérifions ces calculs ci-dessous:

Code

%reset -f
%matplotlib inline


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

from scipy.optimize import fsolve

import sympy as sym
sym.init_printing()

from IPython.display import display, Math

Code

gamma = sym.Symbol('\gamma')
p = 1
bm1 = 0
a0 = gamma**2-sym.Rational(9,2)*gamma+sym.Rational(11,2)
a1 = gamma**2-sym.Rational(7,2)*gamma+sym.Rational(3,2)
b0 = sym.Rational(1,2)*gamma*(gamma-2)
b1 = sym.Rational(1,2)*gamma*(gamma-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)],gamma))

$\displaystyle 2 \gamma^{2} - 8 \gamma + 7$

$\displaystyle \frac{5 \gamma - 9}{6}$

$\displaystyle \left[ \left( 3,\right)\right]$

Dans la suite on considerera $\gamma=3$.

Code

display(Math('a_0='+sym.latex(a0.subs(gamma,3))))
display(Math('a_1='+sym.latex(a1.subs(gamma,3))))
display(Math('b_0='+sym.latex(b0.subs(gamma,3))))
display(Math('b_1='+sym.latex(b1.subs(gamma,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{\gamma=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 \\ \text{multiplicité de $r_j$ supérieur à 1} \implies |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 $\gamma=3$.
La méthode est d’ordre au moins 2 car $\xi(2)=1$: $ \xi(2)= \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{\gamma=3}{=} 1 $
La méthode n’est pas d’ordre 3 car $\xi(3)\neq1$: $ \xi(3)= \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{\gamma=3}{\neq}1 $
On définit la solution exacte (calculée en utilisant le module sympy) pour estimer les erreurs:

Code

t0     = 0
tfinal = 1
y0     = 2

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

Code

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

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}{1 - \frac{e^{- t}}{2}}$

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

Code

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


H = []
err_mp = []
N = 10
for k in range(7):
    N+=20
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu_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))]))

plt.figure(figsize=(8,5))
plt.plot(np.log(H), np.log(err_mp), 'r-o', label='Multipas $\gamma=3$')
plt.xlabel('$\ln(h)$')
plt.ylabel('$\ln(e)$')
plt.legend(bbox_to_anchor=(1.04, 1), loc='upper left')
plt.title(f'Ordre = {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err_mp), <span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:g}')
plt.grid(True)

1.8 Exercice : convergence d’un schéma predictor-corrector basé sur deux schémas 👣

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Soit les méthodes multipas $ \begin{aligned} u_{n+1}&=u_{n}+\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{aligned} $

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

1.8.1 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:

Code

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.

1.8.2 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:

Code

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.

1.8.3 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:

Code

%reset -f
%matplotlib inline

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

Code

import numpy as np
from scipy.optimize import fsolve

def multipas_P(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:3] = sol_exacte(tt[:3])
    for n in range(2,len(tt) - 1):
        k_0 = phi(tt[n],uu[n])
        k_1 = phi(tt[n-1],uu[n-1])
        k_2 = phi(tt[n-2],uu[n-2])
        uu[n+1] = uu[n]+h/12* ( 23*k_0-16*k_1+5*k_2 )
    return uu

def multipas_C(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:3] = sol_exacte(tt[:3])
    for n in range(2,len(tt) - 1):
        k_0 = phi(tt[n],uu[n])
        k_1 = phi(tt[n-1],uu[n-1])
        k_2 = phi(tt[n-2],uu[n-2])
        temp = fsolve ( lambda x : -x+uu[n]+h/24* ( 9*phi(tt[n+1],x) + 19*k_0-5*k_1+k_2 ), uu[n])
        uu[n+1] = temp[0]
    return uu

def multipas_PC(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:3] = sol_exacte(tt[:3])
    for n in range(2,len(tt) - 1):
        k_0 = phi(tt[n],uu[n])
        k_1 = phi(tt[n-1],uu[n-1])
        k_2 = phi(tt[n-2],uu[n-2])
        u_tilde = uu[n]+h/12 * (23*k_0-16*k_1+5*k_2)
        uu[n+1] = uu[n]+h/24* ( 9*phi(tt[n+1],u_tilde) + 19*k_0-5*k_1+k_2 )
    return uu

Code

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


H = []
err_p = []
err_c = []
err_pc = []
N = 10
for k in range(7):
    N += 20
    tt = np.linspace(t0, tfinal, N + 1)
    h = tt[1] - tt[0]
    yy = sol_exacte(tt)
    uu_p = multipas_P(phi, tt, sol_exacte)
    uu_c = multipas_C(phi, tt, sol_exacte)
    uu_pc = multipas_PC(phi, tt, sol_exacte)
    H.append(h)
    err_p.append(np.linalg.norm(uu_p-yy, np.inf))
    err_c.append(np.linalg.norm(uu_c-yy, np.inf))
    err_pc.append(np.linalg.norm(uu_pc-yy, np.inf))
    
print (f'Multipas P ordre  = {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err_p) , <span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:1.2f}' )
print (f'Multipas C ordre  = {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err_c) , <span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:1.2f}' )
print (f'Multipas PC ordre = {</span>np<span class="sc">.</span>polyfit(np.log(H),np.log(err_pc), <span class="dv">1</span>)[<span class="dv">0</span>]<span class="sc">:1.2f}')

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

Multipas P ordre  = 2.99
Multipas C ordre  = 3.86
Multipas PC ordre = 4.00

1.9 Exercice : convergence d’un schéma predictor-corrector basé sur deux schémas 👣

On notera $\varphi_k\stackrel{\text{déf}}{=}\varphi(t_k,u_k)$. Soit les méthodes multipas $ \begin{aligned} u_{n+1}&=\frac{4}{3}u_{n}-\frac{1}{3}u_{n-1}+h\left(b_0\varphi_{n}+b_1\varphi_{n-1}\right) &\text{[P]} \\ u_{n+1}&=\frac{4}{3}u_{n}-\frac{1}{3}u_{n-1}+h\frac{2}{3}\varphi_{n+1} &\text{[C]} \end{aligned} $

Étudier consistance, ordre et zéro-stabilité des deux méthodes.
Écrire et étudier consistance, ordre et zéro-stabilité de la méthode prédictor-corrector associée.
Calculer empiriquement l’ordre de convergence des trois méthodes 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}$.

Les deux schémas multipas linéaires à $q=2$ pas sont de la forme $ u_{n+1} = a_0 u_n + a_1 u_{n-1} + h \left( b_0 \varphi_n + b_1 \varphi_{n-1} \right) + h b_{-1} \varphi_{n+1} $ avec $ \begin{array}{c|cccccc} & a_0 & a_1 & b_0 & b_1 & b_{-1} \\ \hline \text{[P]} & \frac{4}{3} & -\frac{1}{3} & b_0 & b_1 & 0 \\ \text{[C]} & \frac{4}{3} & -\frac{1}{3} & 0 & 0 & \frac{2}{3} \end{array} $

Les deux schémas ont le même premier polynôme caractéristique : $ \varrho(r)=r^2-a_0r-a_1=(r-1)\left(r-\frac{1}{3}\right). $ Ses racines sont $r_0=1$ et $r_1=\frac{1}{3}$. Comme elles appartiennent au disque unité et sont simples, les deux méthodes sont zéro-stables.

D’après la première barrière de Dalquist, un schéma linéaire à $q=2$ pas consistant et zéro-stable est au mieux

d’ordre $\omega = q+2 = 4$ s’il est implicite,
d’ordre $\omega = q = 2$ s’il est explicite.

Étudions l’ordre de convergence en calculant $\xi(i)$ pour $i=0,\dots,4$ (les deux méthodes ayant le même nombre de pas, on a la même expression pour $\xi(i)$) : $ \begin{cases} \xi(0) = a_{0} + a_{1} \\ \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1} \\ \xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right) \\ \xi(3) = - a_{1} + 3 \left(b_{1} + b_{-1}\right) \\ \xi(4) = a_{1} - 4 \left(b_{1} - b_{-1}\right) \end{cases} $

Pour le schéma [C], on a $\xi(0)=\xi(1)=\xi(2)=1$ et $\xi(3)\neq1$, donc la méthode est d’ordre 2 : $ \begin{cases} \xi(0) = a_{0} + a_{1} = \frac{4}{3} -\frac{1}{3}=1\\ \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1} = \frac{1}{3}+\frac{2}{3}=1\\ \xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right) = -\frac{1}{3}+2\frac{2}{3}=1\\ \xi(3) = - a_{1} + 3 \left(b_{1} + b_{-1}\right) = \frac{1}{3}+3\frac{2}{3}=\frac{7}{3}\neq1 \end{cases} $

Pour le schéma [P], nous allons chercher $b_0$ et $b_1$ pour obtenir un schéma d’ordre 2. Pour cela on doit résoudre le système linéaire : $ \begin{cases} \xi(0) = a_{0} + a_{1} = \frac{4}{3} -\frac{1}{3} = 1 \\ \xi(1) = - a_{1} + b_{0} + b_{1} + b_{-1} = \frac{1}{3}+ b_{0} + b_{1}=1\\ \xi(2) = a_{1} - 2 \left(b_{1} - b_{-1}\right) = -\frac{1}{3}-2b_1=1 \end{cases} $ On trouve $b_0=\frac{4}{3}$ et $b_1=-\frac{2}{3}$.

Le schéma predictor-corrector s’écrit $ \begin{cases} u_0 = y_0,\\ u_1 \approx y_1 ,\\ \tilde u = \frac{4}{3}u_{n}-\frac{1}{3}u_{n-1}+h\left(\frac{4}{3}\varphi_{n}-\frac{2}{3}\varphi_{n-1}\right) &\text{[P]} \\ u_{n+1}=\frac{4}{3}u_{n}-\frac{1}{3}u_{n-1}+h\frac{2}{3}\varphi(t_{n+1},\tilde u) \end{cases} $ Comme le prédictor est d’ordre 2 et le corrector d’ordre 2 aussi, le schéma prédictor-corrector est d’ordre 2 (car $\min(\omega_{[C]},\omega_{[P]}+1)=2$).

Code

%reset -f

import numpy as np
from scipy.optimize import fsolve


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


# ==================================================
# Multipas
# ==================================================

q = 2
aa = [4/3,-1/3] # P et C
bb = [4/3,-2/3] # P
bm1 = 2/3 # C

def multipas_P(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:q] = sol_exacte(tt[:q])
    for n in range(q-1,len(tt) - 1):
        k_0 = phi(tt[n],uu[n])
        k_1 = phi(tt[n-1],uu[n-1])
        uu[n+1] = aa[0]*uu[n] + aa[1]*uu[n-1] + h*(bb[0]*k_0+bb[1]*k_1)
    return uu

def multipas_C(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:q] = sol_exacte(tt[:q])
    for n in range(q-1,len(tt) - 1):
        k_0 = phi(tt[n],uu[n])
        k_1 = phi(tt[n-1],uu[n-1])
        temp = fsolve ( lambda x : -x + aa[0]*uu[n] + aa[1]*uu[n-1] + h*bm1*phi(tt[n+1],x) , uu[n])
        uu[n+1] = temp[0]
    return uu

def multipas_PC(phi, tt, sol_exacte):
    h  = tt[1] - tt[0]
    uu = np.zeros_like(tt)
    uu[:q] = sol_exacte(tt[:q])
    for n in range(q-1,len(tt) - 1):
        k_0 = phi(tt[n],uu[n])
        k_1 = phi(tt[n-1],uu[n-1])
        u_tilde = aa[0]*uu[n] + aa[1]*uu[n-1] + h*(bb[0]*k_0+bb[1]*k_1)
        uu[n+1] = aa[0]*uu[n] + aa[1]*uu[n-1] + h*bm1*phi(tt[n+1],u_tilde)
    return uu



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


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

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

2 Schémas 🎨 Runge-Kutta 🎨

Aide-mémoire

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

Soit $\omega$ l’ordre de la méthode.
- Barrières de Butcher Soit une méthode RK à $s$ étapes d’ordre $\omega$.
  - Si la méthode est implicite alors $\color{red}\omega\le 2s$.
  - Si la méthode est explicite et $s<5$ alors $\color{red}\omega\le s$.
  - Si la méthode est explicite et $s\ge5$ alors $\color{red}\omega<s$.
- Critères
  - 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, de plus, $\displaystyle\sum_{j=1}^s b_j c_j=\frac{1}{2}$ alors $\omega\ge2$
  - Si, de plus, $\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, de plus, $\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$
La méthode est zéro-stable
Consistance + zéro-stabilité = convergence
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} $

Code

import sympy 
sympy.init_printing()

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


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

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

    display(Markdown(f"**On a {</span>s<span class="op">+</span><span class="dv">1</span><span class="sc">} conditions pour avoir consistance (= pour être d'ordre 1)**"))
    Eq_1 = ordre_1(s, A, b, c, j)
    Eqs.append(Eq_1)
    display(*Eq_1)
    
    display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
    Eq_2 = ordre_2(s, A, b, c, j)
    Eqs.append([Eq_2])
    display(Eq_2)

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

    return Eqs


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


def ordre_1(s, A, b, c, j):
    texte = "\sum_{j=1}^s b_j =" + f"{</span><span class="bu">sum</span>(b)<span class="sc">.</span>simplify()<span class="sc">}"
    texte += r"\text{ doit être égale à }1"
    # display(Math(texte))
    Eq_1 = [sympy.Eq( sum(b).simplify(), 1 )]
    for i in range(s):
        somma = sympy.summation(A[i,j],(j,0,s-1)).simplify()
        texte = r'\sum_{j=1}^s a_{'</span><span class="op">+</span><span class="bu">str</span>(i)<span class="op">+</span><span class="vs">r'j}=' + sympy.latex( somma )
        texte += r"\text{ doit être égale à }"+sympy.latex(c[i])
        # display(Math( texte ))
        Eq_1.append(sympy.Eq( somma, c[i] ))
    return Eq_1    


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


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


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

    
    
# Exemple de conditions avec s=2 étages
s = 2
type = "semi-implicite" 
display(Markdown(f"### Conditions pour un schéma {</span><span class="bu">type</span><span class="sc">} avec s = {</span>s<span class="sc">} étages "))
Eqs = ordre_RK(s, type)

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

Matrice de Butcher

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

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

$\displaystyle b_{0} + b_{1} = 1$

$\displaystyle a_{00} = c_{0}$

$\displaystyle a_{10} + a_{11} = c_{1}$

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

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

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

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

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

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

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

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

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

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

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

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

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

Code

%reset -f

import sympy 
sympy.init_printing()

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

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

Correction point 1 : schéma

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

Il s’agit d’un schéma à $s=2$ étages explicite.

Correction point 2 : ordre de convergence

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

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

Code

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

    j = sympy.symbols('j')
    
    if A is None:
        A = sympy.MatrixSymbol('a', s, s)
    else:
        A = sympy.Matrix(A)
    if c is None:
        c = sympy.symbols(f'c_0:{</span>s<span class="sc">}')    
    if b is None:
        b = sympy.symbols(f'b_0:{</span>s<span class="sc">}')
    
    
    display(Markdown("**Matrice de Butcher**"))
    But = matrice_Butcher(s, A, b, c)
    
    Eqs = []

    display(Markdown(f"**On a {</span>s<span class="op">+</span><span class="dv">1</span><span class="sc">} conditions pour avoir consistance = pour être d'ordre 1**"))
    Eq_1 = ordre_1(s, A, b, c, j)
    Eqs.append(Eq_1)
    
    display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
    Eq_2 = ordre_2(s, A, b, c, j)
    Eqs.append([Eq_2])

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

    return Eqs


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


def ordre_1(s, A, b, c, j):
    texte = "\sum_{j=1}^s b_j =" + f"{</span><span class="bu">sum</span>(b)<span class="sc">.</span>simplify()<span class="sc">}"
    texte += r"\text{ doit être égale à }1"
    display(Math(texte))
    Eq_1 = [sympy.Eq( sum(b).simplify(), 1 )]
    for i in range(s):
        somma = sympy.summation(A[i,j],(j,0,s-1)).simplify()
        texte = r'\sum_{j=1}^s a_{'</span><span class="op">+</span><span class="bu">str</span>(i)<span class="op">+</span><span class="vs">r'j}=' + sympy.latex( somma )
        texte += r"\text{ doit être égale à }"+sympy.latex(c[i])
        display(Math( texte ))
        Eq_1.append(sympy.Eq( somma, c[i] ))
    return Eq_1    


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


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


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



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

s = 2

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

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

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

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

Matrice de Butcher

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

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

$\displaystyle \sum_{j=1}^s b_j =1\text{ doit être égale à }1$

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

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

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

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

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

Correction point 3 : A-stabilité

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

On obtient alors la relation de récurrence suivante :

Code

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

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

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

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

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

$\displaystyle k_{0} = - \beta u_{n}$

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

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

Code

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

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

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

Code

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

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

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

Code

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

sympy.solve(abs(q)<1)

$\displaystyle x < 2$

Correction point 4 : test de l’ordre

Code

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

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

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

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

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

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

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


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


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

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

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

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

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

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

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

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

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

Code

%reset -f

import sympy 
sympy.init_printing()

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

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

Correction point 1 : schéma

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

Il s’agit d’un schéma à $s=2$ étages implicite.

Correction point 2 : ordre de convergence

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

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

Code

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

    j = sympy.symbols('j')
    
    if A is None:
        A = sympy.MatrixSymbol('a', s, s)
    else:
        A = sympy.Matrix(A)
    if c is None:
        c = sympy.symbols(f'c_0:{</span>s<span class="sc">}')    
    if b is None:
        b = sympy.symbols(f'b_0:{</span>s<span class="sc">}')
    
    
    display(Markdown("**Matrice de Butcher**"))
    But = matrice_Butcher(s, A, b, c)
    
    Eqs = []

    display(Markdown(f"**On a {</span>s<span class="op">+</span><span class="dv">1</span><span class="sc">} conditions pour avoir consistance = pour être d'ordre 1**"))
    Eq_1 = ordre_1(s, A, b, c, j)
    Eqs.append(Eq_1)
    
    display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
    Eq_2 = ordre_2(s, A, b, c, j)
    Eqs.append([Eq_2])

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

    return Eqs


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


def ordre_1(s, A, b, c, j):
    texte = "\sum_{j=1}^s b_j =" + f"{</span><span class="bu">sum</span>(b)<span class="sc">.</span>simplify()<span class="sc">}"
    texte += r"\text{ doit être égale à }1"
    display(Math(texte))
    Eq_1 = [sympy.Eq( sum(b).simplify(), 1 )]
    for i in range(s):
        somma = sympy.summation(A[i,j],(j,0,s-1)).simplify()
        texte = r'\sum_{j=1}^s a_{'</span><span class="op">+</span><span class="bu">str</span>(i)<span class="op">+</span><span class="vs">r'j}=' + sympy.latex( somma )
        texte += r"\text{ doit être égale à }"+sympy.latex(c[i])
        display(Math( texte ))
        Eq_1.append(sympy.Eq( somma, c[i] ))
    return Eq_1    


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


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


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



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

s = 2

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

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

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

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

Matrice de Butcher

$\displaystyle \left[\begin{matrix}0 & \frac{1}{4} & - \frac{1}{4}\\\frac{2}{3} & \frac{1}{4} & \frac{5}{12}\\ & \frac{1}{4} & \frac{3}{4}\end{matrix}\right]$

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

$\displaystyle \sum_{j=1}^s b_j =1\text{ doit être égale à }1$

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

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

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

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

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

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

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

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

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

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

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

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

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

Correction point 3 : A-stabilité

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

On obtient alors la relation de récurrence suivante :

Code

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

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

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

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

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

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

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

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

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

Code

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

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

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

Code

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

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

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

Code

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

sympy.solve(abs(q)<1)

$\displaystyle \text{True}$

Correction point 4 : test de l’ordre

Code

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

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

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

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

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

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

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


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


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

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

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

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

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

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

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

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

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

Code

%reset -f

import sympy 
sympy.init_printing()

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

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

Correction point 1 : schéma

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

Il s’agit d’un schéma à $s=2$ étages implicite.

Correction point 2 : ordre de convergence

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

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

Code

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

    j = sympy.symbols('j')
    
    if A is None:
        A = sympy.MatrixSymbol('a', s, s)
    else:
        A = sympy.Matrix(A)
    if c is None:
        c = sympy.symbols(f'c_0:{</span>s<span class="sc">}')    
    if b is None:
        b = sympy.symbols(f'b_0:{</span>s<span class="sc">}')
    
    
    display(Markdown("**Matrice de Butcher**"))
    But = matrice_Butcher(s, A, b, c)
    
    Eqs = []

    display(Markdown(f"**On a {</span>s<span class="op">+</span><span class="dv">1</span><span class="sc">} conditions pour avoir consistance = pour être d'ordre 1**"))
    Eq_1 = ordre_1(s, A, b, c, j)
    Eqs.append(Eq_1)
    
    display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
    Eq_2 = ordre_2(s, A, b, c, j)
    Eqs.append([Eq_2])

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

    return Eqs


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


def ordre_1(s, A, b, c, j):
    texte = "\sum_{j=1}^s b_j =" + f"{</span><span class="bu">sum</span>(b)<span class="sc">.</span>simplify()<span class="sc">}"
    texte += r"\text{ doit être égale à }1"
    display(Math(texte))
    Eq_1 = [sympy.Eq( sum(b).simplify(), 1 )]
    for i in range(s):
        somma = sympy.summation(A[i,j],(j,0,s-1)).simplify()
        texte = r'\sum_{j=1}^s a_{'</span><span class="op">+</span><span class="bu">str</span>(i)<span class="op">+</span><span class="vs">r'j}=' + sympy.latex( somma )
        texte += r"\text{ doit être égale à }"+sympy.latex(c[i])
        display(Math( texte ))
        Eq_1.append(sympy.Eq( somma, c[i] ))
    return Eq_1    


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


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


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



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

s = 2

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

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

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

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

Matrice de Butcher

$\displaystyle \left[\begin{matrix}\frac{1}{3} & \frac{5}{12} & - \frac{1}{12}\\1 & \frac{3}{4} & \frac{1}{4}\\ & \frac{3}{4} & \frac{1}{4}\end{matrix}\right]$

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

$\displaystyle \sum_{j=1}^s b_j =1\text{ doit être égale à }1$

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

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

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

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

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

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

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

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

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

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

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

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

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

Correction point 3 : A-stabilité

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

On obtient alors la relation de récurrence suivante :

Code

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

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

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

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

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

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

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

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

Code

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

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

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

Code

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

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

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

Code

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

sympy.solve(abs(q)<1)

$\displaystyle \text{True}$

Correction point 4 : test de l’ordre

Code

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

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

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

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

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

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

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


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


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

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

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

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

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

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

2.5 🎨 Exercice ($s=3$, explicite)

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

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

Code

%reset -f

import sympy 
sympy.init_printing()

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

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

Correction point 1 : schéma

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

Il s’agit d’un schéma à $s=3$ étages explicite.

Correction point 2 : ordre de convergence

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

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

Code

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

    j = sympy.symbols('j')
    
    if A is None:
        A = sympy.MatrixSymbol('a', s, s)
    else:
        A = sympy.Matrix(A)
    if c is None:
        c = sympy.symbols(f'c_0:{</span>s<span class="sc">}')    
    if b is None:
        b = sympy.symbols(f'b_0:{</span>s<span class="sc">}')
    
    
    display(Markdown("**Matrice de Butcher**"))
    But = matrice_Butcher(s, A, b, c)
    
    Eqs = []

    display(Markdown(f"**On a {</span>s<span class="op">+</span><span class="dv">1</span><span class="sc">} conditions pour avoir consistance = pour être d'ordre 1**"))
    Eq_1 = ordre_1(s, A, b, c, j)
    Eqs.append(Eq_1)
    
    display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
    Eq_2 = ordre_2(s, A, b, c, j)
    Eqs.append([Eq_2])

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

    return Eqs


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


def ordre_1(s, A, b, c, j):
    texte = "\sum_{j=1}^s b_j =" + f"{</span><span class="bu">sum</span>(b)<span class="sc">.</span>simplify()<span class="sc">}"
    texte += r"\text{ doit être égale à }1"
    display(Math(texte))
    Eq_1 = [sympy.Eq( sum(b).simplify(), 1 )]
    for i in range(s):
        somma = sympy.summation(A[i,j],(j,0,s-1)).simplify()
        texte = r'\sum_{j=1}^s a_{'</span><span class="op">+</span><span class="bu">str</span>(i)<span class="op">+</span><span class="vs">r'j}=' + sympy.latex( somma )
        texte += r"\text{ doit être égale à }"+sympy.latex(c[i])
        display(Math( texte ))
        Eq_1.append(sympy.Eq( somma, c[i] ))
    return Eq_1    


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


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


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



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

s = 3

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

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

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

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

Matrice de Butcher

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

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

$\displaystyle \sum_{j=1}^s b_j =1\text{ doit être égale à }1$

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

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

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

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

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

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

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

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

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

Correction point 3 : A-stabilité

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

On obtient alors la relation de récurrence suivante :

Code

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

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

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

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

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

$\displaystyle k_{0} = - \beta u_{n}$

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

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

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

Code

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

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

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

Code

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

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

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

Code

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

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

$\beta h<2.51274532661833$

Correction point 4 : test de l’ordre

Code

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

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

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

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

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

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

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


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


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

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

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

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

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

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

2.6 🎨 Exercice ($s=3$, implicite)

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

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

Code

%reset -f

import sympy 
sympy.init_printing()


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


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

Correction point 1 : schéma

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

Il s’agit d’un schéma à $s=3$ étages implicite.

Correction point 2 : ordre de convergence

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

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

Code

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

    j = sympy.symbols('j')
    
    if A is None:
        A = sympy.MatrixSymbol('a', s, s)
    else:
        A = sympy.Matrix(A)
    if c is None:
        c = sympy.symbols(f'c_0:{</span>s<span class="sc">}')    
    if b is None:
        b = sympy.symbols(f'b_0:{</span>s<span class="sc">}')
    
    
    display(Markdown("**Matrice de Butcher**"))
    But = matrice_Butcher(s, A, b, c)
    
    Eqs = []

    display(Markdown(f"**On a {</span>s<span class="op">+</span><span class="dv">1</span><span class="sc">} conditions pour avoir consistance = pour être d'ordre 1**"))
    Eq_1 = ordre_1(s, A, b, c, j)
    Eqs.append(Eq_1)
    
    display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
    Eq_2 = ordre_2(s, A, b, c, j)
    Eqs.append([Eq_2])

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

    return Eqs


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


def ordre_1(s, A, b, c, j):
    texte = "\sum_{j=1}^s b_j =" + f"{</span><span class="bu">sum</span>(b)<span class="sc">.</span>simplify()<span class="sc">}"
    texte += r"\text{ doit être égale à }1"
    display(Math(texte))
    Eq_1 = [sympy.Eq( sum(b).simplify(), 1 )]
    for i in range(s):
        somma = sympy.summation(A[i,j],(j,0,s-1)).simplify()
        texte = r'\sum_{j=1}^s a_{'</span><span class="op">+</span><span class="bu">str</span>(i)<span class="op">+</span><span class="vs">r'j}=' + sympy.latex( somma )
        texte += r"\text{ doit être égale à }"+sympy.latex(c[i])
        display(Math( texte ))
        Eq_1.append(sympy.Eq( somma, c[i] ))
    return Eq_1    


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


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


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



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

s = 3

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

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

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

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

Matrice de Butcher

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

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

$\displaystyle \sum_{j=1}^s b_j =1\text{ doit être égale à }1$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Correction point 3 : A-stabilité

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

On obtient alors la relation de récurrence suivante :

Code

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

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

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

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

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

$\displaystyle k_{0} = - \beta u_{n}$

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

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

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

Code

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

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

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

Code

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

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

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

Code

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

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

$\beta h<6$

Correction point 4 : test de l’ordre

Code

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

phi = lambda t,y: -y

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

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

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

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

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


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


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

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

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

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

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

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

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

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

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

Code

%reset -f

import sympy 
sympy.init_printing()

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

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

Correction point 1 : schéma

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

Il s’agit d’un schéma à $s=3$ étages explicite.

Correction point 2 : ordre de convergence

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

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

Code

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

    j = sympy.symbols('j')
    
    if A is None:
        A = sympy.MatrixSymbol('a', s, s)
    else:
        A = sympy.Matrix(A)
    if c is None:
        c = sympy.symbols(f'c_0:{</span>s<span class="sc">}')    
    if b is None:
        b = sympy.symbols(f'b_0:{</span>s<span class="sc">}')
    
    
    display(Markdown("**Matrice de Butcher**"))
    But = matrice_Butcher(s, A, b, c)
    
    Eqs = []

    display(Markdown(f"**On a {</span>s<span class="op">+</span><span class="dv">1</span><span class="sc">} conditions pour avoir consistance = pour être d'ordre 1**"))
    Eq_1 = ordre_1(s, A, b, c, j)
    Eqs.append(Eq_1)
    
    display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
    Eq_2 = ordre_2(s, A, b, c, j)
    Eqs.append([Eq_2])

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

    return Eqs


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


def ordre_1(s, A, b, c, j):
    texte = "\sum_{j=1}^s b_j =" + f"{</span><span class="bu">sum</span>(b)<span class="sc">.</span>simplify()<span class="sc">}"
    texte += r"\text{ doit être égale à }1"
    display(Math(texte))
    Eq_1 = [sympy.Eq( sum(b).simplify(), 1 )]
    for i in range(s):
        somma = sympy.summation(A[i,j],(j,0,s-1)).simplify()
        texte = r'\sum_{j=1}^s a_{'</span><span class="op">+</span><span class="bu">str</span>(i)<span class="op">+</span><span class="vs">r'j}=' + sympy.latex( somma )
        texte += r"\text{ doit être égale à }"+sympy.latex(c[i])
        display(Math( texte ))
        Eq_1.append(sympy.Eq( somma, c[i] ))
    return Eq_1    


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


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


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



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

s = 3

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

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

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

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

Matrice de Butcher

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

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

$\displaystyle \sum_{j=1}^s b_j =b_0 + b_1\text{ doit être égale à }1$

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

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

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

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

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

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

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

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

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

Code

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

Matrice de Butcher

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

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

$\displaystyle \sum_{j=1}^s b_j =1\text{ doit être égale à }1$

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

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

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

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

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

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

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

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

Correction point 3 : A-stabilité

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

On obtient alors la relation de récurrence suivante :

Code

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

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

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

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

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

$\displaystyle k_{0} = - \beta u_{n}$

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

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

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

Code

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

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

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

Code

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

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

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

Code

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

sympy.solve(abs(q)<1)

$\displaystyle x < 2$

Correction point 4 : test de l’ordre

Code

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

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

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

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

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

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

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


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


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

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

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

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

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

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

2.8 Catalogue de schémas 🎨 Runge-Kutta

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 : étudier l’ordre théorique de convergence et le vérifier empiriquement.

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

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

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

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

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

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

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

Retour au sommet

TD2 - Exercices schémas multipas et Runge-Kutta

1 Schémas 👣 multipas 👣 linéaires

1.1 👣 Exercice : consistance d’un schéma à deux pas explicite

1.1.1 Correction

1.2 👣 Exercice : convergence d’un schéma à 2 pas implicite

1.2.1 Correction

1.3 👣 Exercice : convergence d’un schéma à 3 pas explicite

1.3.1 Correction

1.4 👣 Exercice : schéma à deux pas d’ordre 4

1.4.1 Correction

1.5 👣 Exercice : convergence d’un schéma à 2 pas implicite

1.5.1 Correction

1.6 👣 Exercice : convergence d’un schéma à 2 pas explicite

1.6.1 Correction

1.7 👣 Exercice : convergence d’un schéma à 2 pas explicite

1.7.1 Correction

1.8 Exercice : convergence d’un schéma predictor-corrector basé sur deux schémas 👣

1.8.1 Correction schéma P

1.8.2 Correction schéma C

1.8.3 Correction schéma PC

1.9 Exercice : convergence d’un schéma predictor-corrector basé sur deux schémas 👣

2 Schémas 🎨 Runge-Kutta 🎨

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

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

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

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

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

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

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

2.8 Catalogue de schémas 🎨 Runge-Kutta