计算机辅助设计MATLAB试题2013-A卷答案.doc

试 题： 学号： 姓名：注：每题10分，总分70分。一、基本操作类：编写代码实现以下功能，1）如何快速自学命令dsolve？如何查询它属于哪个命令集？ 2）x为1到9的5维等距数组，数组y中的第个元素满足，不使用for循环求数组y；3）用两种方法解方程Ax=b，A为4阶随机矩阵，b为4阶列向量。答：建立M文件编写以下代码：1） help dsolve2） x=linspace(1,9,5);y= x.3 +exp(x)+tan(sqrt(1./ x-1)+sin(75*pi/180)3） x=A-1*b;或x= Ab;二、基本操作、矩阵操作类：卫星姿态动力学方程（欧拉方程）为。某一时刻,控制力矩,已知卫星惯量矩阵为,试求出该时刻的。创建相应矩阵、数列，编写代码（用dw表示）。I=50,-15,20;-15,62,10;20,10,50;T=0 0.1 0;w=4 8 3*1e-2/57.3;dw=I-1*(T-cross(w,I*w);三、数据可视化：已知二阶系统的单位阶跃响应（欠阻尼状态）为，其中，系统固有频率。在不同坐标系内绘制系统阻尼比分别取0.4、0.7时响应图像。要求编写输出为函数，并在另一M文件中调用该函数实现该功能。t= 0:0.1:15;wn=sqrt(2);zita1=0.3;c1=shuchu(wn,zita1,t);zita2=0.7;c2=shuchu(wn,zita2,t);subplot(211);plot(t,c1);title(输出(zeta=0.4);grid on;subplot(212);plot(t,c2);title(输出(zeta=0.7);grid on;xlabel(t/s);function y=shuchu(wn,zita,t)wd=wn*sqrt(1-zita2);y=1-exp(-zita*wn*t).*(cos(wd*t)-zita*wn/wd*sin(wd*t);四、数据可视化类：在同一坐标内绘制图形，实现如下要求。1）用红色的实线绘制曲线；2）用蓝色的点划线绘制曲线,；3）设置x轴注释为time/s，y轴注释为y(t)，并添加图例标注，添加标题名为数据可视化。4）用黑色圆圈标出最大值点。 t=linspace(0,10,200);y1=1-1.02*exp(-0.2*pi*t).*sin(pi*t+1.37);y2=1+1.02*exp(-0.2*pi*t);y3=1-1.02*exp(-0.2*pi*t);plot(t,y1,r,t,y2,b-.,t,y3,b-.);grid on;axis tight;xlabel(time/s)ylabel(y(t)legend(y_1=1-1.02*exp(-0.2*pi*t)*sin(pi*t+1.37),y_2(y_3)=11.02*exp(-0.2*pi*t);text(2.5,y2_max,leftarrow最大值点)title(数据可视化);y1_max,i_max=max(y1)； %查找y1最大值及其下标plot(t(i_max),y1_max,ko);text(t(i_max)+0.3,y1_max,leftarrow最大值点)； %在图上标注最大值点五、动态系统仿真类，带挠性太阳帆的卫星，俯仰回路转动方程和帆板振动方程（只考虑第三阶扭转模态）为，其中为卫星俯仰轴转动惯量，为耦合系数，为控制力矩，为扭转模态坐标，为帆板阻尼系数，为扭转模态频率。初值为，常数如下。其中，。利用matlab的ode45积分器进行求解, 积分时间为100s, 绘制随时间变化曲线。答：主函数main5：clc;clear;tspan = 0 100;y0 = 3*pi/180;0;0;0;t,yout = ode45(spacecraft_ode,tspan,y0);subplot(211);plot(t,yout(:,1)*180/pi);grid on;title(简单挠性体动力学);ylabel(theta(o),rotation,0,position,-5,1);subplot(212);plot(t,yout(:,3);grid on;ylabel(xi,rotation,0,position,-5,-2e-5);xlabel(t(s);M函数，命名为spacecraft_odefunction dydt = spacecraft_ode(t,y)theta = y(1);dtheta = y(2);xi = y(3);dxi = y(4);Kp = 22.5;Kd = 48.75; Iy = 424.85;B=0.5572; k=B2/Iy;zeta =0.0397; Omega =2.4804*pi;Ty = -Kp*theta - Kd*dtheta; ddxi=1/(1-2*k)*(Ty/Iy*B-2*zeta*Omega*dxi-Omega2*xi);ddtheta = 1/Iy*(Ty+2*B*ddxi);dydt = dtheta;ddtheta;dxi;ddxi;六、动态仿真类：已知卫星轨道方程为，其中，(单位为km)，（单位为km/s）；要求采用eular方法设计积分器,编写其ode文件及主函数代码，仿真卫星飞行6小时的飞行轨道。答：主函数main5：clc;clear;y0=6700;0;0;9.7;t,y=euler_ode(funy,0,3600*6,y0,36000);plot(y(:,1),y(:,2),0,0,rx);xlabel(x/km);ylabel(y/km);ode函数，命名为funyfunction dy=funy(y)u=3.986e5;r=y(1);y(2);v=y(3);y(4);dy=v; -u/norm(r)3*r;欧拉积分器，命名为euler_odefunction t,y=euler_ode(funy,t0,tf,y0,N)h=(tf-t0)/N;t=t0+1:N*h;y(1,:)=y0;for k=1:N f=feval(funy,y(k,:); f=f; y(k+1,:)=y(k,:)+h*f;end 七、动态系统仿真类（ode45复杂调用格式）：对第五题采用ode复杂调用方式进行数值计算，系统预定仿真时间为100s，但当时程序自动停止，编写其ode文件及主函数代码。M文件主函数clc;clear;T,Y,TE,YE,IE=ode45(spacecraft_ode ,);编写M函数，命名为odefilefunction varargout=odefile(t,y,flag,tspan,y0)switch flag case , varargout1=f(t,y);case init, varargout1:3 = init(tspan,y0);case events,varargout1:3 = events(t,y);otherwise , error(Unknown flag flag .);endfunction dydt = f(t,y)theta = y(1);dtheta = y(2);xi = y(3);dxi = y(4);Kp = 22.5;Kd = 48.75;Iy = 424.85;B=0.5572; k=B2/Iy;zeta =0.0397; Omega =2.4804*pi;Ty = -Kp*theta - Kd*dtheta; ddxi=1/(1-2*k)*(Ty/Iy*B-2*zeta*Omega*dxi-Omega2*xi);ddtheta = 1/Iy*(Ty+2*B*ddxi);dydt = dtheta;ddtheta;dxi;ddxi;function tspan,y0,options=i