Julia延迟微分方程：指定constant_lags时Kuramoto模型不收敛

我正在尝试模拟带时滞的Kuramoto模型，其中时滞由矩阵D定义——D的第ij个元素代表振荡器i与j之间的时滞，振荡器间通过加权连接矩阵C相连。模型公式如下：

$$\frac{d\theta_i}{dt} = \omega_i + K \sum_j \left[ C_{ij} \cdot \sin(\theta_j(t - D_{ij}) - \theta_i(t)) \right]$$

其中$\omega_i$是固有频率，$\theta$是振荡器相位。

我编写的模型函数如下：

function kuramoto_dde!(dtheta, theta, h, p, t)
    omega, C, K, D = p

    for i in eachindex(omega)
        sum_coup = 0.
        # 计算求和项（模型中的第二项）
        for j in eachindex(omega)
            if D[i, j] != 0 # 仅对存在连接的项应用时滞
                thetaj = h(p, t - D[i, j])[j]
            else
                thetaj = theta[j]
            end
            sum_coup += C[i, j] * sin(thetaj - theta[i])
        end
        # 写入完整方程
        dtheta[i] = omega[i] + K * sum_coup
    end
    nothing
end

仿真运行代码：

using MAT, DifferentialEquations, LinearAlgebra
# 读取包含矩阵C和D的MATLAB文件
C = matread(raw"C:\Users\duodenum\Desktop\brain_stuff\pr_basics\data\chaudhuri_sc.mat")["J"]
D = matread(raw"C:\Users\duodenum\Desktop\brain_stuff\pr_basics\data\chaudhuri_D.mat")["D"]

n = size(C)[1]
dt = 0.002
tspan = (0., 300.)
timevector = 0:dt:tspan[2]

# 定义历史函数
h(p, t) = rand(n) * 2. * pi

# 固有频率来自均值为omega_0、方差为dispersion的高斯分布
omega_0 = 10.
dispersion = 9.
omega = 2. * pi .* (omega_0 .+ sqrt(dispersion) .* (randn((n, 1))))

# 随机初始条件
theta_0 = 2. * pi * rand(n)

K = 1000
p = (omega, C, K, D)

prob = DDEProblem(kuramoto_dde!, theta_0, h, tspan, p; constant_lags = D)

alg = MethodOfSteps(Rosenbrock23())
sol = solve(prob, alg; saveat = timevector)
println("SOLUTION COMPLETE")

运行后出现错误：

Warning: dt(0.0) <= dtmin(2.220446049250313e-16) at t=0.0, and step error estimate = 1.0. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase C:\Users\duodenum\.julia\packages\SciMLBase\GmToj\src\integrator_interface.jl:599

如果不指定constant_lags，代码能运行但速度极慢：

prob = DDEProblem(kuramoto_dde!, theta_0, h, tspan, p)

alg = MethodOfSteps(Rosenbrock23())
sol = solve(prob, alg; saveat = timevector)
println("SOLUTION COMPLETE")

我试过用constant_lags = D[:]或constant_lags = D[:]'，但都没用。请问这种情况下，该如何正确指定constant_lags来解决问题？

感谢开发的优秀工具包（以及整个SciML生态系统！）

解决方法 #

问题出在constant_lags的传入方式——它需要的是所有唯一非零时滞值的集合，而不是整个D矩阵。

具体步骤： #

1. 提取D矩阵中所有非零的时滞值，去重得到唯一的时滞列表
2. 将这个列表作为constant_lags传入DDEProblem

修改后的代码片段：

# 提取所有非零时滞并去重
unique_lags = unique(filter(x -> x != 0, D))
# 传入constant_lags
prob = DDEProblem(kuramoto_dde!, theta_0, h, tspan, p; constant_lags = unique_lags)

为什么这样有效？ #

• constant_lags的作用是告诉求解器哪些固定时滞会被用到，这样求解器可以提前预存这些时滞对应的历史值，大幅提升效率
• 直接传入D矩阵会导致求解器把矩阵里的每个元素都当成独立时滞，包括大量重复值和0值，这会干扰求解器的步长计算，引发你遇到的dt=0错误
• 去重后的非零时滞列表才是求解器真正需要的信息，既能避免步长问题，又能让求解器优化历史值的存储和查询，提升运行速度

另外，建议检查D矩阵中的时滞值是否存在极小的数值误差（比如接近0的非零值），如果有的话可以设置一个阈值过滤掉，比如：

unique_lags = unique(filter(x -> x > 1e-8, D))

这样可以避免因为数值精度问题引入不必要的微小滞，进一步稳定求解过程。

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