Skip to content

Instantly share code, notes, and snippets.

@retifrav
Created Sep 3, 2020
Embed
What would you like to do?
Python matplotlib animation issue on Mac OS
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import math
# define x coordinates
x0, x1 = 0, 1
n = 500
x = np.linspace(x0, x1, n)
dx = x[1] - x[0]
# define time coordinates
t0, t1 = 0, 1
m = 500
t = np.linspace(t0, t1, m)
dt = t[1] - t[0]
# define velocity
uleft, uright = 1.0, 1.0
c = np.linspace(uleft, uright, x.shape[0])
# calculate stability
stab = np.amax(c) * dt / dx
print(f"stability r: {stab}")
# position in x
xs = 0.5
# width of the wavelet
rick = 0.1
xxsq = np.square((x - xs) / rick)
r = (1 - 12 * np.pi * xxsq) * np.exp(-6 * np.pi * xxsq)
# define two subplots
fig, ax = plt.subplots(2)
# plot ricker wavelet
ax[0].plot(x,r)
ax[0].set(
xlabel="x",
ylabel="displacement",
title="initial displacement, Ricker wavelet",
ylim=[-1.0, 1.0]
)
# initial displacement
u0 = r
# make sure boundary conditions are fulfilled
u0[0], u0[-1] = 0, 0
# initial velocity
ut0 = np.zeros(n)
# estimate second space derivative at t = 0
uxx0 = np.concatenate((u0, np.zeros(2))) \
- 2 * np.concatenate((np.zeros(1), u0, np.zeros(1))) \
+ np.concatenate((np.zeros(2), u0))
# estimate the first time step
u1 = u0 + dt * ut0 + dt * dt / 2 / dx / dx * c * c * uxx0[1:-1]
# initialize u with zeros
u = np.zeros(n)
# define matrix to make all spatial derivatives at once
A = np.diag(2 * np.ones(n - 2)) \
+ np.diag(-np.ones(n - 3), -1) \
+ np.diag(-np.ones(n - 3), 1)
A = dt * dt / dx / dx * A
# the free inner points where no boundary conditions are prescribed
free = np.linspace(1, n - 2, n - 2, dtype="int")
# velocity squared
c2 = c * c
c2 = c2[free]
# time loop
tmp0 = u0[free] # first time step
tmp1 = u1[free] # second time step
#prepare the subplot for the displacement
ax[1].set(
ylim=(-1.0, 1.0),
xlabel="x",
ylabel="displacement",
title="explicit fd solution of 1d wave eqn"
)
def animate(it):
global tmp0, tmp1
# predict u in the next time step
u[free] = 2 * tmp1 - tmp0 - c2 * np.dot(A, tmp1)
# update the previous time steps
tmp0 = tmp1
tmp1 = u[free]
line, = ax[1].plot(x, u)
return line,
# total time of animation in ms
totaltime = 1000
tint = totaltime / m
ani = animation.FuncAnimation(
fig,
animate,
np.arange(m),
interval=tint,
blit=True,
repeat=False
)
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment