Références: Benzoni-Gavage 2010 ou Demailly 2006 ou Agarwal and O'Regan 2008 pour un exposé complet; Boularas et al 2009 pour un exposé vivant et accessible.
Considérerons une équation différentielle du premier ordre:
\begin{align*} \dot x(t) = f(t,x(t))\,,\ t\geq 0\,,\ x(0)=x_{0} \tag{EDO} \end{align*}où $x(t)\in\mathbb{R}^d$ et $f:\mathbb{R}_+\times\mathbb{R}^d\to \mathbb{R}^d$. $x_0$ est la condition initiale. Lorsque $f$ ne dépend pas de $t$, l'ED est dite autonome, lorsque $f$ est linéaire ou affine en $x$, l'ED est dite linéaire.
L'existence et l'unicité d'une solution de (EDO sont assurées par le théorème de Cauchy-Lipschitz. (EDO) est un problème local en temps, nous nous contentons ici de donner quelques éléments dans le cas plus simple de solutions globales en supposant que $f$ est globalement lipschitzienne.
Hypothèse: La fonction $f$ est dite globalement lipschitzienne en $x$ uniformément en $t$, s'il existe une constante $L$ telle que:
\begin{align*} |f(t,x)-f(t,y)| \leq L\,|x-y|\,,\ \forall t,x,y\,. \end{align*}Une fonction lipschitzienne est continue mais pas nécessairement dérivable. Lorsque $f$ admet une dérivée uniformément bornée par une constante $C$, i.e. $|\partial f(t,x)/\partial x|\leq C$ pour tout $(t,x)$, alors elle est lipschitzienne avec constante $C$.
On dira que (EDO) admet une solution unique si, étant données deux fonctions $x_i(t)$ vérifiant (EDO), $i=1,2$, alors $x_1(t)=x_2(t)$ pour tout $t\geq 0$.
Proposition (existence et unicité): Supposons que $(t,x)\to f(t,x)$ soit continue et globalement lipschitzienne en $x$ uniformément en $t$, alors (EDO) admet une et une seule solution.
Lemme (inégalité de Gronwall): Soit $I=[a,\infty]$ ou $[a,b]$ ou $[a,b[$ avec $a<b$. Soit $\alpha$, $\beta$, $u$ des fonctions réelles définies sur $I$. On suppose que $\beta$ et $u$ sont continues et que la partie négative de $\alpha$ est intégrable sur tout sous-intervalle borné et fermé de $I$. Si $\beta\geq 0$ et $u$ vérifie: \begin{align*} u(t) \leq \alpha(t) + \int_a^t \beta(s) u(s)\,\mathrm{d}s,\qquad \forall t\in I \end{align*} alors \begin{align*} u(t) \le \alpha(t) + \int_a^t\alpha(s)\beta(s)\exp\biggl(\int_s^t\beta(r)\,\mathrm{d}r\biggr)\mathrm{d}s,\qquad t\in I. \end{align*} Si de plus $\alpha$ est non-décroissante alors: \begin{align*} u(t) \leq \alpha(t)\exp\biggl(\int_a^t\beta(s)\,\mathrm{d}s\biggr),\qquad t\in I. \end{align*}
L'unicité fait appel à l'inégalité de Gronwall. Elle permet de démontrer que si $x_{i}(\cdot)$, $i=1,2$, sont deux solutions de (EDO), alors le fait que $f$ soit globalement lipschitzienne implique que:
$$ |x_{1}(t)-x_{2}(t)| \leq e^{L\,|t|}\,|x_{1}(0)-x_{2}(0)|\,,\forall t\geq 0\,. $$Donc si $x_{1}(0)=x_{2}(0)=x_0$ alors $x_{1}(t)=x_{2}(t)$ pour tout $t$.
Il existe deux méthodes classiques pour exhiber des solutions du problème de Cauchy et donc montrer l'existence d'une solution. Nous verrons plus loin la méthode d'approximation d'Euler, une autre approche classique, puissante mais moins constructive, consiste à faire appel aux approximations successives de Picard. On se donne $x^{(0)}(t)$ continue, par exemple $x^{(0)}(t)=0$ pour tout $t$, et on pose:
$$ x^{(n+1)} = \mathcal{I} x^{(n)} $$où $\mathcal{I}$ est un opérateur qui à une fonction $x:\mathbb R_+\to \mathbb R^d$ asscocie une $\mathcal{I}x:\mathbb R_+\to \mathbb R^d$ définie par
$$ [\mathcal{I} x](t) = x_{0} + \int_{t_{0}}^t f(s,x(s))\,{\rm d} s\,,\ t\geq 0\,. $$La solution de (EDO) apparait alors comme un point fixe de l'opérateur $\mathcal{I}$, i.e. $x$ telle que $x=\mathcal{I} x$. La démonstration de l'existence d'une solution de (EDO) consiste alors à démontrer l'existence d'un point fixe de l'opérateur $\mathcal I$ en démontrant que cet opérateur est contractant. La preuve de ce résultat est donnée dans [Benzoni-Gavage 2010 p. 147].
On considère l'ED autonome:
$$ \dot x(t) = 2\,\sqrt{x(t)}\,,\ x(0)=0\,. $$Il est clair que $x_{1}(t)=0$, $t\geq 0$ est une solution de cette équation ainsi que $x_{2}(t)=t^2$, $t\geq 0$. Par une technique de recollement consistant à mettre bout à bout ces deux solutions, on peut cosntruire une infinité de solutions. En effet:
$$ x_{1}(t) = \left\{ \begin{array}{ll} 0\,,& 0\leq t\leq a \\ (t-a)^2\,, & t\geq a \end{array} \right. $$défini une solution différente pour tout $a\geq 0$. Ici $f(x)=2\,\sqrt{x}$ n'est pas globalement uniformément lipschitzienne à cause de son comportement en $O$.
In [20]:
import numpy as np
import matplotlib.pyplot as plt
a = 3
t0, t1 = 0, 10
le_temps = np.linspace(t0,t1,300, endpoint=True)
courbe = (le_temps-a)**2
courbe[le_temps-a<0] = 0
plt.plot(le_temps, courbe)
plt.plot([t0,a],[0,0])
plt.axvline(a, color='b', linestyle='dashed', linewidth=1)
plt.show()
Considérons un modèle de croissance de population - peu réaliste - où le taux de croissance per capita est linéaire en la taille de population:
$$ \dot x(t)=\lambda\,x(t)^2\,,\ x(0)=x_{0} $$On peut vérifier que la solution de cette équation est:
$$ x_{t} = \frac{x_{0}}{1-\lambda\,x_{0}\,t} $$qui est définie sur $t\in[0,\frac{1}{\lambda\,x_{0}}[$ et on a effectivement $x(t)\uparrow\infty$ lorsque $t\uparrow \frac{1}{\lambda\,x_{0}}$. Il n'est pas possible d'étendre la solution à tout $t>0$ dans la mesure où $f(x)=\lambda\,x^2$ n'est pas à croissance au plus linéaire, donc pas globalement lipschitzienne, à cause de son comportement en $+\infty$.
In [21]:
import numpy as np
import matplotlib.pyplot as plt
lbd, x0 = 2, 1
t0, t1 = 0, 1/(lbd*x0)
le_temps = np.linspace(t0,t1,300, endpoint=False)
courbe = lambda t: x0/(1-lbd*x0*t)
plt.plot(le_temps, courbe(le_temps))
plt.plot([t1,t1],[0,100], linestyle='--')
plt.ylim([0,100])
plt.legend()
plt.show()
On considère l'EDO linéaire vectoriel:
\begin{align*} \dot x(t) = A \,x(t) \,,\hskip 1em \ t\geq 0\,,\ x(0)=x_{0} \tag{EDO L1} \end{align*}où $x(t)$ est de dimension $d$ et $A$ est une matrice constante de dimension $d\times d$. On note que la fonction $f(x)=A\,x$ est globalement lipschitzienne et donc rentre dans le cadre précédent. Pour comprendre comment la solution de cette équation se construit, on peut reprendre le procédé de construction des approximations successives de Picard. On se donne $x^{(0)}(t)=0$ pour tout $t$, et on pose:
\begin{align*} x^{(n)}(t) = x_0 + A \,\int_0^t x^{(n-1)}(s) \,{\rm d}s \end{align*}donc
\begin{align*} x^{(n)}(t) &= x_0 + A \,\int_0^t \big[x_0 + A \,\int_0^{t_1} x^{(n-2)}(t_2) \,{\rm d}t_2\big] \,{\rm d}t_1\\ &= x_0 + t\,A\,x_0 + A^2 \,\int_0^t \int_0^{t_1} x^{(n-2)}(t_2) \,{\rm d}t_2 \,{\rm d}t_1\\ &= x_0 + t\,A\,x_0 + \frac{t^2}{2}\,A^2\,x_0+ A^3 \,\int_0^t \int_0^{t_1} \int_0^{t_2} x^{(n-3)}(t_3) \,{\rm d}t_3 \,{\rm d}t_2 \,{\rm d}t_1\\ &= x_0 + t\,A\,x_0 + \frac{t^2}{2}\,A^2\,x_0+ \frac{t^3}{3!}\,A^3\,x_0+ \cdots+ \frac{t^{n-1}}{(n-1)!}\,A^{n-1}\,x_0 \end{align*}où $A^n=A\times\cdots\times A$ ($n$ fois). Lorsque $n\to \infty$ nous avons déjà indiqué que $x^{(n)}(t)$ converge vers la solution $x(t)$ de (EDO L1) et on obtient:
\begin{align*} x(t) = \big[I + t\,A + \frac{t^2}{2}\,A^2 + \frac{t^3}{3!}\,A^3+ \cdots \big] \,x_0 \end{align*}On définit l'exponentielle de matrice:
\begin{align*} e^{t\,A} = I + t\,A + \frac{t^2}{2}\,A^2 + \frac{t^3}{3!}\,A^3+ \cdots \end{align*}on obtient finalement que:
\begin{align*} x(t) = e^{t\,A} \,x_0 \end{align*}La matrice exponentielle $e^M$ d'une matrice $M$ de dimension $d\times d$ est définie par: \begin{align*} e^{M} = I + M + \frac{1}{2}\,M^2 + \frac{1}{3!}\,M^3+ \cdots \end{align*} et vérifie
In [22]:
from scipy.linalg import expm, norm
"exponentiel d'une matrice, norme d'une matrice"
matrice_nulle = np.zeros((20,20))
matrice_identite = np.eye(20)
norm(expm(matrice_nulle) - matrice_identite)
Out[22]:
La solution de l'EDO linéaire:
\begin{align*} \dot x(t) = A \,x(t) + b(t) \,,\hskip 1em \ t\geq 0\,, \tag{EDO L2} \end{align*}peut s'obtenir par la méthode de la variation de la constante qui consiste à chercher une solution $x(t)$ de (EDO L2) qui soit de la forme $x(t)=e^{t\,A}\,y(t)$, le but est de déterminer $y(t)$ de telle sorte que $x(t)$ soit solution de (EDO L2). La dérivée de $x(t)$ s'écrit:
\begin{align*} { \frac{\rm d}{{\rm d}t}} x(t) = { \frac{\rm d}{{\rm d}t}}[e^{t\,A}\,y(t)] = { \frac{\rm d e^{t\,A}}{{\rm d}t}} \, y(t) + e^{t\,A}\, \dot y(t) = A\, e^{t\,A} \, y(t) + e^{t\,A}\, \dot y(t) \end{align*}mais d'après (EDO L2) $\dot x(t)=A \,x(t) + b(t)$, donc:
\begin{align*} A \,x(t) + b(t) = A\, x(t) + e^{t\,A}\, \dot y(t) \end{align*}et ainsi $\dot y(t) = e^{-t\,A}\,b(t)$ et comme $x(0)=I\,y(0)=x_0$, on obtient $y(t)=x_0+\int_0^t e^{-s\,A}\,b(s){\rm d} s$ et:
\begin{align*} x(t) = e^{t\,A}\,x_0 + e^{t\,A} \int_0^t e^{-s\,A}\,b(s){\rm d} s \end{align*}qui est donc la solution explicite de (EDO L2).
On considère l'EDO linéaire \begin{align*} \dot x(t) = A(t) \,x(t) \,,\hskip 1em \ t\geq 0\,, \tag{EDO L3} \end{align*} où la matrice dépend continuement du temps $t$. L'approche précédente se généralise en introduisant la résolvante qui est une famille de matrice $\Phi(t,s)$ indicée par $t,s\geq 0$ définie par:
\begin{align*} {\textstyle \frac{\rm d}{{\rm d}t}} \Phi(t,s) = A(t) \,\Phi(t,s) \,,\hskip 1em \ t\geq s\,,\, \Phi(s,s)=I \tag{résolvante} \end{align*}On peut vérifier que $\Phi(t,s)=\Phi(t,u)\times \Phi(u,s)$ et que $[\Phi(t,s)]^1=\Phi(s,t)$. Bien qu'on ne dispose pas de solution explicite à l'équation de la résolvante, on peut montrer que la solution $x(t)$ de (EDO L3) s'exprime en fonction de la résolvante de la façon suivant:
\begin{align*} x(t) = \Phi(t,0)\,x_0 \end{align*}et que la solution $x(t)$ de \begin{align*} \dot x(t) = A(t) \,x(t) + b(t) \end{align*} est donnée par \begin{align*} x(t) = \Phi(t,0)\,x_0 + \Phi(t,0) \int_0^t \Phi(0,s)\,b(s){\rm d} s \end{align*}
On considère le problème de Cauchy autonome suivant:
\begin{align*} \dot x(t) = f(x(t))\,,\ t\geq 0\,,\ x(0)=x_{0} \tag{EDO} \end{align*}pour laquelle on suppose qu'il existe une solution globale. On peut par exemple supposé que $f$ est globalement lipschitzienne.
Supposons que la solution $x(t)$ de (EDO) converge vers un point $x^*$ de $\mathbb{R}^d$, alors nécessairement $x(t)$ va cesser d'évoluer au "bout d'un certain temps" de telle sorte que $\dot x(t)=0$ et donc nécessairement $f(x^*)=0$:
Un point $x^*\in\mathbb R^d$ tel que $f(x^*)=0$ est appelé point d'équilibre de (EDO).
Dans un premier temps sans se poser la question du comportement asymptotique de la solution de l'EDO, on peut chercher à déterminer les points d'équilibre de l'EDO.
Dans un second temps, on peut chercher à caractériser la nature de ces points d'équilibre:
Un point d'équilibre $x^*$ sera dit stable si:
\begin{align*} \forall \epsilon>0\,,\ \exists \eta>0\,:\ |x_{0}-x^*| <\eta \Rightarrow |x(t)-x^*|<\epsilon\,,\ \forall t\geq 0\,. \end{align*}Si de plus:
\begin{align*} |x(t)-x^*|\to 0\textrm{ lorsque }t\to \infty \end{align*}alors le point est dit asymptotiquement stable. Un point d'équilibre qui n'est pas stable est dit instable.
Il existe une méthode d'analyse de la stabilité des EDO assez générale dite de Lyapounov, elle est puissante mais assez difficile à utiliser en pratique.
Nous allons présenter ici une étude locale de la stabilité consistant à linéariser (EDO) autour du point d'équilibre $x^*$ considéré. Cela revient donc à étudier un système linéaire. Nous présentons donc la notion de stabilité dans le cadre linéaire et nous reviendrons plus tard sur le cas linéaire.
On considère le problème de Cauchy autonome linéaire suivant:
\begin{align}\label{eqCauchyLin} \dot x(t) = A\,x(t)\,,\ t\geq 0\,,\ x(0)=x_{0} \end{align}Le point $0$ est un point d'équilibre de ce système. Notons que:
\begin{align*} x(t) = e^{t\,A}\,x_0 \textstyle = (I+t\,A+\frac{t^2}{2}\,A^2+\frac{t^3}{3!}\,A^3+\cdots)\,x_0 \end{align*}maintenant supposons que $x_0$ soit tel qu'il existe un réel $\lambda$ (voire un complexe, cf. infra) vérifiant $A\,x_0=\lambda\,x_0$ alors $A^n\,x_0=\lambda^n\,x_0$ et:
\begin{align*} x(t) \textstyle = x_0+t\,\lambda\,x_0+\frac{(t\,\lambda)^2}{2}\,x_0+\frac{(t\,\lambda)^3}{3!}\,x_0+\cdots = e^{\lambda\,t}\,x_0 \end{align*}Donc lorsque $\lambda<0$ (du moins sa partie réelle est strictement négative dans la cas complexe) alors $x(t)$ converge exponentiellement vite vers $0$.
Les $x_0$ tels qu'il existe $A\,x_0=\lambda\,x_0$ sont des vecteurs de $\mathbb R^d$ qui ne sont pas modifiés par la matrice $A$ (à un coefficient de proportionnalité près, éventuellement négatif), ils sont appelés vecteurs propres de la matrices $A$ et les $\lambda$ les valeurs propres de la matrice $A$.
Valeurs propres et vecteurs propres d'une matrice
Remarque: En mathématiques, ces concepts sont traditionnellement introduits dans le cadre des endomorphismes. Les endormorsphismes sur $\mathbb R^d$ sont simplement les applications linéaires de $\mathbb R^d$ dans lui-même, mais ces applications se mettent toutes sous la forme $x\to A\,x$ pour une certaine matrice carrée $A$ de taille $d\times d$; on peut donc identifier ces endomorphismes aux matrices associées et donc introduire ces concepts directement sur les matrices. Toutefois la plupart des références en français sur le sujet refusent ce raccourci, nous proposons donc une référence en anglais Jeffrey 2010.
Etant donnée une matrice réelle $A$ de dimension $n\times n$, un scalaire $\lambda\neq 0$ sera appelé valeur propre de $A$ lorsqu'il existe $x\in\mathbb R^d$ non-nul tel que
\begin{align*} A\,x = \lambda x \end{align*}le vecteur $x$ est alors appelé vecteur propre de $A$ associé à la valeur propre $\lambda$.
Il est possible de trouver une solution $x$ non nulle à $ A\,x = \lambda\,x $ si et seulement si le déterminant de la matrice $A-\lambda\,I$ est nul (lorsque ce déterminant est non nul, alors le seul $x$ vérifiant $ (A-\lambda\,I)\,x = 0 $ est $x=0$). On définit le spectre de la matrice $A$ par:
\begin{align*} Sp(A) := \left\{\lambda\in\mathbb C\,;\, \det(A-\lambda\,I)=0\right\} \end{align*}notez que
\begin{align*} p_A(\lambda) := \det(A-\lambda\,I) \end{align*}est un polynôme dit caractéristique de degré $d$, et les valeurs propres sont les racines de ce polynôme c'est à dire les solution de:
\begin{align*} p_A(\lambda) = 0 \end{align*}appelée équation caractéristique. Le théorème de d'Alembert-Gauss assure que le polynôme caractéristique admet $d$ racines (éventuellement complexes et non nécessairement distinctes).
Notez que de considèrer la solution de \eqref{eqCauchyLin} avec une condition initiale complexe ne pose aucune difficulté.
Le point $0$ est un point d'équilibre de \eqref{eqCauchyLin} et:
Ce résultat s'appuie en fait sur la compréhension de la résolvante du système \eqref{eqCauchyLin} qui est en effet donnée par [Agarwal and O'Regan 2008 p. 133]:
\begin{align*} \Phi(t,0) = \bigl[ e^{\lambda_{1}\,t}\,v_{1},\dots,e^{\lambda_{n}\,t}\,v_{n}\bigr] \end{align*}où $\lambda_{i}$ sont les valeurs propres de $A$ et $v_{i}$ des vecteurs propres associés.
Exemple 1
\begin{align*} \frac{d}{dt}\begin{pmatrix}x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix}0 & 1 \\ -1 & 0 \end{pmatrix}\, \begin{pmatrix}x_1 \\ x_2 \end{pmatrix} \end{align*}alors
\begin{align*} \det(A-\lambda\,I) = \det\begin{pmatrix}-\lambda & 1 \\ -1 & -\lambda \end{pmatrix} = \lambda^2+1 \end{align*}le polynôme caractéristique admet deux racines complexes $\lambda_1 = +i$ et $\lambda_2 = -i$. En fait $0$ n'est pas stable, pour le comprendre il suffit de remarquer que: $A^{4n+1}=A,\, A^{4n+2}=-I,\, A^{4n+3}=-A,\, A^{4n+4}=I$ pour tout $n$, donc: \begin{align*} e^{t\,a} = \begin{pmatrix} \cos(t) & \sin(t) \\ -\sin(t) & \cos(t) \end{pmatrix} \end{align*}
Ainsi
\begin{align*} x(t) = e^{t\,a} x(0) = \begin{pmatrix} \cos(t)\,x_1(0) +\sin(t)\,x_2(0) \\ -\sin(t) \,x_1(0)+ \cos(t)\,x_2(0) \end{pmatrix} \end{align*}et
\begin{align*} |x(t)|^2= x_1(t)^2+x_2(t)^2 = [\cos(t)\,x_1(0) +\sin(t)\,x_2(0)]^2+[ -\sin(t) \,x_1(0)+ \cos(t)\,x_2(0)]^2 = |x(0)|^2 \end{align*}donc les trajectoires de $x(t)$ décrivent des cercles centrés en 0 et de rayon égal à $|x(0)|^2$.
En dimension 2 on peut avoir une classification précises des points d'équilibre en fonction du déterminant de $A$ et de sa trace:
\begin{align*} \textrm{tr}(A) = A_{1,1}+A_{2,2}+\cdots+A_{d,d} \end{align*}On obtient:
[détails ici]
On considère le problème de Cauchy autonome non-linéaire suivant:
\begin{align*} \dot x(t) = f(x(t))\,,\ t\geq 0\,,\ x(0)=x_{0} \tag{EDO NL A} \end{align*}et $x^*$ un point d'équilibre. On suppose $f$ différentiable et à dérivées continues.
On peut supposer que $f(0)=0$ et donc que $x^*=0$ en faisant le changement de variable $x\to x-x^*$ dans (EDO NL A). On se place donc dans le cas $f(0)=0$ et on étudie les propriétés du point d'équilibre 0. On note $J_f(0)$ la matrice jacobienne de $f$ au point d'équilibre~$0$.
Exemple: Les points d'équilibre de l'équation logistique.
Une méthode classique consiste à étudier le système linéarisé en $0$:
\begin{align*} \dot x(t) = J_f(0)\,x(t)\,,\ t\geq 0\,,\ x(0)=x_{0} \end{align*}où $J_f(x)$ est la matrice jacobienne de $f$ en $x$:
\begin{align*} J_f(x) := \left(\begin{matrix} \partial f_{1}(x)/\partial x_{1} & \partial f_{1}(x)/\partial x_{2} & \cdots & \partial f_{1}(x)/\partial x_{d} & \\ \partial f_{2}(x)/\partial x_{1} & \partial f_{2}(x)/\partial x_{2} & \cdots & \partial f_{2}(x)/\partial x_{d} & \\ \vdots & \vdots & & \vdots & \\[0.5em] \partial f_{d}(x)/\partial x_{1} & \partial f_{d}(x)/\partial x_{2} & \cdots & \partial f_{d}(x)/\partial x_{d} & \end{matrix}\right) \end{align*}Alors:
Si toutes les valeurs propres ont des parties réelles nulles alors $0$ est stable pour l'EDO linéarisée mais on ne peut rien dire de $0$ pour l'EDO non linéaire de départ.
Le module integrate de SciPy propose deux outils pour intégrer numériquement des EDO: integrate.odeint et integrate.ode.
integrate.odeint est une fonction implémentant le solveur LSODA d'ODEPACK qui bascule automatiquement entre une méthode de prédiction-correction d'Adam Adams pour des problèmes non-raides à une méthode BDF pour les problèmes raides.
integrate.ode est une classe offrant une interface orientée objet donnant accès à plusieurs solveurs.
integrate.ode est plus flexible mais integrate.odeint est nettement plus accessible. Nous nous limiterons ici à ce dernier solveur.
La fonction odeint comprend 3 arguments obligatoires : une fonction définissant le membre de droit de l'EDO, i.e. la fonction $f$, un array donnant la condition initiale, est un array de valeurs de $t$ où $x(t)$ doit être calculé. La fonction du membre de droite comprends 2 arguments obligatoires, un array $x$ et un scalaire $t$, et des arguments optionnels.
Voir [Johansson 2015 p. 223].
On considère une solution $x(t)$ de (EDO) définie pour $t\in[0,T]$ avec $t_{0}=0$. On se donne $\delta>0$, on pose $t_{k}=k\,\delta$ et on considère le développement de Taylor:
\begin{align*} x(t+\delta) = x(t) + \delta\,x'(t) + \frac{\delta^2}{2}\,x''(t) + \frac{\delta^3}{3!}\,x'''(t)+\cdots \end{align*}À l'ordre $n$, on a la formule:
\begin{align*} x(t+\delta) = x(t) + \delta\,x'(t) + \frac{\delta^2}{2}\,x''(t) + \frac{\delta^3}{3!}\,x'''(t)+\cdots+\frac{\delta^n}{n!}\,x^{(n)}(\tau) \end{align*}pour un certain $\tau\in[t,t+\delta]$
L'approximation d'Euler de $x(t)$ consiste à tronquer à l'ordre 1:
\begin{align*} x(t+\delta) = x(t) + \delta\,x'(t) +\frac{\delta^2}{2}\,x^{''}(\tau) \end{align*}donc
\begin{align*} x(t+\delta) = x(t) + \delta\,f(t,x(t)) + O(\delta^2) \end{align*}(ordre 1: erreur locale en $O(\delta^2)$, erreur globale en $O(\delta)$)
On obtient l'approximaton: \begin{align*} \tilde x_{k+1} = \tilde x_k+\delta\,f(t_{k},\tilde x_k) \end{align*}
$\delta$ est le pas de l'approximation, on définit $\tilde x(t)$ comme étant l'interpolée linéaire des points $\tilde x_k$.
Cette méthode reprend la méthode d'Euler en utilisant une approximation de la dérivée au point milieu. Le développement de Taylor à l'rdore 3 donne:
\begin{align*} x(t+\delta) = x(t) + \delta\,x'(t) + \frac{\delta^2}{2}\,x''(t) + \frac{\delta^3}{3!}\,x'''(\tau) \end{align*}mais
$x'(t+\delta/2) = x'(t)+ \frac{\delta}{2}\,x''(t)+O(\delta^2)$ donc $\frac{\delta^2}{2}\,x''(t)= \delta[x'(t+\delta/2) - x'(t)] +O(\delta^3)$ et
\begin{align*} x(t+\delta) &= x(t) + \delta x'(t) + \frac{\delta^2}{2}\,\frac{x'(t+\delta/2)-x'(t)}{\delta/2} + O(\delta^3) \\ &= x(t) + \delta x'(t+\delta/2) + O(\delta^3) \\ &= x(t) + \delta f(t+\delta/2,x(t+\delta/2)) + O(\delta^3) \\ &= x(t) + \delta f\big(t+\delta/2,x(t)+(\delta/2)\,f(t,x(t))\big) + O(\delta^3) \end{align*}il s'agit donc d'un schéma d'ordre 2:
\begin{align*} \tilde x_{k+1} = \tilde x_k+\delta\,f\Big(t_{k}+\frac{\delta}{2},\tilde x_k+\frac{\delta}{2}\,f(t_k,\tilde x_k)\Big) \end{align*}
In [2]:
import numpy as np
import matplotlib.pyplot as plt
def euler( f, x0, t ):
n = len( t )
x = np.array( [x0] * n )
for i in range(n-1):
x[i+1] = x[i] + ( t[i+1] - t[i] ) * f( x[i], t[i] )
return x
def rk2a( f, x0, t ):
n = len( t )
x = np.array( [ x0 ] * n )
for i in range( n - 1 ):
h = t[i+1] - t[i]
k1 = 0.5 * h * f( x[i], t[i] )
x[i+1] = x[i] + h * f( x[i] + k1, t[i] + 0.5 * h )
return x
def f( x, t ):
return x * np.sin( t )
x0, T = -1., 10.0
n = 51
t = np.linspace( 0, T, n )
t2 = np.linspace( 0, T, 2*n )
x = -np.exp( 1.0 - np.cos( t ) )
x_euler = euler( f, x0, t )
x_euler2 = euler( f, x0, t2 )
x_rk2 = rk2a( f, x0, t )
plt.plot(t,x,label=r'$x(t)$')
plt.plot(t,x_euler,label='Euler',linestyle='--')
plt.plot(t2,x_euler2,label=r'Euler $\delta/2$',linestyle='--')
plt.plot(t,x_rk2,label='RK2',linestyle='--')
plt.legend()
plt.show()
[^](#ref-1) [^](#ref-5) Sylvie Benzoni-Gavage. 2010. Calcul différentiel et équations différentielles.
[^](#ref-2) Demailly, J.P.. 2006. Analyse numérique et équations fifférentielles.
[^](#ref-3) Driss Boularas and Daniel Fredon and Daniel Petit. 2009. Mini Manuel de Mathématiques pour les sciences de la vie et de l'environnement.
[^](#ref-4) [^](#ref-6) [^](#ref-8) R.P. Agarwal and D. O'Regan. 2008. An introduction to ordinary differential equations.
[^](#ref-7) Alan Jeffrey. 2010. Matrix Operations for Engineers and Scientists: An Essential Guide in Linear Algebra.
[^](#ref-9) Robert Johansson. 2015. Numerical Python. A Practical Techniques Approach for Industry.
[^](#ref-10) E. Hairer and S. P. N\orsett and G. Wanner. 1993. Solving Ordinary Differential Equations I : Nonstiff Problems.