如何在DifferentialEquations.jl中用回调模拟弹簧质量系统的锤击冲量？

核心原理 #

你的方程里的狄拉克δ函数10δ(t-1)代表t=1秒时的10牛秒冲量，物理上冲量会瞬间改变系统的速度（冲量=动量变化，若系统质量m=1，速度变化量Δu'=冲量/m=10）。我们不需要在ODE里直接处理δ函数，而是用PresetTimeCallback在指定时间点修改系统状态。

修正后的完整代码 #

1. 定义二阶ODE

# 二阶ODE：u'' = -2u' -10u（冲量通过回调处理，不在方程中体现δ项）
function hammer_time(du, u, p, t)
    # du 对应速度 u'(t)，u 对应位移 u(t)
    return -2*du - 10*u
end

2. 设置初始条件与问题

begin
    du0 = 0.0  # 初始速度 u'(0) = 0
    u0 = 0.0   # 初始位移 u(0) = 0
    tspan = (0.0, 10.0)
end

prob = SecondOrderODEProblem(hammer_time, du0, u0, tspan)

3. 实现冲量回调

begin
    hit_times = [1.0]
    # 冲量作用：仅修改速度（integrator.u 是 [速度, 位移] 的数组）
    # 质量m=1，冲量10牛秒，所以速度瞬间增加10
    affect!(integrator) = integrator.u[1] += 10.0

    cb = PresetTimeCallback(hit_times, affect!)
    # save_positions=(true,true) 确保保存冲量前后的状态，便于观察突变
    sol = solve(prob, Tsit5(), callback=cb, save_positions=(true, true))
end

为什么你的原代码出错？ #

• 你之前直接修改integrator.u，但integrator.u在二阶ODE问题中是包含速度和位移的数组，直接加法会同时改变两者，不符合冲量仅改变速度的物理规律。
• 速度变化量计算错误：根据冲量公式，Δv=冲量/质量，原方程中u''的系数为1（即质量m=1），所以速度应增加10而非0.5。

验证（可选） #

t<1时系统静止，解为0；t≥1时的解析解为：
$$u(t) = \frac{10}{3}e^{-(t-1)}\sin(3(t-1))$$
你可以通过sol(t)获取数值解，和解析解对比验证正确性。

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