algorithm - 提高 n 体运动的准确性

Question

我正在尝试用牛顿运动定律和万有引力来解决 n 体问题，以获得速度分布。

我已经尝试使用 RK4，然后使用 leapfrog 算法，同时软化势能，但我无法获得足够准确的结果，因为如果两个物体靠得太近，数字会上升，能量会变为数量巨大，因此不再保留，使分发无用。

似乎唯一有效的技巧是减少时间步长。然而，这会带来巨大的运行时间代价。

我正在寻找改进代码的方法，以便保留能量并使程序运行时间合理。相关部分是:(请注意:所有质量都相等，G 设置为 1，我需要将较大的时间间隔与尽可能大的时间步相结合以获得良好的结果。代码中的变量 e 是潜在的软化剂，是设置 0.005)

如果您对时间步/e/其他常量的良好值也有建议，我们将不胜感激。

for i=1:(length(t)-1)

for a1=1:n
    for a2=1:n
        if a1~=a2
     r(a1,a2)=sqrt((qstar(a1,i)-qstar(a2,i))^2+(qstar(a1+n,i)-qstar(a2+n,i))^2+(qstar(a1+2*n,i)-qstar(a2+2*n,i))^2); %Distances

     b(a1,a2)=1/(r(a1,a2)+e)^3; %Force component

     PE(1,i)=PE(1,i)+1/r(a1,a2); %Potential Energy

         A2(a1,a1)=A2(a1,a1)-b(a1,a2);      %Forming the ODE Matrix
         A2(a1,a2)=b(a1,a2);
         A3(a1,a1+n)=A3(a1,a1+n)-b(a1,a2);
         A3(a1,a2+n)=b(a1,a2);
         A4(a1,a1+2*n)=A4(a1,a1+2*n)-b(a1,a2);
         A4(a1,a2+2*n)=b(a1,a2);

        end
    end
    V(a1,i)=sqrt((qstar(a1+3*n,i))^2+(qstar(a1+4*n,i))^2+(qstar(a1+5*n,i))^2); %velocity
end

   PE(1,i)=0.5*PE(1,i);

    A=vertcat(A1,A2,A3,A4);

    if i==1 %RK4
         k1 = A*qstar(:,i);  % Approx for y gives approx for deriv
  y1 = qstar(:,i)+k1*h/2;      % Intermediate value (using k1)

  k2 = A*y1;        % Approx deriv at intermediate value.
  y2 = qstar(:,i)+k2*h/2;      % Intermediate value (using k2)

  k3 = A*y2;        % Another approx deriv at intermediate value.
  y3 = qstar(:,i)+k3*h;        % Endpoint value (using k3)

  k4 = A*y3;        % Approx deriv at endpoint value.

  qstar(:,i+1) = qstar(:,i) + (k1+2*k2+2*k3+k4)*h/6; % Approx soln
    else %leapfrog

  acce = A*qstar;  % Approx for accelaration
  qstar(:,i+1) = qstar(:,i-1) + 2*h*acce(:,i); % Approx soln
    end



        A2=zeros(n,6*n);
        A3=zeros(n,6*n);
        A4=zeros(n,6*n);

        end

我构建了一个矩阵，乘以 qstar 将给出运动定律。qstar 存储位置和速度，它是这样的:

x1,x2,...,xn,y1,...,yn,z1,...,zn,px1,...,pxn,py1,...,pyn,pz1,... zn

只有垂直的。 x1 表示第 1 号物体在 x 轴上的位置，例如 py4 表示第 4 号物体沿 y 轴的动量。

Answer 1

你应该使用 symplectic integrator ，因为它会近似保存能量。