0
$\begingroup$

I'm solving a problem that involves the planetary motion of Earth, the Sun and Jupiter, and how the bodies affect each other. The Sun is centered at the origin and Jupiter orbits it at a fixed radius $R$. Let ($xj$, $yj$) = $R(\cos(2\pi*t/Tj, \sin(2\pi*t/Tj)$ be the coordinates of Jupiter at time $t$.

The coordinates of Earth are given by ($x$,$y$) where $x$ and $y$ obey the equations of motion:

$x'' = (-4π^2)(x)/((x^2 + y^2)^{3/2}) - ((4π^2)k(x - xj)/(((x-xj)^2 +(y - yj)^2)^{3/2})$

$y'' = (-4π^2)(y)/((x^2 + y^2)^{3/2}) - ((4*π^2)k(y - yj)/(((y-yj)^2 + (y - yj)^2)^{3/2})$

It is the motion of Earth that I am interested in, so I solved the two second order differential equations above using the fourth order Runge-Kutta routine.

I present my code below:

MATLAB:

clear all
close all

k = 0.3;          %defining the value of k
R = 5;            %value of R, the radius of Jupiter's orbit

t_end =12;      %the time for which the motion will be followed

N = 1000;       %the number of iterations to make
dt = t_end/N;   %the time step

t(1)= 0;        %initialising the time array
xj(1) = R;      %initialising the xj array (x-coordinate of Jupiter)
yj(1) = 0;      %initialising the yj array (y-coordinate of Jupiter)

vj = (2*pi)*sqrt(1/R);  %constant velocity of Jupiter
Tj = (2*pi)*R/vj;       %period of motion of Jupiter

x_1(1) = 1;             %initial x-coordinate of earth
x_2(1) = 0;             %initial x-velocity of earth

y_1(1) = 0;             %initial y-coordinate of earth
y_2(1) = 2*pi$;          %initial y-velocity of earth

for n = 1:N

    t(n+1) = t(n) + dt;             %increase time by dt
    xj(n+1) = R*cos(2*pi*t(n)/Tj);  %calculate the next xj
    yj(n+1) = R*sin(2*pi*t(n)/Tj);  %calculate the next yj

    A1 = dt*xvel(x_1,x_2(n),y_1,y_2,t);     %Begin the RK4 routine
    B1 = dt*xacc(x_1(n),x_2,y_1(n),y_2,t);
    C1 = dt*yvel(x_1,x_2,y_1,y_2(n),t);
    D1 = dt*yacc(x_1(n),x_2,y_1(n),y_2,t);

    A2 = dt*xvel(x_1,x_2(n)+ B1/2,y_1,y_2,t + 0.5*dt);
    B2 = dt*xacc(x_1(n)+ A1/2,x_2,y_1(n)+ C1/2,y_2,t + 0.5*dt);
    C2 = dt*yvel(x_1,x_2,y_1,y_2(n)+D1/2,t + 0.5*dt);
    D2 = dt*yacc(x_1(n)+ A1/2,x_2,y_1(n)+C1/2,y_2,t + 0.5*dt);

    A3 = dt*xvel(x_1,x_2(n)+ B2/2,y_1,y_2,t + 0.5*dt);
    B3 = dt*xacc(x_1(n)+ A2/2,x_2,y_1(n)+ C2/2,y_2,t + 0.5*dt);
    C3 = dt*yvel(x_1,x_2,y_1,y_2(n)+D2/2,t + 0.5*dt);
    D3 = dt*yacc(x_1(n)+ A2/2,x_2,y_1(n)+C2/2,y_2,t + 0.5*dt);

    A4 = dt*xvel(x_1,x_2(n)+ B3,y_1,y_2,t + dt);
    B4 = dt*xacc(x_1(n)+ A3,x_2,y_1(n)+ C3,y_2,t + dt);
    C4 = dt*yvel(x_1,x_2,y_1,y_2(n)+D3,t + dt);
    D4 = dt*yacc(x_1(n)+ A3,x_2,y_1(n)+C3,y_2,t + dt);

    x_1(n+1) = x_1(n) + (1/6)*(A1 + 2*A2 + 2*A3 + A4);
    x_2(n+1) = x_2(n) + (1/6)*(B1 + 2*B2 + 2*B3 + B4);
    y_1(n+1) = y_1(n) + (1/6)*(C1 + 2*C2 + 2*C3 + C4);
    y_2(n+1) = y_2(n) + (1/6)*(D1 + 2*D2 + 2*D3 + D4);

end

plot(x_1,y_1,'b',xj,yj,'r-',0,0,'o')       %plot the orbit of Earth, Jupiter,and Sun

I have defined the functions in the RK-4 routine as follows:

function funcyvel = yvel(x_1,x_2,y_1,y_2,t)

    funcyvel = y_2;


function funcxvel = xvel(x_1,x_2,y_1,y_2,t)

    funcxvel = x_2;


function funcxacc = xacc(x_1,x_2,y_1,y_2,t)

    n = 1;
    R =5;

    k=0.3;
    t(1) = 0;
    t_end = 12;
    N = 1000;
    dt = t_end/N;

    vj = (2*pi)*sqrt(1/R);
    Tj = (2*pi)*R/vj;

    for n = 1:N
        t(n+1) = t(n) + dt;
    end

    function juptx = xj(dummy)  %function returns the x-coordinates of Jupiter
        juptx = R*cos(2*pi*t/Tj);      
    end

    function jupty = yj(dummy)   %function returns the y-coordinates of Jupiter
        jupty = R*sin(2*pi*t/Tj);   
     end

    funcxacc = (-4*pi^2)*(x_1)/((x_1.^2 + y_1.^2).^{3/2}) - (((4*pi^2*k)*(x_1 - xj))/(((x_1-xj).^2 +(y_1 - yj).^2).^{3/2}));

end




function funcyacc = yacc(x_1,x_2,y_1,y_2,t)

    R =5;
    k=0.3;
    t(1) = 0;
    t_end = 12;
    N = 1000;

    dt = t_end/N;

    vj = (2*pi)*sqrt(1/R);
    Tj = (2*pi)*R/vj;

    for n = 1:N
        t(n+1) = t(n) + dt;
    end

    function juptx = xj(dummy)      %function returns the x-coordinates of Jupiter

        juptx = R*cos(2*pi*t/Tj);

    end

     function jupty = yj(dummy)     %function returns the y-coordinates of Jupiter 

        jupty = R*sin(2*pi*t/Tj);

    end

    funcyacc = (-4*pi^2)*(y_1)/((x_1.^2 + y_1.^2).^{3/2}) - (((4*pi^2*k)*(y_1 - yj))/(((x_1-xj).^2 +(y_1 - yj).^2).^{3/2}));

end

PROBLEM

When I run the above script, it plots the orbits of Earth and Jupiter with the Sun at the centre but Earth's orbit seems unstable and it varies slightly with time. If I vary the initial conditions for Earth, then it sometimes has a tendency to fly off into outer space after a certain period of time. I have tried reducing the time step dt, but I keep getting a non-constant orbit.Is this a result of the inaccuracy of the RK_4 method? Or is there an error in the way I handled the coordinates of Jupiter?

As a side note, I would appreciate suggestions on how to make the script run faster, since for large N it takes a very long time.

$\endgroup$
  • 1
    $\begingroup$ You should expect a non-constant orbit from a method like RK4, because it does not build conservation of energy into the method. There are conservative numerical methods out there for these sorts of issues. But the three-body problem is quite hard, even numerically. Anyway, for speed you should not nest your functions so much, and also you can use MATLAB's ode45, which implements RK4. $\endgroup$ – Ian Oct 3 '15 at 15:32
  • 1
    $\begingroup$ The buzzword for the type of method you want is "symplectic". $\endgroup$ – Ian Oct 3 '15 at 15:49
  • $\begingroup$ I see! Thank you Ian! $\endgroup$ – 2good4this Oct 3 '15 at 16:38
  • $\begingroup$ @Ian: ode45 implements the Dormand-Price method, also called dopri5. This has almost nothing to do with the classical Runge-Kutta method, except the generality that it is a Runge-Kutta method (or better 2 RK methods in one). $\endgroup$ – LutzL Oct 3 '15 at 19:17
  • $\begingroup$ @LutzL My impression was that ode45 was classical RK4 with a time step selection algorithm. Apparently (looking at the Butcher tableau for the Dormand-Price method) I was wrong. So thanks for the correction. $\endgroup$ – Ian Oct 3 '15 at 19:25
1
$\begingroup$

Your system is not conservative. Thus there is no implied stability of the orbit. You have an oscillator that is periodically excited. Resonance is one possible result.


Looking at your code, it is one of the finer examples for the necessity of structured programming. There should be only one ODE system function implementing the derivatives as vector, and the application of the integration method then just uses the vectors.


There is a sweet spot of maximal accuracy for numerical integration methods. You can test if your step size dt is above this point by comparing the results for step sizes dt, 2*dt, 4*dt. For the RK4 method the differences of the successive results should be apart by a very clear factor 16.

For tame examples with first, second etc. derivatives in the scale of $1$, the rule is that the sweet spot is at dt=1e-3. Above that the discretization error increases, below that the accumulation of the floating point errors becomes progressively larger.

The point of this reasoning is that you should not have problems with the run-time of large N because there is no reason to use overly large N.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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