I am a student in physics, and out of curiosity I wanted to create a code that makes an animation of a wave propagating in space, eventually reflecting on the boundaries of that space.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.integrate import solve_ivp
from scipy.fftpack import diff as psdiff
#%% Functions
def wave_eq(t, y, c, L, dx):
"""
Wave equation
Parameters
----------
y : list
Initial conditions.
t : numpy.ndarray
Time array.
c : float
Phase velocity.
L : float
Length of the system.
Returns
-------
dydt : numpy.ndarray
"""
N = len(y)//2
ux = np.gradient(y[:N], dx)
uxx = np.gradient(ux, dx)
dudt = y[N:]
du_tdt = c**2*uxx # wave eq
dydt = np.append(dudt, du_tdt)
return dydt
#%% Solve
## Space, discretized grid
L = 50
dx = 0.01
x = np.arange(0, L, dx)
## Time
T = 20
dt = 0.1
t = np.arange(0, T, dt)
## Parameters
c = 4 # phase velocity
A = 1 # amplitude
## Initial conditions
u0 = A*np.exp(-((x-L/2)/2)**2)
# u0 = 2*np.cosh(x)**(-1)
u_t0 = np.zeros(len(x))
y0 = np.append(u0, u_t0)
## Solve
print("Calculating...")
sol = solve_ivp(wave_eq, t_span=[0, T], t_eval=t, y0=y0, args=(c, L, dx))
N = len(sol.y)//2
y = sol.y[:N]
#%% Wave reflection
dc = int(c*dt/dx) # at each iteration, the wave advances c*dt meters
## Reflection at the right hand side
refl_r = np.zeros((len(y), len(y.T)))
for i in range(len(y.T)):
if i > 0:
refl_r.T[i] += np.roll(refl_r.T[i-1], -dc)
refl_r.T[i][-dc:] -= y.T[i-1][:-dc-1:-1]
## Reflection at the left hand side
refl_l = np.zeros((len(y), len(y.T)))
for i in range(len(y.T)):
if i > 0:
refl_l.T[i] += np.roll(refl_l.T[i-1], dc)
refl_l.T[i][:dc] -= y.T[i-1][dc-1::-1]
# y += refl_r + refl_l
#%% Animation
plt.close('all')
fig, ax = plt.subplots()
line1, = ax.plot(x, y.T[0])
# line2, = ax.plot(x, refl_r.T[0])
# line3, = ax.plot(x, refl_l.T[0])
ax.set_ylim([-A, A])
def animate(i):
line1.set_ydata(y.T[i]) # update the data.
# line2.set_ydata(refl_r.T[i])
# line3.set_ydata(refl_l.T[i])
return line1, #line2, line3,
ani = animation.FuncAnimation(fig, animate, frames=int(T/dt), interval=50, blit=True, repeat=False)
print("Movie")
plt.show()
#%% Save animation - Use one of the following methods to save the animation
# ani.save("movie.gif", fps=20)
# writer = animation.FFMpegWriter(fps=15, metadata=dict(artist='Me'), bitrate=1800)
# ani.save("movie.gif", writer=writer)
# FFwriter = animation.FFMpegWriter()
# ani.save('animation.mp4', writer = FFwriter, fps=10)
# ani.save('movie.mp4', writer='ffmpeg')
Here is the entire code.
In the part "Solve", I define the different constants, variables, inital conditions, etc. Then I get the the output y = sol.y[:N] which is the solution of the equation that I wish to display as an animation.
I have also added a part called "Wave reflection" which simulates de reflection of the wave at 0m and at 50m. Then in the part "Animation" I get the animation. It has worked pretty well in my first try:
In my second attempt I increased the time to see what would happen. The second reflection did not work as I expected (sorry I couldn't include the gif file because it was too large). It looked like something was interfering with it and there would not be the usual pi phase shift, but instead no phase shift at all.
I investigated the problem and I had the idea to make an animation of the solution without taking into account the reflections. In this third attempt, I was surprised to find out that a "reflection" would occur - that is, I would see the waves coming back as if they were reflected but with a great delay. This delay is exactly the time it takes for the first reflected waves to travel accross the whole space and hit the boundaries - so the two overlapped in my second attempt.
Do you know why this happens? Does the scipy solve_ivp function work in a way that I am not aware of, that could explain this phenomenon? I would like to solve this problem so that I can create an animation of sinusoidal wave which would be generated on the left hand side of the space and would be reflected on both boundaries, then showcasing the appearance of a standing wave over time.
Jared is absolutely right: you don't have to post-process your results to add wave reflection: you simply have to make sure that you have the correct side boundary condition in the first place.
Since u = 0 at x=0 and x=L for all time, you simply have to set du/dt=0 at these points; i.e. add the line
dudt[0] = 0; dudt[N-1] = 0
immediately after setting dudt[] as an array.
Your code in full (minus lots of commented-out bits) below. The waves reflect beautifully and immediately, with the desired phase change.
Other boundary conditions, e.g. dudx=0 at boundaries are also possible.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy.integrate import solve_ivp
from scipy.fftpack import diff as psdiff
def wave_eq(t, y, c, L, dx):
N = len(y)//2
ux = np.gradient(y[:N], dx)
uxx = np.gradient(ux, dx)
dudt = y[N:]
dudt[0] = 0; dudt[N-1] = 0 ### BOUNDARY CONDITIONS ###
du_tdt = c**2*uxx # wave eq
dydt = np.append(dudt, du_tdt)
return dydt
## Space, discretized grid
L = 50
dx = 0.01
x = np.arange(0, L, dx)
## Time
T = 20
dt = 0.1
t = np.arange(0, T, dt)
## Parameters
c = 4 # phase velocity
A = 1 # amplitude
## Initial conditions
u0 = A*np.exp(-((x-L/2)/2)**2)
u_t0 = np.zeros(len(x))
y0 = np.append(u0, u_t0)
## Solve
print("Calculating...")
sol = solve_ivp(wave_eq, t_span=[0, T], t_eval=t, y0=y0, args=(c, L, dx))
N = len(sol.y)//2
y = sol.y[:N]
#%% Animation
plt.close('all')
fig, ax = plt.subplots()
line1, = ax.plot(x, y.T[0])
ax.set_ylim([-A, A])
def animate(i):
line1.set_ydata(y.T[i]) # update the data.
return line1, #line2, line3,
ani = animation.FuncAnimation(fig, animate, frames=int(T/dt), interval=50, blit=True, repeat=False)
print("Movie")
plt.show()