有限差分法收敛到错误的稳态（python）

Question

我正在尝试重新创建文章中的数值结果。特别是，我想对以下 PDE 系统进行数值求解：

在该系统中，Fn 代表猎物种群密度，Cn 代表捕食者种群密度，En 代表捕食者的能量。 (3.57) 方程模拟了当时间 t 位于两个连续自然数（即 n）之间时这些不同群体如何相互作用

当作者对非空间系统进行数值求解时，即删除 (3.57) 方程中的所有扩散项，他们产生以下结果（他们使用的参数位于图像底部）：

当我用 Python 对系统进行数值求解时，我得到以下结果：

我的结果看起来非常相似，但是，我的解决方案的峰值始终是文章中解决方案的两倍左右。这种差异不会通过采用更小的时间步长而消失，所以我认为我在 Python 中实现前向欧拉方法的方式有问题。任何帮助将不胜感激。

import matplotlib.pyplot as plt

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

h=0.01 #step size
T=15000 #number of steps

D=[] #domain

for i in range(0,T+1):
    D.append(i*h)

F=0.05 #prey initial condition
C=2.5 #predator initial condition
E=0 #energy initial condition

Fval=[F] #list that will store all calculated prey values
Cval=[C] #list that will store all calculated predator values
Eval=[E] #list that will store all calculated energy values

for i in range(1,T+1):
    if (i*h)%1==0: #between year dynamics
        Fnew=ρ*F
        Cnew=E+C
        Enew=0
    else: #within year dynamics
        Fnew=(1+h*r1*(1-F))*F-((h*F*C)/(ω+θ1*F+θ2*C))
        Cnew=(1-h*m1)*C
        Enew=((h*β*F*C)/(ω+θ1*F+θ2*C))+(1-h*δ1)*E
        
    Fval.append(Fnew)
    Cval.append(Cnew)
    Eval.append(Enew)

    F=Fnew
    C=Cnew
    E=Enew

#plotting
plt.plot(D, Fval,color='red',linestyle='dotted')
plt.plot(D, Cval,color='blue',linestyle='dashdot')
plt.plot(D, Eval,color='black')
plt.xlabel('time t')
plt.ylabel('$F_n(t)$, $C_n(t)$, $E_n(t)$')
plt.legend(['$F_n(t)$', '$C_n(t)$','$E_n(t)$']) 
plt.xlim(0,50)
plt.ylim(0)
plt.tick_params(axis="x", direction="in")
plt.tick_params(axis="y", direction="in")
plt.title('Spatial homogeneous impulsive solutions')
plt.show()

plt.plot(D, Fval,color='blue')
plt.xlabel('time t')
plt.ylabel('$F_n(t)$')
plt.xlim(140,150)
plt.ylim(0)
plt.tick_params(axis="x", direction="in")
plt.tick_params(axis="y", direction="in")
plt.title('Spatial impulsive solution')
plt.show()

plt.plot(D, Cval,color='blue')
plt.xlabel('time t')
plt.ylabel('$C_n(t)$')
plt.xlim(140,150)
plt.ylim(0)
plt.tick_params(axis="x", direction="in")
plt.tick_params(axis="y", direction="in")
plt.title('Spatial impulsive solution')
plt.show()

plt.plot(D, Eval,color='blue')
plt.xlabel('time t')
plt.ylabel('$E_n(t)$')
plt.xlim(140,150)
plt.ylim(0)
plt.tick_params(axis="x", direction="in")
plt.tick_params(axis="y", direction="in")
plt.title('Spatial impulsive solution')
plt.show()

Answer 1

我相信非线性就是问题所在。您所采用的前向欧拉形式非常适合线性 ODE，但您的第一个 ODE 是相当非线性的。 （我将这些称为 ODE，因为当您消除空间扩散时，这些只是

t

中的耦合方程。）

似乎存在通过

F_n^2

项的非线性阻尼，以及通过

F_n C_n

项的共振阻尼。这两者对变化都非常敏感，并且太敏感而无法通过缩小网格大小来平滑。

我的建议是看看作者是否确实使用前向欧拉（它并不总是合适！）根据经验，我通常会首先尝试使用 RK4（四阶 Runge-Kutta） 来解决这个问题。

RK4 通常非常擅长求解与此类似的非线性方程，前提是您及时使用精细离散化。

祝你好运！