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

Specific Domains Modelling & Simulations
,
1

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:

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!

2

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

3 Likes
3

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.

1 Like
4

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 :smiley:

2 Likes
5

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

1 Like
6

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!

1 Like