from IPython.display import display, Latex
from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
from math import *
Les schémas de Runge-Kutta sont des méthodes à un pas qui approchent l'intégrale $\int_{t_n}^{t_{n+1}} \varphi(t,y(t))\dt$ par une formule de quadrature en des points donnés toujours à l'intérieur de l'intervalle $[t_{n};t_{n+1}]$: \begin{align*} t_{n,i} & \stackrel{\text{déf}}{=}t_n+c_ih \qquad\text{avec }c_i\in[0;1] \\ \int_{t_n}^{t_{n+1}} \varphi(t,y(t))\dt &\approx h\sum_{i=1}^sb_i \varphi(t_{n,i},y(t_{n,i})). \end{align*} Le problème est que, si $t_{n,i}$ n'est pas un point de la discrétisation, l'on ne connait pas $y(t_{n,i})$. L'évaluation des $\varphi(t_{n,i},y(t_{n,i}))$ en les points intérieurs et alors approchée par 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). $$ Cela donne $$ y_{n+1} \approx y_n + h \sum_{i=1}^s\left[b_i \varphi\left(t_{n,i},y_n+h\sum_{j=1}^{s}a_{ij} \varphi(t_{n,j},y_{n,j})\right)\right] $$
Une méthode de Runge-Kutta à $s\ge1$ étages s'écrit: \begin{cases} u_0=y(t_0)=y_0,\\ u_{n+1}=u_{n}+h \displaystyle\sum_{i=1}^sb_iK_i& n=0,1,\dots N-1\\ K_i=\displaystyle\varphi\left( t_n+hc_i,u_n+h\sum_{j=1}^{s}a_{ij}K_j \right) & i=1,\dots s. \end{cases}
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ée matrice de Butcher de la méthode Runge-Kutta considérée.
Théorème [Ordre 1]
La méthode de Runge-Kutta est consistante (i.e. d'ordre 1) ssi
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)$.
Théorème [Ordre 2, 3 et 4]
La méthode de Runge-Kutta est
Théorème
L'ordre d'une méthode RK **explicite** à $s$ étapes ne peut pas être plus grand que $s$.
De plus, il n’existe pas de méthode à $s$ étapes d’ordre $s$ si $s \ge 5$.
L'ordre d'une méthode RK **implicite** à $s$ étapes ne peut pas être plus grand que $2s$.
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.
Rappelons la définition de schéma A-stable:
Soit $\beta>0$ un nombre réel positif et considérons le problème de Cauchy \begin{cases} y'(t)=-\beta y(t), &\text{pour }t>0,\\ y(0)=y_0 \end{cases} où $y_0\neq0$ est une valeur donnée. Sa solution est $y(t)=y_0e^{-\beta t}$ donc $$\lim_{t\to+\infty}y(t)=0.$$ Soit $h>0$ un pas de temps donné, $t_n=nh$ pour $n\in\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.
Une méthode de Runge-Kutta à $s\ge1$ étages pour $y'(t)=-\beta y(t)$ s'écrit: \begin{cases} u_0=y_0,\\ u_{n+1}=u_{n}+h \displaystyle\sum_{i=1}^sb_iK_i& n=0,1,\dots N-1\\ K_i=\displaystyle -\beta \left( u_n+h\sum_{j=1}^{s}a_{ij}K_j \right) & i=1,\dots s. \end{cases}
On peut monter que, contrairement aux méthodes multi-pas, la taille des régions de stabilité absolue des méthodes RK augmente quand l’ordre augmente.
Étudier les méthodes RK avec $s=1$:
Ces méthodes ont 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 \implies \begin{array}{c|cc} c_1 & c_1\\ \hline & 1 \end{array} \quad \implies \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} $$
Par exemple, si on choisi $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) $$
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 si on remplace le $u_n+\frac{h}{2} K_1=\frac{u_n+\left(u_n+h K_1\right)}{2}$ dans la définition de $K_1$ par $u_{n+1}$ on obtient $$ K_1 = \varphi\left(t_{n}+\frac{h}{2},u_n+\frac{h}{2} K_1\right) = \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(t_{n}+\frac{h}{2},\frac{u_n+u_{n+1}}{2}\right) $$
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.
Si le schéma est explicite alors il est d'ordre 1.
Une méthode de Runge-Kutta à 1 étage pour $y'(t)=-\beta y(t)$ 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:
Pour nos exemples, on obtient:
Étudier les méthodes RK avec $s=2$:
Les méthodes avec $s=2$ ont 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} $$
Voici un exemple, appelée 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.
from sympy import *
gamma=sqrt(3)/6
c=[1/2-gamma,1/2+gamma]
b=[1/2,1/2]
A=[[1/4,1/4-gamma],[1/4+gamma,1/4]]
s=len(c)
print("Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i]).simplify()==c[i])
print("\nOrdre 2")
print(sum([b[i]*c[i] for i in range(s)]).simplify()==1/2)
print("\nOrdre 3")
print(sum([b[i]*c[i]**2 for i in range(s)]).simplify()==1/3)
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==1/6)
print("\nOrdre 4")
print(sum([b[i]*c[i]**3 for i in range(s)]).simplify()==1/4)
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)]).simplify()==1/8)
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)]).simplify()==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()==1/24)
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} $$
Par exemple, si on choisi $c_1=0$, $c_1=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)) $$
Les méthodes avec $s=2$ explicites ont 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} $$
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 deux 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é $$ \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} $$
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$.
Un schéma RK à $3$ étages explicite a pour matrice de Butcher $$ \begin{array}{c|cc} 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 deux exemples:
Étant une méthode explicite à 3 étapes, elle a au mieux ordre 3. Vérifions qu'elle est effectivement d'ordre 3:
c=[0,1/3,2/3]
b=[1/4,0,3/4]
A=[[0,0,0],[1/3,0,0],[0,2/3,0]]
s=len(c)
print("Consistance")
print(sum(b)==1)
for i in range(s):
print(sum(A[i])==c[i])
print("\nOrdre 2")
print(sum([b[i]*c[i] for i in range(s)])==1/2)
print("\nOrdre 3")
print(sum([b[i]*c[i]**2 for i in range(s)])==1/3)
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==1/6)
Étant une méthode explicite à 3 étapes, elle a au mieux ordre 3. On vérifie ci-dessous qu'elle est d'ordre 3:
from fractions import Fraction
c=[0,Fraction(1,2),1]
b=[Fraction(1,6),Fraction(2,3),Fraction(1,6)]
A=[[0,0,0],[Fraction(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("\nOrdre 2")
print(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))
print("\nOrdre 3")
print(sum([b[i]*c[i]**2 for i in range(s)])==Fraction(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,6))
Étant une méthode explicite à 3 étapes, elle a au mieux ordre 3. On vérifie ci-dessous qu'elle n'est que d'ordre 2:
c=[0,Fraction(1,2),1]
b=[-Fraction(1,6),Fraction(4,3),-Fraction(1,6)]
A=[[0,0,0],[Fraction(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("\nOrdre 2")
print(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))
print("\nOrdre 3")
print(sum([b[i]*c[i]**2 for i in range(s)])==Fraction(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,6))
Un schéma RK à $4$ étages explicite a pour matrice de Butcher $$ \begin{array}{c|cc} 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):
c=[0,Fraction(1,2),Fraction(1,2),1]
b=[Fraction(1,6),Fraction(1,3),Fraction(1,3),Fraction(1,6)]
A=[[0,0,0,0],[Fraction(1,2),0,0,0],[0,Fraction(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("\nOrdre 2")
print(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))
print("\nOrdre 3")
print(sum([b[i]*c[i]**2 for i in range(s)])==Fraction(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,6))
print("\nOrdre 4")
print(sum([b[i]*c[i]**3 for i in range(s)])==Fraction(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)])==Fraction(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)])==Fraction(1,24))
c=[0,Fraction(1,4),Fraction(1,2),1]
b=[Fraction(1,6),0,Fraction(2,3),Fraction(1,6)]
A=[[0,0,0,0],[Fraction(1,4),0,0,0],[0,Fraction(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("\nOrdre 2")
print(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))
print("\nOrdre 3")
print(sum([b[i]*c[i]**2 for i in range(s)])==Fraction(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,6))
print("\nOrdre 4")
print(sum([b[i]*c[i]**3 for i in range(s)])==Fraction(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)])==Fraction(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)])==Fraction(1,24))
from fractions import Fraction
c=[0,Fraction(1,3),Fraction(2,3),1]
b=[Fraction(1,8),Fraction(3,8),Fraction(3,8),Fraction(1,8)]
A=[[0,0,0,0],[Fraction(1,3),0,0,0],[-Fraction(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("\nOrdre 2")
print(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))
print("\nOrdre 3")
print(sum([b[i]*c[i]**2 for i in range(s)])==Fraction(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,6))
print("\nOrdre 4")
print(sum([b[i]*c[i]**3 for i in range(s)])==Fraction(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)])==Fraction(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)])==Fraction(1,24))
Un schéma RK à $5$ étages explicite a pour matrice de Butcher $$ \begin{array}{c|cc} 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:
from fractions import Fraction
c=[0,Fraction(1,3),Fraction(1,3),Fraction(1,2),1]
b=[Fraction(1,6),0,0,Fraction(2,3),Fraction(1,6)]
A=[[0,0,0,0,0],[Fraction(1,3),0,0,0,0],[Fraction(1,6),Fraction(1,6),0,0,0],[Fraction(1,8),0,Fraction(3,8),0,0],[Fraction(1,2),0,-Fraction(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("\nOrdre 2")
print(sum([b[i]*c[i] for i in range(s)])==Fraction(1,2))
print("\nOrdre 3")
print(sum([b[i]*c[i]**2 for i in range(s)])==Fraction(1,3))
print(sum([b[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,6))
print("\nOrdre 4")
print(sum([b[i]*c[i]**3 for i in range(s)])==Fraction(1,4))
print(sum([b[i]*c[i]*A[i][j]*c[j] for i in range(s) for j in range(s)])==Fraction(1,8))
print(sum([b[i]*A[i][j]*c[j]**2 for i in range(s) for j in range(s)])==Fraction(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)])==Fraction(1,24))