04_edo.md



Équations différentielles ordinaires

Introduction
Pour illustrer le concept d’équations différentielles, nous allons
considérer pour commencer des systèmes qui évoluent dans le temps
(évolution d’une population, taux d’intérêts, circuits électriques,
...).

Mouvement rectiligne uniforme
Imaginons que nous connaissons la fonction décrivant le vitesse d’une
particule au cours du temps et notons la v(t). Nous savons également
que la vitesse d’une particule est reliée à l’évolution au cours du temps
de sa position. Cette dernière peut être notée, x(t). En particulier,
nous avons que la vitesse n’est rien d’autre que la dérivée de la
position. On peut donc écrire une équation reliant la vitesse à la
position x'(t)=v(t). Cette équation est appelée équation
différentielle, car elle fait intervenir non seulement les fonctions
x(t) et v(t), mais également la dérivée de la fonction x(t). Si
maintenant nous précisons ce que vaut la fonction v(t) nous pourrons
résoudre cette équation. Comme le nom de la sous-section le laisse
entendre, nous nous intéressons à un mouvement rectiligne uniforme, qui décrit

le mouvement d’un objet qui se déplace à
vitesse constante, v(t)=v. Nous cherchons ainsi à résoudre
l’équation différentielle x'(t)=v. Ou en d’autres termes, nous
cherchons la fonction dont la dérivée donne une constante[^3]. Vous savez sans
doute que l’ensemble de fonctions satisfaisant la contrainte précédente
est x(t)=v\cdot t+B, où B est une constante arbitraire. Cette solution
générale  n’est pas
unique, car nous obtenons une infinité de solutions (comme quand nous avons
calculé la primitive d’une fonction au chapitre précédent). Afin de
trouver une solution unique, nous devons imposer une condition, typiquement une “condition initiale”
à notre équation différentielle. En effet, si nous imposons la condition
initiale x(t_0)=x_0, il vient
x(t_0)=x_0=v\cdot t_0+B \Leftrightarrow B=x_0-v\cdot t_0.
Finalement, la solution du problème différentiel est donnée par
x(t)=v\cdot (t-t_0)+x_0.

Remarque {-}
La solution de l’équation différentielle x'(t)=v,\ x(t_0)=x_0,
revient à calculer \begin{aligned}
\int x'(t){\mathrm{d}}t=\int v {\mathrm{d}}t,\\
x(t)=v\cdot t + B.\end{aligned}

Mouvement rectiligne uniformément accéléré
Dans le cas du mouvement rectiligne d’un objet dont on le connaît que
l’accélération, a(t), on peut également écrire une équation
différentielle qui décrirait l’évolution de la position de l’objet en
fonction du temps. En effet, l’accélération d’un objet est la deuxième
dérivée de la position, soit x''(t)=a(t), ou encore la première
dérivée de la vitesse. \begin{aligned}
v'(t)&=a(t),\\
x'(t)&=v(t).\end{aligned}
Par simplicité supposons que l’accélération est constante, a(t)=a, donc que le mouvement est uniformément accéléré.
On
doit  résoudre[^4] x''(t)=a, ou \begin{aligned}
v'(t)&=a,\\
x'(t)&=v(t).\end{aligned}{#eq:xpv} Pour résoudre ce système
d’équations  nous résolvons la première équation
pour v(t) pour trouver v(t)=a\cdot t+C. En substituant ce résultat dans
l’@eq:xpv, on a x'(t)=a\cdot t+C. On peut ainsi
directement intégrer des deux côtés comme vu dans la sous-section
précédente \begin{aligned}
\int x'(t){\mathrm{d}}t&=\int (a\cdot t+C){\mathrm{d}}t,\nonumber\\
x(t)&=\frac{a}{2}\cdot t^2+C\cdot t + D.\end{aligned} On voit  que
la position d’un objet en mouvement rectiligne uniformément accéléré est
donné par une parabole. Cette équation a néanmoins  encore deux
constantes indéterminées. Pour les déterminer, on doit  imposer deux
conditions initiales. Une possibilité est d’imposer une condition
initiale par équation v(t_0)=v_0,\mbox{ et } x(t_0)=x_0. On obtient
v(t_0)=v_0=a\cdot t_0+C \Leftrightarrow C=v_0-a\cdot t_0, et
x(t_0)=x_0=\frac{a}{2}\cdot t_0^2+D \Leftrightarrow D=x_0-\frac{a}{2}\cdot t_0^2.
Finalement la solution est donnée par
x(t)=\frac{a}{2}\cdot (t^2-t_0^2)+v_0\cdot (t-t_0)+x_0.

Remarque {-}
La solution du problème différentiel peut également se calculer de
la façon suivante x''(t)=a,\ x(t_0)=x_0,\ v(t_0)=v_0. revient à
calculer \begin{aligned}
\int \int x''=\int \int a,\\
x(t)=\frac{a}{2}t^2+C\cdot t + D.\end{aligned}

Évolution d’une population
Imaginons une colonie de bactéries dont nous connaissons le taux de
reproduction r. Nous connaissons le nombre de ces bactéries au temps
t, qui est donné par n(t). Nous souhaitons connaître la population
au temps t+\delta t. On a donc
n(t+\delta t)=n(t)+(r\delta t)\cdot n(t)=n(t)(1+r\delta t).{#eq:evolpop}
Imaginons que le taux de reproduction r=1/3600 s^{-1}, que la
population à un temps donné t_0 est de n(t_0)=1000, et qu’on veuille
connaître la population après \delta t=1h=3600s. Il vient alors
n(t_0+3600)=(1+1/3600 \cdot 3600)\cdot n(t_0)=2\cdot1000=2000.
Imaginons maintenant que nous voulions calculer la population après
\delta t=2h=7200s. Nous avons deux façons de faire. Soit nous
utilisons le résultat précédent n(t_1)=2000 avec t_1=t_0+3600 et
évaluons la population après une heure supplémentaire
(\delta t_1=3600s)
n(t_1+3600)=(1+1/3600 \cdot 3600)\cdot n(t_1)=2\cdot 2000=4000.{#eq:comp}
Soit nous reprenons l’équation de départ (voir l'@eq:evolpop) et nous
obtenons
n(t_0+7200)=(1+1/3600 \cdot 7200)\cdot n(t_0)=3\cdot 1000=3000. On
voit que ces deux résultats ne sont pas égaux. Effectuer deux itérations
de notre algorithme discret avec un pas d’itération de \delta t, ne
correspond pas à effectuer une seule itération avec un pas deux fois
plus grand (2\delta t). Néanmoins cela devrait être le cas pour
\delta t\rightarrow 0.
Pour nous en convaincre faisons l’exercice suivant. Reprenons l’@eq:comp que vous pouvons réécrire comme
n(t_0+2\delta t)=n(t_1+\delta t)=(1+r\delta t) n(t_1)=(1+r  \delta t)(1+r  \delta t) n(t_0)=(1+r\delta t)^2 n(t_0).
Si à présent nous comparons les résultats obtenus pour
\delta t_1=2\delta t dans l’@eq:evolpop on a
\begin{aligned}
n_1&=(1+r\delta t)^2 n(t_0)=(1+2r\delta t+(r\delta t)^2) n(t_0),\\
n_2&=(1+2r\delta t) n(t_0).\end{aligned} On trouve donc finalement
que n_2-n_1=(r\delta t)^2n(t_0). On a donc que la différence tend bien
vers 0 quand \delta t tend vers 0.
Afin de voir plus en détail ce qu’il se passe lorsque
\delta t\rightarrow 0, reprenons l’équation de départ
(l'@eq:evolpop), divisons la par \delta t et arrangeons les
termes. Il vient \frac{n(t+\delta t)-n(t)}{\delta t}=r\cdot n(t). En
prenant la limite \delta t\rightarrow 0 on voit apparaître la dérivée
dans le membre de gauche de l’équation ci-dessus
\lim\limits_{\delta t\rightarrow 0} \frac{n(t+\delta t)-n(t)}{\delta t}=n'(t)=r\cdot n(t).{#eq:cont}
On voit qu’on a construit ici une équation différentielle à partir d’un
système discret.
Nous pouvons à présent résoudre l’équation différentielle ci-dessus en
se souvenant que la fonction dont la dérivée est proportionnelle à la
fonction de départ est l’exponentielle. Il vient
n(t)=C\exp(r t),{#eq:sol_pop} où C est une constante. Il est
en effet élémentaire de montrer que cette solution satisfait l’@eq:cont. On voit également qu’il nous manque une condition pour
avoir l’unicité de la solution ci-dessus (on ne connaît toujours pas
C). La constante peut-être obtenue à l’aide d’une condition initiale
(correspondant au n(t_0) de tout à l’heure). Si n(t_0)=n_0, nous trouvons
pour C n(t_0)=C\exp(r t_0)=n_0 \Leftrightarrow C=n_0\exp(-r t_0).
substituant cette relation dans l'@eq:sol_pop, on
obtient n(t)=n_0\exp(r (t-t_0)).

Autres illustrations de l’utilisation des équations différentielles
La plupart des systèmes naturels (ou moins naturels) peuvent être
décrits à l’aide d’équations différentielles. Nous allons en écrire
quelques exemples ci-dessous.

Systèmes proies-prédateurs
Considérons un système où nous avons des prédateurs (des guépards) et
des proies (des antilopes)[^5]. Supposons que les antilopes se
reproduisent exponentiellement vite et que leur seul moyen de mourir est
de se faire manger par les guépards et que la chance de se faire manger
est proportionnelle au nombre de guépards. Les guépards meurent
exponentiellement vite de faim et se reproduisent proportionnellement au
nombre d’antilopes se trouvant dans le système.
Avec ces hypothèses, on peut écrire le système d’équations suivant (a
est le nombre d’antilopes, et g le nombre de guépards)
\begin{aligned}
\frac{{\mathrm{d}}a}{{\mathrm{d}}t}&= \underbrace{k_a a(t)}_{(1)}-\underbrace{k_{g,a}g(t) a(t)}_{(2)},\\
\frac{{\mathrm{d}}g}{{\mathrm{d}}t}&= -\underbrace{k_g g(t)}_{(3)} +\underbrace{k_{a,g} a(t)g(t)}_{(4)}\end{aligned}
Le terme (1) représente la reproduction des antilopes avec taux k_a.
Le terme (2) représente la mort des antilopes qui se font manger par
les guépards avec un taux k_{g,a} (la chance qu’un guépard rencontre
une antilope). Le terme (3) est la mort des guépards avec un taux
k_g. Finalement le terme (4) est la reproduction des guépards
proportionnelle au nombre d’antilopes avec un taux k_{a,g}.
Nous avons à faire ici à un système d’équations différentielles. Nous
n’allons pas nous intéresser aux détails de la résolution de ce système mais
simplement étudier le comportement de la solution (voir la @fig:lkA et @fig:lkB).


Circuits électriques: le circuit RC
Supposons que nous ayons le circuit RC de la Fig. @fig:rc, où nous
avons une résistance (de résistance R) branchée en série avec une
capacité (de capacité électrique C). Sur ce circuit nous avons une
source qui délivre une tension U. Nous avons également un interrupteur
qui quand il est en position (a) relie le circuit RC à la source, ce
qui a pour effet de chargé la capacité. En position (b) la capacité se
décharge et son énergie est dissipée dans la résistance.

Nous souhaitons étudier la variation de la chute de tension dans la
capacité U_c lorsque:


nous mettons l’interrupteur en position (a).


puis lorsque la capacité est chargée, nous mettons l’interrupteur en
position (b).


Les chutes de tension dans la capacité et la résistance sont
respectivement données par U_C=Q/C,\quad U_R=R I, où Q est la
charge de la capacité et I le courant traversant la résistance. Nous
avons par la loi de Kirchoff que U=U_C+U_R.{#eq:tot_tension} De
plus le courant traversant la résistance est donné par I(t)=Q'(t).
En combinant ces équations, nous obtenons
U_C'(t)+\frac{U_C(t)}{RC}=\frac{U}{RC}. Nous avons également la
condition initiale U_C(0)=0 (la tension au moment de la mise de
l’interrupteur en position (a) est nulle).
Lors de la mise de l’interrupteur en position (b) nous avons
simplement que l'@eq:tot_tension devient
0=U_C+U_R.{#eq:tot_tension_0} On a donc que l’équation
différentielle pour l’évolution de la chute de tension dans la capacité
devient U_C'(t)+\frac{U_C(t)}{RC}=0. Et la condition initiale
devient U_C(0)=U.
Pour cette dernière équation nous avons déjà calculé une solution très
similaire et on a U_C(t)=U\exp(-t/(RC)). La tension dans la capacité
va décroître exponentiellement vite. Pour le cas de l’interrupteur en
position (a) la solution est U_C(t)=U(1-\exp(-t/(RC))). La tension
augmente exponentiellement au début, puis au fur et à mesure que la
capacité se charge il devient de plus en plus difficile de la charger.
L’augmentation de la tension se fait donc de plus en plus lentement
jusqu’à devenir une asymptote horizontale en U.

Taux d’intérêts composés
Nous voulons étudier l’augmentation d’un capital c(t) au cours du
temps qui est soumis à un taux d’intérêt annuel r qui est composé
après chaque intervalle \delta t. On peut également inclure des
dépôts/retraits d sur l’intervalle \delta t. La valeur du capital
après un intervalle \delta t est de
c(t+\delta t)=c(t)+(r\delta t )c(t)+d\delta t.{#eq:cap_discr}
Supposons qu’on a un capital de départ 1000 \mathrm{CHF}, un taux
d’intérêts annuel de 1\% et un dépôt annuel de 100\mathrm{CHF}.
Après deux mois (\delta t=2/12=1/6)  le capital devient
c(1/6)=1000+0.01/6\cdot 1000 +100/6=1018.3\mathrm{CHF}. Si
maintenant, nous voulons avoir la valeur du capital à n’importe quel
moment dans le temps, nous allons prendre \delta t\rightarrow 0. En
divisant l'@eq:cap_discr par \delta t, et en
réarrangeant les termes, on obtient c'(t)=rc(t)+d. En supposant que
c(t=0)=c_0 (le capital initial), cette équation différentielle a pour
solution c(t)=\frac{d}{r}(e^{rt}-1)+c_0e^{r t}. Cette solution a
pour les paramètres précédents la forme suivante sur une période de 100
ans.


Définitions et théorèmes principaux

Définition (Équation différentielle ordinaire) {-}
Soit y une fonction dérivable n fois et dépendant d’une seule
variable. Une équation différentielle ordinaire est un équation de
la forme F(x,y,y',y'',...,y^{(n)})=0, où F est une fonction, et
y', y'', ..., y^{(n)} sont les dérivées première, deuxième, ...,
n-ème de y.


Illustation {-}
L’équation suivante est une équation différentielle ordinaire
y''+4y'+8y+3x^2+9=0.

Afin de résoudre cette équation, nous cherchons une solution de la forme
y=f(x). On dit également que nous cherchons à intégrer l’équation
différentielle.
Afin de classifier les équation différentielles, considérons les
définitions suivantes

Définition (Ordre)  {-}
L’ordre d’une équation différentielle est l’ordre le plus haut des
dérivées de y qui y apparaissent. L’ordre de l’équation différentielle
F(x,y,y',y'',...,y^{(n)})=0 est de n, si n\neq 0.

Illustration {-}
L’équation différentielle suivante est d’ordre 3
4y'''+x\cdot y'+4y+6x=0.

Définition (Condition initiale) {-}
Une condition initiale pour une équation différentielle d’ordre n, est
un ensemble de valeurs, y_0, y_1, ..., y_{n-1} donnée telles que
pour une valeur x_0 donnée on a
y(x_0)=y_0,\ y'(x_0)=y_1,\ ...,\ y^{(n-1)}(x_0)=y_{n-1}.
Nous souhaitons maintenant savoir sous quelles conditions une équation
différentielle admet une solution et si elle est unique. Nous n’allons
pas vraiment écrire ni démontrer le théorème d’existence et d’unicité
des équations différentielles ordinaires, mais simplement en donner une
version approximative et la discuter


Théorème (Existence et unicité) {-}
Soit D\subseteq{\real} le domaine de définition de la fonction
y. Soit y:D\rightarrow E\subseteq {\real} une fonction à valeur
réelle continue et dérivable sur D, et
f:D\times E\rightarrow F\subseteq{\real} une fonction à deux variables continue
sur D\times E. Alors, le système suivant (également appelé problème de
Cauchy) \begin{aligned}
&y'=f(y,x),\\
&y(x=x_0)=y_0,
\end{aligned} admet une unique solution y(x).

Ce théorème peut être étendu à une équation d’un ordre arbitraire, n,
possédant n-1 conditions initiales. En effet, n’importe quel équation
différentielle d’un ordre n peut être réécrite sous la forme de n
équations différentielles d’ordre 1. Pour illustrer cette propriété
considérons l’équation différentielle suivante y''+3y'+y+3x=0. Si
nous définissons z=y', nous avons le système suivant à résoudre
\begin{aligned}
y'=z,\\
z'+3y'+y+3x=0.\end{aligned} Nous voyons que ce système est d’ordre 1,
mais que nous avons augmenté le nombre d’équations à résoudre.
Cette propriété peut se généraliser de la façon suivante. Soit une
équation différentielle d’ordre n F(x,y,y',...,y^{(n)})=0. Nous
pouvons définir z_i=y^{(i-1)} et on aura donc que z_{i+1}=z_i'. On
peut ainsi réécrire l’équation différentielle d’ordre n comme étant
\begin{aligned}
&z_{i+1}=z_i',\ i=1,...,n-1\\
F(x,y,y',..,y^{(n)})=0 \Rightarrow &G(x,z_1,z_2,...,z_n)=0.\end{aligned}
Jusqu’ici F peut être totalement arbitraire. Essayons de classifier un
peu les équations différentielles en fonction des propriétés de F.


Définition (Linéarité) {-}
Une équation différentielle ordinaire d’ordre n est dite linéaire si
on peut l’écrire sous la forme
a_0(x)\cdot y(x)+a_1(x)\cdot y'(x)+...+a_n(x)\cdot y^{(n)}(x)=b(x).
Si les coefficients a_i ne dépendent pas de x, alors l’équation est
dite à coefficients constants.

L’équation ci-dessus a les propriétés suivantes


Les a_i ne dépendent que de x (ils ne peuvent pas dépendre de
y).


Les y et toutes leur dérivées ont un degré polynomial de 1.


Illustration {-}
L’équation suivante est linéaire y''+4x\cdot y'=e^x.
L’équation
suivante n’est pas linéaire y\cdot y''+4x\cdot y'=e^x.

Définition (Homogénéité) {-}
Une équation différentielle ordinaire est dite homogène si le terme
dépendant uniquement de x est nul. Dans le cas où nous avons à faire à
une équation différentielle linéaire, cela revient à dire que b(x)=0.

Illustration (Homogénéité) {-}
Les équations suivantes sont homogènes \begin{aligned}
&y''+4x\cdot y\cdot y'+3x^2\cdot y^3=0,\\
&2y'''+5x^2\cdot y'=0.
\end{aligned} Les équations suivantes ne le sont pas
\begin{aligned}
&y''+4x\cdot y\cdot y'+3x^2\cdot y^3=4x+2,\\
&2y'''+5x^2\cdot y'=1.
\end{aligned}


Exercice (Homogénéité) {-}
Pour chacune de ces équations différentielles ordinaires
donner tous les qualificatifs possibles. Si l’équation est inhomogène
donner l’équation homogène associée. \begin{aligned}
&y^{(4)}+4x^2 y=0,\\
&y'+4x^2 y^2=3x+2,\\
&\frac{1}{y+1}y''+4x^2 y^2=0,\\
&y'=y,\\
&4y''+4x y=1.
\end{aligned}

La solution  des équations différencielles inhomogènes se
trouve de la façon suivante.


Trouver la solution générale de l’équation différentielle homogène associée,
notons-la y_h(x).


Trouver une solution particulière à l’équation inhomogène, notons-la
y_0(x).


La solution sera donnée par la somme de ces deux solutions
y=y_h+y_0.

Techniques de résolution d’équations différentielles ordinaires d’ordre 1
Ici nous considérerons uniquement les équations différentielles
ordinaires d’ordre 1. Pour certains types d’équations différentielles,
il existe des techniques standard pour les résoudre. Nous allons en voir
un certain nombre.

Équations à variables séparables


Définition (Équations à variable séparables) {-}
On dit qu’une équation différentielle d’ordre 1 est à variables
séparables, si elle peut s’écrire sous la forme suivante
y' a(y)=b(x).


Illustration {-}
L’équation suivante est à variables séparables
e^{x^2+y^2(x)}y'(x)=1.

Pour ce genre d’équations, la solution se trouve de la façon suivante.
Nous commençons par écrire la dérivée, y'={\mathrm{d}}y/{\mathrm{d}}x
et on obtient \begin{aligned}
\frac{{\mathrm{d}}y}{{\mathrm{d}}x} a(y)=b(x),\\
a(y){\mathrm{d}}y=b(x){\mathrm{d}}x.\end{aligned} On peut maintenant
simplement intégrer des deux côtés et on obtient
\int a(y){\mathrm{d}}y=\int b(x){\mathrm{d}}x. Si nous parvenons à
résoudre les intégrales nous obtenons une solution pour y(x) (cette
solution n’est peut-être pas explicite). Il existe le cas simple où
a(y)=1 et il vient y=\int b(x){\mathrm{d}}x.


Exemple {-}
Résoudre l’équation différentielle suivante n'(t)=r\cdot n(t).

Solution {-}
En
écrivant n'={\mathrm{d}}n /{\mathrm{d}}t, on réécrit l’équation
différentielle sous la forme
\frac{1}{n} {\mathrm{d}}n=r{\mathrm{d}}t, qu’on intègre
\begin{aligned}
\int \frac{1}{n} {\mathrm{d}}n&=\int r{\mathrm{d}}t,\nonumber\\
\ln(n)&=r\cdot t+C,\nonumber\\
n(t)&=e^{r\cdot t+C}=A\cdot e^{r\cdot t},\end{aligned} où A=e^C.


Exercice {-}


Résoudre l’équation différentielle suivante c'(t)=rc(t)+d.


Résoudre l’équation différentielle suivante
x\cdot y(x) \cdot y'(x)=1.


Équations linéaires {#sec:eq_lin}
Pour une équation du type y'(x)=a(x)\cdot y(x)+b(x),{#eq:lin}
on doit résoudre le problème en deux parties.
supposons que nous connaissons une
solution “particulière” à cette équation. Notons la y_p. Si nous
faisons maintenant le changement de variables y=y_h+y_p et remplaçons
ce changement de variables dans l’équation ci-dessus nous obtenons
y_p'(x)+y_h'(x)=a(x)\cdot y_p(x)+a(x)\cdot y_h(x)+b(x).{#eq:lin_hp}
Comme y_p est solution de l'@eq:lin on a
y_p'(x)=a(x)\cdot y_p(x)+b(x). En remplaçant cette relation dans
l'@eq:lin_hp il vient y_h'(x)=a(x)\cdot y_h(x).
Cette équation différentielle n’est rien d’autre que l’équation homogène
correspondant à @eq:lin.
Nous voyons qu’une équation inhomogène se résout en trouvant la
solution générale à l’équation homogène correspondante et en y ajoutant
une solution particulière.
Revenons donc à la résolution de l’équation différentielle linéaire
d’ordre un. La première partie de la résolution consiste à résoudre
l’équation homogène associée à l'@eq:lin
y'(x)=a(x)\cdot y(x). Cette équation se résout par séparation des
variables. La solution est donc y_h(x)=Ce^{\int a(x){\mathrm{d}}x}.
Puis nous devons chercher une solution dite particulière de l’équation
inhomogène. Pour ce faire nous utilisons la méthode de la variation de
la constante. Il s’agit de trouver une solution particulière qui aura la
même forme que la solution de l’équation homogène, où C dépendra de
x (d'où le nom de méthode de variation de la constante)
y_p(x)=C(x)e^{\int a(x){\mathrm{d}}x}. En remplaçant cette équation
dans l'@eq:lin, on obtient \begin{aligned}
C'(x)e^{\int a(x){\mathrm{d}}x}+C(x)\cdot a(x)e^{\int a(x){\mathrm{d}}x}&=a(x)\cdot C(x) e^{\int a(x){\mathrm{d}}x}+b(x),\nonumber\\
C'(x)&=\frac{b(x)}{e^{\int a(x){\mathrm{d}}x}}.
\end{aligned} Il nous reste donc à résoudre cette équation
différentielle pour C(x) qui est une équation à variables séparables où
on aurait un a(c)=1. On intègre donc directement cette équation
pour obtienir
C(x)=\int \frac{b(x)}{e^{\int a(x){\mathrm{d}}x}}{\mathrm{d}}x.
Finalement, on a que la solution de l’équation générale de l’équation
inhomogène est
y=y_p+y_h=\left(\int \frac{b(x)}{e^{\int a(x){\mathrm{d}}x}}{\mathrm{d}}x+C\right)e^{\int a(x){\mathrm{d}}x}.

Exemple {-}
Résoudre l’équation suivante
U_C'(t)+\frac{U_C(t)}{RC}=\frac{U}{RC}.{#eq:rc_inhom}

Solution {-}
On
commence par résoudre l’équation homogène
{U_C}_h'(t)+\frac{{U_C}_h(t)}{RC}=0. D’où on obtient
{U_C}_h=A\cdot e^{-\frac{1}{RC} t}. Puis par variations des
constantes, on essaie de déterminer la solution particulière de la forme
{U_C}_p=B(t)\cdot e^{-\frac{1}{RC} t}. En remplaçant cette forme de
solution dans l'@eq:rc_inhom, on obtient
B'(t)=\frac{U}{RC}\cdot e^{\frac{1}{RC} t}. Qui donne par
intégration B(t)=U e^{\frac{1}{RC} t}+D. Finalement, il vient que
U_c(t)=\left(U e^{\frac{1}{RC} t}+D+A\right)e^{-\frac{1}{RC}t}=U+(D+A)e^{-\frac{1}{RC}t}=U+Ce^{-\frac{1}{RC}t},
où C=D+A. Pour le cas de la charge du condensateur, on a de plus
U_c(0)=0. On peut donc fixer la constante C=-U.
Résoudre les équations différentielles suivantes

Exercice {-}


y'+2y=t^2


y'+y=\frac{1}{1+e^t}.


Équations de Bernouilli
Il existe des équations particulières qui peuvent se ramener à des
équations linéaires via des changements de variables.
Une classe particulière sont les équations de Bernouilli, qui s’écrit
y'(x)+a(x)\cdot y(x)+b(x)\cdot y^n(x)=0,{#eq:bernouilli} où
r\in{\real}.
Cette équation peut être réécrite sous la forme
\frac{y'(x)}{y^n(x)}+\frac{a(x)}{y^{n-1}(x)}+b(x)=0.{#eq:bernouilli_2}
Dans ce cas là, en effectuant le changement de variable suivant
z=y^{1-n}, on obtient (exercice)
z'(x)+(1-n)a(x)\cdot z(x)+(1-n)b(x)=0. On a donc ramené l’équation
de Bernouilli à une équation linéaire que nous savons résoudre à l’aide
de la méthode de la section @sec:eq_lin.


Exemple {-}
Résoudre l’équation de Bernouilli suivante y'-y-x\cdot y^6=0.

Solution {-}
Avec
la substitution z=y^5, on obtient z'-5z+5x=0. Cette équation se
résout en trouvant d’abord la solution de l’équation
homogène z_h'-5z_h=0, qui est donnée par z_h=Ae^{5x}. En
remarquant qu’une solution particulière à z_p'-5z_p+5x=0, peut être de
la forme z_p=x+B (avec B une constante) on obtient \begin{aligned}
1-5(x+B)+5x=0,\nonumber\\
1-5B=0\Rightarrow B=\frac{1}{5}.\end{aligned} Et finalement
z=z_h+z_p=Ae^{5x}+x+\frac{1}{5}. Il nous reste à présent à calculer
y=z^{1/5} et on a y=\left(Ae^{5x}+x+\frac{1}{5}\right)^{1/5}.


Équation de Riccati
L’équation de Riccati qui est de la forme
y'(x)+a(x)+b(x)\cdot y(x)+c(x)\cdot y^2(x)=0,{#eq:riccati} et
est donc quadratique en y. On notera que c’est une équation de
Bernouilli (avec n=2 et qui est inhomogène).
Cette équation a une propriété intéressante. Si nous connaissons une
solution particulière à l’équation inhomogène, notons la y_p, alors la
solution générale peut être trouvée de la façon suivante.
Faisons le changement de variable suivant y=y_h+y_p. L’équation
ce-dessus devient donc
y_p'+y_h'+a(x)+b(x)\cdot y_p+b(x)\cdot y_h+c(x)\cdot (y_p^2+2y_p(x)y_h(x)+y_h^2)=0.
En utilisant que y_p est solution de l’équation de Riccati, on a
y_h'+a(x)+(b(x)+2y_p(x)c(x))\cdot y_h+c(x)\cdot y_h^2=0. Cette
équation est une équation de Bernouilli avec n=2. On sait donc comment
la résoudre.
--

Exercice {-}
Résoudre l’équation de Riccati suivante y'+y^2-\frac{2}{x^2}=0.
Indication: la solution particulière a la forme y=\frac{a}{x}, avec
a une constante.
--
De plus, ce genre d’équation peut-être transformée via un changement de
variables en une équation linéaire d’ordre deux. Si c(x) est
dérivable, alors on peut faire le changement de variables suivant
v=y\cdot c(x), et on a donc que v'=y' c+y c'. En insérant ces
relations dans l'@eq:riccati, il vient
v'(x)+d(x)+e(x)\cdot v(x)+v^2(x)=0,{#eq:riccati_2} où nous
avons nommé d(x)=a(x)\cdot c(x) et e(x)=\frac{c'(x)}{c(x)}+b(x). Si
à présent nous faisons un autre changement de variables
v(x)=-\frac{z'(x)}{z(x)}, on obtient que l’équation ci-dessus peut
se réécrire comme
z''(x)+e(x)\cdot z'(x)+d(x)\cdot z(x)=0.{#eq:riccati_3}
L’équation de Riccati (une équation d’ordre un non-linéaire et
inhomogène) est ainsi transformée en une équation linéaire d’ordre deux.

Equations différentielles ordinaires d’ordre deux
Dans cette section, nous allons étudier des cas particuliers d’équations
différentielles que nous savons intégrer. Cela sera toujours des
équations linéaires.
De façon générale ces équations s’écrivent
a(x)y''(x)+b(x)y'(x)+c(x)y(x)=d(x), où
a,b,c,d:{\real}\rightarrow{\real} sont des fonctions
réelles. Avant de résoudre l’équation générale, nous allons considérer
des plus simples.

EDO d’ordre deux homogène à coefficients constants
Ce genre d’équations s’écrit sous la forme
a y''(x)+by'(x)+cy(x)=0.{#eq:edo2_cch} Voyons maintenant
comment résoudre cette équation.
Ces équations ont des propriétés intéressantes dûes à la linéarité de
l’équation différentielle.


Propriétés {-}
Ces propriétés (qui caractérisent le mot "linéaires") sont à démontrer en exercice.


Soit f(x) une solution de l'@eq:edo2_cch, alors
pour C\in{\real} Cf(x) est également
solution de @eq:edo2_cch.


Soient f(x) et g(x) deux solutions de l’équation
@eq:edo2_cch, alors h(x)=f(x)+g(x)
est également solution de @eq:edo2_cch.


De ces deux propriétés, on déduit la propriété suivante. Soient
f(x) et g(x) deux solutions de l'@eq:edo2_cch,
et C_1,C_2\in{\real},  h(x)=C_1f(x)+C_2g(x)
est aussi  solution de l'@eq:edo2_cch.


Afin de simplifier la discussion prenons une EDO d’ordre deux à
coefficients constants particulière y''+3y'+2y=0.{#eq:edo2_ex}
On va supposer que cette équation a pour solution une fonction de la
forme y(x)=e^{\lambda x}. Substituons cette forme de solution dans
l’équation de départ, on obtient \begin{aligned}
\lambda^2 e^{\lambda x}+3\lambda e^{\lambda x}+2\lambda^2 e^{\lambda x}=0,\nonumber\\
\lambda^2+3\lambda +2=0,\end{aligned}s où on a utilisé que
e^{\lambda x} ne peut jamais s’annuler pour le simplifier entre les
deux lignes. La seconde ligne ci-dessus, s’appelle le polynôme
caractéristique de notre EDO d’ordre 2.
Il nous reste à présent à déterminer \lambda ce qui est un simple
problème d’algèbre. Le polynome ci-dessus se factorise simplement en
(\lambda+1)(\lambda+2)=0, on a donc pour solution \lambda=-1, et
\lambda=-2.
On a donc immédiatement deux solutions à notre équation différentielle
y_1(x)=e^{-x},\quad y_2(x)=e^{-2x}. On vérifie aisément que ces deux
équations vérifient l'@eq:edo2_ex. Précédemment, nous
avons vu que la linéarité de ces équations différentielles, faisait
qu’on pouvait contrsuire des solutions plus générales. En effet, on peut
montrer que la solution la plus générale à cette EDO est
y(x)=C_1 y_1(x)+C_2y_2(x)=C_1e^{-x}+C_2e^{-2x}. On constate qu’il y
a deux constantes à déterminer pour avoir une solution unique. Pour ce
faire il faudra donner deux conditions initiales. Une sur y(x) et une
sur y'(x). Par exemple on pourrait avoir y(0)=1 et y'(0)=0 et on
obtient \begin{aligned}
C_1+C_2&=1,\\
-C_1-2C_2&=0.\end{aligned} Ce système d’équations ordinaires a pour
solution C_1=2,\quad C_2=-1. On a donc finalement
y(x)=2e^{-x}-e^{-2x}.
A présent, nous pouvons généraliser cette méthode pour l’équation
@eq:edo2_cch a y''(x)+by'(x)+cy(x)=0. En faisans la même
subsitution que précédemment, y=e^{\lambda x}, on a \begin{aligned}
&a \lambda^2e^{\lambda x}+b\lambda e^{\lambda x} +ce^{\lambda x}=0,\\
&a \lambda^2+\lambda b+c=0.\end{aligned} L’équation ci-dessus doit
être résolue pour \lambda. Nous savons comment résoudre ce genre
d’équation du second degré. La solution est donnée par
\lambda=\frac{-b\pm\sqrt{\Delta}}{2a}, où \Delta = b^2-4ac. On a
deux solutions \begin{aligned}
\lambda_1=\frac{-b-\sqrt{\Delta}}{2a},\\
\lambda_2=\frac{-b+\sqrt{\Delta}}{2a}.\end{aligned}
Il y a trois cas  possibles: \Delta > 0, \Delta = 0,
\Delta < 0.

Le cas \Delta>0

Dans ce cas, on a que \lambda_1,\lambda_2\in{\real} sont réels.
La solution est donc donnée par (comme on l’a vu au paravant)
y(x)=C_1e^{\lambda_1 x}+C_2e^{\lambda_2 x}.

Le cas \Delta=0

Ici, \lambda_1=\lambda_2=\lambda=-b/(2a) et \lambda est réel.
Dans ce cas-là les choses se compliquent un peu. Si on utilisait
directement la formule ci-dessus, on aurait y(x)=Ce^{\lambda x},
avec C\in{\real}. Par contre, cette solution ne peut pas
satisfaire deux conditions initiales comme nous avons vu précédemment.
Il fau donc travailler un peu plus. Supposons que y(x) est donné par
la fonction suivante y(x)=z(x)e^{\lambda x}, avec z(x) une
fonction réelle. En substituant cela dans l’équation générale, on a
az''+(2\lambda a+b)z'+(a\lambda^2+b\lambda+c)z=0. En utilant que
\lambda=-b/(2a) et \Delta =0 il vient z''=0. La solution de
cette équation est z=C_1+xC_2. On obtient donc comme solution
générale de l’équation différentielle y(x)=(C_1+C_2 x)e^{\lambda x}.

Le cas \Delta<0

Dans ce cas-là, on a deux solutions complexes (la racine d’une nombre
négatif n’est pas réelle). Les racines sont de la forme
\begin{aligned}
\lambda_1=\frac{-b+i\sqrt{|b^2-4ac|}}{2a},
\lambda_2=\frac{-b-i\sqrt{|b^2-4ac|}}{2a},\end{aligned} où i est l'unité
imaginaire. En écrivant u=-b/(2a) et v=\sqrt{|b^2-4ac|}/(2a),
on peut écrire \lambda_1=u+iv et \lambda_2=u-iv. On a donc que
\lambda_2 est le complexe conjugué de \lambda_1, ou
\lambda_1=\bar{\lambda}_2. En utilisant ces notations dans notre
exponentielle, on a \begin{aligned}
y_1&=e^{(u+iv)x}=e^{ux}e^{ivx},\\
y_2&=e^{(u-iv)x}=e^{ux}e^{-ivx}.\end{aligned} En se rappelant de la
linéarité des solutions des EDO linéaires, on peut écrire la forme
générale de la solution comme (C_1,C_2\in {\real})
y=C_1y_1+C_2y_2=C_1e^{ux}e^{ivx}+C_2e^{ux}e^{-ivx}=e^{ux}(C_1e^{ivx}+C_2e^{-ivx}).{#eq:sol2}
En utilisant la formule d’Euler \begin{aligned}
e^{ivx}&=(\cos(vx)+i\sin(vx)),\\
e^{-ivx}&=e^{ux}(\cos(vx)-i\sin(vx)),\end{aligned} on peut réécrire
l'@eq:sol2 comme \begin{aligned}
y&=e^{ux}\left(C_1(\cos(vx)+i\sin(vx))+C_2(\cos(vx)-i\sin(vx))\right),\nonumber\\
&=e^{ux}\left((C_1+C_2)\cos(vx)+i(C_1-C_2)\sin(vx))\right),\nonumber\\
&=e^{ux}\left(C_3\cos(vx)+C_4\sin(vx))\right),\end{aligned} où on a
définit C_3\equiv C_1+C_2 et C_4\equiv i(C_1-C_2).
Résoudre les EDO d’ordre 2 à coefficiens constants suivantes:


y''+y'+y=0,


y''+4y'+5y=0, y(0)=1, y'(0)=0.


y''+5y'+6y=0, y(0)=2, y'(0)=3.


2y''-5y'+2y=0, y(0)=0, y'(0)=1.


Résolution numérique d’équations différentielles ordinaires
Pour la plupart des problèmes d’ingénierie classique, les solutions des
équations différentielles sont trop compliquées à calculer
analytiquement (si elles sont calculables). Il est donc nécessaire d’en
obtenir des solutions approximées numériquement.

Problématique
Le problème à résoudre est une EDO avec condition initiale qui peut
s’écrire de la façon suivante y'=F(t,y),\quad y(t_0)=y_0, où F est
une fonction de y et de t, et où y_0 est la condition initiale.
Nous cherchons donc à connaître l’évolution de y(t) pour t>t_0.

Méthode de résolution: la méthode d’Euler
Afin de résoudre ce genre de problème numériquement il existe une grande
quantité de techniques. Ici nous allons en considérer une relativement
simple, afin d’illustrer la méthodologie (vous en verrez une autre dans
le TP).
Nous cherchons donc à évaluer y(t=t_0+\delta t), étant donné y_0,
\delta t et F(t,y). Intégrons donc simplement notre EDO entre t_0
et t_0+\delta t dans un premier temps et on obtient
\int_{t_0}^{t_0+\delta t} y' {\mathrm{d}}t=\int_{t_0}^{t_0+\delta t} F(t,y){\mathrm{d}}t.
Le théorème fondamental du calcul intégral nous dit que cette équation
peut s’écrire
y(t_0+\delta t)-y(t_0)=\int_{t_0}^{t_0+\delta t} F(t,y){\mathrm{d}}t,
y(t_0+\delta t)-y_0=\int_{t_0}^{t_0+\delta t} F(t,y){\mathrm{d}}t.{#eq:edo_app_gen}
Ont doit donc intégrer le membre de droite de cette équation. Pour ce
faire nous pouvons utiliser une des techniques vues dans le chapitre
précédent. Par exemple, on peut choisir la méthode des rectangle à
gauche. Cette équation devient \begin{aligned}
&y(t_0+\delta t)-y_0=\delta t F(t_0,y(t_0)),\nonumber\\
&y(t_0+\delta t)=y_0+\delta t F(t_0,y(t_0)).\end{aligned} Cette
dernière équation nous permet donc d’évaluer y(t_0+\delta t)
connaissant y_0. Cette méthode s’appelle “méthode d’Euler” et est une
dite explicite, car y(t_0+\delta t) ne dépend que de la valeur de
y évaluée au temps t_0.
Si plutôt que d’utiliser la méthode des rectangle à gauche pour
approximer l’intégrale de l'@eq:edo_app_gen, nous
utilisons la méthodes des rectangles à droite on a
y(t_0+\delta t)=y_0+\delta t F(t_0+\delta t,y(t_0+\delta t)). Dans
ce cas, on voit que la valeur y(t_0+\delta t) est calculée par rapport
à la valeur d’elle même. Dépendant de la forme de F on ne peut pas
résoudre cette équation explicitement. On a donc à faire à une équation
sous forme implicite. Cette façon d’approximer une EDO est dite
méthode d’Euler implicite.
Sans entrer dans les détails, la différence entre une méthode explicite
et une méthode implicite est une question de stabilité numérique. En
effet, les méthodes explicites peuvent devenir numériquement instables
(la solution numérique s’éloigne exponentiellement vite de la solution
de l’EDO) si \delta t devient “trop grand” (la contrainte du la taille
de \delta t s’appelle CFL, pour Courant-Friedrich-Lévy). Les méthodes
implicites ne souffrent pas de ce problème de stabilité, en revanche
elles sont plus coûteuses en temps de calcul et en complexité
algorithmique, étant donné qu’elles requièrent la résolution d’une
équation implicite.
Notre but initial était de connaître l’évolution de y(t) pour t>t_0.
Pour déterminer la valeur de y(t_1) avec t_N=t_0+N\delta t, il
suffit donc d’effectuer N pas de la méthode d’intégration choisie (ici
la méthode d’Euler explicite). On a donc que
y(t_0+N\delta t)=y_0+\delta t\sum_{i=1}^{N}F(t_i,y_i), où
t_i=t_0+i\cdot\delta t et y_i=y(t_i). Le deuxième terme du membre de
droite de cette équation est la même que la formule d’intégration en
plusieurs pas pour la méthode du rectangle (voir l’équation
@eq:rect_gauche). On a vu que cette méthode a une erreur
d’ordre \delta t. On peut en conclure que l’erreur que la précision de
la méthode d’Euler est également d’ordre \mathcal{O}(\delta t).

Méthode de résolution: la méthode de Verlet
Cette méthode d’intégration est utilisée pour l’intégration numérique
d’EDO d’ordre deux avec une forme particulière qui est donnée par
x''(t)=a(x(t)),{#eq:x2} où F est une fonction de x(t). On a
également les conditions initiales x(t_0)=x_0 et x'(t_0)=v_0. Cette
forme d’équation différentielle est bien connue en physique sous la
forme \vec F=m\vec a, qui peut s’écrire \begin{aligned}
&\vec{F}=m \vec a(t)=m \vec x''(t),\nonumber\\
&\frac{\vec{F}}{m}= \vec x''(t),\end{aligned} qui est de la forme de
l’EDO de départ de l'@eq:x2. La force peut avoir
différentes forme. Cela peut être la forme de gravité \vec F=m \vec g,
de frottement \vec F=-\zeta \vec v=-\zeta x'(t), etc ou une
combinaison de toutes ces forces.
Dans la section précédente, nous avons vu l’algorithme d’Euler pour
résoudre des EDO. Cette méthode a pour avantage sa simplicité de codage,
son faible coût de calcul, mais a pour désavantage son manque de
précision. Dans un certain nombres d’applications, telles que les
moteurs physiques pour les graphismes dans les jeux vidéos, ce manque de
précision est inacceptable et une meilleure méthode doit être utilisée.
Dans le TP vous avez vu les méthodes de Runge-Kutta. Ces méthodes
améliorent la précision de façon spectaculaire, mais ont en général un
coû de calcul trop élevé.
La méthode de Verlet qu’on va voir ci-dessous est augmente combine un
faible coût de calcul et une amélioration notable de la précision. Elle
est en effet très répandue dans l’industrie du jeu vidéo pour intégrer
les équations différentielles omniprésentes dans les moteurs physiques.
La méthode de Verlet s’écrit (en utilisant les notations de la section
précédente)
x(t_{n+1})=x(t_n)+\delta t v(t_n)+\frac{1}{2}\delta t^2 a(x(t_n)).{#eq:verlet_gen}
Considérons d’abord le terme v(t_n). Ce terme est approximé ici comme
v(t_n) = \frac{x(t_{n+1})-x(t_{n-1})}{2\delta t}. En remplaçant
cette approximation dans l’équation ci-dessus, il vient
x(t_{n+1})=x(t_n)+\frac{x(t_{n+1})-x(t_{n-1})}{2}+\frac{1}{2}\delta t^2 a(x(t_n)),
2x(t_{n+1})=2x(t_n)+x(t_{n+1})-x(t_{n-1})+\delta t^2 a(x(t_n)),
x(t_{n+1})=2x(t_n)-x(t_{n-1})+\delta t^2 a(x(t_n)).{#eq:verlet_novel}
On voit ici que cette formule est inutilisable pour évaluer x(t_1) (ce
qui veut dire que n=0 dans le cas ce-dessus), car elle fait intervenir
x(t_{-1}) dans le membre de droite. Pour résoudre ce problème il
suffit d’évaluer x(t_1) grâce à l'@eq:verlet_gen où
n=0 \begin{aligned}
x(t_{1})&=x(t_0)+\delta t v(t_0)+\frac{1}{2}\delta t^2 a(x(t_0)),\nonumber\\
x(t_{1})&=x_0+\delta t v_0+\frac{1}{2}\delta t^2 a(x_0),\end{aligned}
où x_0 et v_0 sont les conditions initiales de notre problème.
Esuite les itérations suivantes (n>0) sont calculables directement
avec l'@eq:verlet_novel. Un autre avantage
considérable de ce modèle est qu’il est très simple d’y inclure une
force de frottement proportionnelle à la vitesse. Sans entrer dans les
détails de la dérivation du schéma on a
x(t_{n+1})=(2-\delta t\zeta)x(t_n)-(1-\delta t\zeta)x(t_{n-1})+\delta t^2 a(x(t_n)).