在Julia的DifferentialEquations.jl中如何为非规则时间序列设置时间步长？

研究谐振子时，我们原本使用毫秒级规则时间序列w_i驱动微分方程，示例代码如下：

ζ = 1/(4π)                                           # 阻尼比

function oscillator!(du,u,p,t)
    du[1] = u[2]                                     # y'(t) = z(t)
    du[2] = -2*ζ*p(t)*u[2] - p(t)^2*u[1]             # z'(t) = -2ζw(t)z(t) -w(t)^2y(t)
end

y0 = 0.0                                            # 初始位置
z0 = 0.0002                                         # 初始速度
u0 = [y0, z0]                                       # 初始状态向量
tspan = (0.0,10)                                    # 时间区间

dt = 0.001                                          # 时间步长

w = t -> freq[Int(floor(t/dt))+1]                   # 时间序列

prob = ODEProblem(oscillator!,u0,tspan,w)           # 定义ODEProblem

sol = solve(prob,DP5(),adaptive=false,dt=0.001)

现需解决的问题：当w_i为毫秒级非规则时间序列（如下表所示）时，应如何设置时间步长？

可行方案 #

方案1：以突变点为强制停止步（推荐）

非规则时间序列的核心是w(t)在离散时间点发生突变，因此需要让求解器精准捕捉这些节点，同时兼顾计算效率。最优做法是提取所有非规则时间点作为求解器的tstops参数，配合自适应步长使用：

1. 提取时间序列中的时间值为数组ts，对应的w值为数组ws
2. 根据物理场景选择插值方式（阶跃插值适合离散采样值，线性插值适合连续变化场景），定义w(t)的插值函数
3. 在求解时传入tstops=ts，强制求解器在每个突变点停止计算，避免步长跨过突变点
4. 使用自适应步长求解器（如原代码中的DP5()），让求解器在平缓区间自动放大步长提升效率

示例代码：

using DifferentialEquations, Interpolations

# 从非规则序列中提取的时间点（单位：秒）和对应w值
ts = [0.023, 0.125, 0.383, 0.587, 0.590, 0.802, 0.906]
ws = [4.3354, 2.34225, -2.0312, -0.280142, 6.28319, 9.82271, -5.21289]

# 定义w(t)：这里用阶跃插值（每个时间区间内w保持常数），也可替换为linear_interpolation做线性插值
itp = interpolate(ws, SteppedInterpolation(ts))
w(t) = itp(t)

ζ = 1/(4π)                                           # 阻尼比

function oscillator!(du,u,p,t)
    du[1] = u[2]                                     # y'(t) = z(t)
    du[2] = -2*ζ*w(t)*u[2] - w(t)^2*u[1]             # z'(t) = -2ζw(t)z(t) -w(t)^2y(t)
end

y0 = 0.0                                            # 初始位置
z0 = 0.0002                                         # 初始速度
u0 = [y0, z0]                                       # 初始状态向量
tspan = (0.0, 10.0)                                 # 时间区间

prob = ODEProblem(oscillator!, u0, tspan)

# 传入tstops强制捕捉突变点，启用自适应步长
sol = solve(prob, DP5(), tstops=ts, adaptive=true)

方案2：固定步长为最小时间间隔

如果场景必须使用固定步长，需取所有相邻时间点的最小间隔作为dt，确保每个突变点都被覆盖。这种方法计算量较大，仅适合对步长有严格要求的场景：

1. 计算所有相邻时间点的差值，取最小值作为固定步长dt
2. 定义w(t)函数，通过查找t所在的区间返回对应的w值
3. 求解时关闭自适应步长，使用固定dt

示例代码调整：

# 提取时间点和对应w值
ts = [0.023, 0.125, 0.383, 0.587, 0.590, 0.802, 0.906]
ws = [4.3354, 2.34225, -2.0312, -0.280142, 6.28319, 9.82271, -5.21289]

# 计算最小时间间隔作为固定步长
dt_min = minimum(diff(ts))
dt = dt_min  # 也可选用更小的步长如0.001，确保覆盖所有点

# 定义w(t)：查找t所在的时间区间，返回对应w值
function w(t)
    idx = searchsortedlast(ts, t)
    idx == 0 && return ws[1]
    idx > length(ws) && return ws[end]
    return ws[idx]
end

ζ = 1/(4π)                                           # 阻尼比

function oscillator!(du,u,p,t)
    du[1] = u[2]                                     # y'(t) = z(t)
    du[2] = -2*ζ*w(t)*u[2] - w(t)^2*u[1]             # z'(t) = -2ζw(t)z(t) -w(t)^2y(t)
end

y0 = 0.0                                            # 初始位置
z0 = 0.0002                                         # 初始速度
u0 = [y0, z0]                                       # 初始状态向量
tspan = (0.0, 10.0)                                 # 时间区间

prob = ODEProblem(oscillator!, u0, tspan)

# 关闭自适应步长，使用固定dt
sol = solve(prob, DP5(), adaptive=false, dt=dt)

关键注意事项 #

• 插值方式选择：若w(t)是传感器采样的离散阶跃值，用阶跃插值；若为连续变化量，优先用线性插值
• 避免步长跨过突变点：如果不设置tstops，自适应步长可能跳过w(t)突变的时间点，导致求解结果偏离实际
• 效率优先：方案1的自适应步长+强制停止步组合，在保证精度的同时能最大程度提升计算效率

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