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
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,
and I iterate over all items of vector via for loop
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.
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.)