Intégration symbolique et numérique d'une équation différentielle ordinaire avec Mathematica

Ce qui suit peut être considéré comme un mini-cours sur l'intégration d'une équation différentielle.
Une équation différentielle ordinaire est une relation entre une variable x, une fonction inconnue y(x) et ses dérivées successives jusqu'à un certain ordre.
Dans cette page, on va montrer comment intégrer, avec Mathematica, une équation différentielle sur l'exemple simple de l'équation différentielle linéaire d'ordre 1 inhomogène suivante : y' - 4 x y = x.
On pourra ensuite appliquer la méthode utilisée dans des cas plus compliqués où , par exemple, l'équation différentielle n'est pas linéaire.
Dans Mathematica, on écrira :
equadiff = y'[x] - 4 x y[x] == x
Cherchons d'abord une solution symbolique en exécutant la commande DSolve qui comporte 3 arguments : l'équation à intégrer, ici equadiff, la fonction inconnue, ici y, et la variable, ici x.
sols = DSolve[equadiff,y,x]
La réponse de Mathematica semble au premier abord très compliquée à comprendre.
Attention: le texte de cette page s'applique à la version 4.0 de Mathematica. Dans la nouvelle version, 4.1, on a préféré utiliser plus souvent l'écriture alternative d'une fonction pure. Par exemple, la fonction pure #^2 & qui calcule le carré de son argument peut s'écrire Function[{x}, x^2]. De prime abord, cette façon d'écrire semble plus facile à lire. Pourtant, à l'usage, on se rend compte que l'écriture abrégée est extrêmement pratique…
La solution peut se résumer ainsi :
Si l'on remplace y par l'expression donnée par le membre de droite de la règle de substitution y -> ..., alors l'équation différentielle est satisfaite, c'est-à-dire le membre de gauche est identiquement égal au membre de droite.
Le membre de droite de la règle de substitution contient trois symboles #, C[1] et & que nous allons expliquer:
D'abord C[1] est une constante arbitraire qui a pour but de nous rappeler qu'il n'y a pas qu'une solution, mais toute une famille de solutions dépendant d'un paramètre. Pour une équation d'ordre plus élevé, on aurait plusieurs constantes arbitraires C[1], C[2],...etc.
Les symboles # et & sont utilisés pour écrire ce qu'on appelle une fonction pure, ou fonction anonyme, qui transforme son argument # d'après la loi décrite par l'expression à gauche du symbole & .
Ainsi, 3 + #^2 & est la loi qui ajoute à 3 le carré de son argument #, quel que soit cet argument. On peut l'écrire aussi (3 + #^2) & pour bien montrer sur quoi porte la loi.
Comme pour les constantes arbitraires, Mathematica numérote l'argument #1, car on pourra aussi travailler avec des lois portant sur plusieurs arguments, autrement dit des fonctions pures de plusieurs variables.
Comme nous le verrons dans un instant et comme la pratique le montre dans beaucoup d'autres cas, les fonctions pures sont extrêmement utiles et simplifient de beaucoup l'écriture de commandes aussi bien au niveau symbolique, que numérique ou graphique.
Pour vérifier que la solution symbolique obtenue satisfait effectivement l'équation différentielle proposée, il suffit d'écrire :
verif = equadiff /. sols[[1]] //Simplify
La réponse "True" s'interprète comme suit: c'est vrai que si, dans l'équation différentielle proposée, je remplace la fonction y par le membre de droite de la règle de substitution contenue dans sols, alors le membre de gauche de l'équation est identiquement égal au membre de droite.
Parmi la famille des solutions de notre équation différentielle, une seule satisfait la condition (dite condition initiale ) : y(0) = 3. En effet, cette condition nous permet de calculer la valeur de C[1], comme suit :
Extrayons le membre de droite de la solution générale, évaluons-le en 0 et imposons-lui d'être égal à 3. Ceci se traduit par l'équation :
equaalg = sols[[1,1,2]][0] == 3
qui est une équation algébrique pour C[1]:
res = Solve[equaalg,C[1]]
On trouve 13/4 pour la valeur de C[1]. Si l'on décide d'appeler f cette solution particulière, on pourra la définir par :
f[x_]:= Evaluate[sols[[1,1,2]][x] /. (C[1]-> 13/4)]
et dessiner son graphe dans l'intervalle [0,1] en exécutant :
graphesols = Plot[f[x]//Evaluate,{x,0,1}]
Pour nous préparer au cas où Mathematica ne pourrait pas trouver de solution symbolique, essayons d'intégrer l'équation numériquement. La commande NDSolve fait appel à une panoplie d'algorithmes dont l'explication n'a pas sa place ici. Sachons cependant qu'il faut ici préciser non seulement la condition initiale, mais encore l'intervalle dans lequel on aimerait que la solution soit calculée. On a donc à exécuter une commande du type suivant (on la nomme soln pour solution numérique). Cette commande a aussi trois arguments, la liste formée de l'équation différentielle et de la condition initiale, la fonction cherchée y, la variable indépendante x et le domaine de variation de x :
soln = NDSolve[{equadiff,y[0]==3},y,{x,0,1}]
On obtient
{{y -> InterpolatingFunction[{{0., 1.}}, <>]}}
La réponse est donc une fonction d'interpolation, une sorte de fonction virtuelle que nous allons nommer ynum (pour ynumérique) et définir ainsi :
ynum[x_]:= Evaluate[soln[[1,1,2]][x]]
Si l'on essaye d'évaluer ynum en un x quelconque, c'est-à-dire d'exécuter ynum[x], on obtiendra :
InterpolatingFunction[{{0., 1.}}, <>][x]
Par contre, pour une valeur numérique de x, on aura un résultat numérique.
Par exemple
ynum[.3] donne 3.64097 et ynum[1/Sqrt[2]] donne 8.58448
On peut dessiner cette fonction virtuelle, en calculer les dérivées successives, l'intégrer, et même en trouver une primitive!
Bornous-nous ici à dessiner son graphe :
grapheynum = Plot[ynum[x]//Evaluate,{x,0,1}]
Comparons les graphes de la solution symbolique et de la solution numérique :
Show[graphesols,grapheynum]
A en juger par la figure, les deux graphes coincident parfaitement.
Pour en avoir le coeur net, translatons le graphe de f vers le haut, juste assez pour le séparer de celui de ynum.
graphesolstranslate = Plot[f[x]+ .5 //Evaluate,{x,0,1}]
Reprenons la comparaison du graphe de la solution numérique avec le graphe de la solution symbolique translaté cette fois :
Show[graphesolstranslate,grapheynum]
Visuellement, il y a semble-t-il correspondance parfaite entre la solution symbolique et la solution numérique. Et pourtant...
Portons la solution numérique ynum dans l'équation différentielle qui peut s'écrire y'[x]- 4 x y[x] - x == 0.
Plus précisément, calculons le membre de gauche de cette équation, qui devrait donner identiquement zéro et traçons-en le graphe qui devrait se confondre avec l'axe des x. Le résultat graphique vaut plus que mille mots :
Plot[(y'[x]- 4 x y[x] - x /. {y[x]-> ynum[x],y'[x]-> D[ynum[x],x]})//Evaluate, {x,0,1}]
Pour apprécier la qualité de l'approximation de f par ynum, on "calcule" la "distance" entre la fonction sur l'intervalle [0,1] donnée par le membre de gauche de l'équation différentielle et ynum. Une façon de définir cette distance est de calculer l'intégrale de la valeur absolue de la différence entre les deux fonctions. Dans notre cas, cela se traduirait par :
NIntegrate[Abs[(y'[x]- 4 x y[x] - x /. {y[x]-> ynum[x],y'[x]-> D[ynum[x],x]}) - 0]//Evaluate,{x,0,1}]
Après un message d'avertissement, Mathematica donne comme valeur 0.000196714.
Contentons-nous de ce résultat dans le cadre de cette page, non sans avoir dessiné le graphe de la fonction à intégrer :
Plot[Abs[(y'[x]- 4 x y[x] - x /. {y[x]-> ynum[x],y'[x]-> D[ynum[x],x]})- 0]//Evaluate,{x,0,1},PlotRange->All]
Cette dernière figure nous montre qu'une fonction peut être "très proche" d'une solution d'une équation différentielle sans qu'elle soit elle-même une solution de cette équation. L'étude approfondie de ce phénomène occupe une place non négligeable de l'analyse.

Retour à ma page personnelle



©jmt : création le 001206 - dernière modification le 020723