Comme nous l’avons vu au chapitre précédent, les barrières sur l’ordre de convergence des méthodes multi-pas linéaires sont de plus en plus sévères. Nous allons alors étudier une autre famille de méthodes, à savoir les méthodes de Runge-Kutta. La stratégie est nouvelle par rapport à celle des méthodes multi-pas : en effet, on construit des méthodes à un seul pas, mais nous passons à une stratégie multi-étapes s’appuyant sur l’information en quelques points supplémentaires, situés à l’intérieur de chaque sous-intervalle de la discrétisation du domaine.
Schémas
Les schémas de Runge-Kutta sont des méthodes à un pas permettant d’approcher l’intégrale
\(
\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t
\)
par une formule de quadrature en des points situés à l’intérieur de l’intervalle \([t_n,t_{n+1}]\) :
\(
t_{n,i} \stackrel{\text{déf}}{=} t_n + c_i h, \quad \text{avec } c_i \in [0,1],
\) \(
\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\mathrm{d}t
\approx
h\sum_{i=1}^{s}b_i\varphi(t_{n,i},y(t_{n,i})).
\) Cependant, si \(c_i\neq0\) et \(c_i\neq1\), alors \(t_{n,i}\) ne correspond pas à un point de discrétisation, ce qui signifie que \(y(t_{n,i})\) est inconnu.
Pour contourner ce problème, l’évaluation de \(\varphi(t_{n,i},y(t_{n,i}))\) en ces points intermédiaires est approchée à l’aide d’une autre formule de quadrature :
\(
\varphi(t_{n,i},y(t_{n,i}))
\approx
\varphi\left(t_{n,i},y_n+h\sum_{j=1}^{s}a_{ij}\varphi(t_{n,j},y(t_{n,j}))\right).
\) On en déduit alors l’approximation suivante :
\(
y_{n+1}
\approx
y_n
+
h
\sum_{i=1}^{s}b_i\varphi\left(t_{n,i},y_n+h\sum_{j=1}^{s}a_{ij}\varphi(t_{n,j},y_{n,j})\right).
\)
Une méthode de Runge-Kutta à \(s\ge1\) étages contient \(s^2+2s\) paramètres et s’écrit: \(\begin{cases}
u_0=y(t_0)=y_0,\\
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\\
u_{n+1}=u_{n}+h \displaystyle\sum_{i=1}^sb_iK_i& n=0,1,\dots N-1.
\end{cases}\)
Matrice de Butcher
Les coefficients sont généralement organisés en deux vecteurs \(\mathbf{b}=(b_1,\dots,b_s)^T\), \(\mathbf{c}=(c_1,\dots,c_s)^T\) et une matrice \(\mathbb{A}=(a_{ij})_{1\le i,j\le s}\).
Le tableau
\(
\begin{array}{c|c}
\mathbf{c} & \mathbb{A}\\
\hline
&\mathbf{b}^T
\end{array}
\)
est appelé matrice de Butcher de la méthode de Runge-Kutta considérée.
Méthode explicite. Si \(a_{ij}=0\) pour \(j\ge i\) (i.e. la matrice \(\mathbb{A}\) est triangulaire inférieure stricte), alors chaque \(K_i\) peut être calculé explicitement en fonction des \(i-1\) coefficients \(K_1,\dots,K_{i-1}\) déjà connus.
Méthode semi-implicite (ou diagonalement implicite, DIRK). Si \(a_{ij}=0\) pour \(j> i\) (i.e. la matrice \(\mathbb{A}\) est triangulaire inférieure), alors chaque \(K_i\) est solution de l’équation non linéaire
\(
\boxed{K_i}=\varphi\left(t_n+c_ih, u_n+ha_{ii}\boxed{K_i}+h\sum_{j=1}^{i-1}a_{ij}K_j\right)
\)
Un schéma semi-implicite implique donc la résolution de \(s\) équations non linéaires indépendantes.
Si, de plus, tous les termes diagonaux sont égaux (\(a_{ii}=\gamma\) pour \(i=1,\dots,s\)), on parle de singly diagonal implicit method (SDIRK).
Méthode implicite. Dans tous les autres cas, la méthode est implicite, ce qui signifie qu’il faut résoudre un système non linéaire de dimension \(s\) pour obtenir les \(K_i\).
Bien noter que la méthode est implicite non pas parce que \(u_{n+1}\) dépend de lui-même, mais parce qu’au moins un \(K_i\) dépend de lui-même.
Notons aussi que, si \(c_s=1\) et \(a_{sj}=b_j\) pour tout \(j=1,\dots,s\), alors \(K_s=\varphi(t_{n+1},u_{n+1})\).
L’augmentation du coût de calcul pour les schémas implicites peut rendre leur utilisation onéreuse.
Un compromis souvent adopté consiste à utiliser des méthodes semi-implicites.
Convergence et ordre de convergence
L’étude de l’ordre des méthodes de Runge-Kutta utilise la théorie de Butcher. Elle repose sur des outils qui combinent l’analyse numérique, la théorie des graphes et la différenciation des champs vectoriels. Plus tard, des relations avec la théorie des groupes et la théorie quantique des champs ont été découvertes. Une monographie de référence sur la théorie de l’ordre de Butcher et les problèmes connexes est rédigée par John C. Butcher lui-même. L’outil de base qui relie tous les domaines susmentionnés est la notion d’arbres racinés. Nous n’énoncerons ici que les résultats.
Soit \(\omega\) l’ordre de la méthode (la méthode est dite consistante si \(\omega\ge1\)).
Théorème de Lax-Richtmyer (consistance + zéro-stabilité \(\iff\) convergence)
Le schéma est convergent (d’ordre \(\omega\ge1\)) ssi il est zéro-stable et consistante (d’ordre \(\omega\)).
Proposition (Zéro-stabilité)
Les méthodes de Runge-Kutta sont zéro-stables.
Proposition (Consistance)
La méthode de Runge-Kutta est consistante (i.e. \(\omega\ge1\)) ssi
\(
\begin{cases}
\displaystyle\sum_{j=1}^s b_{j}=1&
\\
\displaystyle c_i=\sum_{j=1}^s a_{ij}, & i=1,\dots,s
\end{cases}
\)
Remarque En particulier, pour une méthode explicite, on aura \(c_1=a_{1,j}=0\) ainsi \(K_1=\varphi(t_n,u_n)\) et \(c_2=a_{21}\) ainsi \(K_2=\varphi(t_n+c_2h,u_n+hc_2K_1)\).
Proposition (Ordre de consistance)
Une méthode de Runge-Kutta est consistante d’ordre \(\omega\ge1\) si elle satisfait les conditions suivantes :
| \(\omega \geq 1\) |
\(\begin{cases} \displaystyle 1 = \sum_{j=1}^s b_j \\ \displaystyle c_i = \sum_{j=1}^s a_{ij}, \quad i=1,\dots,s \end{cases}\) |
| \(\omega \geq 2\) |
Les conditions pour \(\omega \geq 1\) et \(\displaystyle\sum_{j=1}^s b_j c_j = \frac{1}{2}\) |
| \(\omega \geq 3\) |
Les conditions pour \(\omega \geq 2\) et \(\begin{cases} \displaystyle\sum_{j=1}^s b_j c_j^2 = \frac{1}{3} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i a_{ij} c_j = \frac{1}{6} \end{cases}\) |
| \(\omega \geq 4\) |
Les conditions pour \(\omega \geq 3\) et \(\begin{cases} \displaystyle\sum_{j=1}^s b_j c_j^3 = \frac{1}{4} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i c_i a_{ij} c_j = \frac{1}{8} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s b_i a_{ij} c_j^2 = \frac{1}{12} \\ \displaystyle\sum_{i=1}^s\sum_{j=1}^s\sum_{k=1}^s b_i a_{ij}a_{jk} c_k = \frac{1}{24} \end{cases}\) |
(Pour les méthodes ERK cf. page 135 de E. Hairer, S. R Norsett, G. Wanner, Solving Ordinary Differential Equations I. Nonstiff Problems, Second Revised Edition, 2008)
Calcul explicite de l’ordre de consistance avec Sympy
La fonction ordre_RK sert à analyser l’ordre de consistance d’une méthode de Runge-Kutta donnée. Elle prend en entrée :
s : le nombre d’étages de la méthode Runge-Kutta (RK),
A, b, c (optionnels) : les coefficients de la matrice de Butcher. Si ces coefficients ne sont pas fournis, ils sont remplacés par des variables symboliques.
La fonction affiche successivement : la matrice de Butcher, qui résume les coefficients de la méthode RK, et les conditions nécessaires pour que la méthode soit d’un certain ordre (de 1 à 4).
- Exemple d’utilisation sans paramètres optionnels
ordre_RK(s=2, type="explicite")
- Exemple d’utilisation avec des paramètres optionnels (matrice de Butcher donnée)
A = sympy.Matrix([[0, 0], [sympy.Rational(1, 2), 0]])
b = sympy.Matrix([sympy.Rational(1, 2), sympy.Rational(1, 2)])
c = sympy.Matrix([0, 1])
s = len(c)
ordre_RK(s, A, b, c)
Code
# =============================================================================
# Fonctions pour calculer les conditions de consistance pour un schéma RK.
#
# La fonction ordre_RK calcule les conditions pour les ordres 1, 2, 3 et 4 (même si théoriquement cela n'est pas possible)
# - affiche la matrice de Butcher
# - affiche les conditions de consistance obtenues pour chaque ordre par la fonction ordre_i
# - renvoie une liste d'équations pour une éventuelle résolution
#
# Chaque fonction ordre_i calcule les conditions pour être d'ordre i
# - peut afficher les conditions de consistance [actuellement en commentaire]
# - renvoie la liste d'équations
#
# =============================================================================
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:{s}')
if b is None:
b = sympy.symbols(f'b_0:{s}')
display(Markdown("**Matrice de Butcher**"))
But = matrice_Butcher(s, A, b, c)
Eqs = []
display(Markdown(f"**On a {s+1} conditions pour avoir consistance (= pour être d'ordre 1)**"))
Eq_1 = ordre_1(s, A, b, c, j)
Eqs.append(Eq_1)
display(*Eq_1)
display(Markdown("**On doit ajouter 1 condition pour être d'ordre 2**"))
Eq_2 = ordre_2(s, A, b, c, j)
Eqs.append([Eq_2])
display(Eq_2)
display(Markdown("**On doit ajouter 2 conditions pour être d'ordre 3**"))
Eq_3 = ordre_3(s, A, b, c, j)
Eqs.append(Eq_3)
display(*Eq_3)
display(Markdown("**On doit ajouter 4 conditions pour être d'ordre 4**"))
Eq_4 = ordre_4(s, A, b, c, j)
Eqs.append(Eq_4)
display(*Eq_4)
return Eqs
def matrice_Butcher(s, A, b, c):
But = sympy.Matrix(A)
But = But.col_insert(0, sympy.Matrix(c))
last = [sympy.Symbol(" ")]
last.extend(b)
But = But.row_insert(s,sympy.Matrix(last).T)
display(Math(sympy.latex(sympy.Matrix(But))))
return But
def ordre_1(s, A, b, c, j):
texte = r"\sum_{j=1}^s b_j =" + f"{sum(b).simplify()}"
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_{'+str(i)+r'j}=' + sympy.latex( somma )
texte += r"\text{ doit être égale à }"+sympy.latex(c[i])
# display(Math( texte ))
Eq_1.append(sympy.Eq( somma, c[i] ))
return Eq_1
def ordre_2(s, A, b, c, j):
texte = r'\sum_{j=1}^s b_j c_j=' +sympy.latex(sum([b[i]*c[i] for i in range(s)]).simplify())
texte += r"\text{ doit être égale à }\frac{1}{2}"
# display(Math(texte))
Eq_2 = sympy.Eq( sum([b[i]*c[i] for i in range(s)]).simplify(), sympy.Rational(1,2) )
return Eq_2
def ordre_3(s, A, b, c, j):
texte = r'\sum_{j=1}^s b_j c_j^2='
texte += sympy.latex( sum([b[i]*c[i]**2 for i in range(s)]).simplify() )
texte += r"\text{ doit être égale à }\frac{1}{3}"
# display(Math(texte))
Eq_3_1 = sympy.Eq( sum([b[i]*c[i]**2 for i in range(s)]).simplify(), sympy.Rational(1,3) )
texte = r'\sum_{i,j=1}^s b_i a_{ij} c_j='
somma = sum([b[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify()
texte = texte + sympy.latex(somma)
texte += r"\text{ doit être égale à }\frac{1}{6}"
# display(Math(texte))
Eq_3_2 = sympy.Eq( sum([b[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify(), sympy.Rational(1,6) )
return [Eq_3_1, Eq_3_2]
def ordre_4(s, A, b, c, j):
texte = r'\sum_{j=1}^s b_j c_j^3='
texte += sympy.latex( sum([b[i]*c[i]**3 for i in range(s)]).simplify() )
texte += r"\text{ doit être égale à }\frac{1}{4}"
# display(Math(texte))
Eq_4_1 = sympy.Eq( sum([b[i]*c[i]**3 for i in range(s)]).simplify(), sympy.Rational(1,4) )
texte = r'\sum_{i,j=1}^s b_i c_i a_{ij} c_j='
somma = sum([b[i]*c[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify()
texte = texte + sympy.latex(somma)
texte += r"\text{ doit être égale à }\frac{1}{8}"
# display(Math(texte))
Eq_4_2 = sympy.Eq( sum([b[i]*c[i]*A[i,j]*c[j] for j in range(s) for i in range(s)]).simplify(), sympy.Rational(1,8) )
texte = r'\sum_{i,j=1}^s b_i a_{ij} c_j^2='
somma = sum([b[i]*A[i,j]*c[j]**2 for j in range(s) for i in range(s)]).simplify()
texte = texte + sympy.latex(somma)
texte += r"\text{ doit être égale à }\frac{1}{12}"
# display(Math(texte))
Eq_4_3 = sympy.Eq( sum([b[i]*A[i,j]*c[j]**2 for j in range(s) for i in range(s)]).simplify(), sympy.Rational(1,12) )
texte = r'\sum_{i,j,k=1}^s b_i a_{ij} a_{jk} c_k='
somma = sum([b[i]*A[i,j]*A[j,k]*c[k] for k in range(s) for j in range(s) for i in range(s)]).simplify()
texte = texte + sympy.latex(somma)
texte += r"\text{ doit être égale à }\frac{1}{24}"
# display(Math(texte))
Eq_4_4 = sympy.Eq( sum([b[i]*A[i,j]*A[j,k]*c[k] for k in range(s) for j in range(s) for i in range(s)]).simplify(), sympy.Rational(1,24) )
return [Eq_4_1, Eq_4_2, Eq_4_3, Eq_4_4]
# =============================================================================
# Calcul des conditions de consistance pour un schéma RK
# On doit indiquer
# le nombre d'étages s
# le type de schéma parmi implicite, semi-implicite ou explicite
# éventuellement les coefficients A, b et c
# =============================================================================
# Exemple sans paramètres
s = 2
type = "implicite"
display(Markdown(f"<mark>Conditions pour un schéma {type} avec s = {s} étages</mark>"))
ordre_RK(s=s, type="implicite")
# Exemple : méthode à 2 étages, correspondant à Crank-Nicolson
A = [ [sympy.S(0) , sympy.S(0)],
[sympy.Rational(1,2) , sympy.Rational(1,2)]
]
b = [sympy.Rational(1,2), sympy.Rational(1,2)]
c = [sympy.S(0), sympy.S(1)]
s = len(c)
display(Markdown(f"<mark>Exemple de méthode de Runge-Kutta à {s} étages (Crank-Nicolson)</mark>"))
ordre_RK(s, A=A, b=b, c=c);
Conditions pour un schéma implicite avec s = 2 étages
\(\displaystyle \left[\begin{matrix}c_{0} & a_{0, 0} & a_{0, 1}\\c_{1} & a_{1, 0} & a_{1, 1}\\ & 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_{0, 0} + a_{0, 1} = c_{0}\)
\(\displaystyle a_{1, 0} + a_{1, 1} = 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 b_{0} c_{0} a_{0, 0} + b_{0} c_{1} a_{0, 1} + b_{1} c_{0} a_{1, 0} + b_{1} c_{1} a_{1, 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 b_{0} c_{0}^{2} a_{0, 0} + b_{0} c_{0} c_{1} a_{0, 1} + b_{1} c_{0} c_{1} a_{1, 0} + b_{1} c_{1}^{2} a_{1, 1} = \frac{1}{8}\)
\(\displaystyle b_{0} c_{0}^{2} a_{0, 0} + b_{0} c_{1}^{2} a_{0, 1} + b_{1} c_{0}^{2} a_{1, 0} + b_{1} c_{1}^{2} a_{1, 1} = \frac{1}{12}\)
\(\displaystyle b_{0} c_{0} a_{0, 0}^{2} + b_{0} c_{0} a_{0, 1} a_{1, 0} + b_{0} c_{1} a_{0, 0} a_{0, 1} + b_{0} c_{1} a_{0, 1} a_{1, 1} + b_{1} c_{0} a_{0, 0} a_{1, 0} + b_{1} c_{0} a_{1, 0} a_{1, 1} + b_{1} c_{1} a_{0, 1} a_{1, 0} + b_{1} c_{1} a_{1, 1}^{2} = \frac{1}{24}\)
Exemple de méthode de Runge-Kutta à 2 étages (Crank-Nicolson)
\(\displaystyle \left[\begin{matrix}0 & 0 & 0\\1 & \frac{1}{2} & \frac{1}{2}\\ & \frac{1}{2} & \frac{1}{2}\end{matrix}\right]\)
On a 3 conditions pour avoir consistance (= pour être d’ordre 1)
\(\displaystyle \text{True}\)
\(\displaystyle \text{True}\)
\(\displaystyle \text{True}\)
On doit ajouter 1 condition pour être d’ordre 2
\(\displaystyle \text{True}\)
On doit ajouter 2 conditions pour être d’ordre 3
\(\displaystyle \text{False}\)
\(\displaystyle \text{False}\)
On doit ajouter 4 conditions pour être d’ordre 4
\(\displaystyle \text{False}\)
\(\displaystyle \text{False}\)
\(\displaystyle \text{False}\)
\(\displaystyle \text{False}\)
Lien entre l’ordre de convergence \(\omega\) et le nombre d’étages \(s\)
Barrière de Butcher
L’ordre \(\omega \geq 1\) d’une méthode de Runge-Kutta à \(s \geq 1\) étapes convergente satisfait les inégalités suivantes :
\(
\omega \leq
\begin{cases}
2s, & \text{si la méthode est implicite}, \\
s, & \text{si la méthode est explicite et } s < 5, \\
s - 1, & \text{si la méthode est explicite et } s \geq 5.
\end{cases}
\)
Exemples de méthodes
- Schéma d’Euler explicite : méthode RK explicite à 1 étage, d’ordre 1 (\(\omega = s < 5\)).
- Schéma d’Euler implicite : méthode RK implicite à 1 étage, d’ordre 1 (\(\omega \leq 2s\)).
- Schéma d’Euler modifié : méthode RK explicite à 2 étages, d’ordre 2 (\(\omega = s < 5\)).
- Schéma de Heun : méthode RK explicite à 2 étages, d’ordre 2 (\(\omega = s < 5\)).
- Schéma de Crank-Nicholson : méthode RK implicite à 2 étages, d’ordre 2 (\(\omega \leq 2s\)).
Les résultats ci-dessus permettent de comprendre la relation entre \(s\) le nombre d’étages d’une méthode de Runge-Kutta et son ordre de convergence \(\omega\).
Méthode implicite : l’ordre maximal est \(2s\), ce qui permet d’obtenir un meilleur ordre de convergence avec moins d’étages, mais au coût d’une complexité supplémentaire (résolution d’un système non-linéaire).
Méthode explicite : si \(s< 5\), l’ordre maximal est \(s\); si \(s\geq 5\), l’ordre est limité à \(s-1\) (strictement inférieur au nombre d’étages donc). Autrement dit, un schéma RK explicite à \(s\) étages est au mieux d’ordre \(s\), mais cela seulement si \(s < 5\). Si \(s\ge5\), alors il sera au mieux d’ordre \(s-1\), mais ce n’est même pas garanti. En effet, les travaux de Butcher démontrent que pour des méthodes d’ordre 7 et 8, un minimum de 9 et 11 étages, respectivement, est requis. En général, le nombre minimum précis d’étages pour qu’une méthode explicite de Runge-Kutta ait l’ordre \(\omega\) reste un problème ouvert. Certaines valeurs sont connues :
\(
\begin{array}{c|cccccccc}
\omega & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\
\hline
s \text{ minimal} & 1 & 2 & 3 & 4 & 6 & 7 & 9 & 11\\
\end{array}
\)
Théorème de Kuntzmann et Butcher pour des méthodes d’ordre 2s
Soit les conditions suivantes:
- \(B_p\): \(\sum_i b_i c_i^{q-1}=\frac{1}{q}\) pour \(q=1,\dots,p\)
- \(C_\eta\): \(\sum_j a_{ij} c_j^{q-1}=\frac{c_i^q}{q}\) pour \(i=1,\dots,s\) et \(q=1,\dots,\eta\)
- \(D_\zeta\): \(\sum_i b_i c_i^{q-1}a_{ij} =\frac{b_j}{q}(1-c_j^q)\) pour \(j=1,\dots,s\) et \(q=1,\dots,\zeta\)
Si les coefficients de la méthode de RK vérifient les conditions \(B_p\), \(C_\eta\) et \(D_\zeta\) pour \(p\le\eta + \zeta + 1\) et \(p \le 2\eta + 2\) alors la méthode est d’ordre \(p\).
Cf. E. Hairer, S.P. Norsett, G. Wanner, Solving Ordinary Differential Equations I, page 208
Remarque (le cas des systèmes) Une méthode de Runge-Kutta peut être aisément étendue aux systèmes d’équations différentielles ordinaires. Néanmoins, l’ordre d’une méthode RK dans le cas scalaire ne coı̈ncide pas nécessairement avec celui de la méthode analogue dans le cas vectoriel. En particulier, pour \(p \ge 4\), une méthode d’ordre \(p\) dans le cas du système autonome \(\mathbf{y}'(t)=\boldsymbol{\varphi}(t,\mathbf{y})\), avec \(\boldsymbol{\varphi}\colon \mathbb{R}^m\to\mathbb{R}^n\) est encore d’ordre \(p\) quand on l’applique à l’équation scalaire autonome \(\mathbf{y}'(t)=\boldsymbol{\varphi}(\mathbf{y})\), mais la réciproque n’est pas vraie.
A-stabilité
Rappelons la définition de schéma A-stable:
Soit \(\lambda\in\mathbb{C}\) avec \(\Re(\lambda)=-\beta\) où \(\beta\) est un nombre réel positif et considérons le problème de Cauchy \(\begin{cases}
y'(t)=\lambda y(t), &\text{pour }t>0,\\
y(0)=y_0
\end{cases}\) où \(y_0\neq0\) est une valeur donnée. Sa solution est \(y(t)=y_0e^{\lambda t}\) donc \(\lim_{t\to+\infty}|y(t)|=0.\) Soit \(h>0\) un pas de temps donné, \(t_n=nh\) pour \(n\in\mathbb{N}\) et notons \(u_n\approx y(t_n)\) une approximation de la solution \(y\) au temps \(t_n\). Si, sous d’éventuelles conditions sur \(h\), on a \( \lim_{n\to+\infty} u_n =0, \) alors on dit que le schéma est A-stable.
Dans le cas \(\lambda\in\mathbb{R}^-\), i.e. \(-\lambda=\beta>0\), une méthode de Runge-Kutta à \(s\ge1\) étages pour \(y'(t)=-\beta y(t)\) s’écrit: \(\begin{cases}
u_0=y_0,\\
K_i=\displaystyle -\beta \left( u_n+h\sum_{j=1}^{s}a_{ij}K_j \right) & i=1,\dots s\\
u_{n+1}=u_{n}+h \displaystyle\sum_{i=1}^sb_iK_i& n=0,1,\dots N-1.
\end{cases}\)
On peut vérifier que, contrairement aux méthodes multi-pas, la taille des régions de stabilité absolue des méthodes RK augmente quand l’ordre augmente.
Régions de stabilité absolue pour certaines méthodes RK explicites.
Source: A. Quarteroni, F. Saleri, P. Gervasio Calcul Scientifique: Cours, exercices corrigés et illustrations en Matlab et Octave.
Exercice : schémas RK à un étage
Étudier les méthodes RK avec \(s=1\) consistants :
- écrire le schéma implicite générique RK à un étage
- étudier son ordre de convergence
- écrire le schéma explicite générique RK à un étage
- étudier son ordre de convergence
- étudier la A-stabilité
Correction : écriture d’un schéma RK à un étage quelconque
Ces méthodes ont pour matrice de Butcher \(
\begin{array}{c|cc}
c_1 & a_{11}\\
\hline
& b_1
\end{array}
\text{ avec }
\begin{cases}
c_1=a_{11}\\
b_1=1
\end{cases}
\quad
\leadsto
\begin{array}{c|cc}
c_1 & c_1\\
\hline
& 1
\end{array}
\quad
\leadsto
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n+hc_1,u_n+h c_{1}K_1\right)\\
u_{n+1} = u_n + h K_1 & n=0,1,\dots N-1
\end{cases}
\)
⚙️ Exemple de schéma implicite
Par exemple, si on choisit \(c_1=1\) on trouve le schéma implicite
\(
\begin{array}{c|cc}
1 & 1\\
\hline
& 1
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_{n+1},u_n+h K_1\right)\\
u_{n+1} = u_n + h K_1 & n=0,1,\dots N-1
\end{cases}
\)
Notons que ce schéma n’est rien d’autre que le schéma d’Euler implicite car \(u_{n+1} = u_n+h K_1\), donc si on remplace le \(u_n+h K_1\) dans la définition de \(K_1\) par \(u_{n+1}\) on obtient
\(
K_1
= \varphi\left(t_{n+1},u_n+h K_1\right)
= \varphi\left(t_{n+1},u_{n+1}\right)
\quad\implies\quad
u_{n+1} = u_n + h K_1=u_n + h \varphi\left(t_{n+1},u_{n+1}\right)
\)
En effet, \(c_s=1\) et \(a_{sj}=b_j\) pour tout \(j=1,\dots,s\) alors \(K_s=\varphi(t_{n+1},u_{n+1})\).
Correction : étude de l’ordre d’un schéma RK à un étage quelconque
Selon le théorème de Butcher :
- si le schéma est implicite, il est au mieux d’ordre \(\omega\le 2s=2\);
- si le schéma est explicite, il est au mieux d’ordre \(\omega\le s=1\).
Code
ordre_RK(1)
\(\displaystyle \left[\begin{matrix}c_{0} & a_{0, 0}\\ & b_{0}\end{matrix}\right]\)
On a 2 conditions pour avoir consistance (= pour être d’ordre 1)
\(\displaystyle b_{0} = 1\)
\(\displaystyle a_{0, 0} = c_{0}\)
On doit ajouter 1 condition pour être d’ordre 2
\(\displaystyle b_{0} c_{0} = \frac{1}{2}\)
On doit ajouter 2 conditions pour être d’ordre 3
\(\displaystyle b_{0} c_{0}^{2} = \frac{1}{3}\)
\(\displaystyle b_{0} c_{0} a_{0, 0} = \frac{1}{6}\)
On doit ajouter 4 conditions pour être d’ordre 4
\(\displaystyle b_{0} c_{0}^{3} = \frac{1}{4}\)
\(\displaystyle b_{0} c_{0}^{2} a_{0, 0} = \frac{1}{8}\)
\(\displaystyle b_{0} c_{0}^{2} a_{0, 0} = \frac{1}{12}\)
\(\displaystyle b_{0} c_{0} a_{0, 0}^{2} = \frac{1}{24}\)
\(\displaystyle \left[ \left[ b_{0} = 1, \ a_{0, 0} = c_{0}\right], \ \left[ b_{0} c_{0} = \frac{1}{2}\right], \ \left[ b_{0} c_{0}^{2} = \frac{1}{3}, \ b_{0} c_{0} a_{0, 0} = \frac{1}{6}\right], \ \left[ b_{0} c_{0}^{3} = \frac{1}{4}, \ b_{0} c_{0}^{2} a_{0, 0} = \frac{1}{8}, \ b_{0} c_{0}^{2} a_{0, 0} = \frac{1}{12}, \ b_{0} c_{0} a_{0, 0}^{2} = \frac{1}{24}\right]\right]\)
⚙️ Exemple de schéma (implicite) d’ordre maximal
Si on cherche une méthode d’ordre \(2\) alors il faut que \(1\times c_1=\frac{1}{2}\): il existe une seule méthode d’ordre 2, elle est implicite et sa matrice de Butcher est
\(
\begin{array}{c|cc}
\frac{1}{2} & \frac{1}{2}\\
\hline
& 1
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}K_1\right)\\
u_{n+1} = u_n + h K_1 & n=0,1,\dots N-1
\end{cases}
\)
Notons que \(u_{n+1} = u_n+h K_1\), donc \(hK_1=u_{n+1}-u_n\). Si on remplace le \(hK_1\) dans la définition de \(K_1\) ou trouve \(K_1 = \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}K_1\right)=\varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}u_{n+1}-\frac{h}{2}u_n\right)=\varphi\left(t_n+\frac{h}{2},h\frac{(u_{n+1}+u_n)}{2}\right)\): on obtient le schéma RK1_M
\(
K_1
= \varphi\left(t_{n}+\frac{h}{2},\frac{u_n+u_{n+1}}{2}\right)
\quad\implies\quad
u_{n+1} = u_n + h K_1=u_n + h \varphi\left(\frac{t_n+t_{n+1}}{2},\frac{u_n+u_{n+1}}{2}\right)
\)
Correction : écriture d’un schéma RK à un étage explicite
Si le schéma est explicite alors
\(
\begin{array}{c|cc}
0 & 0\\
\hline
& 1
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
u_{n+1} = u_n + h K_1 & n=0,1,\dots N-1
\end{cases}
\)
Il existe un seul schéma explicite RK à un étage, il s’agit du schéma d’Euler explicite.
Correction : étude de l’ordre d’un schéma RK à un étage explicite
Si le schéma est explicite alors il est d’ordre 1.
Correction : étude de la A-stabilité
Une méthode de Runge-Kutta à 1 étage pour \(\varphi(t,y)=-\beta y\) s’écrit :
\(\begin{cases}
u_0=y_0,\\
u_{n+1}=u_{n}+h K_1& n=0,1,\dots N-1\\
K_1=-\beta \left( u_n+hc_1K_1 \right).
\end{cases}\)
On trouve donc \((1+\beta hc_1)K_1=-\beta u_n\) ainsi
\(\begin{cases}
u_0 = y_0 \\
u_{n+1} = \left(1-\frac{(\beta h)}{1+c_1(\beta h)}\right)u_n & n=0,1,\dots N-1
\end{cases}\)
Par induction on obtient
\(
u_{n}=\left(1-\frac{\beta h}{1+c_1\beta h}\right)^nu_0.
\)
Par conséquent, \(\lim\limits_{n\to+\infty}u_n=0\) si et seulement si
\(
\left|1-\frac{\beta h}{1+c_1\beta h}\right|<1.
\)
Notons \(x\) le produit \(\beta h>0\) (donc \(x>0\)) et \(q\) la fonction \(q(x)=1-\frac{x}{1+c_1x}\) (avec \(c_1\in[0;1]\)).
\(\begin{align*}
|q(x)|<1
&\quad\iff\quad
-1<1-\frac{x}{1+c_1x}<1
\quad\iff\quad
\begin{cases}
\frac{x}{1+c_1x}<2\\
\frac{x}{1+c_1x}>0
\end{cases}
\\
&\quad\stackrel{\substack{x>0\\c_1\ge0}}{\iff}\quad
\frac{x}{1+c_1x}<2
\quad\stackrel{c_1x\ge0}{\iff}\quad
x<2(1+c_1x)
\\
&\quad\iff\quad
(1-2c_1)x<2
\quad\iff\quad
\begin{cases}
x<\frac{2}{1-2c_1}&\text{si }1-2c_1>0\\
\forall x&\text{sinon.}
\end{cases}
\end{align*}\) Conclusion:
- si \(c_1\ge \frac{1}{2}\) le schéma est inconditionnellement A-stable,
- si \(c_1<\frac{1}{2}\) le schéma est A-stable ssi \(\beta h <\frac{2}{1-2c_1}\).
Pour nos exemples, on obtient:
- si \(c_1=1\) (Euler Implicite) alors le schéma est inconditionnellement A-stable et d’ordre 1,
- si \(c_1=\frac{1}{2}\) (RK1_M) alors le schéma est inconditionnellement A-stable et d’ordre 2,
- si \(c_1=0\) (Euler Explicite) alors le schéma est A-stable ssi \(h<\frac{2}{\beta}\) et il est d’ordre 1.
Exercice : schémas RK à deux étages
Étudier les méthodes RK avec \(s=2\) consistantes:
- écrire le schéma implicite générique RK à deux étages
- écrire le schéma semi-implicite générique RK à deux étages
- écrire le schéma explicite générique RK à deux étages
- étudier l’ordre de convergences des schémas explicites
- étudier la A-stabilité des schémas explicites
Correction : écriture d’un schéma RK à deux étages quelconque
Les méthodes avec \(s=2\) ont pour matrice de Butcher
\(
\begin{array}{c|cc}
c_1 & a_{11} & a_{12}\\
c_2 & a_{21} & a_{22}\\
\hline
& b_1 & b_2
\end{array}
\text{ avec }
\begin{cases}
c_1=a_{11}+a_{12}\\
c_2=a_{21}+a_{22}\\
b_1+b_2=1
\end{cases}
\quad
\implies
\begin{array}{c|cc}
c_1 & a_{11} & c_1-a_{11}\\
c_2 & a_{21} & c_2-a_{21}\\
\hline
& 1-b_2 & b_2
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n+hc_1,u_n+h(a_{11}K_1+ (c_1-a_{11})K_2)\right)\\
K_2 = \varphi\left(t_n+hc_2,u_n+h(a_{21}K_1+ (c_2-a_{21})K_2)\right)\\
u_{n+1} = u_n + h\left((1-b_2)K_1+b_2K_2\right) & n=0,1,\dots N-1
\end{cases}
\)
Code
def RK(tt,phi,y0):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for i in range(len(tt)-1):
# X = [K1,K2]
sys = lambda X : [ -X[0]+phi( tt[i]+h*c1 , uu[i]+h*(a11*X[0]+a12*X[1]) ) ,
-X[1]+phi( tt[i]+h*c2 , uu[i]+h*(a21*X[0]+a22*X[1]) )
]
K_start = phi(tt[i],uu[i])
K1, K2 = fsolve( sys , [K_start,K_start ] )
uu[i+1] = uu[i]+h*(b1*K1+b2*K2)
return uu
Voici un exemple, appelé méthode de Gauss: \(
\begin{array}{c|cc}
\frac{1}{2}-\gamma & \frac{1}{4} & \frac{1}{4}-\gamma\\
\frac{1}{2}+\gamma & \frac{1}{4}+\gamma & \frac{1}{4}\\
\hline
& \frac{1}{2} & \frac{1}{2}
\end{array}
\qquad\text{avec }\gamma=\frac{\sqrt{3}}{6}
\)
Cette méthode est d’ordre 4 et on peut monter que c’est la seule méthode d’ordre 4 à deux étages.
Code
gamma = sympy.sqrt(3)/6
c = [ sympy.Rational(1,2)-gamma, sympy.Rational(1,2)+gamma ]
b = [ sympy.Rational(1,2), sympy.Rational(1,2) ]
A = [ [ sympy.Rational(1,4), sympy.Rational(1,4)-gamma ] ,
[ sympy.Rational(1,4)+gamma, sympy.Rational(1,4) ] ]
ordre_RK(2,A=A,b=b,c=c)
\(\displaystyle \left[\begin{matrix}\frac{1}{2} - \frac{\sqrt{3}}{6} & \frac{1}{4} & \frac{1}{4} - \frac{\sqrt{3}}{6}\\\frac{\sqrt{3}}{6} + \frac{1}{2} & \frac{1}{4} + \frac{\sqrt{3}}{6} & \frac{1}{4}\\ & \frac{1}{2} & \frac{1}{2}\end{matrix}\right]\)
On a 3 conditions pour avoir consistance (= pour être d’ordre 1)
\(\displaystyle \text{True}\)
\(\displaystyle \text{True}\)
\(\displaystyle \text{True}\)
On doit ajouter 1 condition pour être d’ordre 2
\(\displaystyle \text{True}\)
On doit ajouter 2 conditions pour être d’ordre 3
\(\displaystyle \text{True}\)
\(\displaystyle \text{True}\)
On doit ajouter 4 conditions pour être d’ordre 4
\(\displaystyle \text{True}\)
\(\displaystyle \text{True}\)
\(\displaystyle \text{True}\)
\(\displaystyle \text{True}\)
\(\displaystyle \left[ \left[ \text{True}, \ \text{True}, \ \text{True}\right], \ \left[ \text{True}\right], \ \left[ \text{True}, \ \text{True}\right], \ \left[ \text{True}, \ \text{True}, \ \text{True}, \ \text{True}\right]\right]\)
Code
# gamma = sympy.sqrt(3)/6
# c = [ sympy.Rational(1,2)-gamma, sympy.Rational(1,2)+gamma ]
# b = [ sympy.Rational(1,2), sympy.Rational(1,2) ]
# A = [ [sympy.Rational(1,4), sympy.Rational(1,4)-gamma],
# [sympy.Rational(1,4)+gamma, sympy.Rational(1,4)]
# ]
s = len(c)
print("🌟 Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i]).simplify()==c[i])
print("\n🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:")
print(sum([b[i]*c[i] for i in range(s)]).simplify()==sympy.Rational(1,2))
print("\n🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==sympy.Rational(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,6))
print("\n🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**3 for i in range(s)]).simplify()==sympy.Rational(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,12))
print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]).simplify()==sympy.Rational(1,24))
🌟 Consistance
True
True
True
🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:
True
🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:
True
True
🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:
True
True
True
True
Correction : écriture d’un schéma RK à deux étages semi-implicite
Les méthodes semi-implicites vérifient, de plus, \(c_1=a_{11}\) ainsi la matrice de Butcher est \(
\begin{array}{c|cc}
c_1 & c_1 & 0\\
c_2 & a_{21} & c_2-a_{21}\\
\hline
& 1-b_2 & b_2
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n+hc_1,u_n+hc_{1}K_1\right)\\
K_2 = \varphi\left(t_n+hc_2,u_n+h(a_{21}K_1+ (c_2-a_{21})K_2)\right)\\
u_{n+1} = u_n + h\left((1-b_2)K_1+b_2K_2\right) & n=0,1,\dots N-1
\end{cases}
\)
Code
def RK(tt,phi,y0):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for i in range(len(tt)-1):
# X = [K1,K2]
K_start = phi(tt[i],uu[i])
K1 = fsolve( lambda x : -x+phi( tt[i]+h*c1 , uu[i]+h*(a11*x) ) , K_start ) [0]
K2 = fsolve( lambda x : -x+phi( tt[i]+h*c2 , uu[i]+h*(a21*K1+a22*x) ) , K_start ) [0]
uu[i+1] = uu[i]+h*(b1*K1+b2*K2)
return uu
Par exemple, si on choisit \(c_1=0\), \(c_2=1\) et \(a_{21}=a_{22}=b_1=b_2=\frac{1}{2}\) on trouve le schéma semi-implicite \(
\begin{array}{c|cc}
0 & 0 & 0\\
1 & \frac{1}{2} & \frac{1}{2}\\
\hline
& \frac{1}{2} & \frac{1}{2}
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_{n},u_n\right)\\
K_2 = \varphi\left(t_{n+1},u_n+\frac{h}{2}(K_1+K_2)\right)\\
u_{n+1} = u_n + \frac{h}{2}(K_1+K_2) & n=0,1,\dots N-1
\end{cases}
\)
Notons que ce schéma n’est rien d’autre que le schéma de Crank-Nicolson car \(u_{n+1} = u_n+\frac{h}{2}(K_1+K_2)\), donc si on remplace le \(u_n+\frac{h}{2}(K_1+K_2)\) dans la définition de \(K_2\) par \(u_{n+1}\) on obtient \(
K_2
= \varphi\left(t_{n+1},u_n+\frac{h}{2}(K_1+K_2)\right)
= \varphi\left(t_{n+1},u_{n+1}\right)
\) ainsi \(
u_{n+1} = u_n+\frac{h}{2}(K_1+K_2)=u_n + \frac{h}{2}(\varphi\left(t_{n},u_n\right)+\varphi\left(t_{n+1},u_{n+1}\right))
\) En effet, \(c_s=1\) et \(a_{sj}=b_j\) pour tout \(j=1,\dots,s\) alors \(K_s=\varphi(t_{n+1},u_{n+1})\).
Il s’agit donc d’un schéma d’ordre 2.
Code
c = [0, 1]
b = [sympy.Rational(1,2), sympy.Rational(1,2)]
A = [[0, 0], [sympy.Rational(1,2), sympy.Rational(1,2)]]
s = len(c)
print("🌟 Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i])==c[i])
print("\n🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:")
print(sum([b[i]*c[i] for i in range(s)]).simplify()==sympy.Rational(1,2))
print("\n🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==sympy.Rational(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,6))
🌟 Consistance
True
True
True
🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:
True
🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:
False
False
Un autre exemple est la méthode dite DIRK (Diagonally Implicit Runge-Kutta) qui est d’ordre 3: \(
\begin{array}{c|cc}
\frac{1}{3} & \frac{1}{3} & 0\\
1 & 1 & 0\\
\hline
& \frac{3}{4} & \frac{1}{4}
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_{n}+\frac{1}{3}h,u_n+\frac{1}{3}hK_1\right)\\
K_2 = \varphi\left(t_{n+1},u_n+hK_1\right)\\
u_{n+1} = u_n + \frac{h}{4}(3K_1+K_2) & n=0,1,\dots N-1
\end{cases}
\)
Code
c = [sympy.Rational(1,3), 1]
b = [sympy.Rational(3,4), sympy.Rational(1,4)]
A = [[sympy.Rational(1,3),0],[1,0]]
s = len(c)
print("🌟 Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i])==c[i])
print("\n🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:")
print(sum([b[i]*c[i] for i in range(s)]).simplify()==sympy.Rational(1,2))
print("\n🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==sympy.Rational(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,6))
print("\n🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**3 for i in range(s)]).simplify()==sympy.Rational(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,12))
print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]).simplify()==sympy.Rational(1,24))
🌟 Consistance
True
True
True
🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:
True
🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:
True
True
🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:
False
False
False
False
Correction : écriture d’un schéma RK à deux étages explicite
Les méthodes avec \(s=2\) explicites ont pour matrice de Butcher \(
\begin{array}{c|cc}
0 & 0 & 0\\
c_2 & c_{2} & 0\\
\hline
& 1-b_2 & b_2
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+hc_2,u_n+hc_2K_1\right)\\
u_{n+1} = u_n + h\left((1-b_2)K_1+b_2K_2\right) & n=0,1,\dots N-1
\end{cases}
\)
Code
def RK(tt,phi,y0):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for i in range(len(tt)-1):
K1 = phi( tt[i], uu[i] )
K2 = phi( tt[i]+h*c2 , uu[i]+h*(c2*K1) )
uu[i+1] = uu[i]+h*(b1*K1+b2*K2)
return uu
Correction : étude de l’ordre d’un schéma RK à deux étages explicite
- Un schéma RK à deux étage explicite est au mieux d’ordre 2.
- Pour avoir l’ordre 2 il faut que \(b_2c_2=\frac{1}{2}\).
Conclusion: un schéma RK à deux étages explicite d’ordre 2 s’écrit \(
\begin{array}{c|cc}
0 & 0 & 0\\
\frac{1}{2\alpha} & \frac{1}{2\alpha} & 0\\
\hline
& 1-\alpha & \alpha
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+\frac{h}{2}\alpha,u_n+\frac{h}{2}\alpha K_1\right)\\
u_{n+1} = u_n + h\left((1-\alpha)K_1+\alpha K_2\right) & n=0,1,\dots N-1
\end{cases}
\)
Voici quelques cas particuliers :
\(\alpha=\frac{1}{2}\): schéma de Heun \(
\begin{array}{c|cc}
0 & 0 & 0\\
1 & 1 & 0\\
\hline
& \frac{1}{2} & \frac{1}{2}
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi(t_n,u_n)\\
K_2 = \varphi(t_{n+1},u_{n}+hK_1)\\
u_{n+1} = u_n + h\left(\frac{1}{2}K_1+\frac{1}{2}K_2\right) & n=0,1,\dots N-1
\end{cases}
\)
\(\alpha=1\): schéma d’Euler Modifié (ou du point milieu) \(
\begin{array}{c|cc}
0 & 0 & 0\\
\frac{1}{2} & \frac{1}{2} & 0\\
\hline
& 0 & 1
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi(t_n,u_n)\\
K_2 = \varphi(t_{n}+\frac{1}{2}h,u_{n}+\frac{1}{2}hK_1)\\
u_{n+1} = u_n + h K_2 & n=0,1,\dots N-1
\end{cases}
\)
\(\alpha=\frac{2}{3}\): schéma de Ralston \(
\begin{array}{c|cc}
0 & 0 & 0\\
\frac{3}{4} & \frac{3}{4} & 0\\
\hline
& \frac{1}{3} & \frac{2}{3}
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi(t_n,u_n)\\
K_2 = \varphi(t_{n}+\frac{3}{4}h,u_{n}+\frac{3}{4}hK_1)\\
u_{n+1} = u_n + \frac{h}{3} (K_1+2K_2) & n=0,1,\dots N-1
\end{cases}
\)
Correction : étude de la A-stabilité des schémas explicites
Le schéma donné appliqué à cette équation devient \(\begin{cases}
u_0 = y_0 \\
u_{n+1} = \left(1-(\beta h)+c_2b_2(\beta h)^2\right)u_n & n=0,1,\dots N-1
\end{cases}\) Par induction on obtient \(
u_{n}=\left(1-(\beta h)+c_2b_2(\beta h)^2\right)^nu_0.
\) Par conséquent, \(\lim\limits_{n\to+\infty}u_n=0\) si et seulement si \(
\left|1-(\beta h)+c_2b_2(\beta h)^2\right|<1.
\) Notons \(x\) le produit \(\beta h>0\) (donc \(x>0\)) et \(q\) le polynôme \(q(x)=1-x+c_2b_2x^2\).
Si \(c_2b_2=0\), le schéma se réduit au schéma d’Euler explicite.
Si \(c_2b_2>0\), c’est une parabole convexe et le sommet est situé en \(\left(\frac{1}{2c_2b_2}>0,1-\frac{1}{4c_2b_2}\right)\). \(
|q(x)|<1
\quad\iff\quad
x<\frac{1}{c_2b_2}
\) En particulier, lorsque le schéma est d’ordre 2 on a \(b_2c_2=\frac{1}{2}\) et le schéma est A-stable ssi \(h<\frac{2}{\beta}\).
Si \(c_2b_2<0\), c’est une parabole concave et le sommet est situé en \(\left(\frac{1}{2c_2b_2}<0,1-\frac{1}{4c_2b_2}\right)\). \(
|q(x)|<1
\quad\iff\quad
x<\frac{1-\sqrt{1-8b_2c_2}}{2b_2c_2}
\)
Exemples à 3 étages explicites
Un schéma RK à \(3\) étages explicite a pour matrice de Butcher \(
\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\
c_2 & a_{21} & 0&0\\
c_3 & a_{31} &a_{32}&0\\
\hline
& b_1 & b_2 & b_3
\end{array}
\text{ avec }
\begin{cases}
c_2=a_{21}\\
c_3=a_{31}+a_{32}\\
b_1+b_2+b_3=1
\end{cases}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+hc_2,u_n+hc_2K_1\right)\\
K_3 = \varphi\left(t_n+hc_3,u_n+h(a_{31}K_1+ a_{32}K_2)\right)\\
u_{n+1} = u_n + h\left(b_1K_1+b_2K_2+b_3K_3\right) & n=0,1,\dots N-1
\end{cases}
\)
Voici trois exemples:
- schéma de Heun à trois étages \(
\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}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+\frac{h}{3},u_n+\frac{h}{3}K_1\right)\\
K_3 = \varphi\left(t_n+\frac{2}{3}h,u_n+\frac{2}{3}hK_2\right)\\
u_{n+1} = u_n + \frac{h}{4}\left(K_1+3K_3\right) & n=0,1,\dots N-1
\end{cases}
\)
Étant une méthode explicite à 3 étapes, elle est au mieux d’ordre 3. Vérifions qu’elle est effectivement d’ordre 3:
Code
c = [0, sympy.Rational(1,3), sympy.Rational(2,3)]
b = [sympy.Rational(1,4), 0, sympy.Rational(3,4)]
A = [[0, 0, 0],
[sympy.Rational(1,3), 0, 0],
[0, sympy.Rational(2,3), 0]]
s = len(c)
print("🌟 Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i])==c[i])
print("\n🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:")
print(sum([b[i]*c[i] for i in range(s)]).simplify()==sympy.Rational(1,2))
print("\n🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==sympy.Rational(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,6))
print("\n🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**3 for i in range(s)]).simplify()==sympy.Rational(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,12))
print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]).simplify()==sympy.Rational(1,24))
🌟 Consistance
True
True
True
True
🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:
True
🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:
True
True
🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:
False
False
False
False
- schéma de 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}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}K_1\right)\\
K_3 = \varphi\left(t_{n+1},u_n+h(-K_1+2K_2)\right)\\
u_{n+1} = u_n + \frac{h}{6}\left(K_1+4K_2+K_3\right) & n=0,1,\dots N-1
\end{cases}
\)
Étant une méthode explicite à 3 étapes, elle a au mieux d’ordre 3. On vérifie ci-dessous qu’elle est d’ordre 3:
Code
c = [0, sympy.S(1)/2, 1]
b = [sympy.S(1)/6, sympy.S(2)/3, sympy.S(1)/6]
A = [[0, 0, 0],
[sympy.S(1)/2, 0, 0],
[-1, 2, 0]]
s = len(c)
print("🌟 Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i])==c[i])
print("\n🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:")
print(sum([b[i]*c[i] for i in range(s)]).simplify()==sympy.Rational(1,2))
print("\n🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==sympy.Rational(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,6))
print("\n🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**3 for i in range(s)]).simplify()==sympy.Rational(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,12))
print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]).simplify()==sympy.Rational(1,24))
🌟 Consistance
True
True
True
True
🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:
True
🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:
True
True
🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:
True
False
True
False
- schéma \(
\begin{array}{c|ccc}
0 & 0 & 0 & 0 \\
\frac{1}{2} & \frac{1}{2} & 0&0\\
1 & -1 &2&0\\
\hline
& -\frac{1}{6} & \frac{4}{3} & -\frac{1}{6}
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}K_1\right)\\
K_3 = \varphi\left(t_{n+1},u_n+h(-K_1+2K_2)\right)\\
u_{n+1} = u_n + \frac{h}{6}\left(-K_1+8K_2-K_3\right) & n=0,1,\dots N-1
\end{cases}
\)
Étant une méthode explicite à 3 étapes, elle a au mieux d’ordre 3. On vérifie ci-dessous qu’elle n’est que d’ordre 2:
Code
c = [0, sympy.Rational(1,2), 1]
b = [-sympy.Rational(1,6), sympy.Rational(4,3), -sympy.Rational(1,6)]
A = [[0, 0, 0], [sympy.Rational(1,2), 0, 0], [-1, 2, 0]]
s = len(c)
print("🌟 Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i])==c[i])
print("\n🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:")
print(sum([b[i]*c[i] for i in range(s)]).simplify()==sympy.Rational(1,2))
print("\n🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==sympy.Rational(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,6))
🌟 Consistance
True
True
True
True
🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:
True
🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:
False
False
- schéma Radau-IIa (d’ordre 5) \(
\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}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n+\frac{4-\sqrt{6}}{10}h,u_n+h\left(\frac{88-7\sqrt{6}}{360} K_1+\frac{296-169\sqrt{6}}{1800} K_2+\frac{-2+3\sqrt{6}}{225}K_3\right)\right)\\
K_2 = \varphi\left(t_n+\frac{4+\sqrt{6}}{10}h,u_n+h\left(\frac{296+169\sqrt{6}}{1800}K_1+\frac{88+7\sqrt{6}}{360}K_2+\frac{-2-3\sqrt{6}}{225}K_3\right)\right)\\
K_3 = \varphi\left(t_n+h,u_n+h\left(\frac{16-\sqrt{6}}{36}K_1+\frac{16+\sqrt{6}}{36}K_2+\frac{1}{9}K_3\right)\right)\\
u_{n+1} = u_n + h\left(\frac{16-\sqrt{6}}{36}K_1+\frac{16+\sqrt{6}}{36}K_2+\frac{1}{9}K_3\right) & n=0,1,\dots N-1
\end{cases}
\) Notons que \(c_s=1\) et \(a_{sj}=b_j\) pour tout \(j=1,\dots,s\) ainsi \(K_s=\varphi(t_{n+1},u_{n+1})\).
Étant une méthode implicite à 3 étapes, elle a au mieux d’ordre 6. On vérifie ci-dessous qu’elle est d’ordre 5:
Code
sq6 = sympy.sqrt(6)
c = [(4-sq6)/10, (4+sq6)/10, 1]
b = [(16-sq6)/36, (16+sq6)/36, sympy.Rational(1,9)]
A = [[(88-7*sq6)/360, (296-169*sq6)/1800, (-2+3*sq6)/225],
[(296+169*sq6)/1800, (88+7*sq6)/360, (-2-3*sq6)/225],
[(16-sq6)/36, (16+sq6)/36,sympy.S(1)/9] ]
s = len(c)
print("🌟 Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i])==c[i])
print("\n🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:")
print(sum([b[i]*c[i] for i in range(s)]).simplify()==sympy.Rational(1,2))
print("\n🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==sympy.Rational(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,6))
print("\n🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**3 for i in range(s)]).simplify()==sympy.Rational(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,12))
print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]).simplify()==sympy.Rational(1,24))
🌟 Consistance
True
True
True
True
🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:
True
🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:
True
True
🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:
True
True
True
True
Exemples à 4 étages explicites
Un schéma RK à \(4\) étages explicite a pour matrice de Butcher \(
\begin{array}{c|cccc}
0 & 0 & 0 & 0 & 0\\
c_2 & a_{21} & 0&0&0\\
c_3 & a_{31} &a_{32}&0&0\\
c_4 & a_{41} &a_{42}&a_{43}&0\\
\hline
& b_1 & b_2 & b_3 & b_4
\end{array}
\text{ avec }
\begin{cases}
c_2=a_{21}\\
c_3=a_{31}+a_{32}\\
c_4=a_{41}+a_{42}+a_{43}\\
b_1+b_2+b_3+b_4=1
\end{cases}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+hc_2,u_n+h(a_{21}K_1)\right)\\
K_3 = \varphi\left(t_n+hc_3,u_n+h(a_{31}K_1+ a_{32}K_2)\right)\\
K_4 = \varphi\left(t_n+hc_4,u_n+h(a_{41}K_1+ a_{42}K_2+ a_{43}K_3)\right)\\
u_{n+1} = u_n + h\left(b_1K_1+b_2K_2+b_3K_3+b_4K_4\right) & n=0,1,\dots N-1
\end{cases}
\)
Voici deux exemples (le premier est le plus connus):
- schéma RK4-1 (“La” méthode de Runge-Kutta) \(
\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}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2} K_1)\right)\\
K_3 = \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}K_2\right)\\
K_4 = \varphi\left(t_{n+1},u_n+h K_3\right)\\
u_{n+1} = u_n + \frac{h}{6}\left(K_1+2K_2+2K_3+K_4\right) & n=0,1,\dots N-1
\end{cases}
\) Il est bien d’ordre 4:
Code
c = [0, sympy.Rational(1,2), sympy.Rational(1,2), 1]
b = [sympy.Rational(1,6), sympy.Rational(1,3), sympy.Rational(1,3), sympy.Rational(1,6)]
A = [[0, 0, 0, 0],
[sympy.Rational(1,2), 0, 0, 0],
[0, sympy.Rational(1,2), 0, 0],
[0, 0, 1, 0]]
s = len(c)
print("🌟 Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i])==c[i])
print("\n🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:")
print(sum([b[i]*c[i] for i in range(s)]).simplify()==sympy.Rational(1,2))
print("\n🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==sympy.Rational(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,6))
print("\n🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**3 for i in range(s)]).simplify()==sympy.Rational(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,12))
print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]).simplify()==sympy.Rational(1,24))
🌟 Consistance
True
True
True
True
True
🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:
True
🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:
True
True
🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:
True
True
True
True
- schéma 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}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+\frac{h}{4},u_n+\frac{h}{4} K_1)\right)\\
K_3 = \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}K_2\right)\\
K_4 = \varphi\left(t_{n+1},u_n+h (K_1-2K_2+2K_3)\right)\\
u_{n+1} = u_n + \frac{h}{6}\left(K_1+4K_3+K_4\right) & n=0,1,\dots N-1
\end{cases}
\) Il est bien d’ordre 4:
Code
c = [0, sympy.Rational(1,4), sympy.Rational(1,2), 1]
b = [sympy.Rational(1,6), 0, sympy.Rational(2,3), sympy.Rational(1,6)]
A = [ [0, 0, 0, 0],
[sympy.Rational(1,4), 0, 0, 0],
[0, sympy.Rational(1,2), 0, 0],
[1, -2, 2, 0] ]
s = len(c)
print("🌟 Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i])==c[i])
print("\n🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:")
print(sum([b[i]*c[i] for i in range(s)]).simplify()==sympy.Rational(1,2))
print("\n🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==sympy.Rational(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,6))
print("\n🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**3 for i in range(s)]).simplify()==sympy.Rational(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,12))
print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]).simplify()==sympy.Rational(1,24))
🌟 Consistance
True
True
True
True
True
🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:
True
🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:
True
True
🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:
True
True
True
True
- schéma RK4-3 (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}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+\frac{h}{3},u_n+\frac{h}{3} K_1)\right)\\
K_3 = \varphi\left(t_n+\frac{2}{3}h,u_n+\frac{h}{3}(-K_1+3K_2)\right)\\
K_4 = \varphi\left(t_{n+1},u_n+h (K_1-K_2+K_3)\right)\\
u_{n+1} = u_n + \frac{h}{8}\left(K_1+3K_2+3K_3+K_4\right) & n=0,1,\dots N-1
\end{cases}
\) Il est bien d’ordre 4:
Code
c = [0, sympy.Rational(1,3), sympy.Rational(2,3), 1]
b = [sympy.Rational(1,8), sympy.Rational(3,8), sympy.Rational(3,8), sympy.Rational(1,8)]
A = [ [0, 0, 0, 0],
[sympy.Rational(1,3), 0, 0, 0],
[-sympy.Rational(1,3), 1, 0, 0],
[1, -1, 1, 0] ]
s = len(c)
print("🌟 Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i])==c[i])
print("\n🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:")
print(sum([b[i]*c[i] for i in range(s)]).simplify()==sympy.Rational(1,2))
print("\n🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==sympy.Rational(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,6))
print("\n🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**3 for i in range(s)]).simplify()==sympy.Rational(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,12))
print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]).simplify()==sympy.Rational(1,24))
🌟 Consistance
True
True
True
True
True
🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:
True
🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:
True
True
🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:
True
True
True
True
Exemple à 5 étages explicite
Un schéma RK à \(5\) étages explicite a pour matrice de Butcher \(
\begin{array}{c|ccccc}
0 & 0 & 0 & 0 & 0 &0 \\
c_2 & a_{21} & 0&0&0&0\\
c_3 & a_{31} &a_{32}&0&0&0\\
c_4 & a_{41} &a_{42}&a_{43}&0&0\\
c_5 & a_{51} &a_{52}&a_{53}&a_{54}&0\\
\hline
& b_1 & b_2 & b_3 & b_4 & b_5
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+hc_2,u_n+h(a_{21}K_1)\right)\\
K_3 = \varphi\left(t_n+hc_3,u_n+h(a_{31}K_1+ a_{32}K_2\right)\\
K_4 = \varphi\left(t_n+hc_4,u_n+h(a_{41}K_1+ a_{42}K_2+ a_{43}K_3)\right)\\
K_5 = \varphi\left(t_n+hc_5,u_n+h(a_{51}K_1+ a_{52}K_2+ a_{53}K_3+ a_{54}K_4)\right)\\
u_{n+1} = u_n + h\left(b_1K_1+b_2K_2+b_3K_3+b_4K_4+b_5K_5\right) & n=0,1,\dots N-1
\end{cases}
\) avec \(
\begin{cases}
c_2=a_{21}\\
c_3=a_{31}+a_{32}\\
c_4=a_{41}+a_{42}+a_{43}\\
c_5=a_{51}+a_{52}+a_{53}+a_{54}\\
b_1+b_2+b_3+b_4+b_5=1
\end{cases}
\)
Voici un exemple:
- schéma de 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}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+\frac{h}{3},u_n+\frac{h}{3}K_1\right)\\
K_3 = \varphi\left(t_n+\frac{h}{3},u_n+\frac{h}{6}(K_1+K_2)\right)\\
K_4 = \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{8}(K_1+3K_3)\right)\\
K_5 = \varphi\left(t_n+h ,u_n+\frac{h}{2}(K_1-3K_3+4K_4)\right)\\
u_{n+1} = u_n + \frac{h}{6}\left(K_1+4K_4+K_5\right) & n=0,1,\dots N-1
\end{cases}
\)
Code
c = [0, sympy.Rational(1,3), sympy.Rational(1,3), sympy.Rational(1,2), 1]
b = [sympy.Rational(1,6), 0, 0, sympy.Rational(2,3), sympy.Rational(1,6)]
A = [[0, 0, 0, 0, 0], [sympy.Rational(1,3), 0, 0, 0, 0], [sympy.Rational(1,6), sympy.Rational(1,6), 0, 0, 0], [sympy.Rational(1,8), 0, sympy.Rational(3,8), 0, 0], [sympy.Rational(1,2), 0, -sympy.Rational(3,2), 2, 0]]
s = len(c)
print("🌟 Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i])==c[i])
print("\n🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:")
print(sum([b[i]*c[i] for i in range(s)]).simplify()==sympy.Rational(1,2))
print("\n🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==sympy.Rational(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,6))
print("\n🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**3 for i in range(s)]).simplify()==sympy.Rational(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,12))
print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]).simplify()==sympy.Rational(1,24))
print("\n🌟 Pour être d'ordre 5 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**4 for i in range(s)]) == sympy.Rational(1,5))
print(sum([b[i]*c[i]**2*(1-c[i]) for i in range(s)]) + sum([b[i]*A[i][j]*c[j]*(2*c[i]-c[j]) for i in range(s) for j in range(s)]) == sympy.Rational(1,4))
🌟 Consistance
True
True
True
True
True
True
🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:
True
🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:
True
True
🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:
True
True
True
True
🌟 Pour être d'ordre 5 il faut que les conditions suivantes soient vraies:
False
True
Exemple à 6 étages explicite
Un schéma RK à \(6\) étages explicite a pour matrice de Butcher \(
\begin{array}{c|cccccc}
0 & 0 & 0 & 0 & 0 &0&0 \\
c_2 & a_{21} & 0&0&0&0&0\\
c_3 & a_{31} &a_{32}&0&0&0&0\\
c_4 & a_{41} &a_{42}&a_{43}&0&0&0\\
c_5 & a_{51} &a_{52}&a_{53}&a_{54}&0&0\\
c_6 & a_{61} &a_{62}&a_{63}&a_{64}&a_{65}&0\\
\hline
& b_1 & b_2 & b_3 & b_4 & b_5&b_6
\end{array}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+hc_2,u_n+h(a_{21}K_1)\right)\\
K_3 = \varphi\left(t_n+hc_3,u_n+h(a_{31}K_1+ a_{32}K_2\right)\\
K_4 = \varphi\left(t_n+hc_4,u_n+h(a_{41}K_1+ a_{42}K_2+ a_{43}K_3)\right)\\
K_5 = \varphi\left(t_n+hc_5,u_n+h(a_{51}K_1+ a_{52}K_2+ a_{53}K_3+ a_{54}K_4)\right)\\
K_6 = \varphi\left(t_n+hc_6,u_n+h(a_{61}K_1+ a_{62}K_2+ a_{63}K_3+ a_{64}K_4+ a_{65}K_5)\right)\\
u_{n+1} = u_n + h\left(b_1K_1+b_2K_2+b_3K_3+b_4K_4+b_5K_5+b_6K_6\right) & n=0,1,\dots N-1
\end{cases}
\) avec \(
\begin{cases}
c_2=a_{21}\\
c_3=a_{31}+a_{32}\\
c_4=a_{41}+a_{42}+a_{43}\\
c_5=a_{51}+a_{52}+a_{53}+a_{54}\\
c_6=a_{61}+a_{62}+a_{63}+a_{64}+a_{65}\\
b_1+b_2+b_3+b_4+b_5+b_6=1
\end{cases}
\)
Voici un exemple:
- schéma de Butcher (d’ordre 5) \(
\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}
\quad
\implies
\begin{cases}
u_0 = y_0 \\
K_1 = \varphi\left(t_n,u_n\right)\\
K_2 = \varphi\left(t_n+\frac{h}{4},u_n+\frac{h}{4}K_1\right)\\
K_3 = \varphi\left(t_n+\frac{h}{4},u_n+\frac{h}{8}(K_1+K_2)\right)\\
K_4 = \varphi\left(t_n+\frac{h}{2},u_n+\frac{h}{2}(-K_2+2K_3)\right)\\
K_5 = \varphi\left(t_n+\frac{3h}{4},u_n+\frac{h}{16}(3K_1+9K_4)\right)\\
K_6 = \varphi\left(t_{n+1},u_n+\frac{h}{7}(-3K_1+2K_2+12K_3-12K_4+8K_5)\right)\\
u_{n+1} = u_n + \frac{h}{90}\left(7K_1+32K_3+12K_4+32K_5+7K_6\right) & n=0,1,\dots N-1
\end{cases}
\)
Code
c = [0, sympy.Rational(1,4), sympy.Rational(1,4), sympy.Rational(1,2), sympy.Rational(3,4), 1]
b = [sympy.Rational(7,90), 0, sympy.Rational(32,90), sympy.Rational(12,90), sympy.Rational(32,90), sympy.Rational(7,90)]
A = [[0, 0, 0, 0, 0, 0],
[sympy.Rational(1,4), 0, 0, 0, 0, 0],
[sympy.Rational(1,8), sympy.Rational(1,8), 0, 0, 0, 0],
[0, sympy.Rational(-1,2), 1, 0, 0, 0],
[sympy.Rational(3,16), 0, 0, sympy.Rational(9,16), 0, 0],
[sympy.Rational(-3,7), sympy.Rational(2,7), sympy.Rational(12,7), sympy.Rational(-12,7), sympy.Rational(8,7), 0]]
s = len(c)
print("🌟 Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i])==c[i])
print("\n🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:")
print(sum([b[i]*c[i] for i in range(s)]).simplify()==sympy.Rational(1,2))
print("\n🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==sympy.Rational(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,6))
print("\n🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**3 for i in range(s)]).simplify()==sympy.Rational(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).simplify()==sympy.Rational(1,12))
print(sum([b[i]*A[i][j]*A[j][k]*c[k] for i in range(s) for j in range(s) for k in range(s)]).simplify()==sympy.Rational(1,24))
print("\n🌟 Pour être d'ordre 5 il faut que les conditions suivantes soient vraies:")
print(sum([b[i]*c[i]**4 for i in range(s)]) == sympy.Rational(1,5))
print(sum([b[i]*c[i]**2*(1-c[i]) for i in range(s)]) + sum([b[i]*A[i][j]*c[j]*(2*c[i]-c[j]) for i in range(s) for j in range(s)]) == sympy.Rational(1,4))
🌟 Consistance
True
True
True
True
True
True
True
🌟 Pour être d'ordre 2 il faut que la condition suivante soit vraie:
True
🌟 Pour être d'ordre 3 il faut que les conditions suivantes soient vraies:
True
True
🌟 Pour être d'ordre 4 il faut que les conditions suivantes soient vraies:
True
True
True
True
🌟 Pour être d'ordre 5 il faut que les conditions suivantes soient vraies:
True
True
Méthodes adapatatives
La manière la plus simple pour appliquer une méthode de résolution numérique consiste à utiliser un pas \(h\) constant. La principale difficulté est alors de choisir a priori la valeur optimale de ce pas \(h\) : il faut qu’il soit suffisamment petit pour garantir une précision suffisante, mais pas trop petit pour éviter de surcharger inutilement le calcul. Par ailleurs, la solution que l’on cherche à approcher peut nécessiter un pas \(h\) petit à certains endroits et plus grand à d’autres. Ainsi, si l’on accepte l’idée d’utiliser une discrétisation non uniforme, on peut décider de construire une méthode qui choisit à chaque itération \(n\) la taille du pas d’intégration \(h\), potentiellement différents à chaque itération. Pour ce faire, on a besoin d’estimer l’erreur \(\varepsilon_{n}\) que l’on commet lorsque l’on effectue une itération du schéma, afin de le comparer à une erreur « acceptable » “tol” : si \(\varepsilon_n\le \text{tol}\), on pourra augmenter la taille du pas d’intégration (si possible), et si au contraire \(\varepsilon_n> \text{tol}\), la dernière étape effectuée n’est pas acceptable, la taille du pas est diminuée et on répète les deux calculs.
Ce procédé, appelé adaptation du pas de discrétisation, est efficace mais nécessite d’avoir un bon estimateur d’erreur locale. En général, il s’agit d’un estimateur d’erreur a posteriori, car les estimateurs d’erreur a priori sont trop compliqués en pratique. On peut construire l’estimateur d’erreur de deux manières :
- en utilisant un unique schéma, mais avec deux pas de discrétisation différents ;
- en utilisant deux schémas d’ordre différent.
Choix et ajustement du pas : dans la pratique, on impose un encadrement du pas de temps \(h_{\min} \le h_n \le h_{\max}\), où \(h_{\min}\) est fixé par les contraintes de temps de calcul et par la croissance des erreurs d’arrondi lorsque le pas devient trop petit, et \(h_{\max}\) par des considérations de stabilité et de précision globale. À chaque itération, l’erreur locale est estimée par \(
\varepsilon_{n+1} = \frac{\left| u_{n+1} - \hat u_{n+1} \right|}{h_n}.
\) Étant donnée une tolérance prescrite \(\text{tol} >0\), le pas est ajusté selon les règles suivantes : - si \(\varepsilon_{n+1} > \text{tol}\), le pas est jugé trop grand et l’étape est rejetée ; on diminue alors le pas, par exemple en posant \(
h_{n+1} = 0.8 h_n ;
\)
- si \(\varepsilon_{n+1} \le \text{tol}\), l’étape est acceptée ; si de plus \(\varepsilon_{n+1}\) est nettement plus petite que \(\text{tol}\), on peut augmenter le pas pour la suite, par exemple en posant \(
h_{n+1} = 1.25 h_n ;
\)
- sinon, le pas est conservé inchangé.
CAS 1 : un unique schéma, mais avec deux pas de discrétisation différents
Une méthode courante pour garantir la précision dans la résolution d’un problème de Cauchy est de résoudre le problème deux fois en utilisant des pas \(h\) et \(h/2\), puis de comparer les approximations obtenues aux points de maillage correspondant à la grille la plus grossière. Autrement dit, à chaque étape, deux approximations différentes de la solution sont calculées
- on calcule \(u_{n+1}\) avec un pas \(h\) donné,
- on calcule également \(\hat u_{n+1}\) avec un pas \(h/2\) de la même méthode (ce qui nécessite deux étapes pour que cela corresponde à l’approximation au même point \(t_n+h\));
L’erreur est donc estimée par la différence entre ces deux approximations : \(
\varepsilon_{n+1} = \frac{\left| u_{n+1} - \hat u_{n+1} \right|}{h}.
\)
Selon la valeur de \(\varepsilon_{n+1}\), on peut alors décider d’augmenter \(h\) ou de le diminuer pour la prochaine itération.
Exemple : méthode d’Euler explicite à pas adaptatif
La méthode d’Euler explicite est bien adaptée à un calcul dynamique du pas de discrétisation \(h\) tenant compte des variations de l’inconnue sur l’intervalle d’intégration. Ici, l’estimateur d’erreur peut être construit à l’aide de deux pas de discrétisation (typiquement \(h\) et \(h/2\)).
Voyons un exemple d’implémentation de la méthode d’Euler explicite à pas adaptatif.
Code
import numpy as np
import matplotlib.pyplot as plt
def EE( phi, tt, y0 ):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for n in range(len(tt)-1):
uu[n+1] = uu[n]+h*phi(tt[n],uu[n])
return uu
def Adaptatif_EE(phi, t0, tfinal, y0, tol, hinit, hmin, hmax):
tt = [t0]
uu = [y0]
h = hinit
while tt[-1] < tfinal:
t = tt[-1]
u = uu[-1]
# Un pas d'Euler (h)
u_gr = u + h * phi(t, u)
# Deux pas d'Euler (h/2)
v = u + (h/2) * phi(t, u)
u_fi = v + (h/2) * phi(t + h/2, v)
# Estimation de l'erreur
err = abs(u_gr - u_fi) / h
if err > tol:
# Solution refusée et adaptation du pas
h = max(0.8 * h, hmin)
else:
# Solution acceptée
tt.append(t + h)
uu.append(u_fi)
# Adaptation du pas
if err < tol / 3:
h = min(1.25 * h, hmax)
# sinon h inchangé
return np.array(tt), np.array(uu)
# ================================
# MAIN
# ================================
t0, tfinal = 0, 2
y0 = 0
phi = lambda t,y : 1+y
sol_exacte = lambda t : np.exp(t)-1
# ================================
# Comparaison des méthodes
# ================================
plt.figure(figsize=(15,5))
# exacte
tt_ex = np.linspace(t0,tfinal,101)
yy = sol_exacte(tt_ex)
plt.plot(tt_ex, yy, 'r',label=r'Exacte')
plt.xlabel('t')
plt.ylabel('y(t)')
# Pas fixe h
tt = np.linspace(t0,tfinal,11)
h_init = tt[1]-tt[0]
uu_EE = EE(phi, tt, y0)
plt.plot(tt,uu_EE, 'x', label='EE avec h')
# Pas fixe h/2
tt2 = np.linspace(t0,tfinal,21)
uu_EE2 = EE(phi, tt2, y0)
plt.plot(tt2,uu_EE2, 'd', label='EE avec h/2')
# Pas adaptatif
tol = h_init/2
tt_adapt, uu_adapt = Adaptatif_EE(phi, t0, tfinal, y0, tol=tol, hinit=h_init, hmin=h_init/10, hmax=4*h_init)
plt.plot(tt_adapt,uu_adapt, 'o', label='Adaptatif_EE')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.grid();
# ================================
# Graphe de comparaison des erreurs
# ================================
err_EE = np.abs(uu_EE - sol_exacte(tt))
err_EE2 = np.abs(uu_EE2 - sol_exacte(tt2))
err_adapt = np.abs(uu_adapt - sol_exacte(tt_adapt))
plt.figure(figsize=(15,5))
plt.plot(tt, err_EE, 'x-', label='Erreur EE (h)')
plt.plot(tt2, err_EE2, 'd-', label='Erreur EE (h/2)')
plt.plot(tt_adapt, err_adapt, 'o-', label='Erreur Adaptatif')
plt.xlabel("t_n")
plt.ylabel(r"$|u(t_n)-y(t_n)|$")
plt.title("Comparaison des erreurs en fonction du temps")
plt.grid(True, which="both")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()
# ================================
# Comparaison des grilles
# ================================
plt.figure(figsize=(15,5))
plt.plot(tt,np.ones_like(tt)*3,'x',label='Grille grossière')
plt.plot(tt2,np.ones_like(tt2)*2,'d',label='Grille fine')
plt.plot(tt_adapt,np.ones_like(tt_adapt),'o',label='Grille adaptative')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.ylim(0,4)
plt.xlabel('t')
plt.xticks(tt2)
plt.yticks([1,2,3], ['Adaptative','Fine','Grossière'])
plt.grid()
plt.show()
# ================================================================
# Évolution du pas adaptatif
# ================================================================
plt.figure(figsize=(11,5))
plt.plot(tt_adapt[1:], np.diff(tt_adapt), 'o-')
plt.xlabel("$t$")
plt.ylabel("$h$")
plt.title("Évolution du pas adaptatif")
plt.grid(True)
plt.show()
# ================================================================
# Comparaison des approximations en t=tfinal
# ================================================================
from tabulate import tabulate
tab = []
headers = ["Méthode", "t final", "|u(t final) - sol_exacte(tfinal)|", "Points de grille"]
tab.append(["EE avec h", tt[-1], abs(uu_EE[-1]-sol_exacte(tt[-1])), len(tt)])
tab.append(["EE avec h/2", tt2[-1], abs(uu_EE2[-1]-sol_exacte(tt2[-1])), len(tt2)])
idx = np.argmax(tt_adapt>tfinal)-1
tab.append(["Adaptatif_EE", tt_adapt[idx], abs(uu_adapt[idx]-sol_exacte(tt_adapt[idx])), len(tt_adapt)])
# tab.append(["Exacte", tfinal, sol_exacte(tfinal), ""])
display(tabulate(tab, headers=headers, tablefmt="html", floatfmt=".2e"))
| EE avec h |
2.00e+00 |
1.20e+00 |
11 |
| EE avec h/2 |
2.00e+00 |
6.62e-01 |
21 |
| Adaptatif_EE |
1.98e+00 |
4.65e-01 |
19 |
CAS 2 : deux schémas d’ordre différent (méthode de Runge-Kutta emboîtée)
Tout comme la méthode d’Euler explicite, les méthodes RK étant à un pas se prêtent bien aux techniques d’adaptation. On peut construire leur estimateur d’erreur de deux manières :
- en utilisant un schéma RK du même ordre, mais avec deux pas de discrétisation différents (comme pour la méthode d’Euler ci-dessus) ;
- en utilisant deux schémas RK d’ordre différent, mais avec le même nombre \(s\) d’étapes et les mêmes coefficients \(a_{ij}\) et \(c_i\) (seules les coefficients \(b_i\) changent entre les deux méthodes).
La méthode de Runge-Kutta emboîtée utilise cette deuxième approche. L’intérêt de ces deux formules emboîtées est de permettre une estimation de l’erreur locale facilement calculable : \(
\varepsilon_{n+1} = \left| u_{n+1} - \hat u_{n+1} \right|
= h \left| \sum_{i=1}^s (b_i - \hat b_i) K_i \right|
\)
Méthode de Runge-Kutta adaptative Euler Explicite/Heun (RK12)
La méthode de Runge-Kutta adaptative la plus simple implique de combiner la méthode de Heun, d’ordre 2, avec la méthode d’Euler explicite, d’ordre 1. Son tableau de Butcher étendu est le suivant : \(
\begin{array}{c|cc}
0 & 0 & 0 \\
1 & 1 & 0 \\
\hline
\text{(Heun)}& \frac{1}{2} & \frac{1}{2}\\
\text{(EE)}& 1 & 0
\end{array}
\)
Code
import numpy as np
import matplotlib.pyplot as plt
def EE( phi, tt, y0 ):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for n in range(len(tt)-1):
k0 = phi(tt[n],uu[n])
k1 = phi(tt[n]+h,uu[n]+h*k0)
uu[n+1] = uu[n]+h*(1*k0+0*k1)
return uu
def HEUN( phi, tt, y0 ):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for n in range(len(tt)-1):
k0 = phi(tt[n],uu[n])
k1 = phi(tt[n]+h,uu[n]+h*k0)
uu[n+1] = uu[n]+h*(1/2*k0+1/2*k1)
return uu
def RK12( phi, t0, tfinal, y0, tol, hinit, hmin, hmax ):
tt = np.array([t0])
uu = np.array([y0])
uu[0] = y0
h = hinit
n = 0
while tt[n]<tfinal:
k0 = phi(tt[n],uu[n])
k1 = phi(tt[n]+h,uu[n]+h*k0)
# Un pas de EE
# uu_om = uu[n]+h*(1*k0+0*k1)
# Un pas de HEUN
uu_om_p1 = uu[n]+h*(1/2*k0+1/2*k1)
# choix
err = (-1/2*k0+1/2*k1)*h
if err > tol:
# Solution refusée et adaptation du pas
h = max( 0.8*h, hmin )
else:
# Solution acceptée
tt = np.append(tt,tt[n]+h)
uu = np.append(uu,uu_om_p1)
n = n+1
# Adaptation du pas
if err < tol / 3:
h = min( 1.25*h, hmax )
# sinon h inchangé
if tt[n]+h>tfinal:
h = tfinal-tt[n]
return tt,uu
# ================================
# MAIN
# ================================
t0, tfinal = 0, np.pi/2-0.1
y0 = 0
phi = lambda t,y : 1+y**2
sol_exacte = lambda t : np.tan(t)
# ================================
# Comparaison des méthodes
# ================================
plt.figure(figsize=(15,5))
# exacte
tt_ex = np.linspace(t0,tfinal,101)
yy = sol_exacte(tt_ex)
plt.plot(tt_ex, yy, 'r',label=r'$\tan(t)$')
# Pas fixe h
tt = np.linspace(t0,tfinal,11)
h_init = tt[1]-tt[0]
# EE
uu_EE = EE(phi, tt, y0)
plt.plot(tt,uu_EE, 'x', label='EE')
# HEUN
uu_HEUN = HEUN(phi, tt, y0)
plt.plot(tt,uu_HEUN, 'd', label='HEUN')
# Pas adaptatif
toll = 1.e-1
h_min = h_init/10
h_max = 4*h_init
tt_adapt, uu_adapt = RK12(phi, t0, tfinal, y0, tol=toll, hinit=h_init, hmin=h_min, hmax=h_max)
plt.plot(tt_adapt,uu_adapt, 'o', label='Adaptatif_EE')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.xlabel('$t$')
plt.ylabel('$y(t)$')
plt.grid();
# ================================
# Graphe de comparaison des erreurs
# ================================
err_EE = np.abs(uu_EE - sol_exacte(tt))
err_HEUN = np.abs(uu_HEUN - sol_exacte(tt))
err_adapt = np.abs(uu_adapt - sol_exacte(tt_adapt))
plt.figure(figsize=(15,5))
plt.plot(tt, err_EE, 'x-', label='Erreur EE')
plt.plot(tt, err_HEUN, 'd-', label='Erreur HEUN')
plt.plot(tt_adapt, err_adapt, 'o-', label='Erreur Adaptatif')
plt.xlabel("$t_n$")
plt.ylabel(r"$|u(t_n)-y(t_n)|$")
plt.title("Comparaison des erreurs en fonction du temps")
plt.grid(True, which="both")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()
# ================================
# Comparaison des grilles
# ================================
plt.figure(figsize=(15,5))
plt.plot(tt,np.ones_like(tt)*2,'x',label='Grille fixe')
plt.plot(tt_adapt,np.ones_like(tt_adapt),'o',label='Grille adaptative')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.ylim(0,3)
plt.xticks(tt)
plt.xlabel('$t$')
plt.yticks([1,2], ['Adaptative','Fixe'])
plt.grid()
plt.show()
# ================================================================
# Évolution du pas adaptatif
# ================================================================
plt.figure(figsize=(11,5))
plt.plot(tt_adapt[1:], np.diff(tt_adapt), 'o-')
plt.xlabel("$t$")
plt.ylabel("$h$")
plt.title("Évolution du pas adaptatif")
plt.grid(True)
plt.show()
# ================================================================
# Comparaison des approximations en t=tfinal
# ================================================================
from tabulate import tabulate
tab = []
headers = ["Méthode", "t final", "|u(t final)-exacte(tfinal)|", "Points de grille"]
tab.append(["EE", tt[-1], abs(uu_EE[-1]-sol_exacte(tfinal)), len(tt)])
tab.append(["HEUN", tt[-1], abs(uu_HEUN[-1]-sol_exacte(tfinal)), len(tt)])
idx = np.argmax(tt_adapt>tfinal)-1
tab.append(["RK12", tt_adapt[idx], abs(uu_adapt[idx]-sol_exacte(tfinal)), len(tt_adapt)])
display(tabulate(tab, headers=headers, tablefmt="html", floatfmt=".2e"))
| EE |
1.47e+00 |
6.38e+00 |
11 |
| HEUN |
1.47e+00 |
2.50e+00 |
11 |
| RK12 |
1.47e+00 |
9.59e-01 |
17 |
Méthode de Runge-Kutta adaptative Heun/Simpson (RK23)
Une méthode de Runge-Kutta adaptative simple combine la méthode de Heun, d’ordre 2, avec la méthode de Simpso, d’ordre 3. Son tableau de Butcher étendu est le suivant : \(
\begin{array}{c|cc}
0 & 0 & 0 & 0\\
1 & 1 & 0 & 0 \\
\frac{1}{2} & \frac{1}{4} & \frac{1}{4} & 0\\
\hline
\text{(Heun)}& \frac{1}{2} & \frac{1}{2} & 0\\
\text{(Simpson)}& \frac{1}{6} & \frac{1}{6} & \frac{4}{6}
\end{array}
\)
Code
import numpy as np
import matplotlib.pyplot as plt
def coeff_K(phi, h, t, u):
k0 = phi(t,u)
k1 = phi(t+h,u+h*k0)
k2 = phi(t+h/2,u+h/4*k0+h/4*k1)
return k0, k1, k2
def HEUN( phi, tt, y0 ):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for n in range(len(tt)-1):
k0, k1, k2 = coeff_K(phi, h, tt[n], uu[n])
uu[n+1] = uu[n]+h*(1/2*k0+1/2*k1)
return uu
def SIMPSON( phi, tt, y0 ):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for n in range(len(tt)-1):
k0, k1, k2 = coeff_K(phi, h, tt[n], uu[n])
uu[n+1] = uu[n]+h*(1/6*k0+4/6*k2+1/6*k1)
return uu
def RK23( phi, t0, tfinal, y0, tol, hinit, hmin, hmax ):
tt = np.array([t0])
uu = np.array([y0])
uu[0] = y0
h = hinit
n = 0
while tt[n]<tfinal :
k0, k1, k2 = coeff_K(phi, h, tt[n], uu[n])
# Un pas de Heun
# uu_om = uu[n]+h*(1/2*k0+1/2*k1)
# Un pas de Simpson
uu_om_p1 = uu[n]+h*(1/6*k0+1/6*k1+4/6*k2)
# choix
err = h * ( (1/6-1/2)*k0 + (1/6-1/2)*k1 + (4/6-0)*k2 )
if err > tol:
# Solution refusée et adaptation du pas
h = max( 0.8*h, hmin )
else:
# Solution acceptée
tt = np.append(tt,tt[n]+h)
uu = np.append(uu,uu_om_p1)
n = n+1
# Adaptation du pas
if err < tol / 3:
h = min( 1.25*h, hmax )
# sinon h inchangé
if tt[n]+h>tfinal:
h = tfinal-tt[n]
return tt,uu
# ================================
# MAIN
# ================================
t0, tfinal = 0, np.pi/2-0.1
y0 = 0
phi = lambda t,y : 1+y**2
sol_exacte = lambda t : np.tan(t)
# ================================
# Comparaison des méthodes
# ================================
plt.figure(figsize=(15,5))
# exacte
tt_ex = np.linspace(t0,tfinal,101)
yy = sol_exacte(tt_ex)
plt.plot(tt_ex, yy, 'r',label=r'$\tan(t)$')
# Pas fixe h
tt = np.linspace(t0,tfinal,11)
h_init = tt[1]-tt[0]
# HEUN
uu_HEUN = HEUN(phi, tt, y0)
plt.plot(tt,uu_HEUN, 'x', label='HEUN')
# SIMPSON
uu_SIMPSON = SIMPSON(phi, tt, y0)
plt.plot(tt,uu_SIMPSON, 'd', label='SIMPSON')
# Pas adaptatif
tt_adapt, uu_adapt = RK23(phi, t0, tfinal, y0, tol=1.e-1, hinit=h_init, hmin=1e-3, hmax=4*h_init)
plt.plot(tt_adapt,uu_adapt, 'o', label='Adaptatif')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.grid();
# ================================
# Graphe de comparaison des erreurs
# ================================
err_HEUN = np.abs(uu_HEUN - sol_exacte(tt))
err_SIMPSON = np.abs(uu_SIMPSON - sol_exacte(tt))
err_adapt = np.abs(uu_adapt - sol_exacte(tt_adapt))
plt.figure(figsize=(15,5))
plt.plot(tt, err_HEUN, 'x-', label='Erreur HEUN')
plt.plot(tt, err_SIMPSON, 'd-', label='Erreur SIMPSON')
plt.plot(tt_adapt, err_adapt, 'o-', label='Erreur Adaptatif')
plt.xlabel("$t_n$")
plt.ylabel(r"$|u(t_n)-y(t_n)|$")
plt.title("Comparaison des erreurs en fonction du temps")
plt.grid(True, which="both")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()
# ================================
# Comparaison des grilles
# ================================
plt.figure(figsize=(15,5))
plt.plot(tt,np.ones_like(tt)*2,'x',label='Grille fixe')
plt.plot(tt_adapt,np.ones_like(tt_adapt),'o',label='Grille adaptative')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.ylim(0,3);
plt.xticks(tt)
plt.grid()
plt.show()
# ================================================================
# Évolution du pas adaptatif
# ================================================================
plt.figure(figsize=(11,5))
plt.plot(tt_adapt[1:], np.diff(tt_adapt), 'o-')
plt.xlabel("$t$")
plt.ylabel("$h$")
plt.title("Évolution du pas adaptatif")
plt.grid(True)
plt.show()
# ================================================================
# Comparaison des approximations en t=tfinal
# ================================================================
from tabulate import tabulate
tab = []
headers = ["Méthode", "t final", "|u(t final)-exacte(tfinal)|", "Points de grille"]
tab.append(["HEUN", tt[-1], abs(uu_HEUN[-1]-sol_exacte(tfinal)), len(tt)])
tab.append(["SIMPSON", tt[-1], abs(uu_SIMPSON[-1]-sol_exacte(tfinal)), len(tt)])
idx = np.argmax(tt_adapt>tfinal)-1
tab.append(["RK23", tt_adapt[idx], abs(uu_adapt[idx]-sol_exacte(tt_adapt[idx])), len(tt_adapt)])
display(tabulate(tab, headers=headers, tablefmt="html", floatfmt=".2e"))
| HEUN |
1.47e+00 |
2.50e+00 |
11 |
| SIMPSON |
1.47e+00 |
1.03e+00 |
11 |
| RK23 |
1.47e+00 |
7.55e-01 |
10 |
Méthode de Runge-Kutta adaptative Bogacki/Shampine (ode23)
Une méthode de Runge-Kutta adaptative souvent utilisée combine la méthode de Bogacki avec celle de Shampine. Son tableau de Butcher étendu est le suivant : \(
\begin{array}{c|ccc}
0 & 0 & 0 & 0 & 0\\
\frac{1}{2} & \frac{1}{2} & 0 & 0 & 0\\
\frac{3}{4} & 0 & \frac{3}{4} & 0 & 0\\
1 & \frac{2}{9} & \frac{1}{3} & \frac{4}{9} & 0\\
\hline
& \frac{2}{9} & \frac{1}{3} & \frac{4}{9} & 0 \\
& \frac{7}{24} & \frac{1}{4} & \frac{1}{3} & \frac{1}{8}
\end{array}
\)
Code
import numpy as np
import matplotlib.pyplot as plt
def coeff_K(phi, h, t, u):
k0 = phi(t,u)
k1 = phi(t+h/2,u+h/2*k0)
k2 = phi(t+h*3/4,u+h*3/4*k1)
k3 = phi(t+h,u+h*2/9*k0+h/3*k1+h*4/9*k2)
return k0, k1, k2, k3
def Bogacki( phi, tt, y0 ):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for n in range(len(tt)-1):
k0, k1, k2, k3 = coeff_K(phi, h, tt[n], uu[n])
uu[n+1] = uu[n]+h*(2/9*k0+1/3*k1+4/9*k2)
return uu
def Shampine( phi, tt, y0 ):
h = tt[1]-tt[0]
uu = np.zeros_like(tt)
uu[0] = y0
for n in range(len(tt)-1):
k0, k1, k2, k3 = coeff_K(phi, h, tt[n], uu[n])
uu[n+1] = uu[n]+h*(7/24*k0+1/4*k1+1/3*k2+1/8*k3)
return uu
def edo23( phi, t0, tfinal, y0, tol, hinit, hmin, hmax ):
tt = np.array([t0])
uu = np.array([y0])
uu[0] = y0
h = hinit
n = 0
while tt[n]<tfinal :
k0, k1, k2, k3 = coeff_K(phi, h, tt[n], uu[n])
# Un pas de Bogacki
# uu_om = uu[n]+h*(2/9*k0+1/3*k1+4/9*k2)
# Un pas de Shamphine
uu_om_p1 = uu[n]+h*(7/24*k0+1/4*k1+1/3*k2+1/8*k3)
# choix
err = h * ( (7/24-2/9)*k0 + (1/4-1/3)*k1 + (1/3-4/9)*k2 + (1/8-0)*k3 )
if err > tol:
# Solution refusée et adaptation du pas
h = max( 0.7*h, hmin )
else:
# Solution acceptée
tt = np.append(tt,tt[n]+h)
uu = np.append(uu,uu_om_p1)
n = n+1
# Adaptation du pas
if err < tol/4:
h = min( 1.2*h, hmax )
# sinon h inchangé
if tt[n]+h>tfinal:
h = tfinal-tt[n]
return tt,uu
# ================================
# MAIN
# ================================
t0, tfinal = 0, np.pi/2-0.1
y0 = 0
phi = lambda t,y : 1+y**2
sol_exacte = lambda t : np.tan(t)
# ================================
# Comparaison des méthodes
# ================================
plt.figure(figsize=(15,5))
# exacte
tt_ex = np.linspace(t0,tfinal,21)
yy = sol_exacte(tt_ex)
plt.plot(tt_ex, yy, 'r',label=r'$\tan(t)$')
# Pas fixe h
tt = np.linspace(t0,tfinal,21)
h_init = tt[1]-tt[0]
# Bogacki
uu_Bogacki = Bogacki(phi, tt, y0)
plt.plot(tt,uu_Bogacki, 'x', label='Bogacki')
# Shampine
uu_Shampine = Shampine(phi, tt, y0)
plt.plot(tt,uu_Shampine, 'd', label='Shampine')
# Pas adaptatif
tt_adapt, uu_adapt = edo23(phi, t0, tfinal, y0, tol=1.e-3, hinit=h_init, hmin=1e-3, hmax=4*h_init)
plt.plot(tt_adapt,uu_adapt, 'o', label='Adaptatif')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.grid();
# ================================
# Graphe de comparaison des erreurs
# ================================
err_Bogacki = np.abs(uu_Bogacki - sol_exacte(tt))
err_Shampine = np.abs(uu_Shampine - sol_exacte(tt))
err_adapt = np.abs(uu_adapt - sol_exacte(tt_adapt))
plt.figure(figsize=(15,5))
plt.plot(tt, err_Bogacki, 'x-', label='Erreur Bogacki')
plt.plot(tt, err_Shampine, 'd-', label='Erreur Shampine')
plt.plot(tt_adapt, err_adapt, 'o-', label='Erreur Adaptatif')
plt.xlabel("$t_n$")
plt.ylabel(r"$|u(t_n)-y(t_n)|$")
plt.title("Comparaison des erreurs en fonction du temps")
plt.grid(True, which="both")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()
# ================================
# Comparaison des grilles
# ================================
plt.figure(figsize=(15,5))
plt.plot(tt,np.ones_like(tt)*2,'x',label='Grille fixe')
plt.plot(tt_adapt,np.ones_like(tt_adapt),'o',label='Grille adaptative')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.ylim(0,3);
plt.xticks(tt)
plt.grid()
plt.show()
# ================================================================
# Évolution du pas adaptatif
# ================================================================
plt.figure(figsize=(11,5))
plt.plot(tt_adapt[1:], np.diff(tt_adapt), 'o-')
plt.xlabel("$t$")
plt.ylabel("$h$")
plt.title("Évolution du pas adaptatif")
plt.grid(True)
plt.show()
# ================================================================
# Comparaison des approximations en t=tfinal
# ================================================================
from tabulate import tabulate
tab = []
headers = ["Méthode", "t final", "|u(t final)-exacte(tfinal)|", "Points de grille"]
tab.append(["Bogacki", tt[-1], abs(uu_Bogacki[-1]-sol_exacte(tfinal)), len(tt)])
tab.append(["Shampine", tt[-1], abs(uu_Shampine[-1]-sol_exacte(tfinal)), len(tt)])
idx = np.argmax(tt_adapt>tfinal)-1
tab.append(["edo23", tt_adapt[idx], abs(uu_adapt[idx]-sol_exacte(tt_adapt[idx])), len(tt_adapt)])
display(tabulate(tab, headers=headers, tablefmt="html", floatfmt=".2e"))
| Bogacki |
1.47e+00 |
2.74e-01 |
21 |
| Shampine |
1.47e+00 |
1.74e-01 |
21 |
| edo23 |
1.47e+00 |
1.72e-01 |
28 |
Runge-Kutta-Fehlberg (RKF45) à pas adaptatif
Voici la matrice de Butcher associée à la méthode de Runge-Kutta-Fehlberg (RKF45). Ce tableau représente les coefficients des deux méthodes Runge-Kutta utilisées dans la méthode RK45 pour calculer les coefficients \(k_0\) à \(k_5\), ainsi que les coefficients \(b_0\) à \(b_5\) et \(\hat b_0\) à \(\hat b_5\) utilisés pour combiner les \(K_i\) afin d’obtenir une approximation à l’étape suivante. \(
\begin{array}{c|cccccc}
0 & & & & & & \\
\frac{1}{4} & \frac{1}{4} & & & & & \\
\frac{3}{8} & \frac{3}{32} & \frac{9}{32} & & & & \\
\frac{12}{13} & \frac{1932}{2197} & -\frac{7200}{2197} & \frac{7296}{2197} & & & \\
1 & \frac{439}{216} & -\frac{8}{27} & \frac{3680}{513} & -\frac{845}{4104} & & \\
\frac{1}{2} & -\frac{8}{27} & \frac{2}{9} & -\frac{3544}{2565} & \frac{1859}{4104} & -\frac{11}{40} & \\
\hline
& \frac{16}{135} & 0 & \frac{6656}{12825} & \frac{28561}{56430} & -\frac{9}{50} & \frac{2}{55} \\
& \frac{25}{216} & 0 & \frac{1408}{2565} & \frac{2197}{4104} & -\frac{1}{5} & 0 \\
\end{array}
\)
Voyons un exemple d’implémentation de la méthode de Runge-Kutta-Fehlberg (RKF45) à pas adaptatif.
Code
import numpy as np
import matplotlib.pyplot as plt
def coeff_K(phi, h, t, u):
k = np.zeros(6)
k[0] = phi(t,u)
k[1] = phi(t+h/4,u+h/4*k[0])
k[2] = phi(t+3*h/8,u+3*h/32*k[0]+9*h/32*k[1])
k[3] = phi(t+12*h/13,u+1932*h/2197*k[0]-7200*h/2197*k[1]+7296*h/2197*k[2])
k[4] = phi(t+h,u+439*h/216*k[0]-8*h*k[1]+3680*h/513*k[2]-845*h/4104*k[3])
k[5] = phi(t+h/2,u-8*h/27*k[0]+2*h*k[1]-3544*h/2565*k[2]+1859*h/4104*k[3]-11*h/40*k[4])
return k
def RK4( phi, tt, y0 ):
h = tt[1]-tt[0]
uu = np.empty_like(tt)
uu[0] = y0
for n in range(len(tt)-1):
kk = coeff_K(phi, h, tt[n], uu[n])
uu[n+1] = uu[n] + h*( 25/216*kk[0] + 1408/2565*kk[2] + 2197/4104*kk[3] - 1/5*kk[4] )
return uu
def RK5( phi, tt, y0 ):
h = tt[1]-tt[0]
uu = np.empty_like(tt)
uu[0] = y0
for n in range(len(tt)-1):
kk = coeff_K(phi, h, tt[n], uu[n])
uu[n+1] = uu[n] + h* ( 16/135*kk[0] + 6656/12825*kk[2] + 28561/56430*kk[3] -9/50*kk[4] + 2/55*kk[5] )
return uu
def RK45( phi, t0, tfinal, y0, tol, hinit, hmin, hmax ):
tt = np.array([t0])
uu = np.array([y0])
uu[0] = y0
h = hinit
n = 0
while tt[n]<tfinal:
kk = coeff_K(phi, h, tt[n], uu[n])
# Un pas de RK5
uu_om_p1 = uu[n] + h* ( 16/135*kk[0] + 6656/12825*kk[2] + 28561/56430*kk[3] -9/50*kk[4] + 2/55*kk[5] )
# choix
err = h*( (25/216-16/135)*kk[0] +
(1408/2565-6656/12825)*kk[2] +
(2197/4104-28561/56430)*kk[3] +
(-1/5+9/50)*kk[4] +(0-2/55)*kk[5] )
if err > tol:
# Solution refusée et adaptation du pas
h = max( 0.8*h, hmin )
else:
# Solution acceptée
tt = np.append(tt,tt[n]+h)
uu = np.append(uu,uu_om_p1)
n = n+1
# Adaptation du pas
if err < tol / 3:
h = min( 1.25*h, hmax )
# sinon h inchangé
if tt[n]+h>tfinal:
h = tfinal-tt[n]
return tt,uu
# ================================
# MAIN
# ================================
t0, tfinal = 0, 4
y0 = 1
phi = lambda t,y : y*(2-y)
import sympy as sp
t = sp.symbols('t')
y = sp.Function('y')
sol_par = sp.dsolve( sp.Eq(y(t).diff(t), phi(t,y(t)) ) , ics={y(t0): y0} )
display(sol_par)
sol_exacte = sp.lambdify(t,sol_par.rhs,'numpy')
# ================================
# Comparaison des méthodes
# ================================
plt.figure(figsize=(15,5))
# exacte
tt_ex = np.linspace(t0,tfinal,101)
yy = sol_exacte(tt_ex)
plt.plot(tt_ex, yy, 'r',label=r'Exacte')
# Pas fixe h
tt = np.linspace(t0,tfinal,6)
h_init = tt[1]-tt[0]
# RK4
uu_RK4 = RK4(phi, tt, y0)
plt.plot(tt,uu_RK4, 'x', label='RK4')
# RK5
uu_RK5 = RK5(phi, tt, y0)
plt.plot(tt,uu_RK5, 'd', label='RK5')
# Pas adaptatif
tt_adapt, uu_adapt = RK45(phi, t0, tfinal, y0, tol=1.e-5, hinit=h_init, hmin=5e-3, hmax=4*h_init)
plt.plot(tt_adapt,uu_adapt, 'o', label='Adaptatif')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.grid();
# ================================
# Graphe de comparaison des erreurs
# ================================
err_RK4 = np.abs(uu_RK4 - sol_exacte(tt))
err_RK5 = np.abs(uu_RK5 - sol_exacte(tt))
err_adapt = np.abs(uu_adapt - sol_exacte(tt_adapt))
plt.figure(figsize=(15,5))
plt.plot(tt, err_RK4, 'x-', label='Erreur RK4')
plt.plot(tt, err_RK5, 'd-', label='Erreur RK5')
plt.plot(tt_adapt, err_adapt, 'o-', label='Erreur Adaptatif')
plt.xlabel("$t_n$")
plt.ylabel(r"$|u(t_n)-y(t_n)|$")
plt.title("Comparaison des erreurs en fonction du temps")
plt.grid(True, which="both")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()
# ================================
# Comparaison des grilles
# ================================
plt.figure(figsize=(15,5))
plt.plot(tt,np.ones_like(tt)*2,'x',label='Grille fixe')
plt.plot(tt_adapt,np.ones_like(tt_adapt),'o',label='Grille adaptative')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.ylim(0,3);
plt.xticks(tt)
plt.grid()
plt.show()
# ================================================================
# Évolution du pas adaptatif
# ================================================================
plt.figure(figsize=(11,5))
plt.plot(tt_adapt[1:], np.diff(tt_adapt), 'o-')
plt.xlabel("$t$")
plt.ylabel("$h$")
plt.title("Évolution du pas adaptatif")
plt.grid(True)
plt.show()
# ================================================================
# Comparaison des approximations en t=tfinal
# ================================================================
from tabulate import tabulate
tab = []
headers = ["Méthode", "t final", "|u(t final)-exacte(tfinal)|", "Points de grille"]
tab.append(["RK4", tt[-1], abs(uu_RK4[-1]-sol_exacte(tfinal)), len(tt)])
tab.append(["RK5", tt[-1], abs(uu_RK5[-1]-sol_exacte(tfinal)), len(tt)])
idx = np.argmax(tt_adapt>tfinal)-1
tab.append(["RK45", tt_adapt[idx], abs(uu_adapt[idx]-sol_exacte(tt_adapt[idx])), len(tt_adapt)])
display(tabulate(tab, headers=headers, tablefmt="html", floatfmt=".2e"))
\(\displaystyle y{\left(t \right)} = \frac{2}{1 + e^{- 2 t}}\)
| RK4 |
4.00e+00 |
3.36e-04 |
6 |
| RK5 |
4.00e+00 |
1.70e-04 |
6 |
| RK45 |
4.00e+00 |
1.31e-07 |
12 |
Dormand-Prince (DP45) à pas adaptatif
Voici la matrice de Butcher associée à la méthode de Dormand-Prince (DP45). Ce tableau représente les coefficients des deux méthodes Runge-Kutta utilisées dans la méthode RK45 pour calculer les coefficients \(k_0\) à \(k_5\), ainsi que les coefficients \(b_0\) à \(b_5\) et \(\hat b_0\) à \(\hat b_5\) utilisés pour combiner les \(K_i\) afin d’obtenir une approximation à l’étape suivante. \(
\begin{array}{c|cccccc}
0 & & & & & & \\
\frac{1}{5} & \frac{1}{5} & & & & & \\
\frac{3}{10} & \frac{3}{40} & \frac{9}{40} & & & & \\
\frac{4}{5} & \frac{44}{45} & -\frac{56}{15} & \frac{32}{9} & & & \\
\frac{8}{9} & \frac{19372}{6561} & -\frac{25360}{2187} & \frac{64448}{6561} & -\frac{212}{729} & & \\
1 & \frac{9017}{3168} & -\frac{355}{33} & \frac{46732}{5247} & \frac{49}{176} & -\frac{5103}{18656} & \\
1 & \frac{35}{384} & 0 & \frac{500}{1113} & \frac{125}{192} & -\frac{2187}{6784} & \frac{11}{84} \\
\hline
& \frac{35}{384} & 0 & \frac{500}{1113} & \frac{125}{192} & -\frac{2187}{6784} & \frac{11}{84} & 0\\
& \frac{5179}{57600} & 0 & \frac{7571}{16695} & \frac{393}{640} & -\frac{92097}{339200} & \frac{187}{2100} \\
\end{array}
\)
Voyons un exemple d’implémentation de la méthode de Runge-Kutta-Fehlberg (RKF45) à pas adaptatif.
Code
import numpy as np
import matplotlib.pyplot as plt
def coeff_K(phi, h, t, u):
k = np.zeros(7)
k[0] = phi(t,u)
k[1] = phi(t+h/5,u+h/5*k[0])
k[2] = phi(t+3*h/10,u+3*h/40*k[0]+9*h/40*k[1])
k[3] = phi(t+4*h/5,u+44*h/45*k[0]-56*h/15*k[1]+32*h/9*k[2])
k[4] = phi(t+8*h/9,u+19372*h/6561*k[0]-25360*h/2187*k[1]+64448*h/6561*k[2]-212*h/729*k[3])
k[5] = phi(t+h,u+9017*h/3168*k[0]-355*h/33*k[1]+46732*h/5247*k[2]+49*h/176*k[3]-5103*h/18656*k[4])
k[6] = phi(t+h,u+35*h/384*k[0]+0*k[1]+500*h/1113*k[2]+125*h/192*k[3]-2187*h/6784*k[4]+11*h/84*k[5])
return k
def Dor( phi, tt, y0 ):
h = tt[1]-tt[0]
uu = np.empty_like(tt)
uu[0] = y0
for n in range(len(tt)-1):
kk = coeff_K(phi, h, tt[n], uu[n])
uu[n+1] = uu[n] + h*( 35/384*kk[0] + 0*kk[1] + 500/1113*kk[2] + 125/192*kk[3] - 2187/6784*kk[4] + 11/84*kk[5] )
return uu
def Pri( phi, tt, y0 ):
h = tt[1]-tt[0]
uu = np.empty_like(tt)
uu[0] = y0
for n in range(len(tt)-1):
kk = coeff_K(phi, h, tt[n], uu[n])
uu[n+1] = uu[n] + h* ( 5179/57600*kk[0] + 0*kk[1] + 7571/16695*kk[2] + 393/640*kk[3] - 92097/339200*kk[4] + 187/2100*kk[5] + 1/40*kk[6] )
return uu
def DP45( phi, t0, tfinal, y0, tol, hinit, hmin, hmax ):
tt = np.array([t0])
uu = np.array([y0])
uu[0] = y0
h = hinit
n = 0
while tt[n]<tfinal:
kk = coeff_K(phi, h, tt[n], uu[n])
# Un pas de Pri
# uu_om_p1 = uu[n] + h* ( 5179/57600*k0 + 0*k1 + 7571/16695*k2 + 393/640*k3 - 92097/339200*k4 + 187/2100*k5 + 1/40*k6 )
uu_om_p1 = uu[n] + h*( 35/384*kk[0] + 0*kk[1] + 500/1113*kk[2] + 125/192*kk[3] - 2187/6784*kk[4] + 11/84*kk[5] )
# choix
err = h*( (35/384-5179/57600)*kk[0] +
(0-0)*kk[1] +
(500/1113-7571/16695)*kk[2] +
(125/192-393/640)*kk[3] +
(-2187/6784-92097/339200)*kk[4] +
(11/84-187/2100)*kk[5] +
(0-1/40)*kk[6] )
if err > tol:
# Solution refusée et adaptation du pas
h = max( 0.8*h, hmin )
else:
# Solution acceptée
tt = np.append(tt,tt[n]+h)
uu = np.append(uu,uu_om_p1)
n = n+1
# Adaptation du pas
if err < tol / 3:
h = min( 1.25*h, hmax )
# sinon h inchangé
if tt[n]+h>tfinal:
h = tfinal-tt[n]
return tt,uu
# ================================
# MAIN
# ================================
t0, tfinal = 0, 8
y0 = 1
phi = lambda t,y : y*(2-y)
import sympy as sp
t = sp.symbols('t')
y = sp.Function('y')
sol_par = sp.dsolve( sp.Eq(y(t).diff(t), phi(t,y(t)) ) , ics={y(t0): y0} )
display(sol_par)
sol_exacte = sp.lambdify(t,sol_par.rhs,'numpy')
# ================================
# Comparaison des méthodes
# ================================
plt.figure(figsize=(15,5))
# exacte
tt_ex = np.linspace(t0,tfinal,101)
yy = sol_exacte(tt_ex)
plt.plot(tt_ex, yy, 'r',label=r'Exacte')
# Pas fixe h
tt = np.linspace(t0,tfinal,11)
h_init = tt[1]-tt[0]
# Dorman
uu_Dor = Dor(phi, tt, y0)
plt.plot(tt,uu_Dor, 'v', label='Dorman')
# Prince
uu_Pri = Pri(phi, tt, y0)
plt.plot(tt,uu_Pri, 'd', label='Prince')
# DP45 à pas adaptatif
tt_adapt, uu_adapt = DP45(phi, t0, tfinal, y0, tol=1.e-3, hinit=h_init, hmin=h_init/10, hmax=3*h_init)
plt.plot(tt_adapt,uu_adapt, 'o', label='Adaptatif_EE')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.grid();
# ================================
# Graphe de comparaison des erreurs
# ================================
err_Dor = np.abs(uu_Dor - sol_exacte(tt))
err_Pri = np.abs(uu_Pri - sol_exacte(tt))
err_adapt = np.abs(uu_adapt - sol_exacte(tt_adapt))
plt.figure(figsize=(15,5))
plt.plot(tt, err_Dor, 'v-', label='Erreur Dorman')
plt.plot(tt, err_Pri, 'd-', label='Erreur Prince')
plt.plot(tt_adapt, err_adapt, 'o-', label='Erreur Adaptatif')
plt.xlabel("$t_n$")
plt.ylabel(r"$|u(t_n)-y(t_n)|$")
plt.title("Comparaison des erreurs en fonction du temps")
plt.grid(True, which="both")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()
# ================================
# Comparaison des grilles
# ================================
plt.figure(figsize=(15,5))
plt.plot(tt,np.ones_like(tt)*2,'v',label='Grille grossière')
plt.plot(tt_adapt,np.ones_like(tt_adapt),'o',label='Grille adaptative')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.ylim(0,3);
plt.xticks(tt)
plt.yticks([1,2], ['Adaptative','Fixe'])
plt.grid()
plt.show()
# ================================================================
# Évolution du pas adaptatif
# ================================================================
plt.figure(figsize=(11,5))
plt.plot(tt_adapt[1:], np.diff(tt_adapt), 'o-')
plt.xlabel("$t$")
plt.ylabel("$h$")
plt.title("Évolution du pas adaptatif")
plt.grid(True)
plt.show()
# ================================================================
# Comparaison des approximations en t=tfinal
# ================================================================
from tabulate import tabulate
tab = []
headers = ["Méthode", "t final", "|u(t final)-exacte(tfinal)|", "Points de grille"]
tab.append(["Dorman", tt[-1], abs(uu_Dor[-1]-sol_exacte(tfinal)), len(tt)])
tab.append(["Prince", tt[-1], abs(uu_Pri[-1]-sol_exacte(tfinal)), len(tt)])
idx = np.argmax(tt_adapt>tfinal)-1
tab.append(["DP45", tt_adapt[idx], abs(uu_adapt[idx]-sol_exacte(tt_adapt[idx])), len(tt_adapt)])
display(tabulate(tab, headers=headers, tablefmt="html", floatfmt=".2e"))
\(\displaystyle y{\left(t \right)} = \frac{2}{1 + e^{- 2 t}}\)
| Dorman |
8.00e+00 |
7.86e-08 |
11 |
| Prince |
8.00e+00 |
6.01e-08 |
11 |
| DP45 |
8.00e+00 |
6.76e-07 |
13 |
🔧 Méthode d’Euler explicite adaptative
Lorsqu’on utilise une des méthodes présentées jusqu’ici, une question ouverte est le choix du pas de discrétisation \(h\) pour obtenir une solution précise sans effectuer un nombre excessif de calculs. Une idée naturelle consiste à ajuster dynamiquement le pas de discrétisation en fonction de l’erreur commise à chaque itération. Si l’erreur est inférieure à une tolérance donnée, on peut augmenter le pas de discrétisation pour accélérer le calcul. Inversement, si l’erreur est trop grande, on peut diminuer le pas de discrétisation pour améliorer la précision. Dans ce cas, la discrétisation \(h_n\) n’est plus uniforme mais dépend de l’itération \(n\) : \(t_{n+1} = t_n + h_n\).
Pour estimer l’erreur de discrétisation locale, nous comparons deux approximations : la première calcule \(u_{n+1}\) avec un pas \(h\), la seconde calcule \(u_{n+1}\) avec un pas \(h/2\). L’erreur de discrétisation locale est alors estimée par la différence entre ces deux approximations.
Cette approche est particulièrement adaptée au cas des équations raides, où la solution varie rapidement sur certaines parties de l’intervalle de temps et lentement sur d’autres. Un pas de discrétisation uniforme peut être trop grand pour les parties lentes et trop petit pour les parties rapides. En ajustant le pas de discrétisation en fonction de l’erreur locale, on peut obtenir une solution précise sans effectuer un nombre excessif de calculs.
Étape I : Implémentation de la fonction euler_adaptatif
La méthode d’Euler explicite classique utilise un pas de temps fixe \(h\) : \(u_{n+1} = u_n + h \varphi(t_n, u_n).\) Cependant, pour améliorer la précision tout en minimisant le nombre de calculs, on peut utiliser une version adaptative de la méthode d’Euler. Dans cette version, le pas de temps \(h_n\) est ajusté dynamiquement en fonction d’une estimation de l’erreur locale.
Créez une fonction euler_adaptatif(phi, t0, y0, h0, tol, tfin) qui prend en entrée :
phi : fonction définissant l’équation différentielle \(y'(t) = \varphi(t, y(t))\)t0 : temps initialy0 : condition initiale \(y(t_0)\)h0 : pas de temps initialtol : tolérance pour l’estimation de l’erreur localetfin : temps maximal jusqu’où intégrer l’équation.
La fonction doit renvoyer :
- un tableau des temps \([t_0,t_1,\dots,t_N]\).
- un tableau des valeurs \([u_0,u_1,\dots,u_N]\).
L’algorithme d’adaptation du pas de temps suit les étapes suivantes :
Tant que \(t_n\le t_{\text{fin}}\) :
- Calcul des solutions avec deux pas différents :
- Utiliser la méthode d’Euler explicite avec un pas \(h\) : $ v_{n+1} = u_n + h (t_n, u_n). $
- Utiliser la méthode d’Euler explicite deux fois avec un pas \(h/2\) :
- Calcul intermédiaire : \(\tilde u = u_n + \frac{h}{2} \varphi(t_n, u_n)\)
- Calcul final : $ w_{n+1} = u + (t_n + h/2, u). $
- Ajustement du pas de temps :
- Définir l’erreur locale relative : $ = $
- Si l’erreur est inférieure à la tolérance \(\text{tol}\), ajouter \(t_{n+1}=t_n+h\) et \(u_{n+1}\) aux tableaux de résultats, puis doubler \(h\).
- Si l’erreur est supérieure à la tolérance \(\text{tol}\), diviser \(h\) par 2 (si \(h\) est inférieur à \(h_{\min}\), continuer avec \(h_{\min}\)) et refaire le calcul des solutions avec ce nouveau pas de temps.
Étape 2 : Test de la méthode
Considérons le problème de Cauchy
\(
\begin{cases}
y'(t) = -y(t),\\
y(0) = 1.
\end{cases}
\quad t \in [0, 20].
\)
Nous voulons comparer les solutions obtenues avec la méthode d’Euler explicite et avec la méthode d’Euler adaptative. Nous comparerons ensuite le nombre de points utilisés avec les deux méthodes.
- Utiliser Euler explicite avec un pas de temps fixe \(h = 1\) (valeur maximale qui garantit la monotonie de la solution).
- Utiliser Euler adaptatif avec \(h_0 = 1\) et une tolérance \(\text{tol} = 10^{-3}\).
- Comparer les solutions obtenues en traçant les 2 courbes (et la solution exacte).
- Comparer le nombre de pas effectués par chaque méthode et leurs erreurs respectives.
Correction
On commence par calculer une solution approchée avec la méthode d’Euler explicite. On peut ensuite calculer une solution approchée avec la méthode d’Euler adaptatif. On compare ensuite le nombre de points calculés avec les deux méthodes.
Code
import numpy as np
import matplotlib.pyplot as plt
def EE(phi,tt,y0):
N = len(tt)-1
h = tt[1]-tt[0]
uu = np.zeros(N+1)
uu[0] = y0
for n in range(N):
uu[n+1] = uu[n] + h*phi(tt[n],uu[n])
return uu
def EA(phi,t0,y0,hinit,tol,tfin):
hmin = 1.e-12
tt = np.array([t0])
uu = np.array([y0])
h = hinit
while tt[-1] < tfin:
k1 = phi(tt[-1],uu[-1])
v = uu[-1] + h*k1
u_demi = uu[-1] + h/2*k1
w = u_demi + h/2*phi(tt[-1]+h/2,u_demi)
erreur = np.abs(w - v)/np.linalg.norm(uu,ord=np.inf)
if erreur < tol or h <= hmin:
tt = np.append(tt,tt[-1]+h)
uu = np.append(uu,w)
h = 2*h
else :
h = max(h/2,hmin)
return tt,uu
Code
# TEST
# ========================
phi = lambda t, y : -y
t0, tfin = 0, 40
y0 = 1
# Pour la A-stabilité, il faut que h<2, pour la monotonie h<=1
# h = 1 pour N = (tfin-t0)/h = tfin-t0
N = int(tfin-t0)
tt = np.linspace(t0,tfin,N+1)
# Solution exacte
sol_exacte = lambda t : y0 * np.exp(-t)
plt.plot(tt,sol_exacte(tt),'-',label='exacte', color='green', lw=2)
# EE
uu_EE = EE(phi,tt,y0)
plt.plot(tt,uu_EE,'o-',label='EE', color='blue')
# EA
hinit = tt[1]-tt[0]
tol = 1.e-3
tt_EA, uu_EA = EA(phi,t0,y0,hinit,tol,tfin)
plt.plot(tt_EA,uu_EA,'o-',label='EA', color='red')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');
print(f"Erreur EE : {np.linalg.norm(uu_EE-sol_exacte(tt))}")
print(f"Erreur EA : {np.linalg.norm(uu_EA-sol_exacte(tt_EA))}")
print(f"Nombre de pas EE : {len(tt)-1}")
print(f"Nombre de pas EA : {len(tt_EA)-1}")
# Affichons juste les deux grilles de temps sur une droite horizontale
plt.figure()
plt.plot(tt,np.ones_like(tt),'o-',label='EE',color='blue')
plt.plot(tt_EA,1-np.ones_like(tt_EA),'o-',label='EA',color='red')
plt.grid()
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.title('Grilles de temps');
Erreur EE : 0.3956231069460752
Erreur EA : 0.03012966353642783
Nombre de pas EE : 40
Nombre de pas EA : 47
🔧 Méthode adaptative avec deux schémas distincts
Lorsqu’on utilise une des méthodes présentées jusqu’ici, une question ouverte est le choix du pas de discrétisation \(h\) pour obtenir une solution précise sans effectuer un nombre excessif de calculs. Une idée naturelle consiste à ajuster dynamiquement le pas de discrétisation en fonction de l’erreur commise à chaque itération. Si l’erreur est inférieure à une tolérance donnée, on peut augmenter le pas de discrétisation pour accélérer le calcul. Inversement, si l’erreur est trop grande, on peut diminuer le pas de discrétisation pour améliorer la précision. Dans ce cas, la discrétisation \(h_n\) n’est plus uniforme mais dépend de l’itération \(n\) : \(t_{n+1} = t_n + h_n\).
Pour estimer l’erreur de discrétisation locale, nous comparons deux approximations : pour un même pas \(h\), la première calcule \(u_{n+1}\) avec un schéma, la seconde calcule \(u_{n+1}\) avec un autre schéma. L’erreur de discrétisation locale est alors estimée par la différence entre ces deux approximations.
Cette approche est particulièrement adaptée au cas des équations raides, où la solution varie rapidement sur certaines parties de l’intervalle de temps et lentement sur d’autres. Un pas de discrétisation uniforme peut être trop grand pour les parties lentes et trop petit pour les parties rapides. En ajustant le pas de discrétisation en fonction de l’erreur locale, on peut obtenir une solution précise sans effectuer un nombre excessif de calculs.
Étape I : Implémentation de la fonction euler_heun_adaptatif
Créez une fonction euler_heun_adaptatif(phi, t0, y0, h0, tol, tfin) qui prend en entrée :
phi : fonction définissant l’équation différentielle \(y'(t) = \varphi(t, y(t))\)t0 : temps initialy0 : condition initiale \(y(t_0)\)h0 : pas de temps initialtol : tolérance pour l’estimation de l’erreur localetfin : temps maximal jusqu’où intégrer l’équation.
La fonction doit renvoyer :
- un tableau des temps \([t_0,t_1,\dots,t_N]\).
- un tableau des valeurs \([u_0,u_1,\dots,u_N]\).
L’algorithme d’adaptation du pas de temps suit les étapes suivantes :
Tant que \(t_n\ge t_{\text{fin}}\) :
- Calcul des solutions avec deux pas différents :
- Utiliser la méthode d’Euler explicite avec un pas \(h\) : $ v_{n+1} = u_n + h (t_n, u_n). $
- Utiliser la méthode de Heun avec le même pas \(h\) : $ w_{n+1} = u_n + ((t_n, u_n)+(t_{n+1}, v_{n+1})). $
- Ajustement du pas de temps :
- Définir l’erreur locale relative : $ = $
- Si l’erreur est inférieure à la tolérance \(\text{tol}\), ajouter \(t_{n+1}=t_n+h\) et \(u_{n+1}\) aux tableaux de résultats, puis doubler \(h\).
- Si l’erreur est supérieure à la tolérance \(\text{tol}\), diviser \(h\) par 2 (si \(h\) est inférieur à \(h_{\min}\), continuer avec \(h_{\min}\)) et refaire le calcul des solutions avec ce nouveau pas de temps.
Étape 2 : Test de la méthode
Considérons le problème de Cauchy
\(
\begin{cases}
y'(t) = -y(t),\\
y(0) = 1.
\end{cases}
\quad t \in [0, 20].
\)
Nous voulons comparer les solutions obtenues avec la méthode d’Euler explicite, la méthode de Heun et avec la méthode d’Euler-Heun adaptative. Nous comparerons ensuite le nombre de points utilisés avec les deux méthodes.
- Utiliser Euler explicite avec un pas de temps fixe \(h = 1\) (valeur maximale qui garantit la monotonie de la solution). Même exercice avec Heun.
- Utiliser Euler-Heun adaptatif avec \(h_0 = 1\) et une tolérance \(\text{tol} = 5\times 10^{-3}\).
- Comparer les solutions obtenues en traçant les 3 courbes (et la solution exacte).
- Comparer le nombre de pas effectués par chaque méthode et leurs erreurs respectives.
Correction
On commence par calculer une solution approchée avec la méthode d’Euler explicite. On peut ensuite calculer une solution approchée avec la méthode d’Euler adaptatif. On compare ensuite le nombre de points calculés avec les deux méthodes.
Code
import numpy as np
import matplotlib.pyplot as plt
def EE(phi,tt,y0):
N = len(tt)-1
h = tt[1]-tt[0]
uu = np.zeros(N+1)
uu[0] = y0
for n in range(N):
k_0 = phi(tt[n],uu[n])
uu[n+1] = uu[n] + h*k_0
return uu
def Heun(phi,tt,y0):
N = len(tt)-1
h = tt[1]-tt[0]
uu = np.zeros(N+1)
uu[0] = y0
for n in range(N):
k_0 = phi(tt[n],uu[n])
k_1 = phi(tt[n]+h,uu[n]+h*k_0)
uu[n+1] = uu[n] + h/2*(k_0+k_1)
return uu
def EHA(phi,t0,y0,hinit,tol,tfin):
hmin = 1.e-12
tt = np.array([t0])
uu = np.array([y0])
h = hinit
while tt[-1] < tfin:
k_0 = phi(tt[-1],uu[-1])
v = uu[-1] + h*k_0
k_1 = phi(tt[-1]+h,v)
w = uu[-1] + h/2*(k_0 + k_1)
erreur = np.abs(w - v)/np.linalg.norm(uu,ord=np.inf)
if erreur < tol or h <= hmin:
tt = np.append(tt,tt[-1]+h)
uu = np.append(uu,w)
h = 2*h
else :
h = max(h/2,hmin)
return tt,uu
Code
# TEST
# ========================
phi = lambda t, y : -y
t0, tfin = 0, 40
y0 = 1
# Pour la A-stabilité, il faut que h<2, pour la monotonie h<=1
# h = 1 pour N = (tfin-t0)/h = tfin-t0
N = int(tfin-t0)
tt = np.linspace(t0,tfin,N+1)
# Solution exacte
sol_exacte = lambda t : y0 * np.exp(-t)
plt.plot(tt,sol_exacte(tt),'-',label='exacte', color='green', lw=2)
# EE
uu_EE = EE(phi,tt,y0)
plt.plot(tt,uu_EE,'o-',label='EE', color='blue')
# Heun
uu_Heun = Heun(phi,tt,y0)
plt.plot(tt,uu_Heun,'o-',label='Heun', color='orange')
# EHA
hinit = tt[1]-tt[0]
tol = 5.e-3
tt_EHA, uu_EHA = EHA(phi,t0,y0,hinit,tol,tfin)
plt.plot(tt_EHA,uu_EHA,'o-',label='EA', color='red')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left');
# Affichons juste les deux grilles de temps sur une droite horizontale
plt.figure()
plt.plot(tt,np.ones_like(tt),'o-',label='EE',color='blue')
plt.plot(tt,0.5-np.zeros_like(tt),'o-',label='Heun',color='orange')
plt.plot(tt_EHA,1-np.ones_like(tt_EHA),'o-',label='EHA',color='red')
plt.grid()
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.title('Grilles de temps');
# print(f"Erreur EE : {np.linalg.norm(uu_EE-sol_exacte(tt))}")
# print(f"Erreur Heun : {np.linalg.norm(uu_Heun-sol_exacte(tt))}")
# print(f"Erreur EHA : {np.linalg.norm(uu_EHA-sol_exacte(tt_EHA))}")
# print(f"Nombre de pas EE : {len(tt)-1}")
# print(f"Nombre de pas Heun : {len(tt)-1}")
# print(f"Nombre de pas EHA : {len(tt_EHA)-1}")
# Affichage des résultats avec pandas
# ===================================
import pandas as pd
# Création du DataFrame contenant les résultats
data = {
"Méthode": ["Euler Explicite", "Heun", "Euler-Heun Adaptatif"],
"Erreur": [
np.linalg.norm(uu_EE - sol_exacte(tt)),
np.linalg.norm(uu_Heun - sol_exacte(tt)),
np.linalg.norm(uu_EHA - sol_exacte(tt_EHA))
],
"Nombre de pas": [
len(tt) - 1,
len(tt) - 1,
len(tt_EHA) - 1
]
}
pd.DataFrame(data)
| 0 |
Euler Explicite |
0.395623 |
40 |
| 1 |
Heun |
0.197615 |
40 |
| 2 |
Euler-Heun Adaptatif |
0.009670 |
47 |
