Orbital satellite with Runge Kutta 4th Order Numerical Approximation

Question

Steven Han 2020-11-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/645068-orbital-satellite-with-runge-kutta-4th-order-numerical-approximation

回答： Meysam Mahooti 2021-7-28

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

James Tursa 2020-11-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/645068-orbital-satellite-with-runge-kutta-4th-order-numerical-approximation#answer_541838

编辑：James Tursa 2020-11-14

在 MATLAB Online 中打开

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

James Tursa 2020-11-13

编辑：James Tursa 2020-11-13

在 MATLAB Online 中打开

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 2020-11-13

编辑：James Tursa 2020-11-13

在 MATLAB Online 中打开

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