如何写出三体的MATLAB程序-代码篇

最新推荐文章于 2023-03-31 16:57:23 发布

Hanpu_Liang

最新推荐文章于 2023-03-31 16:57:23 发布

阅读量2.8k

点赞数 13

分类专栏： MATLAB 文章标签： MATLAB 物理模拟物理模型多体运动

本文链接：https://blog.csdn.net/Hanpu_Liang/article/details/88322141

版权

MATLAB 专栏收录该内容

8 篇文章 18 订阅

订阅专栏

如何写出三体的MATLAB程序-代码篇

写在前面

在上文当中我们已经对三个物体之间的受力进行了分析，也说明了在时间 $t$ 下的加速度、速度和位移的计算方式。

本篇中将根据上一篇的公式来写出对应的代码，并且详细说明一下如何去构建一个程序的框架。

本文所有代码均在我的Github中存有备份，可下载后直接运行，点击Github: HanpuLiang/Three-Body-by-MATLAB即可进入。

构建框架

基本变量

我们首先要确定，物体本身具有哪些物理量？

质量、加速度、速度、坐标。

其中加速度和坐标为矢量，当在计算过程当中可以将其正交分解为 $x$ 与 $y$ 轴的标量。

其余的量我们还可以设置一下，诸如：物体运动的时间长度、我们计算所需要的时间间隔、万有引力常数等。

对于我们要做出的图像大小也需要设置一下，比如说 $x$ 轴范围， $y$ 轴范围等。

程序流程

初始化

首先需要初始化，确定三个物体在初始状态的情况:初始坐标、初始速度(大小与方向)。所以一共需要三个量：坐标、速度大小、速度相对 $x$ 轴角度。

物体的加速度可以直接由万有引力公式计算出来，为了计算方便，需要将其正交分解与叠加。

随着时间演变

物体开始运动了，但是因为我们无法给出一段连续的时间，只能计算在极小的时间间隔 $\Delta t$ 之后物体所在的位置。

所以在 $t\to t+Delta t$ 时间，首先计算出当前位置的加速度，然后根据这个加速度算出当前的速度，再根据这个速度算出经历过 $\Delta t$ 时间后的位移变化量，再将这个位移变化量叠加到 $t$ 的位置上。

这样子就完成了一个循环。

可以将其写成一个流程图

作图

之前按道理，我们应该将每一个时间点的值放在一个矩阵内，这样子我们就可以得到随之间变化的所有物理量。

这样子我们就可以直接做出随着时间变化的各个物理量的图，如 $t - v$ 和 $t - a$ 以及 $t-\theta$ 等。

如果我们想要做出小球的运动图，那就需要 $t$ 时刻及其之前(做出尾迹)的数据进行作图。

实际代码

基本变量-代码

首先是初始化

%% 初始条件
% 初始条件为以圆心为(0, 0)半径r的圆上有三个等质量的点
r = 10;
% 坐标(等边三角形)
pos0 = [0, r; r/2*sqrt(3), -r/2; -r/2*sqrt(3), -r/2];
% 速度大小
v0 = [6, 6, 6];
% 速度方向(x轴正方向为参考)
theta0 = [0, 4*pi/3, 2*pi/3];

%% 参数设置
global G dt m x_min x_max y_min y_max time_end isOutVideo;
% 结束时间
time_end = 2;
% 时间间隔
dt = 0.05;
% 万有引力系数，随便设置的
G = 1;
% 质量
m = [1000, 1000, 1000];
% 小球个数
N = length(v0);
% 图像边界
x_min = -25;
x_max = -x_min;
y_min = x_min;
y_max = -y_min;
% 是否输出视频图像
isOutVideo = true;

初始化-代码

然后是将我们的初始值放入矩阵中，因为我们初始值设定的是角度与速度大小，所以首先要把 $v$ 给分解为 $x y$ 轴上的。

%% 初始设置
time = 0:dt:time_end;
% 坐标
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);

迭代开始

迭代的话这一步其实就和我们的逻辑很像了，不过之所以主代码这么简介，是因为我把一大堆复杂的内容全部放到了函数里面，只留个接口调用，这样子可以让主代码更加简洁明了。

%% 迭代开始
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

对于计算加速度的函数，主要的原理还是和上一篇讲的一样，通过万有引力公式求解，然后正交分解并叠加。

% 计算x与y轴加速度变化量da(3x2)
function da = calAcceleration(pos)
    global m;
    % 小球数量
    [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);
    % 因为cos值的两个象限需要区分，所以这里要变换
    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

计算速度的函数，这个就很简单了，简单的速度与加速度公式。

% 计算小球的速度变化
function [vx, vy] = calVelocity(vx_p, vy_p, da)
    global dt;
    vx = vx_p + dt*da(:,1)';
    vy = vy_p + dt*da(:,2)';
end

计算当前坐标，方法同公式。

% 计算小球的位移变化
function pos = calDisplacement(vx, vy, pos_p)
    global dt;
    vx = vx';
    vy = vy';
    % 计算下一时刻的坐标
    pos(:,1) = pos_p(:,1) + (vx(:,1)+vx(:,2))/2*dt;
    pos(:,2) = pos_p(:,2) + (vy(:,1)+vy(:,2))/2*dt;
end

作图-代码

作图的话就没有这么简单了，因为还需要设置一大堆比较麻烦的图像参数。

对于做轨迹图，可以通过以下代码实现

% 做轨迹图像
plotPosition(pos, vx, vy, time)

% 做运动图像并保存视频
function plotPosition(pos, vx, vy, time)
    global isOutVideo;
    figure(1)
    if isOutVideo == true
        % 创建视频句柄，视频名称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;
            % 修改帧参数
            frame.cdata = imresize(frame.cdata, [685, 685]);
            writeVideo(writerObj, frame);
        end
    end
    % 关闭视频句柄
    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')
    % 图像细节调整
    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

对于做时间随速度大小与角度的图像，可以由以下函数实现，这个就比较简单了。

% 做速度随时间图像
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

做出来的图的趋势还是比较有趣的

时间与速度

时间与角度

最后

以上就是全部内容，我将会全部放在我的Github中，地址在文章开头有。

如果你学到了一些东西的话，麻烦点个赞，加个收藏来个关注噢。

Hanpu_Liang

关注

13
点赞
踩
20

收藏

觉得还不错? 一键收藏
0
评论
如何写出三体的MATLAB程序-代码篇

如何写出三体的MATLAB程序-代码篇写在前面在上文当中我们已经对三个物体之间的受力进行了分析，也说明了在时间ttt下的加速度、速度和位移的计算方式。本篇中将根据上一篇的公式来写出对应的代码，并且详细说明一下如何去构建一个程序的框架。本文所有代码均在我的Github中存有备份，可下载后直接运行，点击Github: HanpuLiang/Three-Body-by-MATLAB即可进入。构...
复制链接

扫一扫