试确定系统的稳定性。 %程序如下：
a=[-4,-3,0 ; 1,0,0 ; 0,1,0] ; b=[1;0;0] ; c=[0,1,2] ; d=0 ;
eig(a) %求特征根 rank(ctrb(a,b))
四、实验要求
1．参照已知程序，改动程序中的参数和输入量，验证输出结果。 2．掌握相关命令的用法
3．实验结果写成实验报告五（全部实验完成后交实验报告）。
实验六 连续系统数字仿真的基本算法
一、实验任务
1．理解欧拉法和龙格-库塔法的基本思想；
2．理解数值积分算法的计算精度、速度、稳定性与步长的关系；
14
二、程序举例
1. 取h=0.2，试分别用欧拉法、RK2法和RK4法求解微分方程的数值解，并比较计算精度。 ?2t?(t)?y(t)??yy(t)?
?y(0)?1?
%程序如下 clear
t(1)=0; % 置自变量初值 y(1)=1; y_euler(1)=1; y_rk2(1)=1; y_rk4(1)=1; % 置解析解和数值解的初值 h=0.2; % 步长 % 求解析解 for k=1:5
t(k+1)=t(k)+h;
y(k+1)=sqrt(1+2*t(k+1)); end
% 利用欧拉法求解 for k=1:5
y_euler(k+1)=y_euler(k)+h*(y_euler(k)-2*t(k)/y_euler(k)); end
% 利用RK2法求解 for k=1:5
k1=y_rk2(k)-2*t(k)/y_rk2(k);
k2=(y_rk2(k)+h*k1)-2*(t(k)+h)/(y_rk2(k)+h*k1); y_rk2(k+1)=y_rk2(k)+h*(k1+k2)/2; end
% 利用RK4法求解 for k=1:5
k1=y_rk4(k)-2*t(k)/y_rk4(k);
k2=(y_rk4(k)+h*k1/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k1/2); k3=(y_rk4(k)+h*k2/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k2/2); k4=(y_rk4(k)+h*k3)-2*(t(k)+h)/(y_rk4(k)+h*k3); y_rk4(k+1)=y_rk4(k)+h*(k1+2*k2+2*k3+k4)/6; end % 输出结果
0?t?1 注：解析解：y(t)?1?2t
15
disp(' 时间 解析解 欧拉法 RK2法 RK4法') yt=[t', y', y_euler', y_rk2', y_rk4']; disp(yt)
程序运行结果如下：
时间 解析解 欧拉法 RK2法 RK4法 0 1.0000 1.0000 1.0000 1.0000 0.2000 1.1832 1.2000 1.1867 1.1832 0.4000 1.3416 1.3733 1.3483 1.3417 0.6000 1.4832 1.5315 1.4937 1.4833 0.8000 1.6125 1.6811 1.6279 1.6125 1.0000 1.7321 1.8269 1.7542 1.7321
2. 考虑如下二阶系统：
?(0)?0）?(t)?2Ry(t)?y(t)?0在0?t?10上的数字仿真解（已知：y(0)?100，yy ?，
并将不同步长下的仿真结果与解析解进行精度比较。 说明：
已知该微分方程的解析解分别为：
y ? t ??100e1?t2y?t??100cost1(当R?0)(当R?0.5) 采用RK4法进行计算，选择状态变量： ?x?y1? ??x2?y 则有如下状态空间模型及初值条件
?x??x?12 ??2??x1?2Rx2?x ? y?x1x1(0)?100x2(0)?031003?2t3cost?esint232
采用RK4法进行计算。程序如下： clear
h=input('请输入步长h='); % 输入步长 M=round(10/h); % 置总计算步数 t(1)=0; % 置自变量初值
y_0(1)=100; y_05(1)=100; % 置解析解的初始值(y_0和y_05分别对应
于为R=0和R=0.5)
x1(1)=100; x2(1)=0; % 置状态向量初值
16
y_rk4_0(1)=x1(1); y_rk4_05(1)=x1(1); % 置数值解的初值 % 求解析解 for k=1:M t(k+1)=t(k)+h;
y_0(k+1)=100*cos(t(k+1));
y_05(k+1)=100*exp(-t(k+1)/2).*cos(sqrt(3)/2*t(k+1))+100*sqrt(3)/3*exp(-t(k+1)/2).*sin(sqrt(3)/2*t(k+1)); end
% 利用RK4法求解 % R=0 for k=1:M
k11=x2(k); k12=-x1(k);
k21=x2(k)+h*k12/2; k22=-(x1(k)+h*k11/2); k31=x2(k)+h*k22/2; k32=-(x1(k)+h*k21/2); k41=x2(k)+h*k32; k42=-(x1(k)+h*k31); x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6; x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6; y_rk4_0(k+1)=x1(k+1); end % R=0.5 for k=1:M
k11=x2(k); k12=-x1(k)-x2(k);
k21=x2(k)+h*k12/2; k22=-(x1(k)+h*k11/2)-(x2(k)+h*k12/2); k31=x2(k)+h*k22/2; k32=-(x1(k)+h*k21/2)-(x2(k)+h*k22/2); k41=x2(k)+h*k32; k42=-(x1(k)+h*k31)-(x2(k)+h*k32); x1(k+1)=x1(k)+h*(k11+2*k21+2*k31+k41)/6; x2(k+1)=x2(k)+h*(k12+2*k22+2*k32+k42)/6; y_rk4_05(k+1)=x1(k+1); end
% 求出误差最大值
err_0=max(abs(y_0-y_rk4_0)); err_05=max(abs(y_05-y_rk4_05)); % 输出结果
disp('最大误差（R=0) 最大误差（R=0.5)') err_max=[err_0,err_05]; disp(err_max)
17
步长h取0.0001 到0.5之间变化，并且将数值解直接与解析解比较，求出最大误差数据列入下表中。
步长h 0.0001 0.0005 0.001 0.005 0.01 0.05 0.1 0.5 R=0 5.4330×最大误10-10 差 10-10 10-10 10-9 10-8 10-5 10-4 10-1 1.6969×1.0574×4.1107×6.6029×4.1439×6.6602×4.2988×R=0.5 2.7649×最大误10-11 差 10-12 10-12 10-10 10-9 10-6 10-5 10-2 6.8123×5.3753×4.0902×6.5425×4.1365×6.7152×4.5976×
从上表中可以看出，当步长h=0.001时，总误差最小；当步长h小于0.001时，由于舍入误差变大而使总误差增加；当步长h大于0.001时，则由于截断误差的增加也使得总误差加大。另外，当系统的解变化激烈时（如R=0），误差对步长的变化较为敏感；当系统的解变化平稳时，步长的变化对误差的影响就要缓和得多。数值积分算法确定以后，在选择步长时，需要综合考虑。
三、实验要求
1．参照已知程序，改动程序中的参数和输入量，验证输出结果。 2．分析不同参数取值对输出结果的影响。
3．实验结果写成实验报告六（全部实验完成后交实验报告）。
MATLAB与控制系统仿真
实 验 报 告
18