SOLUTIONSLaboratoire1.ipynb

Kernel: Python 3 (system-wide)

Laboratoire 1 : Premiers pas dans la résolution numérique d'équations différentielles

MAT/GCH217 Automne 2020

Prof : Virginie Charette

Chargés d'exercices : Régis Koch, Kevin Thouin

---

Le but de ce premier labo est de nous familiariser avec l'usage de Python, particulièrement le calcul symbolique (sympy) pour la résolution d'équations différentielles.

Nous nous concentrerons sur deux aspects :

• dessiner le champ de directions associé à une équation différentielle

• résoudre "symboliquement" (vs numériquement) une équation différentielle

---

Commençons par importer quelques bibliothèques :

• matplotlib pour tracer des graphes

• sympy pour pouvoir faire du calcul symbolique : en fait, nous allons plutôt importer des commandes de la bibliothèque sympy, afin d'en alléger l'écriture.

import matplotlib.pyplot as plt
import numpy as np

---

A. Champs de direction

Commençons par tracer le champ de directions d'une équation différentielle. Une commande utile à cette fin est quiver. Exécutez la prochaine cellule pour en savoir plus long sur la syntaxe de cette commande.

? plt.quiver

Comment coder le champ de directions d'une fonction

Supposons qu'on ait l'équation différentielle suivante : $$\frac{dy}{dx}=f(x,y)$$ Rappelons que la pente de la tangente à une courbe $(x,y(x))$ est donnée par $y'$.

Donc le vecteur suivant est tangent à la courbe-solution $(x,y(x))$: $$\vec{v}=\vec{i}+y'(x)\vec{j}$$

Et l'argument pour plt.quiver(X,Y,U,V) sera : X=x, Y=x, U=1, V=f(x,y).

Exemple

Nous allons tracer le champ de directions pour l'équation suivante: $$\frac{dy}{dx}=ay+b$$ où $a,b$ sont des constantes quelconques.

Essayez-le avec différentes valeurs de $a,b$.

x=np.arange(-2,2,.25)
y=np.arange(-2,2,.25)
X,Y=np.meshgrid(x,y)

a=3
b=4

U=1
V=a*Y+b

plt.quiver(X,Y,U,V)
plt.show()

Exercice 1. Tracer le champ de directions de l'équation différentielle suivante :

$$\frac{dy}{dx}=3-\frac{y}{10-x}$$

Choisissez vos bornes pour $x\geq 0,y\geq 0$, de manière à ce qu'on voit bien ce qui se passe.

x=np.arange(0,9.5,.5)
y=np.arange(5,15,.5)
X,Y=np.meshgrid(x,y)

U=1
V=2-Y/(10-X)

plt.quiver(X,Y,U,V)
plt.show()

B. La commande dsolve

La commande dsolve appartient à la bibliothèque sympy. On peut lui donner n'importe quelle équation différentielle pouvant être résolue analytiquement (c'est-à-dire symboliquement, ou à la main) et sympy connaît la solution.

Puisqu'on fait du calcul symbolique, il faut déclarer tous nos symboles et fonctions.

from sympy import Symbol,dsolve,Function,Derivative,Eq,solve

y=Function("y")
x=Symbol("x")

y_=Derivative(y(x),x)

Notons que Derivative retourne une fonction symbolique, tout comme y

La commande dsolve prend une équation différentielle, donnée avec la commande Eq, et retourne une équation $y(x)=$ solution.

Par exemple, utilisons dsolve pour résoudre : $$\frac{dy}{dx}=-5y$$

dsolve(Eq(y_+5*y(x),0),y(x),ics={y(0):1})

$\displaystyle y{\left(x \right)} = e^{- 5 x}$

On utilise .rhs si on veut seulement la valeur de la solution :

dsolve(Eq(y_+5*y(x),0),y(x)).rhs

$\displaystyle C_{1} e^{- 5 x}$

Exercice 2. Trouvez la solution de l'équation différentielle :

$$\frac{dy}{dx}=3-\frac{y}{10-x}$$

dsolve(Eq(y_-3+y(x)/(10-x),0),y(x))

$\displaystyle y{\left(x \right)} = C_{1} x - 10 C_{1} + 3 x \log{\left(x - 10 \right)} - 30 \log{\left(x - 10 \right)}$

C. Conditions initiales

Par défaut, sympy utilisera C1, C2, etc. pour dénoter les constantes arbitraires.

On peut utiliser la commande solve pour trouver la valeur de C1 déterminée par la condition initiale.

Prenons par exemple l'équation différentielle suivante : $$\frac{dy}{dx}+y^2\sin x=0,\quad y(0)=1$$

from sympy import sin
s1=dsolve(Eq(y_+sin(x)*y(x)**2,0),y(x))
s1

$\displaystyle y{\left(x \right)} = - \frac{1}{C_{1} + \cos{\left(x \right)}}$

s1.rhs

$\displaystyle - \frac{1}{C_{1} + \cos{\left(x \right)}}$

Nous utiliserons la commande solve pour trouver C1 satisfaisant la condition initiale $y(0)=1$.

C1=Symbol("C1")

val_init=s1.rhs.subs(x,0)-1
const=solve(val_init,C1)
s1.subs(C1,const[0])

$\displaystyle y{\left(x \right)} = - \frac{1}{\cos{\left(x \right)} - 2}$

const

[-2]

Exercice 3. Trouvez la solution de l'équations différentielle suivante :

$$y'=x^2/y,\quad y(1)=2$$

attention! Il y a deux solutions possibles avant d'imposer la condition initiale!

s1=dsolve(Eq(y_-x**2/y(x),0),y(x),ics={y(1):2})
s1

$\displaystyle y{\left(x \right)} = \frac{\sqrt{6 x^{3} + 30}}{3}$

C1=Symbol("C1")
S1=s1[1]
val_init=S1.rhs.subs(x,1)-2
const=solve(val_init,C1)
S1.subs(C1,const[0])

$\displaystyle y{\left(x \right)} = \frac{\sqrt{6 x^{3} + 30}}{3}$

Exercice 4. Trouvez la solution de l'équation différentielle suivante :

$$\frac{dy}{dx}=3-\frac{y}{10-x},\quad y(0)=1$$

Tracez la courbe solution sur le champ de directions trouvé à l'exercice 1.

s1=dsolve(Eq(y_-3+y(x)/(10-x),0),y(x))
s1

$\displaystyle y{\left(x \right)} = C_{1} x - 10 C_{1} + 3 x \log{\left(x - 10 \right)} - 30 \log{\left(x - 10 \right)}$

10*C1-30*np.log(10)-1

$\displaystyle 10 C_{1} - 70.0775527898214$

x=np.arange(0,9.5,.5)
y=np.arange(0,10,.5)
X,Y=np.meshgrid(x,y)

U=1
V=3-Y/(10-x)

plt.quiver(X,Y,U,V)

plt.plot(x,-7*(x-10)+3*x*np.log(10-x)-30*np.log(10-x))
plt.show()

D. Solutions implicites vs solutions explicites

Nous avons vu plusieurs exemples d'équations différentielles dont la solution est à laisser sous forme implicite.

Prenons par exemple l'équation suivante : $$\frac{dy}{dx}=\frac{y(1-x)}{x(3-y)},\quad y(1)=1$$ Cette équation se résoud par séparation de variables et on a (vérifiez-le!) : $$3\ln(y)-y-\ln(x)+x=0$$

Tracer la courbe-solution se fait avec la commande plot_implicit. Par exemple, pour tracer le cercle $x^2+y^2=3$ :

from sympy import plot_implicit 
u=Symbol("u") 
v=Symbol("v")

plot_implicit(Eq(u**2+v**2,3),(u,-4,4),(v,-4,4))

<sympy.plotting.plot.Plot at 0x7fd8a3045e50>

Exercice 4. Tracez la courbe :

$$3\ln(y)-y-\ln(x)+x=0$$

from sympy import plot_implicit 
from sympy import log
u=Symbol("u") 
v=Symbol("v")

plot_implicit(Eq(3*log(v)-v-log(u)+u,0),(u,0.1,10),(v,0.1,10))

<sympy.plotting.plot.Plot at 0x7fc57bdccf90>

Utiliser dsolve pour des solutions implicites

Maintenant, sympy connaît toutes sortes de fonctions... et aura donc des réponses parfois étonnantes.

Par exemple, si on lui demande de résoudre l'équation : $$\frac{dy}{dx}=\frac{y(1-x)}{x(3-y)}$$

voici ce que dsolve répond :

y=Function("y")
x=Symbol("x")

y_=Derivative(y(x),x)

dsolve(Eq(y_-y(x)*(1-x)/(x*(3-y(x))),0),y(x))

[Eq(y(x), -3*LambertW(-(C1*x*exp(-x))**(1/3)/3)),
 Eq(y(x), -3*LambertW((C1*x*exp(-x))**(1/3)*(1 - sqrt(3)*I)/6)),
 Eq(y(x), -3*LambertW((C1*x*exp(-x))**(1/3)*(1 + sqrt(3)*I)/6))]

La fonction LambertW est définie comme la fonction inverse de $xe^x$. Autrement dit : $$\mbox{LambertW}(xe^x)=x$$