得到结果
num = 0 4.0000 14.0000 22.0000 15.0000 den =1.0000 4.0000 6.2500 5.2500 2.2500
4 s3 + 14 s2 + 22 s + 15 G(s)?4 32s + 4 s+ 6.25 s + 5.25 s + 2.25
(2) 零极点增益模型参数:编写程序>> A=[2.25 -5 -1.25 -0.5
2.25 -4.25 -1.25 -0.25 0.25 -0.5 -1.25 -1
1.25 -1.75 -0.25 -0.75]；
>> B=[4 2 2 0]'; >> C=[0 2 0 2];
>> D=[0];
>> [Z,P,K]=ss2zp(A,B,C,D)
得到结果Z =-1.0000 + 1.2247i -1.0000 - 1.2247i -1.5000
P= -0.5000 + 0.8660i -0.5000 - 0.8660i -1.5000
-1.5000
K = 4.0000
表达式 G(s)?
（3）部分分式形式的模型参数：编写程序>> A=[2.25 -5 -1.25 -0.5
4?s+1-1.2247i??s+1+1.2247i?
?s+0.5-0.866i??s+0.5+0.866i??s+1.5?2.25 -4.25 -1.25 -0.25 0.25 -0.5 -1.25 -1
1.25 -1.75 -0.25 -0.75]；
>> B=[4 2 2 0]'; >> C=[0 2 0 2];
>> D=[0];
>> [num den]=ss2tf(A,B,C,D)
>> [R,P,H]=residue(num,den)
得到结果R = 4.0000 -0.0000 0.0000 - 2.3094i 0.0000 +
2.3094i
P = -1.5000 -1.5000 -0.5000 + 0.8660i -0.5000 -
0.8660i
H =[]
G(s)?
42.3094i2.3094i??
s?1.5s?0.5?0.866is?0.5?0.866i
2-3.用欧拉法求下面系统的输出响应y(t)在0≤t≤1上，h=0.1时的数值。 y'??y,y(0)?1
要求保留4位小数，并将结果与真解y(t)?e?t比较。
?yk?1?yk?h*f(tk,yk)?解：欧拉法?y'?f(tk,yk)（前向欧拉法,可以自启动）其几何意义：把f(t,y)
?y(t)?y0?0在[tk,yk]区间内的曲边面积用矩形面积近似代替。利用matlab提供的m文件编程，得到算法公式。如下所示
（1） m文件程序为 h=0.1;
disp('函数的数值解为'); %显示 ??中间的文字% disp('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 文件求解y?e?t在t?[0,1]的数值 （ %y?e?t是
y'??y,y(0)?1的真解%）
程序为h=0.1;
disp('函数的离散时刻解为'); disp('y='); for t=0:h:1 y=exp(-t); disp(y);
end 保存文件q3.m 在matalb命令行中键入>> q3 函数的离散时刻解为
y= 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679
比较欧拉方法求解与真值的差别 欧1 0.9000 拉 真1 0.9048 值 0.8100 0.8187 0.7290 0.7408 0.6561 0.6703 0.5905 0.6065 0.5314 0.5488 0.4783 0.4966 0.4305 0.4493 0.3874 0.4066 0.3487 0.3679 误0 -0.0048 -0.0007 –0.0118 –0.0142 –0.0160 –0.0174 –0.0183 –0.0188 -0.0192 -0.0192 差 显然误差与h2为同阶无穷小，欧拉法具有一阶计算精度，精度较低，但算法简单。
2-4用二阶龙格库塔法求解2-3的数值解，并于欧拉法求得的结果比较。
h?y?y?(k1?k2)k?k?12??解：我们经常用到 预报-校正法 的二阶龙-格库塔法， ?k1?f(tk,yk)此
?k2?f(tk?h,yk?hk1)?'??f(t,y)?y方法可以自启动，具有二阶计算精度，几何意义：把f(t,y)在[tk,yk]区间内的曲边面积用上下底为fk和fk?1、高为h的梯形面积近似代替。利用matlab提供的m文件编程，得到算法公式。如下所示
（1）m文件程序为 h=0.1;
disp('函数的数值解为'); disp('y='); y=1; for t=0:h:1
disp(y); k1=-y;
k2=-(y+k1*h); y=y+(k1+k2)*h/2;
end
保存文件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
（2） 比较欧拉法与二阶龙格-库塔法求解.（误差为绝对值）
真值 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679
龙库 误差 1 0 0.9050 0.8190 0.7412 0.6708 0.6071 0.5494 0.4972 0.4500 0.4072 0.3685 0.0002 0.0003 0.0004 0.0005 0.0006 0.0006 0.0006 0.0007 0.0006 0.0006 明显误差为h3得同阶无穷小，具有二阶计算精度，而欧拉法具有以阶计算精度，二阶龙格-库塔法比欧拉法计算精度高。
2-5．用四阶龙格-库塔法求解题2-3数值解，并与前两题结果相比较。
h?y?y?(k1?2k2?2k3?k4)k?k?16??k1?f(tk,yk)?hh?解：四阶龙格-库塔法表达式?k2?f(tk?,yk?k1)，其截断误差为h522?hh?k?f(t?,y?k2)kk?322???k4?f(tk?h,yk?hk3)同阶无穷小，当h步距取得较小时，误差是很小的.
(1) 编辑m文件程序
h=0.1;
disp('四阶龙格-库塔方法求解函数数值解为'); disp('y='); y=1; for t=0:h:1 disp(y=); k1=-y;
k2=-(y+k1*h/2); k3=-(y+k2*h/2); k4=-(y+k3*h);
y=y+(k1+2*k2+2*k3+k4)*h/6;
end 保存文件q5.m 在matlab命令行里键入>> q5
得到结果 四阶龙格-库塔方法求解函数数值解为
y= 1 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679
(2)比较这几种方法： 对于四阶龙格-库塔方法 真1 值 龙1 库 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679 0.9048 0.8187 0.7408 0.6703 0.6065 0.5488 0.4966 0.4493 0.4066 0.3679
误0 差 0 0 0 0 0 0 0 0 0 0 显然四阶龙格-库塔法求解精度很高，基本接近真值。三种方法比较可以得到 精度（四阶 ） 〉精度（二阶） 〉精度（欧拉）
?.??x1??a11a12??x1??b1??x1(0)??x10?2-6．已知二阶系统状态方程为?.?????????u;????x?; xbx(0)?x2??a21a22??2??2??2??20???写出取计算步长为h时，该系统状态变量X=[x1,x2]的四阶龙格-库塔法递推关系式。
h?y?y?(k1?2k2?2k3?k4)k?k?16??k1?f(tk,yk)?hh?解：四阶龙格-库塔法表达式?k2?f(tk?,yk?k1)
22?hh?k?f(t?,y?k2)kk?322???k4?f(tk?h,yk?hk3)所以状态变量的递推公式可以写作：
.?a11a12??b1??x1? A=??,B=??,x???可以写成X?AX?Bu
?b2??x2??a21a22?h?X?X?*(k1?2k2?2k3?k4)k?k?16??k1?AXk?Bu?则递推形式?k2?A(Xk?k1*h/2)?Bu
?k?A(X?k*h/2)?Buk2?3?k4?A(Xk?k3*h)?Bu??
2-7单位反馈系统的开环传递函数已知如下
G(s)?5s?100 2s(s?4.6)(s?3.4s?16.35)用matlab语句 、函数求取系统闭环零极点，并求取系统闭环状态方程的可控标准型实现。
解：已知开环传递函数，求得闭环传递函数为
G(s)?5s?100 2s(s?4.6)(s?3.4s?16.35)?5s?100