【运动学】三体运动模拟【含GUI Matlab源码 871期】

获取代码方式1： 完整代码已上传我的资源：【运动学】基于matlab GUI三体运动模拟【含Matlab源码 871期】 点击上面蓝色字体，直接付费下载，即可。

获取代码方式2： 付费专栏Matlab物理应用（初级版）

备注： 点击上面蓝色字体付费专栏Matlab物理应用（初级版），扫描上面二维码，付费29.9元订阅海神之光博客付费专栏Matlab物理应用（初级版），凭支付凭证，私信博主，可免费获得1份本博客上传CSDN资源代码（有效期为订阅日起，三天内有效）； 点击CSDN资源下载链接：1份本博客上传CSDN资源代码

本文基于MATLAB及其GUI界面设计了一个基于经典力学的三体运动数值模拟软件，旨在建立经典力学框架内的空间三质点运动模型，又名为三体运动模型。软件根据当前的质点初始运动参数，运用数值模拟，迭代计算出后续每一时刻各个质点的运动参数，并将计算结果实时显示出来。本软件可用于质点力学与基础天体物理学的自主学习、教学演示和相关领域的科学研究。

基本原理及思路 设三维自由空间中有三个质量、初始位置和初始速度已知的三个质点，分别记作： A(M1,X1,Y1,Z1,U1,V1,W1); B(M2,X2,Y2,Z2,U2,V2,W2); C(M3,X3,Y3,Z3,U3,V3,W3); 其中，X、Y、Z分别表示各个质点在X、Y、Z方向上的坐标，U、V、W分别表示各个质点在X、Y、Z方向上的速度分量。 设fAB表示质点A对质点B的作用力，fBA表示质点B对质点A的反作用力，则三个质点之间的万有引力可表示为f12、f23、f31、f31、f32、f21，根据牛顿第三定律，有： f12=-f21; f23=-f32; f31=-f13; 这样，牛顿第三定律将力的变量由6个减少到3个。 根据牛顿第二定律，质点间相互作用力f12、f23、f31可分别被表达为： f12=gM1M2/R12^3*[X1-X2,Y1-Y2,Z1-Z2]; f23=gM2M3/R23^3*[X2-X3,Y2-Y3,Z2-Z3]; f31=gM3M1/R13^3*[X3-X1,Y3-Y1,Z3-Z1]; 其中g是引力常量，取值为6.67x10-11(N·m2 /kg^2)，R12、R23、R13分别表示两两质点之间的距离，可以表达为： R12=sqrt((X1-X2)2+(Y1-Y2)2+(Z1-Z2)^2); R13=sqrt((X1-X3)2+(Y1-Y3)2+(Z1-Z3)^2); R23=sqrt((X2-X3)2+(Y2-Y3)2+(Z2-Z3)^2); 注意，式中f12、f23、f31表达为三维矢量。 由于三体问题是一个发生在三维自由空间的动力学问题，所以势必要进行矢量分析。本实验建立空间直角坐标系，根据矢量的叠加性原理，将所有三维矢量包括引力、加速度、速度和位置矢量等分别分解到X、Y、Z三个坐标轴上分析，就将矢量运算简化为标量运算。现在，以X轴为例，对三体系统进行经典动力学分析。

设三质点在X轴上的加速度分别为Ax1、Ax2、Ax3，则可以表达为： Ax1=(-f12(1)+f31(1))/M1; Ax2=( f12(1)-f23(1))/M2; Ax3=( f23(1)-f31(1))/M3;

速度可以表达为： U1=U1+Ax1t; U2=U2+Ax2t; U3=U3+Ax3*t;

位置矢量可以表达为： X1=X1+U1t+1/2Ax1t^2; X2=X2+U2t+1/2Ax2t^2; X3=X3+U3t+1/2Ax3*t^2;

其中f12(1)意为三维矢量f12在X轴向上的分量，t为迭代的时间间隔，本实验中t的取值为0.00001。 以此类推，在Y轴和Z轴上也进行类似的分析和计算，就可以在三维自由空间中建立一个简单的三体运动模型。

% %% function varargout = Threebody(varargin) % THREEBODY MATLAB code for Threebody.fig % THREEBODY, by itself, creates a new THREEBODY or raises the existing % singleton*. % % H = THREEBODY returns the handle to a new THREEBODY or the handle to % the existing singleton*. % % THREEBODY(‘CALLBACK’,hObject,eventData,handles,…) calls the local % function named CALLBACK in THREEBODY.M with the given input arguments. % % THREEBODY(‘Property’,‘Value’,…) creates a new THREEBODY or raises the % existing singleton*. Starting from the left, property value pairs are % applied to the GUI before Threebody_OpeningFcn gets called. An % unrecognized property name or invalid value makes property application % pause. All inputs are passed to Threebody_OpeningFcn via varargin. % % *See GUI Options on GUIDE’s Tools menu. Choose “GUI allows only one % instance to run (singleton)”. % % See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to menu_help Threebody

% Last Modified by GUIDE v2.5 28-Oct-2017 15:50:30

% Begin initialization code - DO NOT EDIT gui_Singleton = 1; gui_State = struct(‘gui_Name’, mfilename, … ‘gui_Singleton’, gui_Singleton, … ‘gui_OpeningFcn’, @Threebody_OpeningFcn, … ‘gui_OutputFcn’, @Threebody_OutputFcn, … ‘gui_LayoutFcn’, [] , … ‘gui_Callback’, []); if nargin && ischar(varargin{1}) gui_State.gui_Callback = str2func(varargin{1}); end

if nargout [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); else gui_mainfcn(gui_State, varargin{:}); end % End initialization code - DO NOT EDIT

%% % — Executes just before Threebody is made visible. function Threebody_OpeningFcn(hObject, eventdata, handles, varargin) %% % This function has no output args, see OutputFcn. % hObject handle to figure % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) % varargin command line arguments to Threebody (see VARARGIN)

%初始化天体的质量坐标和速度 handles.Star=Model_Init_1(); handles.type=1; handles.pretype=handles.type; %运行标志位，判断是否在正在运行 global Is_Running; %初始化标志位，判断是否在已初始化 global Is_Init; %默认状态为已初始化且正在运行 Is_Running=true; Is_Init=true; handles.output = hObject; % Update handles structure guidata(hObject, handles);

% UIWAIT makes Threebody wait for user response (see UIRESUME) % uiwait(handles.figure1);

%% % — Outputs from this function are returned to the command line. function varargout = Threebody_OutputFcn(hObject, eventdata, handles) %% % varargout cell array for returning output args (see VARARGOUT); % hObject handle to figure % eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure varargout{1} = handles.output; %持续运行run()函数，等待中断 run();

%% %星体运动函数 function run() %% global Is_Running; global Is_Init; %万有引力常数 G=6.67*10^(-11); %以0.00001的时间步长迭代 T=0.00001; %注意死循环 while (1) handles=guidata(gcf); %退出 if(isempty(handles)) close all; return; end %已初始化意为一个阶段内的运算结束，将结果显示 if(Is_Init==true) %获取参数 star=handles.Star; %--------------------一号星体-------------------- M1=star(1).M; X1=star(1).X; Y1=star(1).Y; Z1=star(1).Z; U1=star(1).U; V1=star(1).V; W1=star(1).W; %--------------------二号星体-------------------- M2=star(2).M; X2=star(2).X; Y2=star(2).Y; Z2=star(2).Z; U2=star(2).U; V2=star(2).V; W2=star(2).W; %--------------------三号星体-------------------- M3=star(3).M; X3=star(3).X; Y3=star(3).Y; Z3=star(3).Z; U3=star(3).U; V3=star(3).V; W3=star(3).W; %清空 cla; axis([-1000 1010 -1000 1000 -1000 1000]); %不采用默认坐标样式，后面自己画线 axis off; %默认视角下的右手系 set(gca,‘XDir’,‘reverse’); set(gca,‘YDir’,‘reverse’); %画网格 x=-1000:125:1000; y=-1000:125:1000; for i=1:length(x) line([x(i) x(i)],[1000 -1000],[0 0],‘Color’,[0.1 0.1 0.1],‘LineWidth’,0.5); line([1000 -1000],[y(i) y(i)],[0 0],‘Color’,[0.1 0.1 0.1],‘LineWidth’,0.5); end %画箭头 line([1000 0],[0 0],[0 0],‘Color’,[1 0 1],‘LineWidth’,2); line([0 0],[1000 0],[0 0],‘Color’,[1 0 1],‘LineWidth’,2); line([0 0],[0 0],[1000 0],‘Color’,[1 0 1],‘LineWidth’,2); %标记坐标轴 text(1050,0,0,‘X’,‘Color’,[1 0 1],‘FontSize’,12); text(0,1050,0,‘Y’,‘Color’,[1 0 1],‘FontSize’,12); text(0,0,1050,‘Z’,‘Color’,[1 0 1],‘FontSize’,12); %标记单位长度 text(125,0,0,‘125’,‘Color’,[1 0 1],‘FontSize’,12); text(0,125,0,‘125’,‘Color’,[1 0 1],‘FontSize’,12);