分别计算一天、一周、半月和一个月的开普勒解析轨道 与积分轨道并进行统计对比的MATLAB代码

由于你没有给定具体的开普勒元素和积分方法，我无法提供完整的代码，但是我可以给你一个大致的框架，你可以根据需要进行调整。

首先，我们需要设置一些初始参数，如引力常数、中心天体的质量、初始时间、时间步长等：

G = 6.6743e-11; % 引力常数
M = 5.97e24; % 地球质量
m = 1000; % 卫星质量
t0 = 0; % 初始时间
dt = 60; % 时间步长，单位为秒

接下来，我们需要计算出初始状态下的位置和速度矢量，以及相应的开普勒元素：

r0 = [7000e3, 0, 0]; % 初始位置矢量，假设在 x 轴上
v0 = [0, sqrt(G*M/norm(r0)), 0]; % 初始速度矢量，假设沿 y 轴正方向
a = G*(M+m)/norm(r0)^2; % 初始加速度大小
e = norm(cross(v0, cross(r0, v0))/G - r0/norm(r0)); % 偏心率
h = cross(r0, v0); % 角动量矢量
i = acos(h(3)/norm(h)); % 轨道倾角
N = cross([0, 0, 1], h); % 升交点矢量
Omega = atan2(N(2), N(1)); % 升交点赤经
w = atan2(cross(h, r0)(3)/norm(h)/norm(r0)*sin(i), dot(N, r0)/norm(N)/norm(r0)*sin(i)); % 近地点幅角
v = acos(dot(h, r0)/norm(h)/norm(r0)); % 真近地点角度
aop = w - v; % 近地点夹角

接下来，我们需要定义一个计算加速度的函数，根据开普勒问题的解析解计算出加速度矢量：

function [a] = kepler_acc(r, m, M)
    G = 6.6743e-11;
    a = -G*M/norm(r)^3*r;
end

然后，我们可以使用 MATLAB 自带的积分函数（如 `ode45`）计算出一天、一周、半月和一个月的积分轨道，并以相同时间步长的方式计算出对应的开普勒解析轨道：

[t1, y1] = ode45(@(t, y) [y(4:6); kepler_acc(y(1:3), m, M)], [t0, t0+86400], [r0, v0]);
[t2, y2] = ode45(@(t, y) [y(4:6); kepler_acc(y(1:3), m, M)], [t0, t0+7*86400], [r0, v0]);
[t3, y3] = ode45(@(t, y) [y(4:6); kepler_acc(y(1:3), m, M)], [t0, t0+15*86400], [r0, v0]);
[t4, y4] = ode45(@(t, y) [y(4:6); kepler_acc(y(1:3), m, M)], [t0, t0+30*86400], [r0, v0]);

a1 = G*(M+m)/2/norm(r0) * (1 - e^2); % 长半轴
T1 = 2*pi*sqrt(a1^3/G/(M+m)); % 周期
t1_kepler = linspace(t0, t0+86400, 1440)';
r1_kepler = zeros(1440, 3);
for i = 1:1440
    E = 2*atan(sqrt((1-e)/(1+e))*tan(v/2));
    M = E - e*sin(E);
    T = sqrt(a1^3/G/(M+m));
    t = T*(i-1)/1440;
    n = 2*pi/T;
    M = n*t + M;
    E = kepler_E(M, e);
    v = 2*atan(sqrt((1+e)/(1-e))*tan(E/2));
    r1_kepler(i,:) = a1*(cos(E)-e)*[cos(Omega)*cos(w+v)-sin(Omega)*sin(w+v)*cos(i), sin(Omega)*cos(w+v)+cos(Omega)*sin(w+v)*cos(i), sin(w+v)*sin(i)];
end

a2 = G*(M+m)/2/norm(r0) * (1 - e^2); % 长半轴
T2 = 2*pi*sqrt(a2^3/G/(M+m)); % 周期
t2_kepler = linspace(t0, t0+7*86400, 10080)';
r2_kepler = zeros(10080, 3);
for i = 1:10080
    E = 2*atan(sqrt((1-e)/(1+e))*tan(v/2));
    M = E - e*sin(E);
    T = sqrt(a2^3/G/(M+m));
    t = T*(i-1)/10080;
    n = 2*pi/T;
    M = n*t + M;
    E = kepler_E(M, e);
    v = 2*atan(sqrt((1+e)/(1-e))*tan(E/2));
    r2_kepler(i,:) = a2*(cos(E)-e)*[cos(Omega)*cos(w+v)-sin(Omega)*sin(w+v)*cos(i), sin(Omega)*cos(w+v)+cos(Omega)*sin(w+v)*cos(i), sin(w+v)*sin(i)];
end

a3 = G*(M+m)/2/norm(r0) * (1 - e^2); % 长半轴
T3 = 2*pi*sqrt(a3^3/G/(M+m)); % 周期
t3_kepler = linspace(t0, t0+15*86400, 21600)';
r3_kepler = zeros(21600, 3);
for i = 1:21600
    E = 2*atan(sqrt((1-e)/(1+e))*tan(v/2));
    M = E - e*sin(E);
    T = sqrt(a3^3/G/(M+m));
    t = T*(i-1)/21600;
    n = 2*pi/T;
    M = n*t + M;
    E = kepler_E(M, e);
    v = 2*atan(sqrt((1+e)/(1-e))*tan(E/2));
    r3_kepler(i,:) = a3*(cos(E)-e)*[cos(Omega)*cos(w+v)-sin(Omega)*sin(w+v)*cos(i), sin(Omega)*cos(w+v)+cos(Omega)*sin(w+v)*cos(i), sin(w+v)*sin(i)];
end

a4 = G*(M+m)/2/norm(r0) * (1 - e^2); % 长半轴
T4 = 2*pi*sqrt(a4^3/G/(M+m)); % 周期
t4_kepler = linspace(t0, t0+30*86400, 43200)';
r4_kepler = zeros(43200, 3);
for i = 1:43200
    E = 2*atan(sqrt((1-e)/(1+e))*tan(v/2));
    M = E - e*sin(E);
    T = sqrt(a4^3/G/(M+m));
    t = T*(i-1)/43200;
    n = 2*pi/T;
    M = n*t + M;
    E = kepler_E(M, e);
    v = 2*atan(sqrt((1+e)/(1-e))*tan(E/2));
    r4_kepler(i,:) = a4*(cos(E)-e)*[cos(Omega)*cos(w+v)-sin(Omega)*sin(w+v)*cos(i), sin(Omega)*cos(w+v)+cos(Omega)*sin(w+v)*cos(i), sin(w+v)*sin(i)];
end

注意，在计算开普勒解析轨道时，我们需要使用 Kepler's equation 来求解真近地点角度和偏心率异常角（`v` 和 `M`）。这里我们可以使用以下函数来求解偏心率异常角：

function [E] = kepler_E(M, e)
    E0 = M;
    E1 = M + e*sin(E0);
    while abs(E1 - E0) > 1e-8
        E0 = E1;
        E1 = M + e*sin(E0);
    end
    E = E1;
end

最后，我们可以绘制出开普勒解析轨道和积分轨道，并进行对比统计：

figure;
subplot(2,2,1);
plot(y1(:,1), y1(:,2), 'r-', r1_kepler(:,1), r1_kepler(:,2), 'b--');
title('One day');
legend('Numerical', 'Analytical');
subplot(2,2,2);
plot(y2(:,1), y2(:,2), 'r-', r2_kepler(:,1), r2_kepler(:,2), 'b--');
title('One week');
legend('Numerical', 'Analytical');
subplot(2,2,3);
plot(y3(:,1), y3(:,2), 'r-', r3_kepler(:,1), r3_kepler(:,2), 'b--');
title('Half a month');
legend('Numerical', 'Analytical');
subplot(2,2,4);
plot(y4(:,1), y4(:,2), 'r-', r4_kepler(:,1), r4_kepler(:,2), 'b--');
title('One month');
legend('Numerical', 'Analytical');

这样，我们就可以得到相应的对比统计图了。请注意，上述代码只是一个大致的框架，你需要根据具体问题进行调整和修改。