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.(adsbygoogle = window.adsbygoogle || []).push({});

My code is given by

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

Please help

**Physics Forums | Science Articles, Homework Help, Discussion**

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Plotting Orbits on MATLAB

Have something to add?

Draft saved
Draft deleted

Loading...

Similar Threads for Plotting Orbits MATLAB |
---|

Binary bitstream from QPSK with MATLAB? |

Surface plot for a 3D function |

Python (VPython) Can someone help me with my orbit simulator? |

Python Superimpose plots |

Need some advice on an appropriate computer language to use |

**Physics Forums | Science Articles, Homework Help, Discussion**