# Solving system of ODE for updating initial conditions

**URL:** <https://discourse.julialang.org/t/solving-system-of-ode-for-updating-initial-conditions/107784>\
**Category:** New to Julia\
**Tags:** question, package, numerics, differentialequation\
**Created:** [December 18, 2023, 6:44pm UTC](https://discourse.julialang.org/t/solving-system-of-ode-for-updating-initial-conditions/107784 "2023-12-18T18:44:47Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![acubed](https://avatars.discourse-cdn.com/v4/letter/a/b19c9b/32.png) [@acubed](https://discourse.julialang.org/u/acubed)\
**Post date:** [December 18, 2023, 6:44pm UTC](https://discourse.julialang.org/t/solving-system-of-ode-for-updating-initial-conditions/107784/1 "2023-12-18T18:44:47Z")

</div>

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

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

```

I also store the solution as

```julia
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,

```julia
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?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 18, 2023, 7:21pm UTC](https://discourse.julialang.org/t/solving-system-of-ode-for-updating-initial-conditions/107784/2 "2023-12-18T19:21:32Z")

</div>

> [@acubed](#):
>
> 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 KKK, so it seems that I should store last elements of `sol.u` in memory and then rewrite `initial_vel` and `initial_coord`

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.

---

<div class="post-metadata">

**Author:** ![acubed](https://avatars.discourse-cdn.com/v4/letter/a/b19c9b/32.png) [@acubed](https://discourse.julialang.org/u/acubed)\
**Post date:** [December 19, 2023, 8:59am UTC](https://discourse.julialang.org/t/solving-system-of-ode-for-updating-initial-conditions/107784/3 "2023-12-19T08:59:10Z")

</div>

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:

```julia
(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.

---

<div class="post-metadata">

**Author:** ![Qfl3x](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/qfl3x/32/16227_2.png) [@Qfl3x](https://discourse.julialang.org/u/Qfl3x)\
**Post date:** [December 19, 2023, 10:18am UTC](https://discourse.julialang.org/t/solving-system-of-ode-for-updating-initial-conditions/107784/4 "2023-12-19T10:18:21Z")

</div>

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

```julia
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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![acubed](https://avatars.discourse-cdn.com/v4/letter/a/b19c9b/32.png) [@acubed](https://discourse.julialang.org/u/acubed)\
**Post date:** [December 19, 2023, 12:54pm UTC](https://discourse.julialang.org/t/solving-system-of-ode-for-updating-initial-conditions/107784/5 "2023-12-19T12:54:47Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Qfl3x](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/qfl3x/32/16227_2.png) [@Qfl3x](https://discourse.julialang.org/u/Qfl3x)\
**Post date:** [December 19, 2023, 1:12pm UTC](https://discourse.julialang.org/t/solving-system-of-ode-for-updating-initial-conditions/107784/6 "2023-12-19T13:12:28Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [December 19, 2023, 1:45pm UTC](https://discourse.julialang.org/t/solving-system-of-ode-for-updating-initial-conditions/107784/7 "2023-12-19T13:45:19Z")

</div>

> [@acubed](#):
>
> indeed, my goal is to change adibatically paramater K and compute some property based on solutions of ODE.

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

---

<div class="post-metadata">

**Author:** ![acubed](https://avatars.discourse-cdn.com/v4/letter/a/b19c9b/32.png) [@acubed](https://discourse.julialang.org/u/acubed)\
**Post date:** [December 19, 2023, 2:38pm UTC](https://discourse.julialang.org/t/solving-system-of-ode-for-updating-initial-conditions/107784/8 "2023-12-19T14:38:07Z")

</div>

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