扩散方程数值求解代码出现IndexError索引错误求助

错误原因 #

报错出现在绘图循环的V[t/dt,:]行。由于浮点数精度误差，t/dt的计算结果可能不是精确的整数（比如0.3 / 0.0005理论值为600，但实际浮点数运算可能得到600.0000000000001），而NumPy数组要求索引必须是整数类型，因此触发IndexError。

解决方案 #

提供两种可行的修复方式：

方式一：强制转换为整数索引 #

直接将t/dt的结果转换为整数，确保索引为整数类型：

# 修改绘图循环部分
for t in plot_times:
    idx = int(t / dt)  # 强制转换为整数
    plt.plot(y, V[idx,:], 'Gray', label='numerical')

方式二：通过时间数组匹配索引（更稳妥） #

先让diffusion函数返回时间数组t，再通过查找最接近目标时间的索引来定位，彻底避免浮点数精度问题：

1. 修改diffusion函数的返回值：

def diffusion(dt,dy,t_max,y_max,viscosity,V0):
    s = viscosity*dt/dy**2 
    y = np.arange(0,y_max+dy,dy) 
    t = np.arange(0,t_max+dt,dt)
    r = len(t)
    c = len(y)
    V = np.zeros([r,c])
    V[:,0] = V0
    for n in range(0,r-1): # time
        for j in range(1,c-1): # space
            V[n+1,j] = V[n,j] + s*(V[n,j-1] - 2*V[n,j] + V[n,j+1]) 
    return y, V, r, s, t  # 新增返回时间数组t

2. 修改绘图部分的代码：

y,V,r,s,t_array = diffusion(dt,dy,t_max,y_max,viscosity,V0)

# plotting:
plt.figure(figsize=(7,5))
plot_times = np.arange(0.2,1.0,0.1)
for t in plot_times:
    # 找到最接近目标时间的索引
    idx = np.argmin(np.abs(t_array - t))
    plt.plot(y,V[idx,:],'Gray',label=f't={t}s')
plt.xlabel('distance from wall (m)',fontsize=12)
plt.ylabel('velocity (m/s)',fontsize=12)
plt.axis([0,y_max,0,V0])
plt.legend()  # 添加图例区分不同时间点

额外提示 #

扩散方程的显式格式要求稳定性条件s = viscosity*dt/dy² ≤ 0.5，当前参数计算得s = 2e-4 * 0.0005 / (0.0005)² = 0.4，满足稳定性要求，无需调整参数。

内容的提问来源于stack exchange，提问作者George.Cito