如何寫出三體的MATLAB程序-代碼篇
寫在前面
在上文當(dāng)中我們已經(jīng)對三個物體之間的受力進(jìn)行了分析，也說明了在時間下的加速度、速度和位移的計算方式。
本篇中將根據(jù)上一篇的公式來寫出對應(yīng)的代碼，并且詳細(xì)說明一下如何去構(gòu)建一個程序的框架。
本文所有代碼均在我的Github中存有備份，可下載后直接運(yùn)行，點(diǎn)擊Github: HanpuLiang/Three-Body-by-MATLAB即可進(jìn)入。
構(gòu)建框架
基本變量
我們首先要確定，物體本身具有哪些物理量？
質(zhì)量、加速度、速度、坐標(biāo)。
其中加速度和坐標(biāo)為矢量，當(dāng)在計算過程當(dāng)中可以將其正交分解為與
軸的標(biāo)量。
其余的量我們還可以設(shè)置一下，諸如：物體運(yùn)動的時間長度、我們計算所需要的時間間隔、萬有引力常數(shù)等。
對于我們要做出的圖像大小也需要設(shè)置一下，比如說軸范圍，
軸范圍等。
程序流程
初始化
首先需要初始化，確定三個物體在初始狀態(tài)的情況:初始坐標(biāo)、初始速度(大小與方向)。所以一共需要三個量：坐標(biāo)、速度大小、速度相對軸角度。
物體的加速度可以直接由萬有引力公式計算出來，為了計算方便，需要將其正交分解與疊加。
隨著時間演變
物體開始運(yùn)動了，但是因?yàn)槲覀儫o法給出一段連續(xù)的時間，只能計算在極小的時間間隔之后物體所在的位置。
所以在時間，首先計算出當(dāng)前位置的加速度，然后根據(jù)這個加速度算出當(dāng)前的速度，再根據(jù)這個速度算出經(jīng)歷過
時間后的位移變化量，再將這個位移變化量疊加到
的位置上。
這樣子就完成了一個循環(huán)。
作圖
之前按道理，我們應(yīng)該將每一個時間點(diǎn)的值放在一個矩陣內(nèi)，這樣子我們就可以得到隨之間變化的所有物理量。
這樣子我們就可以直接做出隨著時間變化的各個物理量的圖，如和
以及
等。
如果我們想要做出小球的運(yùn)動圖，那就需要時刻及其之前(做出尾跡)的數(shù)據(jù)進(jìn)行作圖。
實(shí)際代碼
基本變量-代碼
首先是初始化
%% 初始條件
% 初始條件為以圓心為(0, 0)半徑r的圓上有三個等質(zhì)量的點(diǎn)
r = 10;
% 坐標(biāo)(等邊三角形)
pos0 = [0, r; r/2*sqrt(3), -r/2; -r/2*sqrt(3), -r/2];
% 速度大小
v0 = [6, 6, 6];
% 速度方向(x軸正方向?yàn)閰⒖?
theta0 = [0, 4*pi/3, 2*pi/3];
%% 參數(shù)設(shè)置
global G dt m x_min x_max y_min y_max time_end isOutVideo;
% 結(jié)束時間
time_end = 2;
% 時間間隔
dt = 0.05;
% 萬有引力系數(shù)，隨便設(shè)置的
G = 1;
% 質(zhì)量
m = [1000, 1000, 1000];
% 小球個數(shù)
N = length(v0);
% 圖像邊界
x_min = -25;
x_max = -x_min;
y_min = x_min;
y_max = -y_min;
% 是否輸出視頻圖像
isOutVideo = true;
初始化-代碼
然后是將我們的初始值放入矩陣中，因?yàn)槲覀兂跏贾翟O(shè)定的是角度與速度大小，所以首先要把給分解為
軸上的。
%% 初始設(shè)置
time = 0:dt:time_end;
% 坐標(biāo)
pos = zeros(N, 2, length(time));
pos(:,:,1) = pos0;
% 速度
vx = zeros(length(time), N);
vx(1,:) = v0.*cos(theta0);
vy = zeros(length(time), N);
vy(1,:) = v0.*sin(theta0);
% 加速度大小
a = zeros(length(time), N);
迭代開始
迭代的話這一步其實(shí)就和我們的邏輯很像了，不過之所以主代碼這么簡介，是因?yàn)槲野岩淮蠖褟?fù)雜的內(nèi)容全部放到了函數(shù)里面，只留個接口調(diào)用，這樣子可以讓主代碼更加簡潔明了。
%% 迭代開始
for t = 2:length(time)
% 得到分加速度
da = calAcceleration(pos(:,:,t-1));
% 計算速度
[vx(t,:), vy(t,:)] = calVelocity(vx(t-1,:), vy(t-1,:), da);
% 計算位移
pos(:,:,t) = calDisplacement(vx(t-1:t,:), vy(t-1:t,:), pos(:,:,t-1));
end
對于計算加速度的函數(shù)，主要的原理還是和上一篇講的一樣，通過萬有引力公式求解，然后正交分解并疊加。
% 計算x與y軸加速度變化量da(3x2)
function da = calAcceleration(pos)
global m;
% 小球數(shù)量
[n, ~] = size(pos);
da = zeros(n, 2);
for i = 1:n
dax = 0;
day = 0;
for j = 1:n
if i ~= j
% i小球和j小球相對角度與距離
[theta, r] = calRelatively(pos(i,:), pos(j,:));
% 兩個小球的引力大小
F = calGravitation(r, i, j);
% 第i個小球收到來自j的加速度分量
dax = dax + F/m(i)*cos(theta);
day = day + F/m(i)*sin(theta);
end
end
da(i,:) = [dax, day];
end
end
% 計算兩個小球的相對角度與距離
function [theta, r] = calRelatively(pos1, pos2)
dx = pos2(1) - pos1(1);
dy = pos2(2) - pos1(2);
r = sqrt(dx^2 + dy^2);
theta = acos(dx/r);
% 因?yàn)閏os值的兩個象限需要區(qū)分，所以這里要變換
if dy < 0 && dx >0
theta = -theta;
end
if dx < 0 && dy < 0
theta = (pi-theta)+pi;
end
end
% 計算兩個小球引力大小
function F = calGravitation(r, i, j)
global m G;
F = G*m(i)*m(j)/r^2;
end
計算速度的函數(shù)，這個就很簡單了，簡單的速度與加速度公式。
% 計算小球的速度變化
function [vx, vy] = calVelocity(vx_p, vy_p, da)
global dt;
vx = vx_p + dt*da(:,1)';
vy = vy_p + dt*da(:,2)';
end
計算當(dāng)前坐標(biāo)，方法同公式。
% 計算小球的位移變化
function pos = calDisplacement(vx, vy, pos_p)
global dt;
vx = vx';
vy = vy';
% 計算下一時刻的坐標(biāo)
pos(:,1) = pos_p(:,1) + (vx(:,1)+vx(:,2))/2*dt;
pos(:,2) = pos_p(:,2) + (vy(:,1)+vy(:,2))/2*dt;
end
作圖-代碼
作圖的話就沒有這么簡單了，因?yàn)檫€需要設(shè)置一大堆比較麻煩的圖像參數(shù)。
對于做軌跡圖，可以通過以下代碼實(shí)現(xiàn)
% 做軌跡圖像
plotPosition(pos, vx, vy, time)
% 做運(yùn)動圖像并保存視頻
function plotPosition(pos, vx, vy, time)
global isOutVideo;
figure(1)
if isOutVideo == true
% 創(chuàng)建視頻句柄，視頻名稱three_body.avi
writerObj = VideoWriter('three_body.avi');
open(writerObj);
myMovie(1:length(time)) = struct('cdata', [], 'colormap', []);
end
% 迭代計算得到圖像
for t = 1:length(time)
plotPosVec(pos(:,:,t), vx(t,:), vy(t,:), t, pos)
if isOutVideo == true
frame = getframe;
% 修改幀參數(shù)
frame.cdata = imresize(frame.cdata, [685, 685]);
writeVideo(writerObj, frame);
end
end
% 關(guān)閉視頻句柄
if isOutVideo == true
close(writerObj);
end
end
% 作圖位置+速度矢量
function plotPosVec(pos, vx, vy, t, pos_all)
% 小球
global x_min x_max y_min y_max;
figure(1)
scatter(pos(:,1), pos(:,2), 'ok', 'filled')
% 圖像細(xì)節(jié)調(diào)整
axis equal
box on
grid on
set(gca, 'linewidth', 1.5, 'xtick', floor(linspace(x_min, x_max, 11)), 'ytick', floor(linspace(y_min, y_max, 11)))
hold on
% 三條速度線
for i = 1:length(vx)
line([pos(i,1) pos(i,1)+vx(i)/2], [pos(i,2), pos(i,2)+vy(i)/2], 'linewidth', 1.2)
end
% 添加軌道線
plotCurrentTrace(pos_all, t)
% 添加文本
text(x_max*13/25, y_min*20/25, 'Made By Liang Hanpu', 'horiz', 'center', 'color', 'r')
axis([x_min x_max y_min y_max])
hold off
end
% 輸出軌跡線
function plotCurrentTrace(pos, t)
global x_min x_max y_min y_max;
if t ~= 0
[a, ~, ~] = size(pos);
hold on
axis equal
box on
grid on
set(gca, 'linewidth', 1.5)
axis([x_min x_max y_min y_max])
for i = 1:a
x = zeros(1, t);
y = zeros(1, t);
for j = 1:t
x(j) = pos(i, 1, j);
y(j) = pos(i, 2, j);
end
plot(x, y, 'linewidth', 1.5)
end
end
end
對于做時間隨速度大小與角度的圖像，可以由以下函數(shù)實(shí)現(xiàn)，這個就比較簡單了。
% 做速度隨時間圖像
plotVelocity(vx, vy)
% 輸出三個小球速度變化圖與角度變化圖
function plotVelocity(vx, vy)
global dt time_end;
% 速度
v = sqrt(vx.^2 + vy.^2);
t = (0:dt:time_end)'*ones(1,3);
% 角度
theta = acos(vx./v);
theta(vx>0&vy<0) = 2*pi-theta(vx>0&vy<0);
theta(vx<0&vy<0) = (pi-theta(vx<0&vy<0))+pi;
figure
plot(t, v, 'linewidth', 1.2)
box on
xlabel('Time', 'fontsize', 16)
ylabel('Velocity' , 'fontsize', 16)
set(gca, 'linewidth', 1.2)
figure
plot(t,theta, 'linewidth', 1.2)
xlabel('Time', 'fontsize', 16)
ylabel('Angle \theta' , 'fontsize', 16)
set(gca, 'linewidth', 1.2)
end
做出來的圖形大概可以看成這樣子，可以看出來還是有比較有趣的趨勢的。
最后
以上就是全部內(nèi)容，我將會全部放在我的Github中，地址在文章開頭有。
如果你學(xué)到了一些東西的話，麻煩點(diǎn)個贊，加個收藏來個關(guān)注噢。
