Les équations différentielles ordinaires

Introduction

Les équations différentielles ordinaires (EDO ou ODE en anglais) : un monument de mathématique, le graal d'un physicien ! Pour moi, les équations différentielles et leur prolongement dans l'étude des systèmes dynamiques sont les chapitres les plus beaux des mathématiques, avec la topologie. Mais c'est surtout un outil phénoménal (dans tous les sens du terme) pour le physicien . Réfléchissez : sans équation différentielle, point de mécanique, point d'électricité, et j'en passe et des meilleures... Le boom de la physique au XVII et XVIIIeme siècle fut provoqué par l'introduction du calcul différentiel et de ses filles, les équations différentielles.

Vous trouverez dans cette page qui leur est consacrée de nombreuses références et de multiples renvois aux autres pages de TangenteX.com, car j'ai déjà abordé de nombreuses fois le sujet, tant dans ses applications que dans ses méthodes.

Pourquoi équations différentielles "ordinaires" ? Parce qu'il en existe une autre sorte, les équations différentielles "aux dérivées partielles", ou EDP (PDE en anglais). Elles généralisent le concept à plusieurs variables, alors que nous verrons que les EDO ne portent que sur une seule variable. Le monde des EDP est encore plus merveilleux que les EDO. Si merveilleux qu'un des problèmes du siècle en mathématique porte sur le EDP. Dotation : 1 million de dollards! Qui a dit que les maths ne payez pas !

Allons, plongeons dans ce monde ! Oh nous ne ferons que l'effleurer : d'abord, je vais essayer de vous expliquer ce qu'est une EDO à partir d'une exemple simple. Puis nous examinerons plusieurs EDO que vous avez rencontré sans le savoir dans le programme de physique de lycée. Enfin, parce nous sommes des physiciens numériciens, nous verrons comme résoudre numériquement une EDO avec nos outils, c'est à dire Python ou Maple.

Qu'est-ce qu'une équation différentielle

Observation locale - évolution globale

Une équation différentielle est une porte ouverte sur l'histoire. Elle décrit l'évolution d'un système en permettant de prédire son futur et d'examiner son passé, depuis la naissance (l'état initial) du système. Bien sur, il y a des conditions pour cela : il faut que cette équation ait une solution, et en physique, cette solution doit être unique, sous peine de connaitre le comportement du système sans en connaître l'histoire. Ce sont ces caractéristiques qui font toute la puissance des équations différentielles.

Imaginons un système physique quelconque, un pendule, une population humaine, un amas de noyaux radioactifs. Imaginons que son évolution, son histoire, soit déterminée par une fonction mathématique F qui ne dépendrait que du temps. Nous pouvons observer l'évolution du système en laboratoire ou sur le terrain sur un bref laps de temps. Il ne serait pas très pratique ni très réaliste de se fixer comme objectif d'observer un tas de noyaux d'uranium pendant des millions d'annèes ou une population humaine pendant des centaines d'années pour en connaître l'évolution.
De nos systèmes physiques nous n'avons la plupart du temps qu'une connaissance locale, sur une petit échelle de temps par rapport à la vie du système. Est-il possible, à partir de cette connaissance locale d'en tirer une connaissance globale, de connaître l'histoire du système ? Est-il possible passer du local, du petit ou très petit dans l'échelle du temps du système, au global, c'est à dire à la vie entière du système ?

Cette action de passer du très petit au tout ne vous rappelle-t-elle rien ? C'est ce que vous faites quand vous intégrez une fonction ! A partir d'un très petit élement, je calcule en intégrant un tout, par exemple l'aire sous la courbe représentative de la fonction entre deux bornes. Le calcul différentiel et le calcul intégral sont frères, et d'ailleurs résoudre une équation différentielle se dit aussi "intégrer" cette équation.

Mais à quoi ressemble une équation différentielle ? Le terme "équation" vous est familier. Vous l'avez rencontré lorsque vous avez égalé un polynôme à une seule variable à zéro, comme par exemple avec le plus simple des polynômes ax + b = 0.
Une équation différentielle est une équation dans laquelle la variable n'est plus un réel (ou un entier) mais une fonction. Plus précisément, c'est une équation qui rassemble une fonction inconnue et ses dérivées. Elle peut ne comporter qu'une dérivé première de la fonction inconnue, on dira alors que c'est une équation différentielle de premier ordre. Si elle comporte aussi une dérivée seconde, on dira qu'elle est du second ordre. Et ainsi de suite.
La plus simple des EDO de premier ordre s'écrira F'(t) = k*F(t), k pouvant être un entier ou un réel non nul. Cela suppose bien sur que F soit dérivable sur le domaine de définition.

Est-ce que vous saisissez la beauté de la chose ? Dans cette équation, on lie le local (la dérivée F'(t)) au global, la fonction F(t). C'est quand même extraordinaire que notre monde soit fait de telle sorte que la connaissance du comportement d'un système pendant un bref instant puisse nous permettre de connaître, en théorie, toute l'histoire de ce système !

Conditions initiales et conditions aux limites

Nous voilà en possession d'une équation qui nous permet de déterminer la trajectoire d'un système dans le temps. Remarquons en passant que nous faisons ici de la physique classique, galiléenne, avec un temps absolu dans un espace absolu. Ce qui nous intéresse, en fait, c'est de savoir où va notre système à partir d'un point de départ et où il arrivera, si par hasard il existe des limites à son évolution. Si l'on étudie la chute d'un corps, on laisse tomber ce corps d'une altitude connue. Autre exemple un peu plus compliqué, qui est régi par une équation aux dérivées partielles : lorsque vous chauffez une barre de fer depuis une de ses extrémités et que vous voulez suivre l'évolution de la température dans la barre, cette évolution dépend des deux paramétres : la température de la source de chauffage et la température de l'autre extrémité.

Ces contraintes qui fixent les solutions d'EDO pour t = 0 et pour t tendant vers la limite du domaine de résolution s'appellent respectivement les conditions initiales et les conditions aux limites. Dans la suite, et pour les équations que nous allons rencontrer dans cette page, nous ne nous préoccuperons que des conditions initiales du système. Une règle d'or (qui se démontre bien sur) à retenir : votre système présente autant de conditions initiales que l'ordre de l'équation différentielle qui le décrit. Une EDO de 1er ordre nécessite une condition initiale, une EDO de 2eme ordre, deux conditions initiales, etc.

Vous verrez plus tard en math que chaque EDO appartient à une classe que l'on définit par ses conditions initiales et aux limites. On appelle ces classes des "problèmes". Il existe par exemple le problème de Cauchy qui définit une EDO et ses conditions initiales, par exemple y'(t) = f(y(t), t) et y(0) = y0. C'est le type de problème que nous rencontrerons dans cette page. Il existe aussi d'autres problèmes plus complexes mélant conditions initiales et conditions aux limites, comme les problèmes de Dirichlet ou Neumann.

Un exemple simple pour un sujet complexe

Pour illustrer ce que je viens de dire, je vais aborder l'évolution d'un système physique que vous connaissez bien. Nous l'avons déjà rencontré dans une autre page de TangenteX. Il s'agit de l'évolution d'une population de noyaux radioactifs. Je vais appeler N(t) la fonction qui indique le nombre de noyaux radioactifs à un instant donné de la vie de la population.

Pour comprendre l'évolution du système sur une courte période de temps, je vais mesurer le nombre de noyaux qui "disparaissent" par unité de temps, cette "disparition" étant caractérisée par la désintégration du noyau, qui sera détectée par un compteur Geiger par exemple. Je vais aussi estimer le nombre de noyaux présents dans mon échantillon, en le pesant et en connaissant la masse molaire du produit radioactif, du radon par exemple.

A l'issue de mes différentes mesures, j'ai pu établir un fait important : le nombre de noyaux que se désintègrent est proportionnel au nombre de noyaux présents dans l'échantillon. La variation de la population pendant un instant qui sera considéré comme bref est proportionnelle à la population présente à cet instant. Reste à traduire cela mathématiquement ...

Appelons \( \Delta N(t) \) la perte de noyaux et \( \Delta t \) l'intervalle de temps pendant lequel j'ai fait ma mesure. La variation de la population par unité de temps s'écrira donc \( \dfrac{\Delta N(t)}{ \Delta t} \). Et je vais faire là ma première approximation : je vais faire tendre l'intervalle de temps vers 0 et je vais chercher la limite de mon rapport. Vous me voyez venir j'imagine : je veux arriver à définir la dérivée de la fonction N(t) qui décrit l'évolution de la population en l'identifiant à la limite quand t tend vers 0 du rapport \( \dfrac{\Delta N(t)}{ \Delta t} \). Sur le plan du formalisme mathématique, c'est assez discutable, car jamais je ne pourrais vraiment faire tendre t vers 0, et encore moins faire des mesures dans ces conditions. Les mathématiciens ont des infiniments petits, les physiciens des très petits !

Ainsi donc, je peux traduire mon observation physique par l'équation \( N'(t) = k*N(t) \) ce que je peux écrire aussi \( \dfrac{N'(t)}{N(t)} = k \), avec N'(t) la dérivée de la fonction N(t). Je vais faire maintenant appel à vos souvenirs sur votre cours d'intégration : intégrons les deux membres de cette équation, on obtient \( \int(\dfrac{N'(t)}{N(t)}) = \int(k) \) soit \( \ln(N(t))) = kt + C \), en ramenant les deux constantes d'intégration à une constante que j'appelle C.

Poursuivons notre calcul : pour faire disparaitre ce logarithme pas très parlant, je vais calculer l'exponentielle des deux membres. Eleves de TS, j'imagine que cela vous est familier. J'obtiens donc \( N(t) = \exp(kt + C) \) ou encore \( N(t) = \exp(kt)\exp(C) \). Que peut bien valoir cette constante ? Pour le savoir, nous allons examiner les conditions initiales, c'est à dire la valeur de N(t) quand t = 0. Cette manipulation est une des pratiques les plus fondamentales de la résolution d'une EDO.
Quand t = 0, j'ai \( N(0) = \exp(0)\exp(C) \) soit \( N(0) = \exp(C) \). Ma constante est le nombre de noyaux qui existaient au début de l'expérience, que je noterai classiquement \(N_0 \), ce qui me donne finalement l'équation \( N(t) = N_0\exp(kt) \).

Muni de cette formule, je peux calculer pour tout instant de l'histoire du système, le nombre de noyaux restants dans une population d'atomes radioactifs. Bien sur, pour ça, il faut que je détermine k, ce qui a été fait expérimentalement de façon très précise (c'est la constante de désintégration...).

Par l'observation du comportement local d'un système, j'ai pu calculer le comportement du système durant toute son histoire: merveilleux non ? Et c'est vrai pour beaucoup de systèmes mécaniques, du moins si on ne va pas trop chercher dans les détails...

Solution existante et unique ?

Qui dit équation dit solution(s) à l'équation ! Vous avez appris au lycée qu'une équation avait en principe autant de solutions que son degré, solutions dans le corps des réels ou des complexes. Une équation du premier degré possède une solution, une équation du seconde degrés en possèdent deux, réelles ou complexes, etc.
Dans les équations polynômiales que vous connaissez, une équation a toujours une ou plusieurs solutions, solutions qui sont des nombres réels ou complexes. Mais qu'en est-il d'une EDO ?

Une différence fondamentale est d'abord que la solution est une fonction et pas un nombre. Une autre est qu'il existe une infinité de solutions pour chaque équation ! Reprenez notre exemple précédent : la solution à notre équation \( N'(t) = k*N(t) \) est de la forme \( N(t) = \exp(kt)\exp(C) \) où apparait une constante C qui est quelconque, soit donc une infinité de solutions ! Très gênant quand on prétend décrire l'histoire d'un système physique grâce à la solution de cette équation ! Cela signifie pratiquement que l'on peut imaginer une infinité d'histoires à notre système !
Il existe une autre situation gênante, que ne l'on trouve pas avec les équations polynômiales : l'inexistence d'une solution ! Une EDO sans solution, c'est à dire un système sans histoire. Physiquement, cela semble impossible, mais il existe des EDO sans solution. On peut simplement imaginer qu'aucun système physique n'est décrit par ces équations. Heureusement pour nous, on démontre en mathématique que le problème de Cauchy, encore lui, c'est à dire les EDO avec conditions initales, possède toujours une solution locale unique (sous certaines conditions qui sont remplies par nos EDO).

En résumé, les EDO que nous rencontrerons dans nos excursions physiques seront de la forme y'(t) = F(y(t),t) avec y(t=0) = y0, c'est à dire que nous connaitrons une condition initiale (ou plusieurs si l'équation est de degré supérieur à 1). Et il existe le merveilleux théorème de Cauchy-Lipschitz qui démontre que ces équations possèdent une solution unique, sous certaines conditions toujours remplies dans nos cas. Fin du questionnement ! Temporairement ...

Quelques exemples d'équations différentielles ordinaires

Dans les pages de TangenteX, vous trouverez beaucoup d'exemples de systèmes physiques décrits par des EDO de premier ordre, mais plus généralement de second ordre. Je vais en rappeler quelques uns ici en vous indiquant les pages dans lesquelles ils sont décrits de façon plus détaillée.

La seconde loi de Newton

Nous pouvons l'écrire sous deux formes : une EDO de premier ordre en écrivant \( \dfrac{d\vec{p}}{dt} = \vec{F} \) avec \( \vec{p} \) la quantité de mouvement, et sous la forme d'une EDO de deuxième ordre en écrivant \( \dfrac{d^2\vec{r}}{dt^2} = \vec{F} \), avec \( \vec{r} \) le vecteur position du point dont on étudie le mouvement. On peut même l'écrire sous la forme d'un système de deux équations de premier ordre : { \( \dfrac{d\vec{v}}{dt} = \vec{F} \) , \( \dfrac{d\vec{r}}{dt} = \vec{v} \) }, ce qui nous sera très utile pour résoudre numériquement cette EDO de second ordre.

Cette équation permet, à partir de la connaissance d'une variation locale du mouvement d'un point sous l'effet d'une force, de connaitre la trajectoire de ce point, c'est à sa position pendant toute l'histoire du mouvement. Laplace, dans son "Essai philosophique sur les probabilités" pensait même que si nous connaissions les conditions initiales (position et vitesse) de toutes les particules de l'univers, nous saurions, du moins "une intelligence", capables de calculer son avenir ! Heureusement (c'est ma façon de voir !), un petit problème vint contrecarrer ce mirage !

Cette équation présente une caractéristique : elle est linéaire. Pratiquement, cela signifie qu'une variation faible des conditions initiales provoque une variation proportionnelle dans les solutions. En fait, les modèles physiques donnent rarement naissance à des EDO linéaires, les équations obtenues étant le plus souvent non linéaires. Il se trouve qu'il est bien plus difficile de résoudre des EDO non linéaires, alors les physiciens numériciens utilisent toute une gamme d'artifices pour linéariser les EDO. Nous verrons un exemple dans le paragraphe suivant. Nous verrons aussi que les EDO non linéaires ont une sensibilité aux conditions initiales qui contrecarrent le rêve de Laplace.

L'oscillation du pendule simple

Restons dans le domaine de la mécanique classique pour étudier le mouvement d'un point massif qui oscille au bout d'un fil rigide de longueur constante. Nous allons nous pencher sur la variation de l'angle que forme le fil avec la verticale. C'est le cas classique, j'aurais très bien pu étudier la variation d'altitude du point ! Restons classiques et appelons \( \theta \) cet angle.

Comment obtenir l'EDO qui décrit ce mouvement ? Nous l'avons déjà fait : voyez ici. Je ne reviendrais pas sur cette étude, sauf pour retenir l'EDO obtenue \( \dfrac{d^2\theta}{dt^2} + \omega^2 \sin(\theta) = 0 \). Nous allons plutôt nous attarder sur la forme de l'équation plutôt que sur les détails. Par exemple, je poserai la constante \( \omega \) égale à 1, pour obtenir finalement l'équation \( \dfrac{d^2\theta}{dt^2} + \sin(\theta) = 0 \).

La grande différence entre cette équation et la précédente, c'est le sinus qui traine dans le membre de gauche. Vous n'avez pas l'habitude de le voir, l'équation que vous rencontrerez le plus souvent au lycée sera plutôt \( \dfrac{d^2\theta}{dt^2} + \theta = 0 \). Pourquoi : parce que la fonction sinus() n'est pas réputée pour sa linéarité ! Et que la solution de l'équation non linéaire (avec le sinus) n'est pas simple du tout ! Comment éliminer cet encombrant sinus : simplement en appliquant l'approximation dite "des petits angles", qui nous permet de d'approximer le sinus de l'angle par la mesure de l'angle, exprimée en radian évidemment. On obtient donc une équation linéaire dont la solution est bien connue. Ironie de l'histoire, la fonction sinus est une solution !

Seulement il y a un petit problème : la solution de cette équation linéaire n'est valable que pour les variations faibles de l'angle, afin de respecter ladite approximation. Pas question d'étudier le comportement du pendule dans tout le domaine de variation de \( \theta \), ce qui est très ennuyeux et nous prive de grandes découvertes... Heureusement, la résolution numérique de l'équation non linéaire est tout à fait possible et nous permet de découvrir les secrets du pendule simple. Pour les redécouvrir, je vous invite ici ici, en Scilab ou là, en Maple.

Pourquoi cette fixation sur les EDO non linéaires ? Elles sont pourtant le résultat de presque tous les modèles physiques un tant soit peu réalistes ! Parce que leur mathématique est compliquée, parce qu'il n'existe pas de solution analytique simple. Pourtant leur étude est passionnante, elle nous fait découvrir le monde du chaos et de la dynamique des systèmes non linéaires. Pourtant, leur résolution numérique est relativement simple et nous possédons de bons outils pour y aller.

Nous avons déjà abordé, de très loin, la dynamique des systèmes non linéaires à travers le calcul du célèbre attracteur de Lorenz, vous savez l'effet papillon,ici.

La charge d'un condensateur

Pour la physique et l'établissement de l'EDO, je vous renvoie à la page traitant le sujet. Nous ne retiendrons que l'EDO, dans sa forme dépouillée, c'est à dire en prenant la constante RC égale à 1. Cette EDO est de la forme \( \dfrac{dV}{dt} + V = Vc \), avec \( Vc \) la tension de charge constante.

C'est une EDO de premier ordre, elle ne fait intervenir qu'une dérivée de première ordre, avec second membre, ici la constante \( Vc \). Nous avons déjà rencontré une équation semblable dans l'étude de la désintégration d'une population de noyaux. Mais dans ce cas, l'équation n'avait pas de second membre. Sur le plan mathématique, cela ne fait pas de grosse différence. Pour résoudre une EDO à second membre, quelque soit son ordre, on cherche une solution générale, en considérant l'équation sans second membre, puis on cherche une solution particulière avec le second membre.

Les méthodes de résolution numérique

Résolution analytique ou numérique

La plupart des EDO que vous rencontrerez sur TangenteX possèdent des solutions analytiques plus ou moins simples. Comme son nom l'indique une solution analytique est donnée par l'analyse. Certaines sont très simples, comme celles faisant appel à des exponentielles. D'autres sont beaucoup plus compliquées, comme les solutions des équations non linéaires comme l'équation du pendule.

Certaines EDO n'ont pas de solutions analytiques connues ou bien les solutions sont tellement compliquées que l'on renonce à les utiliser. Dans ces cas, comme dans les cas où il n'est pas possible, pour des raisons pratiques, d'utiliser les solutions analytiques, on utilise la résolution numérique.

TangenteX étant consacré à la physique numérique, nous allons nous intéresser à la résolution numérique des différentes EDO que vous rencontrerez.

Les méthodes de résolution numérique

La méthode d'Euler

C'est une des méthodes les plus simples d'intégration d'une EDO, qui peut être utilisée dans beaucoup de cas. Elle n'est pas très précise mais son algorithme est facile à assimiler. Cette méthode est décrite dans la page, avec des exemples d'application, en C, pour des EDO de premier ordre.

Les algorithmes que nous verrons dans la suite n'utilisent pas cette méthode, mais je vous recommande quand même de l'étudier de près, car tous les algorithmes plus compliqués que nous serons amené à utiliser sont plus ou moins inspirés de la méthode d'Euler (voir par exemple les méthodes aux différences finies).

La méthode de Runge-Kutta

Il existe plusieurs méthodes de Runge-Kutta, dénommées en fonction de leur précision. Ici, nous utilisons la méthode d'ordre 4 dite "RK4". Nous l'avons étudié dans cette page. Son algorithme est inspiré de la méthode d'Euler et du principe plus général de développement limité d'une fonction, pour lequel je vous renvoie à votre cours d'analyse ou à la page sur les différences finies de TangenteX.

La méthode RK4 est assez robuste, relativement précise et permet d'aborder des calculs sérieux. Elle est utilisée, avec quelques améliorations, dans les librairies de résolution d'EDO que l'on trouve dans Maple, Python, Scilab et d'autres. Attention toutefois, elle montre ses limites dans la résolution des EDO dites "raides", c'est à dire dont la solution risque de varier très fortement dans un intervalle restreint du domaine de résolution. Elle peut aussi être consommatrice de temps de calcul, mais c'est tout relatif avec nos PC modernes !

Résoudre une EDO avec Python

Equation 1er degré sans second membre

Je choisis bien sur l'exemple de la décroissance radioactive, sachant que le programme proposé convient à n'importe quelle EDO. Le principe est toujours le même : définir l'EDO dans une routine particulière pour pouvoir la modifier facilement; définir le domaine de résolution; définir les conditions initiales; appeler la routine de résolution de Python; et enfin tracer la courbe représentative de la fonction solution correspondante aux conditions initiales.

Pour illustrer l'exemple, j'ai choisi de tracer la courbe de décroissance d'une population de noyaux de radon 220. Sa demi-vie est de 55,6 s et sa constante de désintégration est de 0.0125 s-1. Cette dernière est la seule constante physique que j'utiliserai dans mon script Python. Je saisirai la population initiale de noyaux (10 000 dans la courbe tracée ci-dessous). Je rappelle que l'EDO qui décrit le comportement de la population est \( dN(t) = -\lambda N(t) dt \) ou encore \( \dfrac{dN(t)}{dt} = -\lambda N(t) \). Je fixe la durée de l'expérience à 400 s, compte tenu de la demi-vie de 55,6 s.

Voici la courbe représentative de la fonction solution de l'EDO. On reconnait la courbe d'une exponentielle décroissante:

Détaillons un peu le script ODEPython1.py qui résoud l'EDO et trace la courbe. Il débute par l'appel aux différentes librairies qui seront utilisées par le script. On note la fonction odeint du package scipy.integrate, qui contient le code de résolution numérique des EDO. Le package pylab contient lui les fonctions graphiques, entres autres...

from scipy.integrate import odeint

from pylab import *

Puis je définie les paramètres physiques de la simulation, la valeur de \( \lambda \) dans notre cas.

Lambda = 0.0125

Vient ensuite la définition de l'EDO à résoudre, que je place dans la fonction EDO, ce qui me permettra de modifier facilement le programme pour l'adapter à d'autres besoins.

def EDO(x,t):

x_point = -Lambda*x

return x_point

La fonction EDO a deux arguments: x et t. L'argument t n'est pas utilisé par la fonction EDO mais est indispensable à la fonction odeint() dont nous parlerons tout à l'heure. Notre fonction EDO calcule la valeur de \( \dfrac{dx}{dt} \) pour la valeur de t. Et comme l'énonce notre équation différentielle, nous savons que \( \dfrac{dx}{dt} = -\lambda x \), ce que nous écrivons précisement. C'est cette valeur qui est retournée par EDO pour chaque valeur de t.

Nous verrons que cette forme est très générale, et qu'il suffit d'adapter cette fonction pour résoudre une autre équation différentielle.

Le script permet ensuite la saisie de la condition initale de l'équation différentielle, ici le nombre de noyaux dans la population au début de l'expérience. Cette valeur est chargée dans la variable N0

Puis vient la saisie des données temporelles de l'experience : l'instant initial (t0), sa durée et l'incrément temporel. Pour la courbe ci-dessous, t0 = 0, la durée est de 400 s et l'incréement (ou pas) temporel est de 1 s. Ces valeurs varieront en fonction de l'expérience et doivent être adpatées à sa nature. Ces données servent à construire le vecteur "time" de l'expérience, qui détermine les bornes de résolution de l'équation différentielle et le nombre de pas à calculer. Attention à ne pas choisir un incrément trop petit qui chargerait inutilement votre ordinateur, sans forcément ajouter de la précision aux résultats ! Le vecteur "time" est construit avec l'instruction suivante:

time = arange(t0, tmax, pastemps)

Tout est prêt, les paramètres sont initialisés, le vecteur temps chargé, nous allons maintenant résoudre notre EDO. Cela se fait le plus simplement du monde, en appelant la fonction odeint() :

N = odeint(EDO,N0,time)

La syntaxe d'appel est des plus simples : on retrouve le nom de la fonction qui définit notre équation différentielle, EDO, la valeur de la condition initiale N0, et le vecteur time, puisque notre équation différentielle dépend du temps.

Les résultats sont retournés dans le vecteur N, qui est de dimension identique au vecteur time. Chaque élements de N donne la population de noyaux non désintégrés au temps t correspondant.

Je ne peux que vous encourager à lire la documentation Python correspondante à la fonction odeint() pour en découvrir sa richesse. Point particulièrement intéressant : les méthodes utilisées pour résoudre les EDO en fonction de leurs caractéristiques.

Il est possible d'exploiter les résultats sous deux formes : les lister sur l'écran ou l'imprimante ou tracer une courbe. J'ai choisi la seconde solution, plus distrayante et plus parlante. Les dernières instructions servent à initialiser la fenêtre graphique (figure()) légender la courbe (title(), xlabel() et ylabel()) puis la tracer (plot()) et l'afficher (show()).

Vous le constatez, résoudre une EDO de premier ordre sans second membre avec Python, ce n'est vraiment pas sorcier !

Equation 1er degré avec seconde membre

Nous allons maintenant examiner la résolution d'une équation différentielle de premier ordre avec second membre constant, par exemple la charge d'un condensateur. L'équation différentielle est \( \dfrac{dV}{dt} + V = Vc \) ou encore pour reproduire le schéma de résolution de odeint \( \dfrac{dV}{dt} = -V + Vc \).

Partons du script Python précédant et modifions-le pour l'adapter à notre nouveau besoin en créant le script ODEPython1.py. D'abord la fonction EDO qui décrit l'équation différentielle à résoudre. Cette fonction devient :

def EDO(x,t):

x_point = x_point = -x + Vc

return x_point

Vc est la tension de charge du condensateur, que j'ai fixé arbitrairement à 5 V. J'ai également choisi de fixer la constante de temps RC à 1 pour faciliter les calculs, mais vous pouvez tout à fait la modifier en écrivant:

x_point = x_point = (1/RC)(-x + Vc)

Voici la courbe représentative de la fonction solution de l'EDO. On reconnait la courbe d'une exponentielle décroissante:

Cette courbe est obtenue pour une durée d'expérience de 10 s, avec un incrément de temps de 0.1 seconde. La tension initiale est égale à 0 V.

Equation 2eme degré linéaire

Pour illustrer le sujet, penchons nous sur un grand classique de mécanique de TS : le mouvement d'un solide ponctuel sous l'effet de la pesanteur, doté d'une altitude initiale z0 et d'une vitesse initiale v0. le but de la manoeuvre est bien sur de calculer la trajectoire dudit solide.

Commençons par déterminer l'équation différentielle du mouvement. Nous allons appliquer la deuxième loi de Newton, qui nous dit que \( \dfrac{d^2z}{dt^2} = -g \), équation que j'obtiens en projetant sur l'axe Oz de mon référentiel, en remarquant que je choisis de calculer la trajectoire dans le plan Oz. Pour me simplifier la vie et parce que ce n'est pas le but de cette page, je vais négliger la résistance de l'air, bien qu'elle produise une EDO intéressante. Je suppose que cette partie ne vous pose pas de problème ...

J'ai donc une EDO du 2eme ordre, avec 2 conditions initiales, v0 et z0. Comment organiser notre programme pour résoudre cette EDO ?

Comme vous l'avez remarqué précédement, la fonction odeint() ne sait traiter que les EDO de premier ordre. Pour contourner cette limitation apparente, je vais transformer mon EDO de second ordre en un système de deux EDO de premier ordre que odeint() sait traiter. J'écris donc que \( \dfrac{d^2z}{dt^2} = -g \) devient le système différentiel :
\( \left\{ \begin{array}{lr} \dfrac{dv}{dt} = -g\\ \dfrac{dz}{dt} = v \end{array} \right. \)

Voyons pratiquement ce que cela donne en examinant le code de la fonction EDO(Z,t), qui comme dans l'exemple précédent, décrit notre système différentiel :

def EDO(Z,t):

z = Z[0]

v = Z[1]

dz = v

dv = -g

return array([dz, dv])

Nous retrouvons les deux arguments de la fonction EDO, le temps t et ici une variable Z, qui contient les résultats des calculs successifs de la solution à notre équation différentielle. Dans l'exemple précédent, il s'agissait un vecteur simple. Ici, Z est un tableau de deux vecteurs, car nous résolvons deux équations différentielles simultanément.

Vous le constatez dans le code: le premier vecteur, stocké dans la colonne d'indice 0, Z[0], contient les altitudes calculées pour chaque incrément de temps; le second vecteur contient lui la vitesse du solide pour chaque incrément de temps, il est stocké dans la colonne 1, Z[1].

Dans les deux instructions suivantes, vour retrouvez notre système différentiel, tel que nous l'avons défini : dz = v et dv = -g. Le "dt" est implicite, représenté par le paramètre t de la fonction. odeint() s'en débrouillera.

Puis la fonction retourne le résultat du calcul dans un tableau de deux variables dz et dv, que odeint() intègrera dans le tableau de vecteurs Z. Le mot clé "array" indique le retour des valeurs dans un tableau.

Voyons les autres spécificités du script ODE2Python1.py. Nous avons affaire à une EDO de 2eme ordre, qui nécessite donc deux conditions initiales. Elles sont saisies par l'utilisateur et stockées dans un tableau, afin d'être utilisées par la fonction odeint(). C'est l'objet de l'instruction :

C0 = array([z0, v0])

La définition du vecteur temps est identique à l'exemple précédent, ça de moins à réécrire ! L'appel de la fonction odeint() se fait aussi de la même manière que dans l'exemple précédent, ce qui est assez rassurant !

Pour l'affichage, j'utilise aussi le même morceau de code, en changeant les libellés tout de même ! Ah si, une petite modification dans l'instruction de plot de la courbe. Dans l'exemple précédent, j'affichais un vecteur simple. Dans notre cas, je veux afficher un vecteur contenu dans un tableau de vecteurs, plus précisément la variation d'altitude en fonction du temps, dont les données sont contenues dans le premier vecteur du tableau, la colonne 0. D'où l'instruction :

plot(time, Z[:,0])

Voici la courbe solution de l'EDO pour z0 = 10 m, v0 = 1 m.s-1, une durée d'expérience de 1.5 seconde et un incrément de temps de 0.01 seconde.

Vous reconnaissez la classique portion de parabole, qui est la forme de solution attendue pour cette équation différentielle.

Avouez que la résolution numérique d'une EDO de 2eme degré linéaire n'est pas bien compliquée avec Python. Il suffit de bien comprendre le mécanisme de décomposition de l'EDO du 2eme ordre en deux EDO de premier ordre et de savoir le coder.

Equation 2eme degré non linéaire

Nous allons voir un autre exemple de la technique de décomposition d'une ODE de 2eme ordre en système de premier ordre en résolvant numériquement l'équation du pendule simple.

Je vous rappelle l'équation différentielle du pendule, qui nous donne la variation de l'angle du pendule avec la verticale en fonction du temps. Il s'agit de \( \dfrac{d^2\theta}{dt^2} + \omega^2 \sin(\theta) = 0 \) ou encore \( \dfrac{d^2\theta}{dt^2} = -\omega^2 \sin(\theta) \). La variable \( \omega \) représente la pulsation propre du pendule. Pour simplifier les choses, je la poserai égale à 1.

Comme dans l'exemple précédent, je vais décomposer cette équation en deux EDO de premier ordre qui forment le système différentiel :
\( \left\{ \begin{array}{lr} \dfrac{d\dot{\theta}}{dt} = -\sin(\theta)\\ \dfrac{d\theta}{dt} = \dot{\theta} \end{array} \right. \)
Je note, conformément aux habitudes des physiciens, \( \dot{\theta} \) la dérivée de l'angle \( \theta \) par rapport au temps.

Le script ODE2Python2.py qui résoud cette EDO est très semblable au programme précédent. Seuls diffèrent la définition de la fonction ODE() et le nom des variables :

def EDO(Theta,t):

theta = Theta[0]

theta_point = Theta[1]

dtheta = theta_point

dtheta_point = -sin(theta)

return array([dtheta, dtheta_point])

Nous retrouvons le même schéma que dans le script précédent, avec la définition du système différentiel que nous avons déterminé plus haut. De manière générale, vous pouvez adapter ce script pour résoudre toute EDO de ce type, à condition de bien identifier le système différentiel.

Voici la courbe solution de l'EDO pour theta0 = 0.5 rd, theta_point0 = 0 rd.s-1, une durée d'expérience de 10 seconde et un incrément de temps de 0.01 seconde.

Vous reconnaissez la classique sinusoïde, qui est la forme de solution attendue pour cette équation différentielle. Vous noterez toutefois que j'ai ici choisi un angle initial petit, et donc que cette courbe est très proche d'une sinusoïde parfaite. Essayez des valeurs plus grandes pour theta0 et vous verrez une nette déformation de la sinusoïde, puis en allant vers des angles de plus en plus grands, des résultats surprenants, mais que vous saurez interpréter...

Résolution d'une équation en calcul symbolique avec Python

Pour l'instant, nous avons résolu numériquement nos EDO. Mais il est possible de les résoudre littéralement, c'est à dire de demander à Python la forme générale et particulière des solutions. Nous allons utiliser pour cela le module de calcul symbolique que nous avons abordé dans cette page. Je vous renvoie à cette page pour la mise en oeuvre de ce module et la description de ses capacités.

Pour illustrer les capacités de sympy, nous prendrons comme exemple l'EDO \( \dfrac{dN(t)}{dt} + \lambda N(t) = 0 \), avec les conditions initiales \( N(t=0) = N_0\)

Commençons par importer sympy et activer son impression écran, en tapant les commandes suivantes au prompt de votre console Python:

from sympy import *

init_printing()

Puis déclarons les variables et la fonction symboliques :

t,Lambda,N0 = symbols("t,Lambda, N0")

N = Function("N")

Je déclare ensuite mon EDO en respectant le formalisme de sympy:

edo = N(t).diff(t) + Lambda*N(t)

Je vais maintenant chercher la solution générale à notre EDO en utilisant la fonction dsolve() de sympy :

dsolve(edo)

sympy me retourne la solution générale à l'équation sous la forme d'une exponentielle :
\( N(t) = c1 \exp(- \lambda t) \)

Reste maintenant à déterminer la solution particulière pour N(t=0). Cela n'est pas vraiment simple sous sympy. On commence d'abord par lister les conditions initiales dans un dictionnaire Python, que j'ai nommé ici condinit :

condinit = {N(0) : N0}

Puis je demande à sympy de me résoudre l'équation avec les conditions initiales appliquées aux membres de droite et de gauche (c'est la forme générale...) :

solpart = Eq(solgen.lhs.subs(t,0).subs(condinit), solgen.rhs.subs(t,0))

Et sympy me dit gentillement que c1 = N0...

Pas vraiment très pratique, surtout quand on le compare à Maple, qui va suivre ...
Notons aussi que sympy ne sait pas résoudre toutes les EDO, en particulier les EDO non linéaires. Par exemple, si vous tentez de résoudre l'équation non linéaire du pendule, il vous dira "NotImplementedError: solve: Cannot solve omega2*sin(theta(t)) + Derivative(theta(t), t, t)". Pour information, comme nous le verrons ci-dessous, Maple fait ça sans aucun problème ! Bref, le module de résolution des EDO en calcul symbolique de sympy n'est pas encore parfait.

Résoudre une EDO avec Maple

Maple est un outil très différent de Python. Ce dernier est un langage de programmation interprété, qui permet de réaliser des programmes (des scripts) de toute nature, du calcul scientifique à la gestion de facturation. Il s'appuit sur une collection considérable de librairies, écrites en Python ou même en C ou d'autres langages (il en existe en FORTRAN, si si !).

Maple n'est rien de cela: c'est un outil de calcul formel. Un programme spécialisé qui permet d'effectuer des calculs numériques et littéraux. Par exemple, vous pouvez lui demander de vous calculer la primitive de sin(x) et l'intégrale de sin(x) entre 0 et 1 : Maple vous donnera le résultat littéral de la primitive de sin() et la valeur de son intégrale entre 0 et 1. Et ce, pour des fonctions qui peuvent être complexes. Alors bien sur, les programmeurs ont ajouté plein de choses autour de Maple pour qu'il puisse faire autre chose que du calcul, mais ce n'est pas son objet. Il serait difficile, mais pas impossible, de gérer une facturation avec Maple !

Nous allons retrouver cette caractéristique en abordant la résolution d'EDO avec Maple. Je vous propose de reprendre les quatres types d'équations que nous avons résolues avec Python pour les résoudre avec Maple et comparer les efforts à produire dans les deux cas.

Je suppose que vous avez Maple sous la main, et je vous propose de manipuler directement dans Maple.

Equation 1er degré sans second membre

Reprenons notre équation de la décroissance radioactive \( \dfrac{dN}{dt} + \lambda N = 0 \) et voyons comment la résoudre avec Maple. Pour faciliter la comparaison, je choisis la même condition initiale (10 000 noyaux) et la même durée (400 s).

Ouvrez donc Maple dans le worksheet courant et tapez la commande suivante après le prompt Maple :

dsolve({diff(N(t), t)+0.125e-1*N(t), N(0) = 10000},N(t))

Examinons cette commande avant de l'exécuter. Nous utilisons ici la fonction dsolve() de Maple qui permet la résolution des équations différentielles. Il s'agit d'une fonction très puissante, que je vous invite à découvrir dans la documentation Maple ou dans votre cours Maple.

La fonction dsolve() nécessite plusieurs paramètres. Ici, nous définissons d'abord notre équation différentielle par diff(N(t), t)+0.125e-1*N(t). diff(N(t), t) désigne la dérivée première de notre fonction N(t) par rapport à t. Notez que nous n'avons fait que recopier l'énoncé de l'EDO ! Puis nous fixons la condition initiale : N(0) = 10000. Remarquez les accolades qui entourent ces deux paramètres. Enfin, nous nommons notre fonction : N(t).

Appuyez sur Return (ou Enter)

Maple vous retourne immédiatement la fonction solution de l'EDO, soit : s := N(t) = 10000*exp(-(1/80)*t). Notez que Maple indique la solution particulière de l'EDO qui tient compte de la condition initiale, mais pas de données numériques qui nous permettraient de tracer la courbe représentative de cette solution. Il s'agit là d'un calcul formel.

Pour obtenir les valeurs numériques qui nous permettrons de tracer une courbe, il faut faire quelques petites modifications à notre commande. Les voici:

s := dsolve({diff(N(t), t)+0.125e-1*N(t), N(0) = 10000},N(t), numeric)

Nous avons ajouté l'argument numeric dans la fonction dsolve pour lui demander d'exécuter non plus un calcul formel mais plutôt une résolution numérique. Puis j'ai stocké les résultats de la fonction dsolve() dans une variable s, un vecteur, pour pouvoir utiliser ces résultats plus tard.

Introduisez cette commande et appuyez sur Return (ou Enter)

Maple vous retourne : s:=proc(x_rkf45) ... end proc, qui vous indique que le calcul est fait et vous dit même que la méthode numérique employée est de la famille des méthodes Runge-Kutta (RK-Felhberg 4.5 précisément).

Voyons maintenant comment tracer la courbe représentative de cette solution. Entrez les commandes suivantes:

with(plots)

odeplot(s, [t, N(t)], 0 .. 400)

La première commande charge le package graphique de Maple. La seconde demande le tracé de la courbe à partir du vecteur s, pour un temps variant entre 0 et 400 secondes. Et vous obtiendrez la courbe suivante :

Sans surprise, elle est identique à la courbe obtenue avec Python. Vous pouvez dans la fonction odeplot() ajouter des options pour placer vos propres titres et légendes. Pour cela, consultez la documentation Maple.

Vous appréciez la puissance de Maple ? En trois commandes, vous pouvez résoudre et tracer la solution d'une EDO ! Le plus difficile est d'appréhender la syntaxe de Maple et sa richesse. Je ne peux que vous recommander une lecture approfondie de la documentation Maple, qui est disponible en ligne.

Equation 1er degré avec seconde membre

Pour rappel, l'équation différentielle est \( \dfrac{dV}{dt} + V = Vc \), avec V(t=0) = 0 et Vc, tension de charge = 5 V. je veux tracer la courbe de charge entre 0 et 10 secondes.

Je vais utiliser les mêmes fonctions, en les adaptant bien sur ! Dans la fonction dsolve(), j'utilise la définition de mon EDO : diff(V(t), t)+V(t)-5, et la condition initiale V(0) = 0. Dans la fonction odeplot(), je précise les limites du tracé, entre 0 et 10 secondes. Voilà donc les lignes de commandes à rentrer comme précédement. Si vous n'avez pas fermé votre worksheet depuis le dernier exemple, il n'est pas utile de recharger le package plots et donc vous pouvez supprimer la commande with(plots).

s := dsolve({diff(V(t), t)+V(t)-5, V(0) = 0}, V(t),numeric)

with(plots)

odeplot(s, [t, V(t)], 0 .. 10)

Vous obtiendrez la courbe suivante :

Là aussi, vous constatez la même efficacité de Maple.

Equation 2eme degré linéaire

Changeons maintenant d'ordre d'EDO pour examiner comment Maple traite les équations différentielles de 2eme ordre, en reprenant notre exemple \( \dfrac{d^2z}{dt^2} + g = 0 \). Je reprends les mêmes conditions initiales : z0 = 10 m et vz0 = 1 m.s-1.

Traduisons cela en commandes Maple. Je suppose que vous n'avez pas fermé votre worksheet, et donc je ne recharge pas le package plots. Voilà ce que cela donne :

s := dsolve({diff(Z(t), t, t)+9.81, Z(0) = 10, D(Z)(0) = 1}, Z(t), numeric)

odeplot(s, [t, Z(t)], 0 .. 1.5)

La différence est dans la définition de la fonction dsolve(). Notez d'abord la présence de deux fois la variable t dans la fonction diff(): cela définit la dérivée seconde par rapport au temps. Vous observez aussi que j'ai défini deux conditions initiales : l'altitude initiale, pour la condition affectant la dérivée première, et la vitesse initiale pour la condition affectant la dérivée seconde: c'est le D(Z)(0) = 1, D désignant l'opérateur différentiel. Pour le reste, rien de changé.

Après avoir lancé ces deux commandes, vous obtiendrez la courbe suivante :

Pas de surprise, nous obtenons bien la même courbe qu'avec Python, mais avec des efforts moindres. Pas besoin de s'embêter à décomposer notre EDO du 2eme ordre en deux EDO de premier ordre, ce qui perturbe souvent les débutants...

Equation 2eme degré non linéaire

Pour finir, voyons comment se débrouille maple avec une équation de 2eme ordre non linéaire, notre \( \dfrac{d^2\theta}{dt^2} + \sin(\theta) = 0 \). Nous reprenons les mêmes conditions initiales, soit un angle initial de 0.5 rd et une vitesse angulaire initiale nulle.

La définition de notre EDO dans la fonction dsolve() ne présente pas de difficulté, nous retrouvons la même forme que dans l'exemple précédent.

s := dsolve({diff(theta(t), t, t)+sin(theta(t)), theta(0) = .5, (D(theta))(0) = 0}, theta(t), numeric)

J'ai un peu enrichi le tracé de la courbe, en mettant un titre et des légendes sur les deux axes:

odeplot(s, [t, theta(t)], 0 .. 10, title = "Oscillation d'un pendule simple", labels = [Temps, Angle])

Après avoir lancé ces deux commandes, vous obtiendrez la courbe suivante :

Pour le fun, voilà la courbe obtenue avec un angle initial de 3 rd. On est loin de la belle sinusoïde du modèle linéaire du pendule simple:

Les scripts Python

Les scripts Python utilisés pour tracer les courbes et faire les simulations décrites ci-dessus sont contenus dans le fichier d'archives EDOPython.rar. Il contient les scripts suivants :

EDOPython1.py
EDOPython2.py
EDO2Python1.py
EDO2Python2.py

Pour conclure

Les équations différentielles ordinaires méritent toute votre attention ! Ce sont des outils incontournables de la physique et même de bien d'autres sciences. Certaines sont devenues des modèles, comme Lorenz (en physique de l'atmosphère), Lokta-Volterra (biologie et dynamique des populations), Van der Pol (oscillateurs, systèmes physiques et biologiques) pour ne citer que celles qui me passent par la tête à l'instant.

Je ne peux que vous encourager à les étudier en profondeur, tant sur le plan mathématique, que je n'ai pas abordé, que sur le plan numérique, que j'ai à peine abordé. Les outils comme Python ou Maple vous permettent d'expérimenter sans difficulté, et de ne consacrer votre attention qu'à la signification physique des EDO. Vous pouvez aussi utiliser Scilab ou MatLab, qui possèdent les mêmes fonctions de type odeint(), et les mêmes principes d'utilisation.

Nous aborderons un autre type d'équation différentielle, les équations aux dérivées partielles ou EDP sur cette page. Elles sont encore plus importantes en physique ! Savez-vous que l'électromagnétisme se résume en 4 EDP qui tiennent sur un post-it ? La physique des ondes, la dynamique des fluides et d'autres se résument même à une équation (chacune quand même ...).