Orbital satellite with Runge Kutta 4th Order Numerical Approximation
24 次查看(过去 30 天)
显示 更早的评论
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 个)
James Tursa
2020-11-12
编辑:James Tursa
2020-11-14
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 个评论
James Tursa
2020-11-13
编辑:James Tursa
2020-11-13
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
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/413408/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/413413/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/413418/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/413448/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/413428/image.jpeg)
Meysam Mahooti
2021-7-28
https://www.mathworks.com/matlabcentral/fileexchange/55430-runge-kutta-4th-order?s_tid=srchtitle
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Reference Applications 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!