import numpy as np
import matplotlib.pyplot as plt
from math import *
def f(y):
return -10*y
def solexacte(t0,y0,t): #t0 temps initial, y0 valeur initiale, t instant t
return y0*e**(-10*(t - t0))
def EulerExplicite(t0, y0, N, tf):
h = tf / N
Y = [y0]
T = [t0]
for i in range (N):
y = Y[-1]
Y.append(y + h * f(y))
T.append(T[-1] + h)
return np.array(T), np.array(Y)
for N in [30,60,100,300]:
T, Y = EulerExplicite(0,1,N,5)
plt.plot(T, Y, label = "Euler Explicite n = "+str(N))
T = np.linspace(0,5,300)
plt.plot(T, solexacte(0,1,T), label = "solution exacte")
plt.legend(loc = "best")
plt.show()
N = 10**(np.arange(1,6))
sup = []
h = []
for n in N:
T, Y = EulerExplicite(0,1,n,5)
sup.append(np.max(np.abs(Y - solexacte(0,1,T))))
h.append(np.max(T[1:] - T[:-1]))
plt.loglog(N, sup)
plt.loglog(N, h)
plt.show()
def eqdiff(U,S,dx):
n, = np.shape(U)
A = 1/(dx**2) * (2 * np.diag(np.ones(n),0) - np.diag(np.ones(n-1),1) - np.diag(np.ones(n-1),-1))
return - A.dot(U) + S
def S(x):
return 100 * np.sin(5 * pi * x)
def U0(x):
return np.sin(pi * x) + 1 / 2 * np.sin(5 * pi * x)
def EulerExpliciteChaleur(t, U0, S, h, Nx, tf):
dx = 1 / (Nx + 1)
X = dx * np.arange(1, Nx + 1)
nbr_pas = int(tf / h)
U = np.zeros((Nx, nbr_pas + 1))
U[:,0] = U0(X)
plt.ylim(-3,3)
line, = plt.plot(X,U0(X))#initialise l'abscisse et l'ordonnée
for i in range (1, nbr_pas + 1):
U[:,i] = U[:,i-1] + h * eqdiff(U[:,i-1],S(X),dx)
line.set_ydata(U[:,i])#mets à jour la donnée de l'ordonnée
plt.pause(0.1)#fais une pause pour suivre l'évolution
plt.show()
EulerExpliciteChaleur(0, U0, S, 0.01, 10, 1)
def S(x):
return 100 * np.sin(5 * pi * x)
def U0(x):
return np.sin(pi * x) + 1 / 2 * np.sin(5 * pi * x)
def EulerImpliciteChaleur(t, U0, S, h, Nx, tf):
dx = 1 / (Nx + 1)
X = dx * np.arange(1, Nx + 1)
A = 1/(dx**2) * (2 * np.diag(np.ones(Nx),0) - np.diag(np.ones(Nx-1),1) - np.diag(np.ones(Nx-1),-1))
nbr_pas = int(tf / h)
U = np.zeros((Nx, nbr_pas + 1))
U[:,0] = U0(X)
plt.ylim(-3,3)
line, = plt.plot(X,U0(X))#initialise l'abscisse et l'ordonnée
for i in range (1, nbr_pas + 1):
U[:,i] = np.linalg.solve(np.eye(Nx) + h * A,U[:,i-1] + h * S(X))
line.set_ydata(U[:,i])#mets à jour la donnée de l'ordonnée
plt.pause(0.1)#fais une pause pour suivre l'évolution
plt.show()
EulerImpliciteChaleur(0, U0, S, 0.01, 10, 1)
On dérive la quantité $H(p,q)$ en fonction du temps, on trouve $$\dfrac{\partial}{\partial t} H(p(t),q(t)) = p'(t) \nabla_p H(p,q) + q'(t) \nabla_q H(p,q) =0$$
Donc l'Hamiltonien est constant au cours du temps.
On peut utiliser le théorème de sortie de tout compact : comme on suppose que les courbes de niveaux de $H$ sont compactes, on sait qu'il n'y a pas d'explosion en temps fini et donc que les solutions maximales sont globales.
Euler Explicite : $$\begin{pmatrix} p_{i+1} \\ q_{i+1} \end{pmatrix} =\begin{pmatrix} p_{i} \\ q_{i} \end{pmatrix} + h \begin{pmatrix} 0 & -2 \\ 2 & 0\end{pmatrix}\begin{pmatrix} p_{i} \\ q_{i} \end{pmatrix} $$
def EulerExpliciteHamiltonien(t0, p0, q0, N, tf):
h = tf / N
P = [p0]
Q = [q0]
T = [t0]
for i in range (N):
p, q = P[-1], Q[-1]
P.append(p - 2 * h * q)
Q.append(q + 2 * h * p)
T.append(T[-1] + h)
return np.array(T), np.array(P), np.array(Q)
T, P, Q = EulerExpliciteHamiltonien(0, 1, 1, 300, 5)
plt.plot(T, P**2 + Q**2)
plt.show()
Euler Implicite : $$\begin{pmatrix} p_{i+1} \\ q_{i+1} \end{pmatrix} =\begin{pmatrix} p_{i} \\ q_{i} \end{pmatrix} + h \begin{pmatrix} 0 & -2 \\ 2 & 0\end{pmatrix}\begin{pmatrix} p_{i+1} \\ q_{i+1} \end{pmatrix} $$ $$\begin{pmatrix} 1 & 2h \\ -2h & 1\end{pmatrix}\begin{pmatrix} p_{i+1} \\ q_{i+1} \end{pmatrix} =\begin{pmatrix} p_{i} \\ q_{i} \end{pmatrix} $$ Et on a $$\begin{pmatrix} 1 & 2h \\ -2h & 1\end{pmatrix}^{-1} = \dfrac{1}{1+4h^2} \begin{pmatrix} 1 & -2h \\ 2h & 1\end{pmatrix}$$
def EulerImpliciteHamiltonien(t0, p0, q0, N, tf):
h = tf / N
coeff = 1 / (1 + 4 * h * h)
P = [p0]
Q = [q0]
T = [t0]
for i in range (N):
p, q = P[-1], Q[-1]
P.append(coeff * (p - 2 * h * q))
Q.append(coeff * (2 * h * p + q))
T.append(T[-1] + h)
return np.array(T), np.array(P), np.array(Q)
T, P, Q = EulerImpliciteHamiltonien(0, 1, 1, 300, 5)
plt.plot(T, P**2 + Q**2)
plt.show()
Euler Symplectique : $$\begin{pmatrix} p_{i+1} \\ q_{i+1} \end{pmatrix} = \begin{pmatrix} p_{i} \\ q_{i} \end{pmatrix} + h \begin{pmatrix} 0 & -2 \\ 2 & 0\end{pmatrix}\begin{pmatrix} p_{i} \\ q_{i+1} \end{pmatrix}$$
def EulerSimplectiqueHamiltonien(t, p0, q0, N, tf):
h = tf / N
coeff = 1 / (1 + 4 * h * h)
P = [p0]
Q = [q0]
T = [t]
for i in range (N):
p, q = P[-1], Q[-1]
Q.append(q + 2 * h * p)
P.append(p - 2 * h * Q[-1])#Q[-1] correspond à q_{i+1} alors que q correspond à q_i
T.append(T[-1] + h)
return np.array(T), np.array(P), np.array(Q)
T, P, Q = EulerSimplectiqueHamiltonien(0, 1, 1, 300, 5)
plt.plot(T, P**2 + Q**2)
plt.show()