# Plotting Orbits on MATLAB

Tags:
1. Aug 16, 2016

### 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
Code (Text):
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.

2. Aug 16, 2016

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

3. Aug 16, 2016

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

4. Aug 16, 2016

### Dr Transport

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

5. Aug 16, 2016

### Student100

Where are you defining t?