多体问题

问题及分析

Matlab代码

给出了核心逻辑的注释，作图暂时还不太了解。

function SunEarthMoon   % M函数文件

load planets;   % 将planets.mat中的变量mass、position、velocity加载过来

[sun, earth, moon] = deal(18, 3, 25);   % sun、earth、moon分别是18、3、25行
list = [sun, earth, moon];  % 1行3列矩阵
G = 6.67e-11; % gravitational constant

dt = 24*3600;   % 作图的时间间隔为一天，每天有24*3600秒
N = length(list);   % N=3，三个天体
mass = mass(list);  % N行1列矩阵，N个天体的质量
position = position(list,:);    % N行3列矩阵，N个天体的坐标，坐标是1行3列的行向量，三个方向的分量
velocity = velocity(list,:);    % N行3列矩阵，N个天体的速度，速度是1行3列的行向量，三个方向的分量
h = plotplanets(position);  %作图函数

for t = 1:365   % 图中总时间为一年，一年365天
    plotplanets(position,h);    % 
    force = zeros(N,3); % N行3列零矩阵，一行表示某个天体在三个方向上的受力
    for i = 1 : N   % 遍历计算各天体间的万有引力。组合数C（3，2）
        Pi = position(i,:); % 某天体坐标
        Mi = mass(i);   % 某天体质量
        for j = (i+1):N     %the i+1 is to create diagonal 
            Mj = mass(j);   % 另一天体质量
            Pj = position(j,:); % 另一天体坐标
            dr =  Pj - Pi;  % 两天体的相对，1行3列矩阵
            forceij = G*Mi*Mj./(norm(dr).^3).*dr;   % 两天体之间的力，1行3列的向量
            force(i,:) = force(i,:) + forceij;  % 规定正方向，将力计算进矩阵
            force(j,:) = force(j,:) - forceij;  % 反作用力与作用力方向相反，将力计算进矩阵
            % 上两行可替换为force([i,j],:) = force([i,j],:)+[forceij; -forceij];
        end
    end
    velocity = velocity + force ./ repmat(mass,1,3)*dt;   % v=v+a*dt a=F/m
    position = position + velocity*dt;  % r=r+v*dt
end 

% -------------------------------------------------------------------------

function h = plotplanets(pos,h) % 
scale = 50;
total_planets = size(pos,1);
[sun, earth, moon] = deal(1, 2, 3);
radius = [50, 30, 20];
marker = {'.r', 'b.','m.'};
pos(moon,:) = pos(earth,:) + scale*(pos(moon,:)-pos(earth,:));
if nargin==1
    hold on; axis image
    axis( [-2 2 -2 2]*1e11 );
    for i = 1:total_planets
        if any(i == [sun, earth, moon])
            h(i) = plot(pos(i,1),pos(i,2),marker{i},'markersize',radius(i));
            plot(pos(i,1), pos(i,2), marker{i}, 'markersize',5);
        else
            h(i) = plot(pos(i,1), pos(i,2), 'k.', 'markersize', 20);
            plot(pos(i,1), pos(i,2), 'k.', 'markersize',5);
        end
    end
else
    for i = 1:total_planets
        set(h(i), 'Xdata', pos(i,1)  , 'Ydata', pos(i,2)  )
        if any(i == [sun, earth, moon])
            plot(pos(i,1), pos(i,2), marker{i}, 'markersize',5);
        else
            plot(pos(i,1), pos(i,2), 'k.', 'markersize',5);
        end
    end
    drawnow
end

结果

作者：@臭咸鱼

本文为作者原创，转载请注明出处：https://chouxianyu.github.io

欢迎讨论交流!