求解 IVP 波动方程 - 不需要的反射

Question

我是一名物理专业的学生，出于好奇，我想创建一个代码来制作波在空间中传播的动画，最终反映在该空间的边界上。

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')

这是完整的代码。

在“求解”部分，我定义了不同的常量、变量、初始条件等。然后我得到输出

y = sol.y[:N]

，这是我希望以动画形式显示的方程的解。

我还添加了一个名为“波浪反射”的部分，它模拟波浪在 0 米和 50 米处的反射。然后在“动画”部分我得到动画。在我的第一次尝试中效果非常好：

在第二次尝试中，我增加了时间来看看会发生什么。第二次反射没有按我的预期工作（抱歉，我无法包含 gif 文件，因为它太大了）。看起来好像有什么东西干扰了它，并且不会出现通常的 pi 相移，而是根本没有相移。

我调查了这个问题，并有想法制作解决方案的动画，而不考虑反射。在第三次尝试中，我惊讶地发现会发生“反射” - 也就是说，我会看到波像被反射一样返回，但有很大的延迟。这个延迟正是第一个反射波穿过整个空间并到达边界所需的时间 - 所以两者在我的第二次尝试中重叠。

你知道为什么会出现这种情况吗？ scipysolve_ivp 函数是否以我不知道的方式工作，可以解释这种现象？我想解决这个问题，以便我可以创建一个正弦波动画，该动画将在空间的左侧生成并在两个边界上反射，然后随着时间的推移展示驻波的外观。

Answer 1

Jared 是绝对正确的：您不必对结果进行后处理来添加波反射：您只需确保首先具有正确的侧面边界条件即可。

由于在 x=0 和 x=L 时 u = 0，因此您只需在这些点设置 du/dt=0 即可；即添加行

dudt[0] = 0;    dudt[N-1] = 0

将 dudt[] 设置为数组后。

下面是您的完整代码（减去大量注释掉的位）。波浪立即美丽地反射，并具有所需的相变。

其他边界条件，例如边界处 dudx=0 也是可能的。

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()