# Different solution if solving a differential equation in chunks or for full time at once

**URL:** <https://discourse.julialang.org/t/different-solution-if-solving-a-differential-equation-in-chunks-or-for-full-time-at-once/121351>\
**Category:** Modelling & Simulations\
**Tags:** diffeq, differentialequation\
**Created:** [October 15, 2024, 9:32pm UTC](https://discourse.julialang.org/t/different-solution-if-solving-a-differential-equation-in-chunks-or-for-full-time-at-once/121351 "2024-10-15T21:32:21Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![shaurya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shaurya/32/30468_2.png) [@shaurya](https://discourse.julialang.org/u/shaurya)\
**Post date:** [October 15, 2024, 9:32pm UTC](https://discourse.julialang.org/t/different-solution-if-solving-a-differential-equation-in-chunks-or-for-full-time-at-once/121351/1 "2024-10-15T21:32:22Z")

</div>

Hi! I’m trying a simulate an orbit in the chaotic 3 body dynamics, and perform two different simulations: one for full time (0, nT) where n is an integer and T is the time period of the orbit. On the other hand, I also perform it by breaking it into chunks: (0,\delta t), (\delta t, 2\delta t) and so on until it reaches nT. The final solution from the previous step is the initial condition for the next step. Theoretically I assume this should mean similar solution (except maybe some numerical inconsistencies.) However, the results are very different for long simulations: such as 50 orbits. I understand the dynamics I’m simulating are chaotic, but I need to work with the second way: integrating by chunks and I’m not sure if that is the correct way here. I’ve attached the code below:

```julia
using DifferentialEquations, LinearAlgebra
const lstar = 238529; # in kms
const tstar = 18913; # in seconds
const μ = 1.901109735892602e-7;

# ODE function
function cr3bp_propagator!(dxdt,x,p,t)
    d=sqrt((x[1]+μ)^2 + x[2]^2+ x[3]^2);
    r=sqrt((x[1]-1+μ)^2 + x[2]^2 + x[3]^2);

    dxdt[1]= x[4];
    dxdt[2]= x[5];
    dxdt[3]= x[6];
    dxdt[4]= 2*x[5] + x[1] - (1-μ)*(x[1]+μ)/(d^3) - μ*(x[1]-1+μ)/(r^3);
    dxdt[5]= -2*x[4] +x[2] - (1-μ)*(x[2])/(d^3) - μ*(x[2])/(r^3);
    dxdt[6]= -(1-μ)*(x[3])/(d^3) - μ*(x[3])/(r^3);
end

function compare_full_chunks(N,nT)

    nrho = [ 1.0025067861477979
                    0
                    4.8938916019621225e-3
                    0
                    -5.3520383247404298e-3
                    0 ];
    x0 = -[50/lstar;0/lstar;50/lstar;0*tstar/lstar;0*tstar/lstar;1e-6*tstar/lstar];
    TP = 2.2460348880483676;
    δt= TP/N; t=0;
    ic = nrho+x0; 
    traj=[]; time=[];tol=1e-12;

    # solving in full at once
    probref=ODEProblem(cr3bp_propagator!,nrho,(0.0,nT*TP));
    solref=solve(probref,TsitPap8(),reltol=1e-12,abstol=1e-12);

    # solving in chunks of δt time steps
    for i in 1:nT*N
        tspan = (t,t+δt);
        prob=ODEProblem(cr3bp_propagator!,ic,tspan);
        sol =solve(prob,TsitPap8(),reltol=tol,abstol=tol);
        append!(traj,sol.u);
        append!(time,sol.t);
        ic = sol.u[end];
        t+=δt;
    end

    return solref.u[end]-traj[end]
end

```

When I call` compare_full_chunks(50,50)`, the returned value which is the difference between solution end points for the two methods is not close to zero and I’m confused why that is the case. Thanks for your help!

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [October 16, 2024, 12:59am UTC](https://discourse.julialang.org/t/different-solution-if-solving-a-differential-equation-in-chunks-or-for-full-time-at-once/121351/2 "2024-10-16T00:59:36Z")

</div>

when simulating chaotic systems for large amounts of time you will get diverging results.

---

<div class="post-metadata">

**Author:** ![Sevi](https://avatars.discourse-cdn.com/v4/letter/s/c67d28/32.png) [@Sevi](https://discourse.julialang.org/u/Sevi)\
**Post date:** [October 16, 2024, 6:24am UTC](https://discourse.julialang.org/t/different-solution-if-solving-a-differential-equation-in-chunks-or-for-full-time-at-once/121351/3 "2024-10-16T06:24:22Z")

</div>

You could try to fix the times the solver steps to in the first case (I believe with the keyword `saveat`? EDIT: Nope, it’s `tstops`, see below) and see if it improves the agreement. But the problem might still persist.

You essentially force the solver to step to predefined times in the second case, but it might choose very different times in the first case. Small rounding errors/differences in each step can definitely lead to arbitrarily far-apart solutions in a chaotic system as @Oscar_Smith mentioned – it’s more or less the definition of a chaotic system.

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [October 16, 2024, 8:40am UTC](https://discourse.julialang.org/t/different-solution-if-solving-a-differential-equation-in-chunks-or-for-full-time-at-once/121351/4 "2024-10-16T08:40:49Z")

</div>

> [@Sevi](#):
>
> You could try to fix the times the solver steps to in the first case (I believe with the keyword `saveat`?) and see if it improves the agreement. But the problem might still persist.

`saveat` is not what you want here. It is `tstops`. Or even better, use a solver with a non-adaptive step and set your large Δt to be a multiple of the small solver step δt.

But honestly, this is all a fools errand. It is very hard to get “exactly same results” for chaotic systems due to even the tiniest change in the order of operations at any point in the simulation.

More importantly, you shouldn’t care about having identical results for chaotic systems. There is nothing confusing about not getting identical results. The confusing thing would be if you _were_ getting identical results with different ways of simulation 😃

---

<div class="post-metadata">

**Author:** ![Datseris](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/datseris/32/13406_2.png) [@Datseris](https://discourse.julialang.org/u/Datseris)\
**Post date:** [October 16, 2024, 8:42am UTC](https://discourse.julialang.org/t/different-solution-if-solving-a-differential-equation-in-chunks-or-for-full-time-at-once/121351/5 "2024-10-16T08:42:05Z")

</div>

> [@shaurya](#):
>
> but I need to work with the second way: integrating by chunks and I’m not sure if that is the correct way here. I

“Correct” here does not have an objective meaning. It is what you want to do. Solving a chaotic system numerically is always an approximation and you always get a fake (formally called “shadow”) orbit, no matter what way you simulate it. If I were you I wouldn’t worry about this at all and proceed simulating in chunks.

---

<div class="post-metadata">

**Author:** ![shaurya](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/shaurya/32/30468_2.png) [@shaurya](https://discourse.julialang.org/u/shaurya)\
**Post date:** [October 16, 2024, 9:56pm UTC](https://discourse.julialang.org/t/different-solution-if-solving-a-differential-equation-in-chunks-or-for-full-time-at-once/121351/6 "2024-10-16T21:56:30Z")

</div>

Thank you! I think this answers my confusion and this is more of a chaotic system property and possibly error accumulation in the long term simulations rather than anything else. I will just proceed with the simulation in chunks. Thanks for you reply, appreciate it!