求解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米处的反射。接着在“动画”部分，我得到了动画。在我的第一次尝试中，它运行得相当不错：

高斯波在中间生成，并在边界上反射。时间设置为T=10秒

在我的第二次尝试中，我增加了时间，想看看会发生什么。第二次反射没有按我预期的那样工作（抱歉，我无法包含gif文件，因为它太大了）。看起来好像有什么东西干扰了它，导致没有通常的π相位偏移，而是完全没有相位偏移。

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

高斯波在中间生成，没有考虑反射。然而，过了一段时间后，我们可以看到波像是被我定义的空间外的物体反射回来一样。这里的时间设置为T=20秒

你知道为什么会这样吗？scipy的solve_ivp函数是否以我不知道的方式工作，能解释这个现象？我想解决这个问题，这样我就可以创建一个正弦波的动画，它会在空间的左侧生成，并在两个边界上反射，然后展示出随时间变化而出现的驻波。

数值解法波动方程初始值问题动画模拟波反射高斯波相位偏移驻波

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

1 个回答

撰写回答