matlab求解二体问题。设质量为m1=0.3,m2=0.03.画出初始条件(x1(0),y1(0))=(2,2), , (x2(0),y2(0))=(0,0), 下两天体的运动轨迹。

在MATLAB中，要模拟二体问题（如行星与太阳系统），我们可以使用经典力学中的牛顿运动定律和欧拉法（Euler's method）来近似解决微分方程。对于给定的质量m1=0.3（主天体）和m2=0.03（卫星），以及初始位置 x1(0)=2, y1(0)=2 和 x2(0)=0, y2(0)=0，你可以按照以下步骤操作：

1. 定义常量：

   G = 6.6743e-11; % 引力常数
   m1 = 0.3;
   m2 = 0.03;
   r01 = [2; 2]; % 初始位置向量 for m1
   r02 = [0; 0]; % 初始位置向量 for m2

2. 定义时间步长和总时间范围：

   dt = 0.01; % 时间间隔
   tspan = [0, 2*86400]; % 2天的时间范围 (秒)

3. 创建向量表示动量和速度：

   p1_0 = zeros(2, 1); % 初始动量 m1
   v1_0 = [2; 2]; % 初始速度 m1 (假设初速度与位置同向)

   p2_0 = zeros(2, 1); % 初始动量 m2
   v2_0 = [0; 0]; % 初始速度 m2 (同样假设)

4. 定义函数来计算作用力：

   function force = gravitational_force(r1, r2)
       r = r1 - r2; % 相对位移
       force = G * m1 * m2 / norm(r)^3 * r;
   end

5. 使用Euler方法迭代更新位置：

   Nsteps = floor(tspan(2) / dt);
   trajectory = cell(2, Nsteps+1);
   trajectory{:, 1} = [r1_0; r2_0];

   for i = 2:Nsteps+1
       % 更新速度
       v1_i = v1_{i-1} + dt * gravitational_force(trajectory{i-1, :}, trajectory{i-1, :}, trajectory{i-1, :}) ./ m2;

       % 更新位置
       r1_i = r1_{i-1} + dt * v1_i;
       r2_i = r2_{i-1} + dt * v2_i;

       % 存储结果
       trajectory{:, i} = [r1_i; r2_i];
   end

6. 绘制轨迹：

   plot(trajectory{1, :}, trajectory{2, :}, 'o-');
   hold on;
   plot(trajectory{1, 2:end}, trajectory{2, 2:end});
   legend('m1', 'm2');
   xlabel('X-axis');
   ylabel('Y-axis');
   title('Two-body problem motion trajectory over 2 days');
   hold off;

这将生成两个质点在2天内的运动轨迹。注意，在实际应用中，你可能需要使用更精确的数值积分方法（如四阶龙格-库塔法RKF45），特别是当时间步长较大时。