How to update the initial condition in DifferentialEquations?


I try to solve a nonlinear differential equation using a split-step method.
The right hand side of my equation has two terms:

\frac{du}{dt} = f(u,t) + g(u,t)

According to the split-step method, I can split every integration step (t_n, t_n+dt) into two (t_n, t_n+dt/2) and (t_n+dt/2,t_n+dt), and then replace the solution of the whole equation by a consecutive solution of the following equations:

\begin{align} & \frac{du}{dt} = f(u,t), \qquad t\in\left(t_n,t_n+\frac{dt}{2}\right)\\ & \frac{du}{dt} = g(u,t), \qquad t\in\left(t_n+\frac{dt}{2},t_n+dt\right) \end{align}

The first equation takes u(t_n) as the initial condition and calculates u(t_n+dt/2). In turn, the second equation takes u(t_n+dt/2) as the initial condition and calculates u(t_n+dt). To solve the first equation I use my own method, but for the second one I want to use DIfferentialEquations. However, in the split-step method I have to change the initial condition of the problem on each step. But redefining the whole problem on each step is too costly. Therefore, I search for a method where I can define the problem once and then update the initial condition on demand.

So far, I think to use the integrators interface. Something like this:

using OrdinaryDiffEq

function g(du, u, p, t)
    @inbounds @. du = u
    return nothing

# Define the problem:
tspan = (0., 10.)
u0 = ones(100)
prob = ODEProblem(g, u0, tspan)
integrator = OrdinaryDiffEq.init(prob, BS3(), dense=false)

# Main loop:
Nt = 100
dt = (tspan[2] - tspan[1]) / Nt
u = copy(u0)   # solution which is updated at each step
for i=1:Nt
    # First half-step of the split-step method:
    @. u = u + 0.01   # dummy example

    # Second half-step of the split-step method:
    @. integrator.u = u   # set the new initial condition?
    # @. integrator.uprev = u   # for some reason, this does not work
    step!(integrator, dt/2, true)
    @. u = integrator.u   # update the solution

However, are there any canonical method to update the initial condition?
Or may be there is some in-place solve method, which can use the original initial condition array u0 and update it with the solution on each step?

Thank you.

Look at Callbacks in the DiffEq manuals, maybe the Discrete Time callback can be adapted for this purpose.

I was thinking of PeriodicCallback but there may be several ways to do what you want.

A DiscreteCallback where the condition is fixed to true will fire after each step, and you can use this to do what you need to do

But what should be the affect function?
Like this?

affect!(integrator) = @. integrator.u = u

where u is the current solution which I want to be the initial condition.

Yup that’s it.

Thank you.

Just the last question. The code with the discrete callback will be equivalent to the code from the first post:

for i=1:Nt
    # modify u
    @. integrator.u = u
    step!(integrator, dt, true)

Is it correct?

No, that’ll actually use the stepping, i.e. be calling the ODE steps. I think you just want to update u and then integrator.t += integrator.dt manually. Then let it do its next step.

Just for completeness, an identical question from a previous discussion: