Solving system of ODE for updating initial conditions

I am dealing with the following task. I have system of 2nd order ODE which depend on parameter K. I solve this system using DifferentialEquations on time interval [0,t_0], where the value t_0 is defined. I define tspan and timevector as

t0 = 1000.0;
dt = 0.05;
tspan = (0., t0);
timevector = 0:dt:tspan[2]

I also store the solution as

prob = SecondOrderODEProblem(my_ODEs, initial_vel, initial_coord, tspan, K)
sol = solve(prob; saveat = timevector)

My goal is:

  1. I would like to solve my system of 2nd order ODE for different values of K from interval from 0 to K_{\max} with step \delta K. To do this, I create vector,
K_max = 2.0;
K_step = 0.1;
K_values = range(0, K_max, step=K_step)

and I iterate over all items of vector via for loop

  1. But I want to use last values of my solutions, i.e. last elemens of sol.u as new initial conditions for next value of K, so it seems that I should store last elements of sol.u in memory and then rewrite initial_vel and initial_coord

For me it seems that more efficient way to do this exist. I have tried to think about callbacks, but I am not sure. Could anyone give advice?

So basically you are solving a system of ODEs where the value of some parameter K changes over time at intervals t_0?

The fact that you are doing so discretely with some step \delta K makes me think that you are actually trying to approximate some other problem — what determines t_0 and \delta K? Are you trying to approximate the limit as K changes arbitrarily slowly?

Some more information about the context here might be helpful.

Roughly speaking yes. But let me try explain in details. First of all, parameters t_0 and \delta K are chosen by hand: their values should be “good enough” to observe the corresponding physical effect. It seems that you are right about arbitrary slow changing of K: indeed, my goal is to change adibatically paramater K and compute some property based on solutions of ODE.

The physical effect that I try to capture is hysteresis. You can think that each ODE of my system descibes micromagnet dynamics and I am interested in how the total time-averaged magnetization depends on interaction K between micromagnets.

To investigate this, I set up initial conditions initial_vel and initial_coord that corresponds to zero total magnetization and then perform following procedure:

(pseudocode)
while current_K <= K_max do:
    find sol of ODE with given current_K on interval [0,t0]
    compute total magnetization from sol
    set new initial_vel and initial_coord
    current_K = current_K + delta_K

Here sol corresponds to the solution of ODE system and new initial_vel and initial_coord corresponds to velocities and coordinates obtained from solution of ODE at moment t_0.

In such set up, the total magnetization is zero and jumps to non-zero value at a critical paramater K_c. The value K_c is known analytically, so choice of t0, dt, and K_step is dictated from good enough convergence between analytical K_c and observed jump of total magnetization.

You can look into remake; There you can define new initial conditions and parameter values like so:

remake(prob; p=new_ps, u0=new_u0)

Just be careful of the format of u0 and p. For p I think a vector of pairs works fine. However, for u0 it may need to be flat.

For Example in my project I use something like:

function resolve(ps, prev_sol=nothing;prob, onlylast=true)
    if prev_sol===nothing 
        updated_prob = remake(prob, p = ps)
    else
        updated_prob = remake(prob, p = ps, u0=vcat([val[end,:] for val in values(prev_sol.u)]...))
    end
    # Solution of the ODE system
    if onlylast
        @suppress return solve(updated_prob, TRBDF2(), save_on=false, save_start=false)   
    else
        return solve(updated_prob, TRBDF2())
    end
end

to “wrap” the previous sol into the new problem. However I’m using MoL/MTK, rather than DE.jl directly so the syntax may be “less forgiving” for me.

Also look into the save_on=false, save_start=false parameters if you only need the last step, may give a performance boost if you model is large enough.

How does it differ from direct substitution of new initial conditions and parameters?

It’s not as costly computationally AFAIK. At least it is for ModelingToolkit.

Then why not change K smoothly in time, rather than in jumps? e.g. let K(t) = \frac{t}{T} K_f and integrate the ODE from t=0 to t=T, looking at larger and larger T? Smoothness is almost always more efficient for ODE solvers, and such a K(t) is a lot easier for you to code as well.

A good way to look at the adiabatic limit T \to \infty of the solution u(T) might be to do Richardson extrapolation via Richardson.jl, e.g. do Richardson.extrapolate(u, T₀, x0=Inf, contract=0.5), which will compute u(T) (where this is a function that calls the ODE solver with K(t) for a given T as above) and repeatedly double T starting at T₀, doing polynomial extrapolation in 1/T at increasingly high orders until the result is converged to the available precision.

(The other thing you could do is to study the steady states of your ODE \frac{du}{dt} = f(u,K) by directly searching for the roots f(u,K) = 0, using a root-finding algorithm rather than an ODE solver. Hysteresis corresponds to the existence of multiple stable roots u for a given K, where stability is determined by the eigenvalues of the Jacobian of f at the roots, i.e. linear-stability analysis.)

1 Like

Will check this and compare with remake. I use the described above scheme only because this approach is well-known in community.