• Exécuter un e cellule : Shift+Entrée
  • Palette de commandes : Ctrl+Shift+p
In [12]:
%xmode Plain

import matplotlib.pyplot as plt
import numpy as np
import math

%matplotlib inline
import matplotlib as mpl
mpl.style.use('fivethirtyeight')

from scipy.integrate import odeint

Exception reporting mode: Plain

Méthode d'Euler

In [2]:
def euler(F, y, temps):
    y = np.array(y)
    res = [y]
    for k in range(len(temps)-1):
        t, h = temps[k], temps[k+1]-temps[k]
        y = y + h * np.array(F(y, t))
        res.append(y)
    return np.array(res)

Exemple basique

EquaDiff : $y'=y$

In [3]:
def F(y, t):
    """ Équation différentielle y'=y """
    return y

x = np.linspace(0, 1, 10)

y_exp = np.exp(x)
plt.plot(x, y_exp, label="exp")

y_euler = euler(F, 1, x)
plt.plot(x, y_euler, '+', label='euler')

plt.legend()
Out[3]:
<matplotlib.legend.Legend at 0x7fcf800d0c88>

Pendule amorti

  • EquaDiff : $\theta'' = -sin(\theta)$
  • Vectorialisation : on pose $V = \begin{pmatrix}\theta \\ \theta'\end{pmatrix}$, et on a $\begin{pmatrix}\theta' \\ \theta''\end{pmatrix} = F\left(\begin{pmatrix}\theta \\ \theta'\end{pmatrix}, t\right) = \begin{pmatrix}\theta' \\ -sin(\theta)\end{pmatrix}$
In [4]:
def F(V, t):
    """ Équation différentielle y''=-sin(y) """
    theta, theta_prime = V
    return (theta_prime, -math.sin(theta))
In [5]:
x = np.linspace(0, 10, 100)
res = euler(F, (0, 1), x)
plt.plot(x, res, '+', label='euler')
plt.legend()
Out[5]:
<matplotlib.legend.Legend at 0x7fcf7e014b00>
In [6]:
x = np.linspace(0, 10, 100)
res = odeint(F, (0, 1), x)
plt.plot(x, res, '+', label='euler')
plt.legend()
Out[6]:
<matplotlib.legend.Legend at 0x7fcf7d959e80>
In [7]:
x = np.linspace(0, 10, 100)
res = euler(F, (0, 1), x)
theta, theta_prime = res[:, 0], res[:, 1]
plt.plot(theta, theta_prime)
Out[7]:
[<matplotlib.lines.Line2D at 0x7fcf7d876550>]
In [8]:
x = np.linspace(0, 10, 100)
res = odeint(F, (0, 1), x)
theta, theta_prime = res[:, 0], res[:, 1]
plt.plot(theta, theta_prime)
Out[8]:
[<matplotlib.lines.Line2D at 0x7fcf7d8560f0>]
In [9]:
def euler_pendule(theta, theta_prime, temps):
    res = [(theta, theta_prime)]
    for k in range(len(temps)-1):
        t, h = temps[k], temps[k+1]-temps[k]
        theta = theta + h * theta_prime
        theta_prime = theta_prime + h * (-sin(theta)) # ATTENTION : pas la bonne formule
        res.append((theta, theta_prime))
    return np.array(res)
In [10]:
def euler_pendule(theta, theta_prime, temps):
    res = [(theta, theta_prime)]
    for k in range(len(temps)-1):
        t, h = temps[k], temps[k+1]-temps[k]
        ancien_theta = theta
        theta = theta + h * theta_prime
        theta_prime = theta_prime + h * (-sin(ancien_theta))
        res.append((theta, theta_prime))
    return np.array(res)
In [11]:
def euler_pendule(theta, theta_prime, temps):
    res = [(theta, theta_prime)]
    for k in range(len(temps)-1):
        t, h = temps[k], temps[k+1]-temps[k]
        theta = res[k][0] + h * res[k][1]
        theta_prime = res[k][1] + h * (-sin(res[k][0]))
        res.append((theta, theta_prime))
    return np.array(res)