如何使用 DifferentialEquations.jl 中的回调来模拟锤子对弹簧质量系统中质量的打击？

Question

我正在教微分方程，并使用带有 DifferentialEquations.jl 的 Pluto 笔记本来补充课程。我们目前正在尝试模拟一个弹簧质量系统，该系统在特定时刻受到 10 牛顿秒（锤击）的冲击。我相信我需要使用回调来模拟这个，但不知道如何实现它。

该系统的微分方程为：

u''(t) + 2u'(t) + 10u(t) = 10*delta(t-1)

，其中

10*delta(t-1)

是 1 秒时 10 牛顿秒脉冲的 delta 函数。初始条件是

u(0) = u'(0) = 0

.

到目前为止我所做的是：

function hammer_time(du, u, p, t)
    -2*du - 10*u
end

begin
    du0 = 0.0
    u0 = 0.0
    tspan = (0.0, 10.0)
end

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

begin
    hit_times = [1.0]
    affect!(integrator) = integrator.u += 0.5

    cb = PresetTimeCallback(hit_times, affect!)
    sol = solve(prob, Tsit5(), callback = cb)
end

在最后一个代码块中我知道第二行是错误的，但我不知道应该用什么来模拟锤击脉冲。我假设它改变了系统的能量状态，但我不知道如何建模。

Answer 1

教会，好问题！

据我所知，SecondOrderODEProblem 不能很好地处理回调变异状态 - 它源于构建这个漂亮界面时的一些细节，虽然对您的用例并不重要。

但是，如果我们将二阶颂歌写成一阶颂歌的系统，我们当然可以完全避免这种情况，让

u = [x; v]

，所以

u' = [v; -2v-10x]

然后我们可以使用更常见的 ODEProblem 接口定义我们的系统

using DifferentialEquations, Plots;

#x = u[1], v = u[2]
function hammer_time(u, p, t)
    return [u[2];
            -2*u[2] - 10*u[1]];
end

u0 = [0.0;
      0.0];
tspan = (0.0, 10.0);

condition(u, t, integrator) = (t == 1);
affect!(integrator) = (integrator.u[2] += 10);
cb = DiscreteCallback(condition, affect!);

sol = solve(prob, Tsit5(), callback=cb, tstops=[1.0])
plot(sol, labels=["x" "v"])

这有这样的效果，即在时间 1，我们将速度 (

u[2]

) 增加 10 个单位。

上面的代码产生

我相信你已经找到了回调的文档，但这里又是为了后代。 https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#Using-Callbacks

一切顺利！