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

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.