Orbital satellite with Runge Kutta 4th Order Numerical Approximation

50 views (last 30 days)
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)');

Answers (2)

James Tursa
James Tursa on 12 Nov 2020
Edited: James Tursa on 14 Nov 2020
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:
  3 Comments
James Tursa
James Tursa on 13 Nov 2020
Edited: James Tursa on 13 Nov 2020
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
James Tursa on 13 Nov 2020
Edited: James Tursa on 13 Nov 2020
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:

Sign in to comment.


Meysam Mahooti
Meysam Mahooti on 28 Jul 2021
https://www.mathworks.com/matlabcentral/fileexchange/55430-runge-kutta-4th-order?s_tid=srchtitle

Categories

Find more on Reference Applications in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!