MATLAB期末大作业.docx
- 文档编号:18240442
- 上传时间:2023-08-14
- 格式:DOCX
- 页数:14
- 大小:687.21KB
MATLAB期末大作业.docx
《MATLAB期末大作业.docx》由会员分享,可在线阅读,更多相关《MATLAB期末大作业.docx(14页珍藏版)》请在冰点文库上搜索。
MATLAB期末大作业
学号:
姓名:
《Matlab/Simulink在数学计算与仿真中的应用》大作业
1.假设地球和火星绕太阳运转的半径分别为
和
,利用comet指令动画显示从地球到火星的转移轨迹(
可以任意取值,要求实时显示探测器、太阳、地球和火星的位置)。
解函数functioncomet(varargin)
[ax,args,nargs]=axescheck(varargin{:
});
error(nargchk(1,3,nargs,'struct'));
%Parsetherestoftheinputs
ifnargs<2,x=args{1};y=x;x=1:
length(y);end
ifnargs==2,[x,y]=deal(args{:
});end
ifnargs<3,p=0.10;end
ifnargs==3,[x,y,p]=deal(args{:
});end
if~isscalar(p)||~isreal(p)||p<0||p>=1
error('MATLAB:
comet:
InvalidP',...
'Theinput''p''mustbearealscalarbetween0and1.');
End
指令%particle_motion
t=0:
5:
16013;
r1=6.7e6;%随便给定参数
%---------------------------
r2=2*r1;
g=9.8;
R=6.378e6;
m=g*R^2;
%内轨道
v_inner=sqrt(m/r1);
w_inner=v_inner/r1;
x_inter=r1*cos(w_inner*t);
y_inter=r1*sin(w_inner*t);
%外轨道
v_outer=sqrt(m/r2);
w_outer=v_outer/r2;
x_outer=r2*cos(w_outer*t);
y_outer=r2*sin(w_outer*t);
%控制器转移轨道
a=(r1+r2)/2;
E=-m/(2*a);
V_near=sqrt(m*(2/r1-2/(r1+r2)));%转移轨道在近地点的速度V_far=sqrt(m*(2/r2-2/(r1+r2)));%转移轨道在远地点的速度
h=r1*V_near;%由于在近地点的速度垂直于位置失量,h是转移轨道的比动量矩
e=sqrt(1+2*E*h^2/m^2);%e为椭圆轨迹的偏心率
TOF=pi*sqrt(a^3/m);%转移轨道是椭圆的一半及飞行时间是周期的一半(开普勒第三定律)
w=pi/TOF;%椭圆轨迹的角速度
c=a*e;
b=sqrt(a^2-c^2);
x_ellipse=a*cos(w*t)-0.5*r1;
y_ellipse=b*sin(w*t);
%动画显示运动轨迹
x=[x_inter;x_outer;x_ellipse]';
y=[y_inter;y_outer;y_ellipse]';
comet(x,y)
%---------------------------
动态图像如下:
2.利用两种不同途径求解边值问题
的解.途径一:
指令symsfg
[f,g]=dsolve('Df=3*f+4*g,Dg=-4*f+3*g','f(0)=0,g(0)=1');
disp('f=');disp(f)
disp('g=');disp(g)
结果(Matlab7.8版本)f=
i/(2*exp(t*(4*i-3)))-(i*exp(t*(4*i+3)))/2
g=exp(t*(4*i+3))/2+1/(2*exp(t*(4*i-3)))
(Matlab6.5版本)
f=exp(3*t)*sin(4*t)
g=exp(3*t)*cos(4*t)
>>
途径二:
%problem2
functiondy=problem2(t,y)
dy=zeros(2,1);
dy
(1)=3*y
(1)+4*y
(2);
dy
(2)=-4*y
(1)+3*y
(2);
[t,y]=ode45('problem2',[02],[01]);
plot(t,y(:
1),'r',t,y(:
2),'b');
图2
3.假设著名的Lorenz模型的状态方程表示为
其中,设
。
若令其初值为
,而
为机器上可以识别的小常数,如取一个很小的正数
,试用尽可能多的方法来求解该微分方程组。
解:
方法一:
新建文件lorenzeq.m
源程序:
functionxdot=lorenzeq(t,x)
xdot=[-8/3*x
(1)+x
(2)*x(3);...
-10*x
(2)+10*x(3);...
-x
(1)*x
(2)+28*x
(2)-x(3)];
再输入命令:
t_final=100;
x0=[0;0;1e-10];
[t,x]=ode45('lorenzeq',[0,t_final],x0);
plot(t,x);
figure;
plot3(x(:
1),x(:
2),x(:
3));
axis([1042-2020-2025]);
得到图形解
:
三维:
图3-1
二维:
图3-2
方法二通过simulink来实现
图1-3
输入命令:
beta=8/3;rho=10;sigma=28;
[t,x]=sim('sumlik_3',[0,100]);plot(t,x)
figure;plot3(x(:
1),x(:
2),x(:
3))
得到和上面一样的图形。
4.已知Apollo卫星的运动轨迹
满足下面的方程
其中,
,试在初值
下分别用ode45命令和Simulink进行求解,并绘制出Apollo位置的
轨迹。
解:
ode45命令:
新建文件aplloeq.m
functiondx=apolloeq(t,x)
mu=1/82.45;mu1=1-mu;
r1=sqrt((x
(1)+mu)^2+x(3)^2);
r2=sqrt((x
(1)-mu1)^2+x(3)^2);
dx=[x
(2);
2*x(4)+x
(1)-mu1*(x
(1)+mu)/r1^3-mu*(x
(1)-mu1)/r2^3;
x(4);
-2*x
(2)+x(3)-mu1*x(3)/r1^3-mu*x(3)/r2^3];
输入命令:
x0=[1.2;0;0;-1.04935751];
options=odeset;
options.RelTol=1e-6;
tic,[t1,y1]=ode45('apolloeq',[0,20],x0,options);
toc,length(t1),plot(y1(:
1),y1(:
3))
得到结果:
Elapsedtimeis0.139725seconds.
ans=
1873
图4
Simulink也能实现上述结果:
5.利用Matlab优化工具箱或该工具箱中的函数求下面的最大值。
解:
源程序:
f=[1-21-3]';%目标函数的负数表示,求此目标函数的最小最优解
A=[0-211;0-16-1];%不等式约束条件,变量前面的系数矩阵
B=[3;4];%不等式约束条件,条件小于的数值矩阵
Aeq=[1131];%等式约束条件,变量前的系数矩阵
Beq=[6];%等式约束条件,等式等于的常数矩阵
xm=[0000];%x向量解得下界
xM=[infinfinfinf];%x向量解的上界
[xf_optkeyc]=linprog(f,A,B,Aeq,Beq,xm,xM);%用linpro函数求解最优值
x’%输出最优解时的x向量解
max_opt=-f_opt%将求得的最小最优解取负数,得到要求目标函数的最大最优解
计算结果:
>>lin3%保存的文件名
Optimizationterminated.
ans=
0.00001.00000.00005.0000
max_opt=
17.0000
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- MATLAB 期末 作业