A list,
In this paper, a three-body motion numerical simulation software based on classical mechanics is designed based on MATLAB and its GUI interface, aiming at establishing the three-particle motion model of space within the framework of classical mechanics, also known as the three-body motion model. According to the initial motion parameters of the particle, the software iteratively calculates the motion parameters of each particle at each subsequent moment by numerical simulation, and displays the calculation results in real time. The software can be used for self-study, teaching demonstration and scientific research in particle mechanics and basic astrophysics.
Suppose there are three particles in three-dimensional free space whose mass, initial position and initial velocity are known, and they are denoted as A(M1,X1,Y1,Z1,U1,V1,W1); B(M2,X2,Y2,Z2,U2,V2,W2); C(M3,X3,Y3,Z3,U3,V3,W3);
Where, X, Y, and Z respectively represent the coordinates of each particle in the X, Y, and Z directions, and U, V, and W respectively represent the velocity components of each particle in the X, Y, and Z directions.
Let fAB represent the force of particle A on particle B, and fBA represent the reaction force of particle B on particle A, then the gravitational attraction between the three particles can be expressed as F12, F23, F31, F31, F32, and F21. According to Newton’s third law, it can be written:
f12=-f21;
f23=-f32;
f31=-f13;
Thus, Newton’s third law reduces the number of variables of force from six to three. According to Newton’s second law, the interaction forces between particles F12, F23 and F31 can be expressed as:
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];
Where, g is the gravitational constant, with a value of 6.67×10-11(N·m2 /kg^2), and R12, R23 and R13 respectively represent the distances between two particles, which can be expressed as:
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);
Note that f12, F23 and F31 are expressed as three-dimensional vectors.
Since the three-body problem is a dynamic problem occurring in three dimensional free space, it is necessary to carry out vector analysis. In this experiment, the space cartesian coordinate system is established. According to the principle of superposition of vectors, all three-dimensional vectors including gravity, acceleration, velocity and position vectors are decomposed into X, Y and Z coordinate axes for analysis, and the vector operation is simplified to scalar operation. Now, taking the X axis as an example, the classical dynamics analysis of the three-body system is carried out.
Suppose that the acceleration of the three particles on the X-axis is Ax1, Ax2 and Ax3 respectively, then it can be expressed as:
Ax1=(-f12(1)+f31(1))/M1; Ax2=( f12(1)-f23(1))/M2; Ax3=( f23(1)-f31(1))/M3; Speed can be expressed as:
U1=U1+Ax1t; U2=U2+Ax2t; U3=U3+Ax3*t; The position vector can be expressed as:
X1=X1+U1t+1/2Ax1t^2;
X2=X2+U2t+1/2Ax2t^2;
X3=X3+U3t+1/2Ax3*t^2;
Where, F12 (1) refers to the component of three-dimensional vector F12 upward on the X-axis, and t is the time interval of iteration. The value of t in this experiment is 0.00001. By analogy, a similar analysis and calculation on the Y and Z axes can build a simple three-body motion model in three-dimensional free space.
Ii. Source code
%
%%
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 anew 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'[]);ifnargin & & 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)% initializes the mass coordinates of the object and the velocity handles.star=Model_Init_1();
handles.type=1; handles.pretype=handles.type; % global Is_Running; % global Is_Init; % Default state is initialized and running 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; % keep running the run() function, waiting for an interrupt run(); %% % star motion functionfunction run(a)
%%
global Is_Running; global Is_Init; % Universal gravitational constant G=6.67*10^ (- 11); % to0.00001Time step iteration T=0.00001; % Note the dead loopwhile (1) handles=guidata(gcf); % exitif(isempty(handles))
close all;
return; End % is initialized to indicate the end of an operation within a phase and display the resultif(Is_Init==trueStar =handles.Star; % -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- 1 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; % -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- 2 stars -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- - 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; % -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- no. 3 stars -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- the 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; % to empty cla. axis([- 1000. 1010 - 1000. 1000 - 1000. 1000]); Do not use the default coordinate style, draw your own line axis off; % Right hand line in default viewset(gca,'XDir'.'reverse');
set(gca,'YDir'.'reverse'); % draw grid 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 % draw arrow 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); % mark the coordinate axis 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); % mark unit length text(125.0.0.'125'.'Color'[1 0 1].'FontSize'.12);
text(0.125.0.'125'.'Color'[1 0 1].'FontSize'.12); % Initial setting of celestial color, dot type, size and other parameters, each cycle reset the entire screenif(M1~=0)
h=line('Color'[1 0 0].'Marker'.'. '.'MarkerSize'.25.'erasemode'.'normal');
line(X1,Y1,Z1,'Color'[1 0 0].'Marker'.'. '.'MarkerSize'.10);
end
if(M2~=0)
i=line('Color'[0 1 0].'Marker'.'. '.'MarkerSize'.25.'erasemode'.'normal');
line(X2,Y2,Z2,'Color'[0 1 0].'Marker'.'. '.'MarkerSize'.10);
end
if(M3~=0)
l=line('Color'[0 0 1].'Marker'.'. '.'MarkerSize'.25.'erasemode'.'normal');
line(X3,Y3,Z3,'Color'[0 0 1].'Marker'.'. '.'MarkerSize'.10); End % refresh screen drawnow; Is_Init=false;
else
[y,fs] = audioread('Windows Background.wav'); % of each calculation200Redraw it twicefor k=1:200; % Calculate distance 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); First determine whether the planet will collide and explode, relative distance to8A limitif((R12< =8)& & (M1~=0)& & (M2~=0))
M1=0;
M2=0;
delete(h);
delete(i); % collision sound(y,fs);end
if ((R23< =8)& & (M3~=0)& & (M3~=0))
M2=0;
M3=0;
delete(i);
delete(l);
sound(y,fs);
end
if ((R13< =8)& & (M1~=0)& & (M3~=0))
M1=0;
M3=0;
delete(h);
delete(l); sound(y,fs); End % Law of gravity f12=G*M1*M2/R12^3*[X1-X2,Y1-Y2,Z1-Z2];
f23=G*M2*M3/R23^3*[X2-X3,Y2-Y3,Z2-Z3];
f31=G*M3*M1/R13^3*[X3-X1,Y3-Y1,Z3-Z1];
if(Is_Running==true) t=T; % Stop motionelse
t=0; End % stars1
if (M1~=0)
Ax1=(f31(1)-f12(1))/M1;
Ay1=(f31(2)-f12(2))/M1;
Az1=(f31(3)-f12(3))/M1; % Note that the coordinates need to be calculated first, and then the speed X1=X1+U1*t+1/2*Ax1*t^2;
Y1=Y1+V1*t+1/2*Ay1*t^2;
Z1=Z1+W1*t+1/2*Az1*t^2; U1=U1+Ax1*t; V1=V1+Ay1*t; W1=W1+Az1*t; End % stars2
if (M2~=0)
Ax2=(f12(1)-f23(1))/M2;
Ay2=(f12(2)-f23(2))/M2;
Az2=(f12(3)-f23(3))/M2;
X2=X2+U2*t+1/2*Ax2*t^2;
Y2=Y2+V2*t+1/2*Ay2*t^2;
Z2=Z2+W2*t+1/2*Az2*t^2; U2=U2+Ax2*t; V2=V2+Ay2*t; W2=W2+Az2*t; End % stars3
if (M3~=0)
Ax3=(f23(1)-f31(1))/M3;
Ay3=(f23(2)-f31(2))/M3;
Az3=(f23(3)-f31(3))/M3;
X3=X3+U3*t+1/2*Ax3*t^2;
Y3=Y3+V3*t+1/2*Ay3*t^2;
Z3=Z3+W3*t+1/2*Az3*t^2; U3=U3+Ax3*t; V3=V3+Ay3*t; W3=W3+Az3*t; End end % resets the position of the starif(M1~=0)
set(h,'XData',X1,'YData',Y1,'ZData',Z1);
line(X1,Y1,Z1,'Color'[1 0 0].'Marker'.'. '.'MarkerSize'.10);
end
if(M2~=0)
set(i,'XData',X2,'YData',Y2,'ZData',Z2);
line(X2,Y2,Z2,'Color'[0 1 0].'Marker'.'. '.'MarkerSize'.10);
end
if(M3~=0)
set(l,'XData',X3,'YData',Y3,'ZData',Z3);
line(X3,Y3,Z3,'Color'[0 0 1].'Marker'.'. '.'MarkerSize'.10); End % Displays data % displays data for the first starset(handles.Text_m1,'String',num2str(M1*10^ (- 19)));
set(handles.Text_x1,'String',num2str(X1));
set(handles.Text_y1,'String',num2str(Y1));
set(handles.Text_z1,'String',num2str(Z1));
set(handles.Text_v1,'String',num2str(V1*10^ (- 3)));
set(handles.Text_u1,'String',num2str(U1*10^ (- 3)));
set(handles.Text_w1,'String',num2str(W1*10^ (- 3))); % displays data for the second starset(handles.Text_m2,'String',num2str(M2*10^ (- 19)));
set(handles.Text_x2,'String',num2str(X2));
set(handles.Text_y2,'String',num2str(Y2));
set(handles.Text_z2,'String',num2str(Z2));
set(handles.Text_u2,'String',num2str(U2*10^ (- 3)));
set(handles.Text_v2,'String',num2str(V2*10^ (- 3)));
set(handles.Text_w2,'String',num2str(W2*10^ (- 3))); % displays data for the third starset(handles.Text_m3,'String',num2str(M3*10^ (- 19)));
set(handles.Text_x3,'String',num2str(X3));
set(handles.Text_y3,'String',num2str(Y3));
set(handles.Text_z3,'String',num2str(Z3));
set(handles.Text_u3,'String',num2str(U3*10^ (- 3)));
set(handles.Text_v3,'String',num2str(V3*10^ (- 3)));
set(handles.Text_w3,'String',num2str(W3*10^ (- 3))); % shows the distance between starsset(handles.value_r12,'String',num2str(R12));
set(handles.value_r23,'String',num2str(R23));
set(handles.value_r13,'String',num2str(R13)); % Set font color % The first star font corresponds to redset(handles.Text_m1,'ForegroundColor'[1 0 0]);
set(handles.Text_x1,'ForegroundColor'[1 0 0]);
set(handles.Text_y1,'ForegroundColor'[1 0 0]);
set(handles.Text_z1,'ForegroundColor'[1 0 0]);
set(handles.Text_u1,'ForegroundColor'[1 0 0]);
set(handles.Text_v1,'ForegroundColor'[1 0 0]);
set(handles.Text_w1,'ForegroundColor'[1 0 0]); % The second star font corresponds to greenset(handles.Text_m2,'ForegroundColor'[0 1 0]);
set(handles.Text_x2,'ForegroundColor'[0 1 0]);
set(handles.Text_y2,'ForegroundColor'[0 1 0]);
set(handles.Text_z2,'ForegroundColor'[0 1 0]);
set(handles.Text_u2,'ForegroundColor'[0 1 0]);
set(handles.Text_v2,'ForegroundColor'[0 1 0]);
set(handles.Text_w2,'ForegroundColor'[0 1 0]);
Copy the code3. Operation results
Fourth, note
Version: 2014 complete code or write to add QQ 1564658423