How to update the initial condition in DifferentialEquations?

General Usage
1

Hello!

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
end

# 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
end

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.

2

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

3

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

4

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

5

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.

6

Yup that’s it.

7

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)
end

Is it correct?

8

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.

9

Just for completeness, an identical question from a previous discussion: https://discourse.julialang.org/t/efficient-way-to-handle-operator-splitting-periodically-updating-solution-mid-solve-in-juliadiffeq/30757/1