有限差分法Python实现问题：PDE系统数值模拟结果不符求助

问题背景 #

正在复现文献中包含猎物（$F_n$）、捕食者（$C_n$）、捕食者能量（$E_n$）的PDE系统数值结果：

• 年内动力学（$n<t<n+1$，$n$为自然数）：由式(3.57)描述种群相互作用
• 年际动力学（$t=n$）：由式(3.4)描述种群更新
• 边界条件：所有变量施加零通量边界条件

采用有限差分法求解，设置扩散系数$D_1=10$时，模拟结果与文献中$D_1=5$的结果相似，怀疑有限差分实现或边界条件有误。

代码实现 #

from matplotlib.animation import FuncAnimation
import numpy as np
import matplotlib.pyplot as plt
import random

ρ=0.55
r1=10
θ1=0.17
θ2=0.1
m1=2.1
β=7
ω=0.01
δ1=0.05
D1=10

k=0.01 #time step
K=5000 #number of time steps

h=0.8 #space step
L=20 #length of space domain
H=int(L/h) #number of space steps

print('r =',(D1*k)/(h**2))

X=np.linspace(0,L,H+1) #space domain

T=np.linspace(0,K*k,K+1) #time domain


F=[] #prey initial conditions
for i in range(0,H+1):
    F.append(0.03+(0.01*(random.random())))

C=[] #predator initial conditions
for i in range(0,H+1):
    C.append(2.4+(0.05*(random.random())))

E=[] #energy initial conditions
for i in range(0,H+1):
    E.append(0.02+(0.05*(random.random())))


Fval=[] #prey values
for i in range(0,K+1):
    Fval.append([])

Cval=[] #predator values
for i in range(0,K+1):
    Cval.append([])

Eval=[] #energy values
for i in range(0,K+1):
    Eval.append([])

Fval[0]=F
Cval[0]=C
Eval[0]=E

for i in range(1,K+1):
    if (i*k)%1==0: #between year dynamics
        Fval[i]=[x*ρ for x in Fval[i-1]]
        Cval[i]=[x+y for x,y in zip(Eval[i-1],Cval[i-1])]
        Eval[i]=[x*0 for x in Eval[i-1]]
        F=Fval[i]
        C=Cval[i]
        E=Eval[i]
        continue
    for j in range(0,H+1): #within year dynamics
        if j==0: #zero flux boundary condition 
            Fnew=(1+k*r1*(1-F[j]))*F[j]-((k*F[j]*C[j])/(ω+θ1*F[j]+θ2*C[j]))+((2*k)/(h**2))*(F[j+1]-F[j])
            Cnew=(1-k*m1)*C[j]+D1*((2*k)/(h**2))*(C[j+1]-C[j])
            Enew=((k*β*F[j]*C[j])/(ω+θ1*F[j]+θ2*C[j]))+(1-k*δ1)*E[j]+D1*((2*k)/(h**2))*(E[j+1]-E[j])
        elif j==H: #zero flux boundary condition
            Fnew=(1+k*r1*(1-F[j]))*F[j]-((k*F[j]*C[j])/(ω+θ1*F[j]+θ2*C[j]))+((2*k)/(h**2))*(F[j-1]-F[j])
            Cnew=(1-k*m1)*C[j]+D1*((2*k)/(h**2))*(C[j-1]-C[j])
            Enew=((k*β*F[j]*C[j])/(ω+θ1*F[j]+θ2*C[j]))+(1-k*δ1)*E[j]+D1*((2*k)/(h**2))*(E[j-1]-E[j])
        else:
            Fnew=(1+k*r1*(1-F[j]))*F[j]-((k*F[j]*C[j])/(ω+θ1*F[j]+θ2*C[j]))+(k/(h**2))*(F[j+1]-2*F[j]+F[j-1])
            Cnew=(1-k*m1)*C[j]+D1*(k/(h**2))*(C[j+1]-2*C[j]+C[j-1])
            Enew=((k*β*F[j]*C[j])/(ω+θ1*F[j]+θ2*C[j]))+(1-k*δ1)*E[j]+D1*(k/(h**2))*(E[j+1]-2*E[j]+E[j-1])
        Fval[i].append(Fnew)
        Cval[i].append(Cnew)
        Eval[i].append(Enew)
    F=Fval[i]
    C=Cval[i]
    E=Eval[i]

#plotting
norm = plt.Normalize(np.min(Cval), np.max(Cval))

ax = plt.axes(projection='3d')
for i in range(0,K+1):
    ax.scatter3D(X,[i*k],Cval[i],c=Cval[i],cmap='inferno',norm=norm);
ax.set_xlabel('space x')
ax.set_ylabel('time t')
ax.set_zlabel('$C_n(x,t)$')
ax.xaxis.set_rotate_label(False)
ax.yaxis.set_rotate_label(False)
ax.zaxis.set_rotate_label(False)
ax.view_init(azim=-120)
plt.show()

关键问题分析 #

1. 零通量边界条件的扩散项实现错误 #

零通量边界条件要求$\frac{\partial u}{\partial x}=0$（$u$为$F,C,E$），即边界外的点值等于边界点值（$u_{-1}=u_0$，$u_{H+1}=u_H$）。代入二阶空间差分公式：
$$\frac{u_{j+1}-2u_j+u_{j-1}}{h^2}$$

• 对于左边界$j=0$，$u_{-1}=u_0$，差分变为$\frac{u_{1}-2u_0+u_0}{h^2}=\frac{u_1 - u_0}{h^2}$
• 对于右边界$j=H$，$u_{H+1}=u_H$，差分变为$\frac{u_H-2u_H+u_{H-1}}{h^2}=\frac{u_{H-1} - u_H}{h^2}$

但代码中，边界处的扩散项多乘了一个因子2：

# 左边界示例
Fnew=... + ((2*k)/(h**2))*(F[j+1]-F[j])
Cnew=... + D1*((2*k)/(h**2))*(C[j+1]-C[j])

正确的写法应该去掉因子2，改为：

# 左边界修正后
Fnew=... + ((k)/(h**2))*(F[j+1]-F[j])
Cnew=... + D1*((k)/(h**2))*(C[j+1]-C[j])

右边界同理，将2*k改为k。这个错误会导致边界处的扩散强度是内部的2倍，抵消了$D_1$参数的变化效果，可能是设置$D_1=10$却得到类似$D_1=5$结果的核心原因。

2. 猎物$F$的扩散系数缺失（若文献中$F$包含扩散项） #

观察代码中$F$的扩散项没有乘任何系数，而$C$和$E$的扩散项乘了$D_1$。需要确认文献中的年内动力学方程是否给$F$也设置了扩散系数：

• 如果文献中$F$的扩散系数为$D_1$，则$F$的扩散项需要补充乘$D_1$
• 如果文献中$F$没有扩散项，则需要删除$F$所有的扩散项代码

3. 年际动力学的时间判断存在浮点数精度风险 #

代码中用(i*k)%1==0判断是否到达年际时间点，由于浮点数精度问题，多次累积后可能出现判断失效（比如$i*k$理论上等于整数，但实际存储为类似1.0000000001或0.9999999999的数值）。建议改为基于整数步长的判断：

# 因为k=0.01，每100步对应时间1单位
if i % 100 == 0:

这样可以完全避免浮点数精度问题。

修正建议 #

1. 先修正边界处的扩散项，删除多余的因子2，重新运行$D_1=5$和$D_1=10$的模拟，对比结果差异
2. 核对文献中$F$的动力学方程，确认是否需要添加扩散系数$D_1$或删除现有扩散项
3. 替换年际时间判断的条件，避免浮点数精度问题
4. 若结果仍不符，可固定初始条件（去掉随机扰动），与文献的初始条件对齐，排除随机因素影响

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