Runge-Kutta - Orbital mechanics application
3 次查看(过去 30 天)
显示 更早的评论
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.
0 个评论
回答(1 个)
James Tursa
2015-10-21
编辑:James Tursa
2015-10-21
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(:);
2 个评论
James Tursa
2015-10-21
编辑:James Tursa
2015-10-21
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Control System Toolbox 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!