1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
| #!/usr/env python3
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import numpy as np
fig, ax = plt.subplots()
ax.invert_yaxis()
##############################################################################
# Construction des variables d'espace
xmax = 10
dx = 1
ymax = 2
dy = 0.1
X = np.arange(0, xmax, dx)
Y = np.arange(0, ymax, dy)
Vx = np.zeros(shape=(len(Y), len(X)))
Vy = np.zeros(shape=(len(Y), len(X)))
##############################################################################
# Paramètres physiques
vhaut = dx/xmax
eta = 1
dt = 1 #ms
##############################################################################
# Calculs
# Initialisation
Vx[0, :] = vhaut
cross, = plt.plot([], [], linestyle='None', marker='+')
# Fonctions utiles
def step():
xx, yy = np.random.randint(0, len(X)-1), np.random.randint(1, len(Y)-1)
cross.set_data([X[xx]], [Y[yy]])
coefx = eta * dx / (dy**2) * dt*1e-3
coefy = eta * dy / (dx**2) * dt*1e-3
Vx[yy, xx] += coefx * (Vx[yy-1, xx] - 2*Vx[yy,xx] + Vx[yy+1, xx])
Vy[yy, xx] += coefy * (Vy[yy, xx-1] - 2*Vy[yy,xx] + Vy[yy, xx+1])
Vx[0, :] = vhaut
Vx[-1, :] = 0
Vy[0, :] = 0
Vy[0, :] = 0
Vx[:, -1] = Vx[:, 0]
q = plt.quiver(X, Y, Vx, Vy, scale=1)
ttxt = plt.text(0.05, 0.95, r'i = 0', transform=plt.gca().transAxes, va='top')
def animate(tt):
for i in range(80):
step()
ttxt.set_text(r'$i = %s$' % tt)
q.set_UVC(Vx, Vy)
anim = FuncAnimation(fig, animate, interval=dt)
plt.show()
|