MATLAB期末大作业.doc

学号：060910305 姓名：杜丹丹Matlab/Simulink在数学计算与仿真中的应用大作业1. 假设地球和火星绕太阳运转的半径分别为和，利用comet指令动画显示从地球到火星的转移轨迹（可以任意取值，要求实时显示探测器、太阳、地球和火星的位置）。解 函数 function comet(varargin)ax,args,nargs = axescheck(varargin:);error(nargchk(1,3,nargs,struct);% Parse the rest of the inputsif nargs 2, x = args1; y = x; x = 1:length(y); endif nargs = 2, x,y = deal(args:); endif nargs 3, p = 0.10; endif nargs = 3, x,y,p = deal(args:); end if isscalar(p) | isreal(p) | p = 1 error(MATLAB:comet:InvalidP, . The input p must be a real scalar between 0 and 1.);End 指令 %particle_motiont = 0:5:16013;r1=6.7e6;%随便给定参数%-r2=2*r1;g=9.8;R=6.378e6;m=g*R2;%内轨道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*h2/m2);%e为椭圆轨迹的偏心率TOF=pi*sqrt(a3/m);%转移轨道是椭圆的一半及飞行时间是周期的一半(开普勒第三定律)w=pi/TOF;%椭圆轨迹的角速度c=a*e;b=sqrt(a2-c2);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利用两种不同途径求解边值问题的解. 途径一：指令syms f gf,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)结果（Matlab 7.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) （Matlab 6.5版本）f=exp(3*t)*sin(4*t) g=exp(3*t)*cos(4*t) 途径二：%problem2function dy=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,0 2,0 1); plot(t,y(:,1),r,t,y(:,2),b);图23假设著名的Lorenz模型的状态方程表示为其中，设。若令其初值为，而为机器上可以识别的小常数，如取一个很小的正数，试用尽可能多的方法来求解该微分方程组。解：方法一：新建文件lorenzeq.m源程序： function xdot=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(10 42 -20 20 -20 25);得到图形解：三维：图3- 1二维：图3-2方法二 通过simulink来实现图1- 2输入命令：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.mfunction dx=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)/r13-mu*(x(1)-mu1)/r23; x(4); -2*x(2)+x(3)-mu1*x(3)/r13-mu*x(3)/r23;输入命令：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)得到结果：Elapsed time is 0.139725 seconds.ans = 1873图 3Simulink也能实现上述结果：5. 利用Matlab优化工具箱或该工具箱中的函数求下面的最大值。解：源程序：f=1 -2 1 -3; %目标函数的负数表示，求此目标函数的最小最优解A=0 -2 1 1;0 -1 6 -1; %不等式约束条件，变量前面的系数矩阵B=3;4; %不等式约束条件，条件小于的数值矩阵Aeq=1 1 3 1; %等式约束条件，变量前的系数矩阵Beq=6; %等式约束条件，等式等于的常数矩阵xm=0 0 0 0; %x向量解得下界xM=inf inf inf inf; %x向量解的上界x f_opt key c=linprog(f,A,B,Aeq,Beq,xm,xM); %用linpro函数求