1、）; %显示 中间的文字%y= %同上%y=1;for t=0:h:1 m=y; disp（y）; %显示y的当前值% y=m-m*h;end保存文件q2.m 在matalb命令行中键入 q2 得到结果 函数的数值解为y= 1 0.9000 0.8100 0.7290 0.6561 0.5905 0.5314 0.4783 0.4305 0.3874 0.3487（2）另建一个m 文件求解在t 0,1的数值 （ %是的真解%）程序为h=0.1;函数的离散时刻解为 y=exp（-t）;end 保存文件q3.m在matlab命令行中键入 q3 函数的离散时刻解为y= 1 0.9048 0.8187
2、 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679比较欧拉方法求解与真值的差别欧拉0.90000.81000.72900.65610.59050.53140.47830.43050.38740.3487真值0.90480.81870.74080.67030.60650.54880.49660.44930.40660.3679误差-0.0048-0.0070.01180.01420.01600.01740.01830.0188-0.0192显然误差与为同阶无穷小，欧拉法具有一阶计算精度，精度较低，但算法简单。我们经常用到 预报-校正法
3、 的二阶龙-格库塔法，此方法可以自启动，具有二阶计算精度，几何意义：区间内的曲边面积用上下底为和、高为h的梯形面积近似代替。 （1）m文件程序为 h=0.1; k1=-y; k2=-（y+k1*h）; y=y+（k1+k2）*h/2; 保存文件q4.m在matlab的命令行中键入 q4 显示结果为 函数的数值解为y= 1 0.9050 0.8190 0.7412 0.6708 0.6071 0.5494 0.4972 0.4500 0.4072 0.3685（1） 比较欧拉法与二阶龙格-库塔法求解.（误差为绝对值）龙库0.90500.81900.74120.67080.60710.54940.
4、49720.45000.40720.36850.00020.00030.00040.00050.00060.0007明显误差为的同阶无穷小，具有二阶计算精度，而欧拉法具有一阶计算精度，二阶龙格-库塔法比欧拉法计算精度高。三、 （20分）分别使用解微分方程方法、控制工具箱、simulink求解具有如下闭环传递函数的系统的阶跃响应。（1）用解微分方程方法：将转化为状态方程，利用matlab语句 num=10; den=1 8 36 40 10; A B C D=tf2ss（num,den）得到结果：A = -8 -36 -40 -10 1 0 0 0 0 1 0 0 0 0 1 0B = 1 0C
5、 = 0 0 0 10D =0 得到状态方程编写m文件求解微分方程组function dx=wffc（t,x）u=1; %阶跃响应，输入为1%dx=-8*x（1）-36*x（2）-40*x（3）-10*x（4）+u;x（1）;x（2）;x（3）;保存文件 wffc.m %注意：保存文件的名字与函数名一致！%在命令行键入 t,x=ode45（wffc,0,8,0;0;0）; y=10*x（:,4）; plot（t,y）; grid 得到结果为下图所示： （2）控制工具箱: sys=tf（num,den）; step（sys）; grid得到阶跃响应结果如图所示：（3）simulink求解:在si
6、mulink模型窗口中建立如下模型，键入该题的传递函数。start后，观察scope中的仿真波形如下：四、 （20分）单位反馈系统的开环传递函数已知如下：用matlab语句 、函数求取系统闭环零极点，并求取系统闭环状态方程的可控标准型实现。用单变量系统四阶龙格-库塔法求解输入阶跃函数r（t）=1（t）;T=1s时的输出量y（t）的动态响应数值解。已知开环传递函数，求得闭环传递函数为 在matlab命令行里键入 a=1 0; b=1 4.6; c=1 3.4 16.35; d=conv（a,b）; e=conv（d,c）e = 1.0000 8.0000 31.9900 75.2100 0 f=
7、0 0 0 5 100; g=e+fg = 1.0000 8.0000 31.9900 80.2100 100.0000%以上是计算闭环传递函数的特征多项式% p=roots（g） %计算特征多项式的根，就是闭环传递函数的极点%p =-0.9987 + 3.0091i -0.9987 - 3.0091i -3.0013 + 0.9697i -3.0013 - 0.9697i m=5 100; z=roots（m）z = -20 %计算零点% 综上：当闭环传函形如时，可控标准型为：； 所以可控标准型是m文件为：function y=hs（A,B,C,D,R,T,h） %T为观测时间,h为计算步长
8、，R为输入信号幅值%数值解为y=0;r=R;x=0;0;N=T/h;for t=1:N; k1=A*x+B*R; k2=A*（x+h*k1/2）+B*R;k3=A*（x+h*k2/2）+B*R;k4=A*（x+h*k3）+B*R;x=x+h*（k1+2*k2+2*k3+k4）/6;y（t）=C*x+D*R;在命令行里键入A= B= C= D= R= T= h=y=hs（A,B,C,D,R,T,h） 得到结果。五、 （20分）系统结构图如图所示，（1） 写出该系统的联结矩阵，并写出联结矩阵非零元素阵（2）上图中，若各环节的传递函数已知为：但重新列写联接矩阵和非零元素矩阵，将程序sp4_2.m完善
9、后，应用sp4_2.m求输出的响应曲线。 （1） 根据图中,拓扑连结关系，可写出每个环节输入受哪些环节输出的影响， 现列如入下: 把环节之间的关系和环节与参考输入的关系分别用矩阵表示出来，即=， =（2）程序为：%输入数据%系统中不能出现纯比例、纯微分环节，含有微分项系数的环节不直接与外加参考输入联接P=input（请按各环节输入参数矩阵P：Wij=input（请输入非零元素连接阵Wij:n=input（请输入环节个数（系统阶次）n=Y0=input（请输入阶跃输入幅值Y0=Yt0=input（请输入各环节初值Yt0=h=input（请输入计算步长h=L1=input（请输入打印间隔点数（每隔
10、点l1输出一次）：T0=input（请输入起始时间T0=Tf=input（请输入终止时间Tf=nout=input（请输入输出环节编号：%形成闭环各系数阵%A=diag（P（:,1）;B=diag（P（:,2）;C=diag（P（:,3）;D=diag（P（:,4）;m=length（Wij（:%求非零元素的个数W0=zeros（n,1）;W=zeros（n,n）;%建立初始W（n*n型方阵）、W0阵（n维列向量）for k=1:m if （Wij（k,2）=0）;W0（Wij（k,1）=Wij（k,3）; else W（Wij（k,1）,Wij（k,2）=Wij（k,3）; endQ=B-D
11、*W;Qn=inv（Q）;%求Q和Q的逆R=C*W-A;V1=C*W0;%求R和V1阵Ab=Qn*R;b1=Qn*V1;%形成闭环系数阵%数值积分求解%Y=Yt0;y=Y（nout）;t=T0;%置初值，做好求解准备N=round（Tf-T0）/（h*L1）;for i=1:%每循环一次，输出一点数据 for j=1:L1;%每输出点之间计算L1次 K1=Ab*Y+b1*Y0; K2=Ab*（Y+h*K1/2）+b1*Y0; K3=Ab*（Y+h*K2/2）+b1*Y0; K4=Ab*（Y+h*K3）+b1*Y0; Y=Y+h*（K1+2*K2+2*K3+K4）/6; y=y,Y（nout）;
12、 t=t,t（i）+h*L1;shuchu=t,y;plot（t,y）执行后在命令窗口输入：1 0.01 1 0;0 0.085 1 0.17;1 0.01 1 0;0 0.051 1 0.15;1 0.0067 70 0;1 0.15 0.21 0;0 1 130 0;1 0.01 0.1 0;1 0.01 0.0044 01 0 1;2 1 1;2 9 -1;3 2 1;4 3 1;4 8 -1;5 4 1;6 5 1;6 7 0.212;7 6 1;8 6 1;9 7 1请输入环节个数（系统阶次）n=9请输入阶跃输入幅值Y0=1请输入各环节初值Yt0=0 0 0 0 0 0 0 0 0请输入计算步长h=0.13请输入起始时间T0=0请输入终止时间Tf=107执行结果为：
