Runge-Kutta - Orbital mechanics application

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
end

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

    end

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.

Answer 1

[James Tursa]

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:

https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

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]

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)];
end

subplot(2,1,1),plot(t,position)
title('position');
xlabel('Time');
ylabel('Displacement');
axis([1 5400 0 10000]);
subplot(2,1,2),plot(t,velocity)
title('velocity');
xlabel('Time');
ylabel('Displacement');
axis([1 5400 0 5]);

(this is the only part changed from the previous code)

[James Tursa]

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.