# How to perform Runge-Kutta 4th order with functions of functions

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. • 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. – Ian Oct 3 '15 at 15:32 • The buzzword for the type of method you want is "symplectic". – Ian Oct 3 '15 at 15:49 • I see! Thank you Ian! – 2good4this Oct 3 '15 at 16:38 • @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). – LutzL Oct 3 '15 at 19:17 • @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. – Ian Oct 3 '15 at 19:25 ## 1 Answer 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.