在 OrdinaryDiffEq 时间积分期间通过回调更改问题和向量大小会保留原始解向量

Question

我正在使用

OrdinaryDiffEq

来求解一组常微分方程。在积分过程中，我想扩展解向量并改变问题。为此，我使用

DiscreteCallback

类型的回调，并使用新问题重新初始化积分器：

using OrdinaryDiffEq

function modify_integrator!(integrator)
    # Define the extended ODE.
    function growth!(du, u, p, t)
        println("growth!")
        du[1] = 0.1*u[1]
        du[2] = 0.2*u[2]
        du[3] = 0.3*u[3]
    end
    u0 = [integrator.u[1]; integrator.u[2]; 1.0]
    tspan = (integrator.t, 10.0)

    # Define the modified problem.
    prob = ODEProblem(growth!, u0, tspan, callback=callbacks)

    # Set up the integrator for the modified problem.
    iter_old = integrator.iter
    naccept_old = integrator.stats.naccept
    integrator = init(prob, integrator.alg, dt=integrator.dt, callback=integrator.opts.callback)
    integrator.iter = iter_old
    integrator.stats.naccept = naccept_old
end


# Define the ODE.
function growth!(du, u, p, t)
    du[1] = 0.1*u[1]
    du[2] = 0.2*u[2]
end

# Define the initial conditions and time span.
u0 = [1.0; 1.0]
tspan = (0.0, 10.0)

# Define the callbacks.
condition(u, t, integrator) = t > 9
cb = DiscreteCallback(condition, modify_integrator!, save_positions=(false,false))
callbacks = CallbackSet(cb)

prob = ODEProblem(growth!, u0, tspan, callback=callbacks)
sol = solve(prob, Tsit5())

integrator

以具有 2 个分量的 u 向量开始，最终结果应该是具有 3 个分量的向量。但我的最终结果就好像修改后的积分器被完全忽略并且原始问题正在得到解决。在我的调试尝试中，我可以验证修改后的

growth!

函数确实被调用，但对最终解决方案没有影响。就好像

ODEProblem

以某种方式保留了原来的问题。

函数

reinit!

似乎几乎可以帮助我，但我不知道如何给它修改问题。

Answer 1

这对我有用。我修改了回调函数，以便它针对新问题创建一个新的积分器。然后我覆盖积分器的问题并调整问题的大小。

function modify_integrator!(integrator)
    # Define the extended ODE.
    function growth!(du, u, p, t)
        println("growth!")
        du[1] = 0.1*u[1]
        du[2] = 0.4*u[2]
        du[3] = 0.3*u[3]
    end
    u0 = [integrator.u[1]; integrator.u[2]; 1.0]
    tspan = (integrator.t, 10.0)

    # Define the modified problem.
    prob = ODEProblem(growth!, u0, tspan, callback=callbacks)

    # Set up the integrator for the modified problem.
    iter_old = integrator.iter
    naccept_old = integrator.stats.naccept
    integrator2 = init(prob, integrator.alg, dt=integrator.dt, callback=integrator.opts.callback)
    integrator.sol.prob.f = integrator2.sol.prob.f
    resize!(integrator, 3)
    integrator.u = integrator2.u
    integrator.du = integrator2.du
    integrator.k = integrator2.k
    integrator.cache = integrator2.cache
    integrator.iter = iter_old
    integrator.stats.naccept = naccept_old
end

这可能不是最干净的解决方案，但它有效，并且可以用作更优雅的起点。