Orbital satellite with Runge Kutta 4th Order Numerical Approximation

I am writing a script for an orbiting satelite around Earth. Currently, when I try to run the code, the code is showing NaN. I am unsure where the mistake was made such that my code is dividing my 0. The equations of motion is correct. Thanks in advance!

%% RK4 MATLAB SCRIPT 

% Parameters
Me = 5.972e24;         % mass of body a
Mq = 1000;         % mass of particle q
G = 6.67408e-11;    % gravitational constant
tf = 30;        % final time
h = 0.01;        % step size
r = 100000;      % distance of satelite from earth


%equations of motion written in state space form
f1 = @(t, r, rdot, theta, thetad) rdot; 
f2 = @(t, r, rdot, theta, thetad) -G*Me/r^2 + r*thetad^2;
f3 = @(t, r, rdot, theta, thetad) thetad;
f4 = @(t, r, rdot, theta, thetad) 2*rdot*thetad/r;

% Initial Conditions (testing random ICs)

r(1) = 1000; %inital position of satellite
rdot(1) = 50; %intial translational rate of satellite
theta(1) = 0; %inital angular position of satellite
thetad(1) = 0; % initial angular rate of satellite

% initial time
t(1) = 0;

% RK4 Loop
for i = 1:tf/h
    
    t(i+1) = t(i) + h;
    
    k1r = f1(t(i), r(i), rdot(i), theta(i), thetad(i));
    k1rdot = f2(t(i), r(i), rdot(i), theta(i), thetad(i));
    k1theta = f3(t(i), r(i), rdot(i), theta(i), thetad(i));
    k1thetad = f4(t(i), r(i), rdot(i), theta(i), thetad(i));
     
    k2r = f1(t(i)+h/2, r(i)+k1r*h/2, rdot(i)+k1rdot*h/2, theta(i)+k1theta*h/2, thetad(i)+k1thetad*h/2);
    k2rdot = f2(t(i)+h/2, r(i)+k1r*h/2, rdot(i)+k1rdot*h/2, theta(i)+k1theta*h/2, thetad(i)+k1thetad*h/2);
    k2theta = f3(t(i)+h/2, r(i)+k1r*h/2, rdot(i)+k1rdot*h/2, theta(i)+k1theta*h/2, thetad(i)+k1thetad*h/2);
    k2thetad = f4(t(i)+h/2, r(i)+k1r*h/2, rdot(i)+k1rdot*h/2, theta(i)+k1theta*h/2, thetad(i)+k1thetad*h/2);
    
    k3r = f1(t(i)+h/2, r(i)+k2r*h/2, rdot(i)+k2rdot*h/2, theta(i)+k2theta*h/2, thetad(i)+k2thetad*h/2);
    k3rdot = f2(t(i)+h/2, r(i)+k2r*h/2, rdot(i)+k2rdot*h/2, theta(i)+k2theta*h/2, thetad(i)+k2thetad*h/2);
    k3theta = f3(t(i)+h/2, r(i)+k2r*h/2, rdot(i)+k2rdot*h/2, theta(i)+k2theta*h/2, thetad(i)+k2thetad*h/2);
    k3thetad = f4(t(i)+h/2, r(i)+k2r*h/2, rdot(i)+k2rdot*h/2, theta(i)+k2theta*h/2, thetad(i)+k2thetad*h/2);
    
    k4r = f1(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3r*h);
    k4rdot = f2(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3r*h);
    k4theta = f3(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3r*h);
    k4thetad = f4(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3r*h);
    
    %updating position and speed
    r(i+1,1) = r(i,1) + (h/6)*(k1r+2*k2r+2*k3r+k4r);
    rdot(i+1,1) = rdot(i,1) + (h/6)*(k1rdot+2*k2rdot+2*k3rdot+k4rdot);
    theta(i+1,1) = theta(i,1) + (h/6)*(k1theta+2*k2theta+2*k3theta+k4theta);
    thetad(i+1,1) = thetad(i,1) + (h/6)*(k1thetad+2*k2thetad+2*k3thetad+k4thetad);
end

% Plots
figure(1)
plot(t,r,'b');
title('Position');
xlabel('time (s)');
ylabel('position (m)');

Answer 1

[James Tursa]

This code

    k4r = f1(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3r*h);
    k4rdot = f2(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3r*h);
    k4theta = f3(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3r*h);
    k4thetad = f4(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3r*h);

should be this instead

    k4r = f1(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3thetad*h);
    k4rdot = f2(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3thetad*h);
    k4theta = f3(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3thetad*h);
    k4thetad = f4(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3thetad*h);

Having said that, your code is WAY too complicated for what it needs to be. You basically have to write four separate RK4 sets of code for four scalar equations. It would be much better if you were to define a single 4-element state vector and then write one set of RK4 equations to handle the propagation. That means 1/4 of the code you have currently written and less chance to make these copy & paste errors like you have done.

Also, these initial conditions

r(1) = 1000; %inital position of satellite
rdot(1) = 50; %intial translational rate of satellite
theta(1) = 0; %inital angular position of satellite
thetad(1) = 0; % initial angular rate of satellite

put the satellite inside the surface of the Earth to begin with, and result in rectilinear motion. I don't think that is what you want.

Additionally,

Q1: Why isn't f2 -G*Me/r^2 - r*thetad^2

EDIT:

Wrong sign correction suggested above. What I should have suggested is leaving f2 as is and making this change

Why isn't f4 -2*rdot*thetad/r

Also, you can see this related thread:

https://www.mathworks.com/matlabcentral/answers/508668-ode45-dynamics-rocket-around-earth-equation-of-motion?s_tid=srchtitle

[Steven Han]

Hi James, thanks a ton for replying. Ideally, I want to plot this such that it output x and y coordinates. I am not too fluent with Matlab (programming in general). with that being said, I have changed the initial conditions to below and yet still receiving NaN errors and rectilinear plot. The equations of motion I used had only a sign error.

r(1) = 42164000; %inital position of satellite
rdot(1) = 50; %intial translational rate of satellite
theta(1) = pi; %inital angular position of satellite
thetad(1) = pi/4; % initial angular rate of satellite

[James Tursa]

Initial theta dot is pi/4 radians per second???

Seems like you should be starting with something close to this (near circular orbit)

thetad(1) = sqrt(G*Me/r(1)^3);

[James Tursa]

If I do the following:

• Make the corrections to the k4 terms noted above
• Make the correction to the sign of the f4 function handle
• Use a reasonable initial condition

Then I get reasonable results. E.g.,

    k4r = f1(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3thetad*h);
    k4rdot = f2(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3thetad*h);
    k4theta = f3(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3thetad*h);
    k4thetad = f4(t(i)+h, r(i)+k3r*h, rdot(i)+k3rdot*h, theta(i)+k3theta*h, thetad(i)+k3thetad*h);

and

    f4 = @(t, r, rdot, theta, thetad) -2*rdot*thetad/r;

and

h = 1;
tf = 2*86400; % Two days

and

r(1) = 42164000; %inital position of satellite
rdot(1) = 0; % circular orbit
theta(1) = 0; % circular orbit
thetad(1) = sqrt(G*Me/r(1)^3); % circular orbit

and

th = t/3600;

% Plots
figure
plot(th,r,'b');
title('R');
xlabel('time (hours)');
ylabel('R (m)');
grid on

figure
plot(th,rdot,'b');
title('Rdot');
xlabel('time (hours)');
ylabel('Rdot (m/s)');
grid on

figure
plot(th,theta*(180/pi),'b');
title('Theta');
xlabel('time (hours)');
ylabel('Theta (deg)');
grid on

figure
plot(th,thetad*(180/pi)*86400,'b');
title('Thetadot');
xlabel('time (hours)');
ylabel('Thetadot (deg/day)');
grid on

x = r.*cos(theta);
y = r.*sin(theta);
figure
plot(x,y,'b');
title('Orbit');
xlabel('X (m)');
ylabel('Y (m)');
axis square
grid on

Produce the following, which all seem reasonable for a circular orbit:

Answer 2

[Meysam Mahooti]

https://www.mathworks.com/matlabcentral/fileexchange/55430-runge-kutta-4th-order?s_tid=srchtitle