有限差分法求解波动方程时动画无变化问题排查求助

问题描述 #

使用有限差分法求解波动方程并尝试动画展示结果，但所有时间步的解完全相同，导致动画始终静止。以下是核心求解代码与动画实现代码：

核心求解代码 #

import numpy as np
from scipy import signal

Tension=100 #N
mu = .1 #

v= (Tension/mu)**0.5
L=1 #string length
Nx=200
SimulationTime=.2


x = np.linspace(0, L, Nx+1) # Mesh points in space
 
dx = x[1] - x[0]
Nt = 1000
dt = SimulationTime/Nt
t = np.linspace(0, Nt*dt, Nt+1)
V = v*dt/dx
V2 = V**2

Lambda=.1

u = np.zeros(Nx+1) # t
u_1 = np.zeros(Nx+1) # t-1dt
u_2 = np.zeros(Nx+1) # t-2dt

# 初始条件函数，边界值需为0
I = lambda x: 1e-3*x*(1-x)  
#I = lambda x: 1.e-3*np.sin(2*np.pi*x/Lambda) 
#I = lambda x: 1.e-3*signal.gausspulse(x-Nx*dx/2, fc=10)

for i in range(0,int(Nx)):
     u_1[i] = I(x[i]) # u_1 存储t=0时刻的解

# 计算第一个时间步（t=dt）
n=0  
for i in range(1, Nx):
     u[i] = u_1[i] -0.5*V2*(u_1[i-1] - 2*u_1[i] + u_1[i+1])  

# 存储所有时间步的解
user_data = {}
user_data['x'] = x
user_data['u'] = []

u_2[:] = u_1
u_1[:] = u
for n in range(1, Nt):
     # 计算当前时间步的所有节点值
     for i in range(1, Nx):
          u[i] = -u_2[i]+2*u_1[i]+V2*(u_1[i+1]-2*u_1[i]+u_1[i-1])
     u[0]=0
     u[Nx]=0
     u_2[:] = u_1
     u_1[:] = u

     user_data['u'].append(u)

动画实现代码 #

import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML

# 处理数据
data = user_data['u']
data.pop(0)

# 创建画布和轴
fig, ax = plt.subplots()
x_data = np.arange(len(data[0]))

# 初始绘图
line, = ax.plot(x_data, data[0], label='弦振动')

# 自定义绘图样式
ax.set_xlim(0, len(data[0]) - 1)
#ax.set_ylim(min(min(data)), max(max(data)))
ax.set_xlabel('位置')
ax.set_ylabel('位移')
ax.legend()

# 动画更新函数
def update(frame):
    line.set_ydata(data[frame])
    return line,

# 创建动画
animation = FuncAnimation(fig, update, frames=len(data), interval=1000, blit=True)

# 显示动画
plt.show()
HTML(animation.to_jshtml())

问题原因 #

1. 数组引用赋值导致数据覆盖：求解循环中user_data['u'].append(u)存储的是numpy数组u的引用，而非独立副本。每次循环更新u时，之前存入列表的所有元素都会同步被修改，最终所有时间步的解都变成最后一步的结果，导致动画静止。
2. 初始帧存储缺失：第一个时间步的解未存入user_data['u']，后续data.pop(0)会移除有效数据，进一步导致数据异常。
3. 初始时间步公式符号错误：波动方程显式格式中，初始速度为0时，第一个时间步的正确公式应为u[i] = u_1[i] + 0.5*V2*(u_1[i-1] - 2*u_1[i] + u_1[i+1])，原代码使用减号会导致初始演化方向错误。
4. 动画x轴坐标错误：动画中使用np.arange(len(data[0]))作为x坐标，而非实际空间位置数组x，导致坐标与物理意义不符。

修复方案 #

1. 存储数组副本 #

修改求解循环中的存储语句，使用u.copy()创建独立副本：

user_data['u'].append(u.copy())

2. 补充初始帧存储 #

在进入循环前，先将第一个时间步的解存入列表：

# 存储所有时间步的解
user_data = {}
user_data['x'] = x
user_data['u'] = []
user_data['u'].append(u.copy())  # 存入t=dt的初始帧

u_2[:] = u_1
u_1[:] = u
for n in range(1, Nt):
    # ... 原有计算代码 ...
    user_data['u'].append(u.copy())

同时移除动画代码中的data.pop(0)，避免删除有效数据。

3. 修正初始时间步公式 #

将第一个时间步的计算改为：

# 计算第一个时间步（t=dt）
n=0  
for i in range(1, Nx):
     u[i] = u_1[i] + 0.5*V2*(u_1[i-1] - 2*u_1[i] + u_1[i+1])  

4. 使用实际空间坐标 #

动画中替换x_data为实际的空间数组：

x_data = user_data['x']
ax.set_xlim(0, L)  # 对应实际长度范围

完整修复后的核心求解代码片段 #

# 存储所有时间步的解
user_data = {}
user_data['x'] = x
user_data['u'] = []
user_data['u'].append(u.copy())  # 存入初始时间步

u_2[:] = u_1
u_1[:] = u
for n in range(1, Nt):
     # 计算当前时间步的所有节点值
     for i in range(1, Nx):
          u[i] = -u_2[i]+2*u_1[i]+V2*(u_1[i+1]-2*u_1[i]+u_1[i-1])
     u[0]=0
     u[Nx]=0
     u_2[:] = u_1
     u_1[:] = u

     user_data['u'].append(u.copy())

完整修复后的动画代码片段 #

# 处理数据
data = user_data['u']

# 创建画布和轴
fig, ax = plt.subplots()
x_data = user_data['x']

# 初始绘图
line, = ax.plot(x_data, data[0], label='弦振动')

# 自定义绘图样式
ax.set_xlim(0, L)
ax.set_ylim(-2e-3, 2e-3)  # 根据初始条件设置合理的y范围
ax.set_xlabel('位置')
ax.set_ylabel('位移')
ax.legend()

# 动画更新函数
def update(frame):
    line.set_ydata(data[frame])
    return line,

# 创建动画
animation = FuncAnimation(fig, update, frames=len(data), interval=50, blit=True)  # 减小interval让动画更流畅

# 显示动画
plt.show()
HTML(animation.to_jshtml())

内容的提问来源于stack exchange，提问作者user23546599