26 views (last 30 days)

Show older comments

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)');

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:

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:

Meysam Mahooti
on 28 Jul 2021

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

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

Start Hunting!