利用8阶RK积分公式计算开普勒卫星积分轨道的MATLAB程序

以下是一个使用 8 阶 RK 积分公式计算开普勒卫星轨道的 MATLAB 程序：

function [t, r, v] = kepler_orbit(r0, v0, t0, tf, h)
% Inputs:
% r0: 初始位置向量
% v0: 初始速度向量
% t0: 初始时间
% tf: 结束时间
% h: 步长
% Outputs:
% t: 时间向量
% r: 位置矩阵
% v: 速度矩阵

% 常数
G = 6.67428e-11;
M = 5.97e24;
mu = G*M;

% 初始条件
r1 = r0;
v1 = v0;
t1 = t0;

% 初始化输出变量
t = t0:h:tf;
n = length(t);
r = zeros(3, n);
v = zeros(3, n);
r(:, 1) = r1;
v(:, 1) = v1;

% 迭代
for i = 1:n-1
    k1_r = v1;
    k1_v = -mu/norm(r1)^3*r1;
    k2_r = v1 + h/2*k1_v;
    k2_v = -mu/norm(r1+h/2*k1_r)^3*(r1+h/2*k1_r);
    k3_r = v1 + h/2*k2_v;
    k3_v = -mu/norm(r1+h/2*k2_r)^3*(r1+h/2*k2_r);
    k4_r = v1 + h*k3_v;
    k4_v = -mu/norm(r1+h*k3_r)^3*(r1+h*k3_r);
    k5_r = v1 + h*(2*k1_v+3*k2_v+4*k3_v)/9;
    k5_v = -mu/norm(r1+h*(2*k1_r+3*k2_r+4*k3_r)/9)^3*(r1+h*(2*k1_r+3*k2_r+4*k3_r)/9);
    k6_r = v1 + h*(7*k1_v+6*k2_v+8*k3_v+3*k4_v)/24;
    k6_v = -mu/norm(r1+h*(7*k1_r+6*k2_r+8*k3_r+3*k4_r)/24)^3*(r1+h*(7*k1_r+6*k2_r+8*k3_r+3*k4_r)/24);
    k7_r = v1 + h*(10*k1_v-8*k3_v+8*k4_v-10*k5_v)/27;
    k7_v = -mu/norm(r1+h*(10*k1_r-8*k3_r+8*k4_r-10*k5_r)/27)^3*(r1+h*(10*k1_r-8*k3_r+8*k4_r-10*k5_r)/27);
    k8_r = v1 + h*(28*k1_v-28*k3_v+12*k4_v+48*k5_v-36*k6_v)/(35);
    k8_v = -mu/norm(r1+h*(28*k1_r-28*k3_r+12*k4_r+48*k5_r-36*k6_r)/(35))^3*(r1+h*(28*k1_r-28*k3_r+12*k4_r+48*k5_r-36*k6_r)/(35));
    r1 = r1 + h/90*(7*k1_r+32*k3_r+12*k4_r+32*k5_r+7*k6_r+12*k7_r+32*k8_r);
    v1 = v1 + h/90*(7*k1_v+32*k3_v+12*k4_v+32*k5_v+7*k6_v+12*k7_v+32*k8_v);
    r(:, i+1) = r1;
    v(:, i+1) = v1;
    t1 = t1 + h;
end

end

这个函数与前面的函数很相似，但使用了更高阶的 RK 积分公式。它也需要初始位置和速度向量，初始和结束时间以及步长作为输入，并返回轨道的时间、位置和速度向量作为输出。你可以像这样调用这个函数：

r0 = [10000; 0; 0];
v0 = [0; 8000; 0];
t0 = 0;
tf = 10000;
h = 10;
[t, r, v] = kepler_orbit(r0, v0, t0, tf, h);

请注意，这只是一个简单的例子，不包括任何误差控制或其他特殊功能。