0
$\begingroup$

I'm trying to evaluate the perturbations magnitude between 2 body orbiting a central one in three dimensions. In order to do this I need to have an estimate of the error, which I did using Richardson extrapolation as described here. I'm using Runge-Kutta 4 and also velocity Verlet. Things seems good in the unperturbed problem, which means I obtain an error $O(10^{-7})$ and $O(10^{-12})$ for Verlet and Runge-Kutta respectively by using a time step of $0.001$. When I consider the perturbations, only velocity Verlet still give me the expected error, while Runge-Kutta give me an error of the same order of Verlet, using the same previous time step. The same happens even if I reduce or increase the time step. Actually for the perturbed problem I already expected something bad related to the error, since that's what I obtained for the Richardson's fraction for Runge-Kutta

enter image description here

As it can be seen, the level of the figure should be about $2^p = 16$, but it's $2$ actually. That sort of spikes are present also in the same plot for velocity Verlet, but for the latter their magnitudes are much less than that in the plot above and the $F_h$ values are about $4$ anyway, as expected for Verlet. I checked all the Runge-Kutta code and the ode system code but I didn't find any error. Moreover I tested both methods with a simpler function obtaining correct results, so I can't understand the reason for such wrong error estimates for Runge-Kutta.

Edit 1

The function code is the following

vector<double> F(double t, vector<double> x, vector<CB> objs, vector<double> s) {
    //objs store the masses of the objects
    //s store the positions of the objects in the order [x1,y1,z1,x2,y2,z2,...];
    //x store the positions and velocity of the target (perturbed) object in the order [x,Vx,y,Vy,z,Vz]
    vector<double> p(x.size(),0);
    double m, d;
    vector<double> y(x.size()/2,0); 

    for (int i = 0; i < objs.size(); i++) {
        //temporary store the mass of the i-th object
        m = objs[i].CB::getmass(); 
        //temporary store the position of the i-th object from the s vector
        y = {s[i*x.size()/2], s[i*x.size()/2 + 1], s[i*x.size()/2 + 2]}; 

        d = sqrt(pow((x[0] - y[0]),2) + pow((x[2] - y[1]),2) + pow((x[4] - y[2]),2));

        //x system
        p[0] = x[1];
        p[1] += -G*m*(x[0] - y[0])/pow(d,3); 
        //y system                           
        p[2] = x[3];
        p[3] += -G*m*(x[2] - y[1])/pow(d,3);
        //z system
        p[4] = x[5];
        p[5] += -G*m*(x[4] - y[2])/pow(d,3);
    }

    return p;
}
Improve this question
$\endgroup$
10
  • $\begingroup$ What is your setup approximately? A starting configuration using real planetary data, or some artificial initial condition like one planet on a circular orbit of radius 1 and period 1, the other at radius 5, with a fixed proportion of masses, so that setting the mass factor turns the interaction on or off? $\endgroup$
    – Lutz Lehmann
    Commented Mar 2, 2021 at 10:59
  • $\begingroup$ @LutzLehmann I use real planetary data taken from JPL Horizon system, cause I need to analyze different planetary configurations, of which the simpler one involve Sun, Jupiter and Saturn, with Jupiter perturbed by Saturn, which corrispondes to the perturbed case of my question. The motion of every planet involved in the simulation start from perihelion. The refernce frame is centered on the Sun. $\endgroup$
    – Zebx
    Commented Mar 2, 2021 at 12:44
  • $\begingroup$ The errors oscillate like the periodic orbits, close to the zero-crossings you get the spikes. This can be prevented by taking the norm over all coordinate differences of the body. There must be something wrong with the method step that reduces the order to 1. How did you implement the perturbation? In the acceleration computation, or as something extra that is added to the method step? $\endgroup$
    – Lutz Lehmann
    Commented Mar 2, 2021 at 12:55
  • $\begingroup$ @LutzLehmann The perturbations are implemented in the acceleration computation. They're "turned on" if more than one object is considered in the computations. $\endgroup$
    – Zebx
    Commented Mar 2, 2021 at 13:06
  • $\begingroup$ Some hour ago I also compared my results to the one obtained using another Runge-Kutta 4 implementaion I found, but I got the same results for $F_h$. $\endgroup$
    – Zebx
    Commented Mar 2, 2021 at 13:32

1 Answer 1

Reset to default
1
$\begingroup$

To explain further what I wrote in the comments: Your function has the form F(t,x,co,xo) where x is the state, co the constants for the other objects and xo the dynamical data for the other objects.

Variant 1: You will get an order 1 approximation out of RK4 if you use the external dynamic in the RK4 loop as (in python notation)

dx1 = dt*F(t[k]       , x[k]        , co, xo[k])
dx2 = dt*F(t[k]+0.5*dt, x[k]+0.5*dx1, co, xo[k])
dx3 = dt*F(t[k]+0.5*dt, x[k]+0.5*dx2, co, xo[k])
dx4 = dt*F(t[k]+    dt, x[k]+    dx3, co, xo[k])
x[k+1]=x[k]+(dx1+2*dx2+2*dx3*dx4)/6

Variant 2: You will get an order 2 approximation if you use linear interpolation as in

dx1 = dt*F(t[k]       , x[k]        , co, xo[k])
dx2 = dt*F(t[k]+0.5*dt, x[k]+0.5*dx1, co, 0.5*(xo[k]+xo[k+1]))
dx3 = dt*F(t[k]+0.5*dt, x[k]+0.5*dx2, co, 0.5*(xo[k]+xo[k+1]))
dx4 = dt*F(t[k]+    dt, x[k]+    dx3, co, xo[k+1])
x[k+1]=x[k]+(dx1+2*dx2+2*dx3*dx4)/6

Variant 3: To get an error that goes down in 4th order with the step size, you need to integrate the other objects with double the rate and half the step size (or find an order 4 interpolation for the data) so that the loop then reads as

dx1 = dt*F(t[k]       , x[k]        , co, xo[2*k])
dx2 = dt*F(t[k]+0.5*dt, x[k]+0.5*dx1, co, xo[2*k+1])
dx3 = dt*F(t[k]+0.5*dt, x[k]+0.5*dx2, co, xo[2*k+1])
dx4 = dt*F(t[k]+    dt, x[k]+    dx3, co, xo[2*k+2])
x[k+1]=x[k]+(dx1+2*dx2+2*dx3*dx4)/6

Variant 4: Of course, having a piecewise polynomial interpolation xo_interp(t) of satisfying accuracy for the external data and using that same interpolation function for all step sizes in the RK4 experiment should also restore the 4th order

dx1 = dt*F(t[k]       , x[k]        , co, xo_interp(t[k]))
dx2 = dt*F(t[k]+0.5*dt, x[k]+0.5*dx1, co, xo_interp(t[k]+0.5*dt))
dx3 = dt*F(t[k]+0.5*dt, x[k]+0.5*dx2, co, xo_interp(t[k]+0.5*dt))
dx4 = dt*F(t[k]+    dt, x[k]+    dx3, co, xo_interp(t[k+1]))
x[k+1]=x[k]+(dx1+2*dx2+2*dx3*dx4)/6
Improve this answer
$\endgroup$
3
  • $\begingroup$ Ok, so it seems I have at least two problems actually: the one which causes the order 1 problem when I consider a perturbing object and the one that causes the spikes. I modified the code in order to solve the first one as explained in your answer and now $F_h$ is better, meaning that at least the flat parts of the $F_h$ plot are around $16$, using variant 3. But why xo is "seen" as order $1$ if I evaluated it with RK4? In other words, why does variant $1$ give me an order $1$ x approsimation if xo was evaluated using RK4? $\endgroup$
    – Zebx
    Commented Mar 3, 2021 at 13:58
  • 1
    $\begingroup$ Because the step function $x_o([t/dt]·dt)$ has error $O(dt)$ against $x_o(t)$, the difference is a saw-tooth function with height $\simeq x_o'(t)dt$. This gives an error $O(dt^2)$ in each step, which accumulates to an error of $O(t·dt)$ over the course of the integration. Similarly one knows that $\frac{x_o(t+dt)+x_o(t)}2-x_o(t+\frac12 dt)\simeq \frac12x_o''(t)dt^2$, so that variant 2 gives an error of $O(dt^2)$ in each step against a step with higher-order exact $x_o$ values. $\endgroup$
    – Lutz Lehmann
    Commented Mar 3, 2021 at 14:05
  • $\begingroup$ Oh ok, I wasn't considering the error accumulation over time. Thank you very much! $\endgroup$
    – Zebx
    Commented Mar 3, 2021 at 20:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.