SOLUTIONS_Laboratoire2.ipynb

Kernel: Python 3

Laboratoire 2 : les équations différentielles d'ordre 2

GCH-MAT217, Automne 2020, V.Charette

L'objectif de ce laboratoire est de regarder un peu comment résoudre des équations différentielles d'ordre 2 avec python.

On a deux grandes catégories d'outils : les outils symboliques, nous permettant de trouver des solutions de type "analytique", et les outils numériques, lorsque les solutions analytiques n'existent pas ou sont au-delà des méthodes habituelles de résolution d'EDOs.

Commençons par les outils symboliques. Plus précisément, nous allons regarder un peu ce qu'on peut faire avec la bibliothèque sympy.

A. Résolution symbolique

Lorsqu'on travaille avec sympy, il faut déclarer les symboles qu'on va utiliser.

J'ai inclus ici quelques commandes, pour montrer un peu comment se comportent les symboles.

from sympy import *

x,t=symbols('x t')

x

$\displaystyle x$

x+1,x+t

(x + 1, t + x)

x.subs(x,1)+1

$\displaystyle 2$

x

$\displaystyle x$

solve(x**2-1)

[-1, 1]

solve(x**2+1)

[-I, I]

Nous avons déjà vu la commande dsolve pour résoudre une équation différentielle d'ordre 1. Ici "y(x).diff(x)" signifie qu'on prend la dérivée une fois par rapport à x. (C'est une syntaxe de type "orienté-objet".)

y=symbols('y',cls=Function)

# résoudre y'-2y=0 
dsolve(y(x).diff(x)-2*y(x),y(x))

$\displaystyle y{\left(x \right)} = C_{1} e^{2 x}$

Remarquez la structure : dsolve(expr, f(x)). Ici, "expr" signifie "expr=0".

Alternativement, on pourrait fournir à dsolve une équation à résoudre. Il faut alors utiliser Eq.

# résoudre y'= 2y 
dsolve(Eq(y(x).diff(x),2*y(x)),y(x))

$\displaystyle y{\left(x \right)} = C_{1} e^{2 x}$

On peut aussi mettre des conditions initiales. (Je ne suis pas certaine que ça marche dans toutes les situations, par contre, par exemple si l'équation à résoudre admet une solution compliquée.)

# résoudre y'-2y=0, y(0)=2 
dsolve(Eq(y(x).diff(x),2*y(x)),y(x),ics={y(0):2})

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

On peut aussi s'en servir pour les équations d'ordre 2. Pour prendre la dérivée seconde, la syntaxe est :

y(x).diff(x,x)

$\displaystyle \frac{d^{2}}{d x^{2}} y{\left(x \right)}$

Par exemple, résolvons : $$\frac{d^2 y}{dx^2}-2\frac{d y}{dx}+3y=0$$

dsolve(y(x).diff(x,x)-2*y(x).diff(x)+3*y(x),y(x))

$\displaystyle y{\left(x \right)} = \left(C_{1} \sin{\left(\sqrt{2} x \right)} + C_{2} \cos{\left(\sqrt{2} x \right)}\right) e^{x}$

Ici aussi, on peut mettre des conditions initiales.

# résoudre y''-2y=0, y(0)=2 , y'(0)=3
dsolve(y(x).diff(x,x)-2*y(x).diff(x)+3*y(x),y(x),ics={y(0):2,y(x).diff(x).subs(x,0):3})

$\displaystyle y{\left(x \right)} = \left(\frac{\sqrt{2} \sin{\left(\sqrt{2} x \right)}}{2} + 2 \cos{\left(\sqrt{2} x \right)}\right) e^{x}$

Exercice 1. Utilisez dsolve pour résoudre les équations différentielles suivantes :

a) $\frac{d^2y}{dx^2}+2\frac{dy}{dx}-y=x$, $y(0)=1$, $y'(0)=2$

b) $\frac{d^2y}{dx^2}+x\frac{dy}{dx}=0$, $y(0)=7$, $y'(0)=5$

c) $\frac{d^2y}{dx^2}+x\frac{dy}{dx}-y=0$

(Remarque : pour c), sympy va tenter de donner une solution sous forme de série.)

dsolve(y(x).diff(x,x)+2*y(x).diff(x)-y(x)-x,y(x),ics={y(0):1,y(x).diff(x).subs(x,0):2})

$\displaystyle y{\left(x \right)} = - x + \left(\frac{3}{2} + \frac{3 \sqrt{2}}{2}\right) e^{x \left(-1 + \sqrt{2}\right)} + \left(\frac{3}{2} - \frac{3 \sqrt{2}}{2}\right) e^{x \left(- \sqrt{2} - 1\right)} - 2$

dsolve(y(x).diff(x,x)+x*y(x).diff(x),y(x),ics={y(0):7,y(x).diff(x).subs(x,0):5})

$\displaystyle y{\left(x \right)} = \frac{5 \sqrt{2} \sqrt{\pi} \operatorname{erf}{\left(\frac{\sqrt{2} x}{2} \right)}}{2} + 7$

dsolve(y(x).diff(x,x)+x*y(x).diff(x)-y(x),y(x))

$\displaystyle y{\left(x \right)} = C_{2} \left(- \frac{x^{4}}{24} + \frac{x^{2}}{2} + 1\right) + C_{1} x + O\left(x^{6}\right)$

Si on veut tracer le graphe du résultat, sympy a une fonction plot.

sol=dsolve(y(x).diff(x,x)-2*y(x).diff(x)+3*y(x),y(x),ics={y(0):2,y(x).diff(x).subs(x,0):3})

plot(sol.rhs,(x,0,5))

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

Exercice 2. Tracez les graphes de chacune des solutions de l'exercice 1a) et 1b).

sol1=dsolve(y(x).diff(x,x)+2*y(x).diff(x)-y(x)-x,y(x),ics={y(0):1,y(x).diff(x).subs(x,0):2})
plot(sol1.rhs,(x,0,5))

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

sol2=dsolve(y(x).diff(x,x)+x*y(x).diff(x),y(x),ics={y(0):7,y(x).diff(x).subs(x,0):5})
plot(sol2.rhs,(x,0,5))

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

B. Résolution numérique

Quand trouver une solution analytique est au-delà de nos forces -- à nous et à sympy -- on doit recourir à une solution numérique.

À la base, les solveurs numériques travaillent avec un système de la forme suivante : $$\frac{dy}{dx}=f(x,y), \quad y(x0)=y0$$ Autrement dit, une équation différentielle du premier ordre avec condition initiale. Grosso modo, le solveur numérique va "suivre la direction tangente", direction donnée par $f(x,y)$, sur de petits intervalles de $x$ : de $x=x_0$ on se rend à la position $x=x_1$ en suivant la direction tangente, on reprend avec la direction tangente pour se déplacer de $x=x_1$ à $x=x_2$, etc.

(Une des difficultés est de décider de la longueur des intervalles de $x$, pour obtenir une solution suffisamment proche de la bonne solution, qui est inconnue.)

Nous on va aussi vouloir tracer le graphe de la solution numérique. Donc on importe les bibliothèques numpy et pyplot. On va aussi importer la commande pour résoudre numériquement un ED à condition initiale.

# je réinitialise x,y,z, au cas où ils serviraient de nouveau...
del x
del y
del z

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-23-bfd7aa83dc6d> in <module>
      2 del x
      3 del y
----> 4 del z

NameError: name 'z' is not defined

import numpy as np
import matplotlib.pyplot as plt

from scipy.integrate import solve_ivp

Afin d'apprivoiser la syntaxe, considérons l'exemple suivant.

def exponential_decay(t, y): return -0.5 * y

sol = solve_ivp(exponential_decay, [0, 10], [2])

plt.plot(sol.t,sol.y[0])

[<matplotlib.lines.Line2D at 0x7fed1b7c1d90>]

la syntaxe pour solve_ivp est : solve_ivp(f(t,y),[t0,tf],[y0]) où :

1- f(t,y) est une fonction "appelable" : par exemple, f peut être une fonction définie par un module ("def")

2- t0 est la valeur initiale de la variable indépendante t

3- tf est la valeur de t où l'on souhaite estimer la valeur de y(tf)

4- y0=y(t0)

La valeur retournée sol contient plusieurs éléments : sol.t retourne les valeurs où la fonction a été approximée et sol.y retourne les valeurs à ces points.

Alternative : solve_ivp(f(t,y),[t0,tf],[y0],dense_output=True) Le dernier argument est optionnel, mais deviendra nécessaire plus loin.

sol.t

array([ 0.        ,  0.11488132,  1.26369452,  3.06074656,  4.81637262,
        6.57504937,  8.33467262, 10.        ])

sol.y

array([[2.        , 1.88835583, 1.0632438 , 0.43316531, 0.18014905,
        0.07480841, 0.03105019, 0.01350782]])

Remarquez que sol.y est en fait une matrice, ayant autant de rangées que la dimension du système à résoudre. Dans le cas d'une seule ED, il y a une seule rangée mais c'est pour cette raison qu'on prend sol.y[0] plutôt que sol.y dans les entrées pour plt.plot.

Si vous trouvez comme moi que ça manque de points, on peut spécifier à quels points de l'intervalle [t0,tf] on doit évaluer.

sol = solve_ivp(exponential_decay, [0, 10], [2], t_eval=np.arange(0, 10, 0.1))

plt.plot(sol.t,sol.y[0])

[<matplotlib.lines.Line2D at 0x7fed1b7b01f0>]

Quoi? solve_ivp prend des systèmes? ben oui...

def lotkavolterra(t, z, a, b, c, d):
    x, y = z
    return [a*x - b*x*y, -c*y + d*x*y]

sol1 = solve_ivp(lotkavolterra, [0, 15], [10, 5], args=(1.5, 1, 3, 1),t_eval=np.arange(0, 10, 0.1))

#t = np.linspace(0, 15, 300)
#z = sol1.sol(t)

#sol = solve_ivp(exponential_decay, [0, 10], [2], t_eval=np.arange(0, 10, 0.1))

plt.plot(sol1.t,sol1.y[0])

#plt.plot(t, z.T)
#plt.xlabel('t')
#plt.legend(['x', 'y'], shadow=True)
#plt.title('Lotka-Volterra System')
plt.show()

Pour nous, l'application qu'on en fera aujourd'hui c'est pour traiter une équation du second ordre : $$\frac{d^2 y}{dx^2} = f(x,y,y')$$ On va transformer cette équation du deuxième ordre en système de DEUX équations du PREMIER ordre.

Pour ce faire, on va remplacer $\frac{dy}{dx}$ par une nouvelle variable, $z$. Alors $$\frac{d^2y}{dx^2}= \frac{dz}{dx}$$ et notre équation devient : $$\begin{align*} \frac{dy}{dx} & = z \\ \frac{dz}{dx}& = f(x,y,z) \end{align*}$$

Considérons par exemple ceci, vu plus haut : $$\frac{d^2 y}{dx^2}-2\frac{dy}{dx}+3y=0$$

Cette équation devient : $$\begin{align*} \frac{dy}{dx} & = z \\ \frac{dz}{dx}& = 2z-3y \end{align*}$$

On va s'inspirer de l'exemple "Lotka-Volterra" pour deviner la syntaxe...

def monsysteme(x, w):
    y,z =w
    return [z,-3*y+2*z]

#y(0)=2 , z(0)=y'(0)=3

sol2 = solve_ivp(monsysteme, [0, 10], [2,3], t_eval=np.arange(0, 10, 0.1))

#x = np.arange(0,10,.1)
#w = sol2.sol(x)

#sol1 = solve_ivp(exponential_decay, [0, 10], [2], t_eval=np.arange(0, 10, 0.1))

plt.plot(sol2.t,sol2.y[0])


#plt.plot(x, w.T) # w.T est la transposée, l'opération est nécessaire pour avoir les valeurs présentées correctement
#plt.xlabel('x')
#plt.legend(['y', 'z=dy/dx'], shadow=True)
#plt.title('Ma solution, et sa dérivée')
plt.show()

Notez que ce qui nous intéresse principalement, c'est le graphe de $y$ et non celui de sa dérivée...

plt.plot(x, w[0]) # Ici on va plutôt aller chercher la première rangée de w, celle contenant les valeurs de y
plt.xlabel('x')
plt.ylabel('y')
plt.title('Ma solution')
plt.show()

Exercice 3. Résolvez numériquement les équations différentielles suivantes. Tracez la solution sur l'intervalle indiqué.

a) $\frac{d^2y}{dx^2}+2\frac{dy}{dx}-y=x$, $y(0)=1$, $y'(0)=2$, $0\leq x\leq 2$

b) $\frac{d^2y}{dx^2}+x\frac{dy}{dx}=0$, $y(0)=7$, $y'(0)=5$, $0\leq x\leq 5$

c) $\frac{d^2y}{dx^2}+x\frac{dy}{dx}-y=x$, $y(0)=2$, $y'(0)=1$, $0\leq x\leq 5$

def equation_a(x, w):
    y,z =w
    return [z,-2*z+y+x]

sola = solve_ivp(equation_a, [0, 2], [1,2], dense_output=True)

x = np.arange(0,2,.1)
w = sola.sol(x)

plt.plot(x, w[0]) 
plt.xlabel('x')
plt.ylabel('y')
plt.title('Ma solution, 3a)')
plt.show()

def equation_b(x, w):
    y,z =w
    return [z,-x*z]

solb = solve_ivp(equation_b, [0, 5], [7,5], dense_output=True)

x = np.arange(0,5,.1)
w = solb.sol(x)

plt.plot(x, w[0]) 
plt.xlabel('x')
plt.ylabel('y')
plt.title('Ma solution, 3b)')
plt.show()

def equation_c(x, w):
    y,z =w
    return [z,-x*z+y+x]

solc = solve_ivp(equation_c, [0, 5], [2,1], dense_output=True)

x = np.arange(0,5,.1)
w = solc.sol(x)

plt.plot(x, w[0]) 
plt.xlabel('x')
plt.ylabel('y')
plt.title('Ma solution, 3c)')
plt.show()