如何在DifferentialEquations.jl中实现积分终止回调求解ODE

问题原因 #

你的回调未正确触发，是因为默认的ContinuousCallback未指定触发方向，积分器可能在单步内直接将h[1]从正值计算为负值，导致在回调触发前就进入了sqrt的无效定义域。

解决方案 #

1. 修正回调的触发方向 #

使用ContinuousCallback时，通过direction=-1指定仅当条件从正值变为负值（即h[1]从正降到0）时触发终止操作，确保在h[1]刚到达0时就停止积分，避免进入负数区域。

2. 修复未定义的参数p #

原代码中ODEProblem的p参数未定义，需移除该参数或赋值（如p=nothing）。

3. 可选：方程内加入数值保护 #

为了双重保险，可在ODE定义中加入判断，当h[1]<=0时强制导数为0，避免意外触发sqrt的定义域错误。

修改后的完整代码 #

function height(dh, h, p, t)
    # 加入数值保护：h<=0时导数设为0
    dh[1] = h[1] > 0 ? -sqrt(h[1]) : 0.0
end

h0 = [14.0]
tspan = (0.0, 10.0)

# 移除未定义的p参数，或设置p=nothing
prob = ODEProblem(height, h0, tspan)

# 定义终止回调：h[1]降至0时停止积分
condition(h, t, integrator) = h[1]
affect!(integrator) = terminate!(integrator)
# direction=-1：仅当条件从正变负时触发
cb = ContinuousCallback(condition, affect!, direction=-1)

# 使用自适应步长算法求解（如Tsit5，默认算法也可）
sol = solve(prob, Tsit5(), callback=cb)

# 查看结果
println("积分终止于t=$(sol.t[end])，此时h=$(sol.u[end])")

关键说明 #

• direction=-1：告诉回调仅在condition(h,t,integrator)的数值从正向下穿越0时触发，精准捕捉h[1]降到0的时刻。
• 方程内的数值保护：即使回调因极端步长未及时触发，也能避免sqrt接收负数参数。
• 使用自适应步长算法（如Tsit5()）：能自动调整步长，减少跨越0点的概率。

内容的提问来源于stack exchange，提问作者G. Church