MATLAB Answers

Runge-Kutta - Orbital mechanics application

22 views (last 30 days)
Sarah on 21 Oct 2015
Edited: James Tursa on 21 Oct 2015
I am attempting to write a code to numerically integrate the equations of motion for 5400 seconds, in increments of 10 seconds using the Runge-Kutta method. I was given 6 orbital elements and was able to find my initial R and V vectors. Below is my script(I have also attached the .m files).
clear all
a = 6652.555663;
e = 0.075;
i = deg2rad(28.5);
OMEGA = deg2rad(40);
omega = deg2rad(30);
nu = deg2rad(0);
u = 3.986e5;
r = (a*(1-e^2))/(1+e*cos(nu));
rp = [r*cos(nu) r*sin(nu)];
vp = [-sqrt(u/(a*(1-e^2)))*sin(nu) sqrt(u/(a*(1-e^2)))*(e+cos(nu))];
Q11 = -sin(OMEGA)*cos(i)*sin(omega) + cos(OMEGA)*cos(omega);
Q12 = cos(OMEGA)*cos(i)*sin(omega) + sin(OMEGA)*cos(omega);
Q13 = sin(i)*sin(omega);
Q21 = -sin(OMEGA)*cos(i)*cos(omega) - cos(OMEGA)*sin(omega);
Q22 = cos(OMEGA)*cos(i)*cos(omega) - sin(OMEGA)*sin(omega);
Q23 = sin(i)*cos(omega);
Q31 = sin(OMEGA)*sin(i);
Q32 = -cos(OMEGA)*sin(i);
Q33 = cos(i);
Q = [Q11 Q12 Q13;Q21 Q22 Q23;Q31 Q32 Q33]';
rp=[rp 0]';
vp = [vp 0]';
r = Q*rp;
v = Q*vp;
N = 540;
a = 0;
b = 5400;
h = (b-a)/N;
t = a;
w = [r(1);r(2);r(3);v(1);v(2);v(3)];
for i = 1:10:5400
k1 = h*X(t,w)
k2 = h*X(t+h/2,w+k1/2);
k3 = h*X(t+h/2,w+k2/2);
k4 = h*X(t+h,w+k3);
w = w + (1/6)*(k1+k2+k3+k4)
t = a+i*h
Here is my function script:
function [ X ] = rky( t, x )
X(1) = x(4);
X(2) = x(5);
X(3) = x(6);
X(4) = (-u*x(1))/((x(1)^2+x(2)^2+x(3)^2))^(3/2);
X(5) = (-u*x(2))/((x(1)^2+x(2)^2+x(3)^2))^(3/2);
X(6) = (-u*x(3))/((x(1)^2+x(2)^2+x(3)^2))^(3/2);
X=[X(1); X(2);X(3);X(4);X(5);X(6)];
I am getting errors in both my function and also where the for loop starts. Everything up until that point is fine. I've never used MATLAB before so an help is greatly appreciated.


Sign in to comment.

Answers (1)

James Tursa
James Tursa on 21 Oct 2015
Edited: James Tursa on 21 Oct 2015
I haven't checked all of your code in detail (and I haven't checked any of your orbital mechanics equations at the front), but here are some initial observations:
1) You wrote this nice derivative function called rky ... and then you never call it. Inside your RK loop you need to call this function. E.g.,
k1 = h * rky(t,w)
k2 = h * rky(t+h/2,w+k1/2);
k3 = h * rky(t+h/2,w+k2/2);
k4 = h * rky(t+h,w+k3);
2) You use the gravitational parameter u inside rky, but it is never set prior to using it. The simplest fix is probably just to add this line at the front of rky:
u = 3.986e5;
To make it more generic, you could pass u into the function (e.g. either directly or via an anonymous function).
3) You should double check your formula for combining the k's. E.g., double check this line:
w = w + (1/6)*(k1+k2+k3+k4)
E.g., check your above line against the RK4 formula from this link:
4) A simple way to turn a variable into a column vector is to use a (:) subscript. E.g., this line:
X = [X(1); X(2);X(3);X(4);X(5);X(6)];
can be simplified to this:
X = X(:);


Sarah on 21 Oct 2015
Thanks, that was a big help! I am trying to plot the solution, but am running into trouble again. The plot runs but it is blank. Any ideas where I'm going wrong?
for t = 1:10:5400
k1 = h*rky(t,w);
k2 = h*rky(t+h/2,w+k1/2);
k3 = h*rky(t+h/2,w+k2/2);
k4 = h*rky(t+h,w+k3);
w = w + (h/6)*(k1+2*k2+2*k3+k4);
t = t+h;
position = [w(1);w(2);w(3)];
velocity = [w(4);w(5);w(6)];
axis([1 5400 0 10000]);
axis([1 5400 0 5]);
(this is the only part changed from the previous code)
James Tursa
James Tursa on 21 Oct 2015
You are not saving the results of each iteration ... you are simply overwriting the "position" and "velocity" variables at each iteration with the current results. To save them, I would suggest putting them each into one large matrix. E.g.,
Just prior to the for loop:
n = numel(1:10:5400) + 1;
position = zeros(3,n);
velocity = zeros(3,n);
m = 1;
position(:,1) = w(1:3);
velocity(:,1) = w(4:6);
Then inside the loop:
m = m + 1;
position(:,m) = w(1:3);
velocity(:,m) = w(4:6);
Then adjust your plots accordingly since you will have a matrix of positions (each column being at one point in time) and another matrix of velocities.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by