Solving Orbital Motion Equations with RK4 Using MATLAB

[Jordan_Tusc]

I'm trying to plot the solutions of the second order differential equation d^2R/dt^2 = GM/R^2 + Lz^2/R^3. I'm reducing this to a system of first order ODEs and then using RK4 to solve this system.

My code is given by

function RK4system()                                      
Tsim   = 10;           % simulate over Tsim seconds
h      = 0.001;         % step size
N      = Tsim/h;       % number of steps to simulate
x4      = ones(2,N);
t4      = zeros(1,N);
theta  = ones(1,N);
G      = 6.67*power(10,-11);
M      = 1.5*1.989*power(10,41);
Lz     = 2.623*power(10,2);
R_1    = 2.623*power(10,20);
U_1    = 0.5;
x(1,1) = U_1;          %IC_ODE_1
x(2,1) = (G*M/(power(R_1,2))+power(Lz,2)/(power(R_1,3)));  %IC_ODE_2

x4 = ones(2,N);
t4 = zeros(1,N);

for n = 1:N-1
   t4(n+1)      = t4(n) + h;
   k1_4         = Orbit(t4(n), x4(:,n));
   k2_4         = Orbit(t4(n)+0.5*h, x4(:,n)+k1_4*0.5*h);
   k3_4         = Orbit(t4(n)+0.5*h, x4(:,n)+0.5*k2_4*h);
   k4_4         = Orbit(t4(n)+h, x4(:,n)+k3_4*h);
   phi4         = (1/6)*(k1_4 + 2*k2_4 + 2*k3_4 + k4_4);
   x4(:,n+1)    = x4(:,n) + phi4*h;
   theta(:,n+1)  = theta(:,n) + h*Angle(t(n), theta(:,n));
   xvar4(:,n+1) = x2(:,n+1)*cos(t4(n));
   yvar4(:,n+1) = x2(:,n+1)*sin(t4(n));
end

figure()
plot(xvar4,yvar4, 'b-', 'LineWidth', 1)
end

function dtheta = Angle(t,R)
Lz = 100000*2.623*power(10,20);
dtheta(1) = Lz/(power(R(1),2));
end

function dx = Orbit(t,R)
G           = 6.67*power(10,-11);
M           = 1.5*1.989*power(10,41);
Lz          = 100000*2.623*power(10,20);
% Equations of motion
dx          = zeros(2,1);
dx(1)       = R(2);
dx(2)       = (G*M/(power(R(1),2))+power(Lz,2)/(power(R(1),3)));
end

I'm not getting the plot I want, If I set V = sqrt(GM/R) and Lz = R*sqrt(GM/R) I should get a circular orbit. If V > sqrt(GM/R) I should get a parabolic orbit. This is not what I'm getting.

Please help

 

[JorisL]

Can you attach some plots for those cases?
Images can help us spot some common errors more quickly than reading code.

One obvious question, what does Orbit() mean?
Finally I usually refrain from using functions like power and for example dot.
Once you start running longer simulations this will slow down your program significantly due to checks built into those functions.
In one case I was able to decrease the runtime of a simulation by a factor 10 by using $a^\prime.a$ for finding the length of a vector.

 

[Dr Transport]

did you try running this rk4 routine against something like a harmonic oscillator to check and see if it gives correct answers. I also suspect that the constants are giving it problems. R1 is on the order of $10^{20}$ and when you square and cube it then divide, your essentially getting zero for some of the terms.

 

[Dr Transport]

$T_{sim}$ is only 10 seconds... it takes a year to orbit the sun..., you barely moved.

 

[Student100]

Where are you defining t?