代碼：

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 constantdt = 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];endendvelocity = velocity + force ./ repmat(mass,1,3)*dt; % v=v+a*dt a=F/mposition = 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==1hold on; axis imageaxis( [-2 2 -2 2]*1e11 );for i = 1:total_planetsif 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);elseh(i) = plot(pos(i,1), pos(i,2), 'k.', 'markersize', 20);plot(pos(i,1), pos(i,2), 'k.', 'markersize',5);endend elsefor i = 1:total_planetsset(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);elseplot(pos(i,1), pos(i,2), 'k.', 'markersize',5);endenddrawnow end

效果圖

總結

以上是生活随笔為你收集整理的多体问题的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。