Nous voulons obtenir la solution d'une équation différentielle du second ordre par la méthode d'Euler.

Nous allons utiliser Python.

On se place dans le cas d'une équation différentielle du second ordre.

$$a\cdot y'' + b \cdot y' + c\cdot y = f(t)$$

Nous allons reformuler cette équation sous la forme : $$y'' = \frac{1}{a}\left(f(t) - b \cdot y' - c\cdot y\right)$$

À un moment $t$, on connaît $y(t)$ et $y'(t)$ (au début on connaît $y(0)$ et $y'(0)$) :

• on peut donc, avec la formule donnée juste avant, calculer $y''(t)$
• on dira que $y'(t + dt) = y'(t) + y''(t)\cdot dt$
• on fait comme si $y'(t)$ changeait de façon affine de $y'(t)$ à $y'(t+dt)$. Cela revient à écrire (je vous épargne les détails) : $y(t + dt) = y(t) + \frac{y'(t) + y'(t+dt)}{2}\cdot dt$

On voit que pour $y(t)$, $y'(t)$ et $dt$ donnés, ainsi que $a$, $b$, $c$ et $f$, on peut calculer $y(t+dt)$ et $y'(t+dt)$.

Recopiez et complétez :

def avance_de_dt(yt, ypt, dt, t, f, a, b, c):
    """
    yt: valeur de y(t)
    ypt: valeur de y'(t)
    dt: temps dt
    t: temps en cours
    f: fonction du second membre
    a,b,c: coefficients de l'équation
    Renvoie le triplet t+dt, y(t+dt), y'(t+dt)
    """
    # à vous

On souhaite tracer la courbe de la solution $y$. Pour cela il nous faut calculer les valeurs de $y(t)$ pour tous les $t$ sur un certain intervalle $[0;T]$.

Bien sûr, on ne peut pas calculer pour toutes les valeurs de $t$ car elles sont en nombre infini. On va donc calculer pour $N$ valeurs, avec $N$ assez grand (1000, 10 000…)

Exemple :

• Si on choisit $T = 10$ (disons que le temps est en secondes)
• Choisissons $N = 100$, on découpe donc $T$ en 100 intervalles de 100 ms = 0.1s. On aura donc $dt = 0.1$
• On va constituer une liste de valeurs y.
  Quand on écrira y[0], cela désignera la première valeur de la liste, donc $y(0)$.
  Quand on écrira y[1], il s'agira de la valeur de $y$ au bout de 1 temps $dt$.
  Quand on écrira y[50], il s'agira de la valeur de $y$ pour $t = 50 dt$.

def historique(y0, yp0, f, a, b, c, N, T):

  """
  y0: valeur de y(0)
  yp0: valeur de y'(0)
  f, a, b, c: éléments de l'équation différentielle
  N: nombre d'intervalles considérés
  T: on veut connaître y(t) pour t dans [0;T]
  renvoie la paire de listes t, y
  """
  # initialisation des listes, d'abord pleines de 0 :
  y = [0]*(N+1)
  yp = [0]*(N+1)
  tps = [0]*(N+1)
  # valeurs initiales :
  y[0] = y0
  yp[0] = yp0
  
  dt = # complétez
  for i in range(N): # répète N fois
      tps[i+1], y[i+1], yp[i+1] = avance_de_dt(... # complétez)
  return tps, y

</code>

On veut utiliser nos fonction avec un exemple :

$$0.1 y''+ 0.12 y' + y = 1 ; y(0)=0$$

Je vous propose le code. Pour comparaison, j'ai fait le calcul de la solution exacte.

from math import exp, sin, cos, sqrt
a = # compléter
b = # compléter
c = # compléter
f = lambda t:30 # fonction t->30, une constante donc
T = 8
N = 100

t, y = historique(#compléter)


w = sqrt(9.64)
y_exact = [1-exp(-0.6*x)*(cos(w*x)+0.6/w*sin(w*x)) for x in t]

# graphique
import matplotlib.pyplot as plt
plt.plot(t, y, label="euler")
plt.plot(t, y_exact, label="exact")
plt.show()

Vous pouvez essayer avec une autre fonction. Dans ce cas bien sûr, la solution exacte que j'ai proposée ne conviendra plus.

On peut par exemple essayer $f(t) = \cos(3t)$.

On peut modifier :

...
f = lambda t:cos(3*t)
...
y_entree = [f(x) for x in t]
# graphique
import matplotlib.pyplot as plt
plt.plot(t, y, label="euler")
plt.plot(t, y_entree, label="entrée")
plt.show()