FitzHugh-Nagumo模型实现无预期振荡输出问题求助

问题背景 #

我是生物领域研究者，未接受过微分方程或动力系统相关训练，正在通过实践学习神经科学中的动力系统直观知识。尝试按照维基百科内容实现FitzHugh-Nagumo模型（Hodgkin & Huxley模型的简化版），编写了动态变量v'和w'的计算代码：

def dv(V, W, R, I):
    return V - ((V**3)/3) - W + (R*I)

def dw(V, W, a, b, tau):
    return (V+a-(b*W))/tau

固定参数与维基百科前两幅图一致：

## fixed params
a = 0.7
b = 0.8
tau = 12.5
R = 0.1
I = 0.5

初始条件选在V零cline右侧、W零cline下方，按相图预期应产生振荡尖峰：

## initial conditions
W = -0.5
V = 0

用简单迭代积分循环（未用专业ODE求解器）实现：

## timescale
T = 200
dt = 0.01

# Initializing arrays to save the solution
fV = np.zeros(int(T/dt))
fW = np.zeros(int(T/dt))
fW[0] = W
fV[0] = V

# Integrate numerically via a for loop
for i in range(int((T/dt)-1)):
    i = i+1

    W+=dw(V, W, a, b, tau) * dt
    V+=dv(V, W, R, I) * dt

    fV[i]=V
    fW[i]=W

但实际运行后无振荡行为，怀疑积分方法有误，请求排查。

问题原因与解决方案 #

问题出在变量更新顺序：当前代码先更新了W，再用更新后的W去计算V的变化量。而显式欧拉法要求基于**同一时刻的状态变量（V和W）**计算导数，再同步更新两者，否则会引入额外误差，破坏系统的振荡特性。

修正后的积分循环代码如下：

## timescale
T = 200
dt = 0.01

# Initializing arrays to save the solution
import numpy as np
fV = np.zeros(int(T/dt))
fW = np.zeros(int(T/dt))
fW[0] = W
fV[0] = V

# Integrate numerically via a for loop
for i in range(int((T/dt)-1)):
    # 先计算当前时刻的导数
    dv_val = dv(V, W, R, I)
    dw_val = dw(V, W, a, b, tau)
    
    # 同步更新V和W
    new_V = V + dv_val * dt
    new_W = W + dw_val * dt
    
    # 更新状态变量
    V, W = new_V, new_W
    
    # 保存结果
    fV[i+1] = V
    fW[i+1] = W

关键说明 #

1. 先计算当前时刻的导数：dv_val和dw_val都基于循环开始时的V和W，确保使用同一时刻的系统状态。
2. 同步更新变量：用同一时刻的导数计算新的V和W，避免部分更新后的状态干扰另一变量的计算逻辑。
3. 统一保存结果：确保下一次循环使用的是同步更新后的最新状态。

运行修正后的代码，就能得到预期的振荡尖峰行为。

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