代碼：

function SunEarthMoon % M函數(shù)文件load planets; % 將planets.mat中的變量mass、position、velocity加載過(guò)來(lái)[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; % 作圖的時(shí)間間隔為一天，每天有24*3600秒 N = length(list); % N=3，三個(gè)天體 mass = mass(list); % N行1列矩陣，N個(gè)天體的質(zhì)量 position = position(list,:); % N行3列矩陣，N個(gè)天體的坐標(biāo)，坐標(biāo)是1行3列的行向量，三個(gè)方向的分量 velocity = velocity(list,:); % N行3列矩陣，N個(gè)天體的速度，速度是1行3列的行向量，三個(gè)方向的分量 h = plotplanets(position); %作圖函數(shù)for t = 1:365 % 圖中總時(shí)間為一年，一年365天plotplanets(position,h); % force = zeros(N,3); % N行3列零矩陣，一行表示某個(gè)天體在三個(gè)方向上的受力for i = 1 : N % 遍歷計(jì)算各天體間的萬(wàn)有引力。組合數(shù)C（3，2）Pi = position(i,:); % 某天體坐標(biāo)Mi = mass(i); % 某天體質(zhì)量for j = (i+1):N %the i+1 is to create diagonal Mj = mass(j); % 另一天體質(zhì)量Pj = position(j,:); % 另一天體坐標(biāo)dr = Pj - Pi; % 兩天體的相對(duì)，1行3列矩陣forceij = G*Mi*Mj./(norm(dr).^3).*dr; % 兩天體之間的力，1行3列的向量force(i,:) = force(i,:) + forceij; % 規(guī)定正方向，將力計(jì)算進(jìn)矩陣force(j,:) = force(j,:) - forceij; % 反作用力與作用力方向相反，將力計(jì)算進(jìn)矩陣% 上兩行可替換為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

效果圖

總結(jié)

以上是生活随笔為你收集整理的多体问题的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。