【运动学】基于matlab GUI三体运动模拟【含Matlab源码 871期】
【摘要】
一、获取代码方式
获取代码方式1： 完整代码已上传我的资源：【运动学】基于matlab GUI三体运动模拟【含Matlab源码 871期】
获取代码方式2： 通过订阅紫极神光博客付费专栏，凭支付凭证，...
一、获取代码方式
获取代码方式1：
完整代码已上传我的资源：【运动学】基于matlab GUI三体运动模拟【含Matlab源码 871期】
获取代码方式2：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。
备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；
二、案例简介
本文基于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);
%初始设置天体颜色、点型、大小等参数，每次循环都重新设置整个画面
if(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
%刷新画面
drawnow;
Is_Init=false;
else
[y,fs] = audioread('Windows Background.wav');
%每计算200次重绘一次
for k=1:200;
%计算距离
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);
%先判断星球会不会碰撞爆炸，相对距离以8位限
if ((R12<=8)&&(M1~=0)&&(M2~=0))
M1=0;
M2=0;
delete(h);
delete(i);
%碰撞音效
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
%万有引力定律
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;
%停止运动
else
t=0;
end
%星体1
if (M1~=0)
Ax1=(f31(1)-f12(1))/M1;
Ay1=(f31(2)-f12(2))/M1;
Az1=(f31(3)-f12(3))/M1;
%注意，需先计算坐标，后计算速度
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
%星体2
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
%星体3
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
%重置星体的位置
if(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
%显示数据
%显示第一颗星体的数据
set(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)));
%显示第二颗星体的数据
set(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)));
%显示第三颗星体的数据
set(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)));
%显示星体之间的距离
set(handles.value_r12,'String',num2str(R12));
set(handles.value_r23,'String',num2str(R23));
set(handles.value_r13,'String',num2str(R13));
%设置字体颜色
%第一颗星体字体对应为红色
set(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]);
%第二颗星体字体对应为绿色
set(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]);
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 79
- 80
- 81
- 82
- 83
- 84
- 85
- 86
- 87
- 88
- 89
- 90
- 91
- 92
- 93
- 94
- 95
- 96
- 97
- 98
- 99
- 100
- 101
- 102
- 103
- 104
- 105
- 106
- 107
- 108
- 109
- 110
- 111
- 112
- 113
- 114
- 115
- 116
- 117
- 118
- 119
- 120
- 121
- 122
- 123
- 124
- 125
- 126
- 127
- 128
- 129
- 130
- 131
- 132
- 133
- 134
- 135
- 136
- 137
- 138
- 139
- 140
- 141
- 142
- 143
- 144
- 145
- 146
- 147
- 148
- 149
- 150
- 151
- 152
- 153
- 154
- 155
- 156
- 157
- 158
- 159
- 160
- 161
- 162
- 163
- 164
- 165
- 166
- 167
- 168
- 169
- 170
- 171
- 172
- 173
- 174
- 175
- 176
- 177
- 178
- 179
- 180
- 181
- 182
- 183
- 184
- 185
- 186
- 187
- 188
- 189
- 190
- 191
- 192
- 193
- 194
- 195
- 196
- 197
- 198
- 199
- 200
- 201
- 202
- 203
- 204
- 205
- 206
- 207
- 208
- 209
- 210
- 211
- 212
- 213
- 214
- 215
- 216
- 217
- 218
- 219
- 220
- 221
- 222
- 223
- 224
- 225
- 226
- 227
- 228
- 229
- 230
- 231
- 232
- 233
- 234
- 235
- 236
- 237
- 238
- 239
- 240
- 241
- 242
- 243
- 244
- 245
- 246
- 247
- 248
- 249
- 250
- 251
- 252
- 253
- 254
- 255
- 256
- 257
- 258
- 259
- 260
- 261
- 262
- 263
- 264
- 265
- 266
- 267
- 268
- 269
- 270
- 271
- 272
- 273
- 274
- 275
- 276
- 277
- 278
- 279
- 280
- 281
- 282
- 283
- 284
- 285
- 286
- 287
- 288
- 289
- 290
- 291
- 292
- 293
- 294
- 295
- 296
- 297
- 298
- 299
- 300
- 301
- 302
- 303
- 304
- 305
- 306
- 307
- 308
- 309
- 310
- 311
- 312
- 313
- 314
- 315
- 316
- 317
- 318
- 319
- 320
四、运行结果
五、matlab版本及参考文献
1 matlab版本
2014a
2 参考文献
[1] 门云阁.MATLAB物理计算与可视化[M].清华大学出版社，2013.
文章来源: qq912100926.blog.csdn.net，作者：海神之光，版权归原作者所有，如需转载，请联系作者。
原文链接：qq912100926.blog.csdn.net/article/details/116355415
【版权声明】本文为华为云社区用户转载文章，如果您发现本社区中有涉嫌抄袭的内容，欢迎发送邮件进行举报，并提供相关证据，一经查实，本社区将立刻删除涉嫌侵权内容，举报邮箱：
cloudbbs@huaweicloud.com
- 点赞
- 收藏
- 关注作者
评论（0）