数值积分法的matlab程序

k_1=a*y1(:,i)+b*u; k_2=a*y1(:,i-1)+b*u; k_3=a*y1(:,i-2)+b*u; k_4=a*y1(:,i-3)+b*u; x1=x+h/24*(55*k_1-59*k_2+37*k_3-9*k_4); k_5=a*x1+b*u; x=x+h/24*(9*k_5+19*k_1-5*k_2+k_3); y=[y;x']; y1=y'; end t1=0:0.01:10; y2=1.84-4.95*t1.*exp(-1.88*t1)-1.5*exp(-1.88*t1)-0.34*exp(-6.24*t 1); subplot(2,2,k) plot(t,y) hold on plot(t1,y2,'k-.','linewidth',2) xlabel('time(sec)') ylabel('·ùÖµy') legend('y','doty','ddoty','y-shuzhi') if (k==1) title('admas h=0.01') axis([0,10,-1,3]) elseif (k==2) title('admas h=0.1') axis([0,10,-1,3]) elseif (k==3) title('admas h=0.105 axis([0,10,-1,3]) else (k==4) title('admas h=0.21 临界') axis([0,20,-3,3]) end end 0.5*临界')
h P y y (55 f n 59 f n 1 37 f n 2 9 f n 3 ) n 1 n 24 四阶亚当姆斯预估校正法： y C y h 9 f P 19 f 5 f f n 1 n n 1 n n 1 n2 24
6、收获。
（1）这次试验我了解到连续系统仿真不同的数值积分法在精度上，步长稳定区以及稳 定性上有了较为直观的理解， 可以基本上掌握了不同数值积分法的选取在稳态精度上的不同。 （2）本次试验在 matlab 编程上有了较高的提高，在编程过程中，对 matlab 编程会出现 的问题以及解决方法有了基本上的掌握，收获较大。
P yn yn f (tn , yn )h 1 改进欧拉法： ； h C P P yn 1 yn [ f (tn , yn ) f (tn 1 , yn 1 )] 2
k1 f (tn , yn ) k2 f (tn h , yn h k1 ) 2 2 h h 四阶龙哥库塔法： k3 f (tn , yn k2 ) ； 2 2 k4 f (tn h, yn hk3 ) yn 1 yn h (k1 2k2 2k3 k4 ) 6
（4）四阶亚当姆斯预估校正法：
clear all close all a=[0,1,0;0,0,1;-22.06,-27,-10]; b=[0;0;40.6]; u=1; for k=1:4 if (k==1) h=0.01; elseif (k==2) h=0.1; elseif (k==3) h=0.105; else (k==4) h=0.21; end x=[0;0;0]; y=x'; t(1)=0; for i=1:3 t(i+1)=h*i; k1=a*x+b*u; k2=a*(x+h/2*k1)+b*u; k3=a*(x+h/2*k2)+b*u; k4=a*(x+h*k3)+b*u; x=x+h/6*(k1+2*k2+2*k3+k4); y=[y;x']; y1=y'; end for i=4:1000 t(i+1)=h*i;
4、实验步骤：
（1）编写 matlab 程序，每种方法选取 0.01、0.1、0.5*h 临界、h 临界四种步长求数值 解，分别与解析解比较。画出 matlab 程序的流程图。
（2）结果输出：生成数据文件，将四种数值方法的结果分别与解析法的计算结果存在 下 面的 ASCII 码 数据 文件 中，分别 命名为： eulerdata.dat ， eulerpdata.dat ， rkdata.dat ， adamsdata.dat。 已生成，例如下图为龙哥库塔法、h=0.01 时的部分数据文件，从数据可知，该方法的 仿真结果较好，与数值解的差别不大。
（2）改进欧拉法：
clear all close all
a=[0,1,0;0,0,1;-22.06,-27,-10]; b=[0;0;40.6]; u=1; for k=1:4 if (k==1) h=0.01; elseif (k==2) h=0.1; elseif (k==3) h=0.16; else (k==4) h=0.321; end x=[0;0;0]; y=x'; t(1)=0; subplot(2,2,k) for i=1:1000 t(i+1)=h*i; x1=x+(a*x+b*u)*h; x=x+[(a*x+b*u)+(a*x1+b*u)]*h/2; y=[y;x']; end t1=0:0.01:10; y1=1.84-4.95*t1.*exp(-1.88*t1)-1.5*exp(-1.88*t1)-0.34*exp(-6.24*t 1); plot(t,y) hold on plot(t1,y1,'k-.','linewidth',2) xlabel('time(sec)') ylabel('·ùÖµy') legend('y','doty','ddoty','y-shuzhi') if (k==1) title('eulerp h=0.01') axis([0,10,-1.5,3.5]) elseif (k==2) title('eulerp h=0.1') axis([0,10,-1.5,3]) elseif (k==3) title('eulerp h=0.16 axis([0,10,-1.5,2]) else (k==4) title('eulerp h=0.321 临界') axis([0,50,-20,5]) 0.5*临界')
控制系统仿真实验
姓名： 王体强 班 级 ： 120324 学 号 ： 12031190 指导教师： 王卫红 时 间 ： 2015 年 4 月 13 日
控制系统仿真实验
1、实验目的：
进一步掌握数值积分法； 进一步掌握 MATLAB 软件的使用方法。
2、百度文库验设备：
数字计算机，MATLAB 软件
3、实验预备：
end end
（3）四阶龙哥库塔法：
clear all close all a=[0,1,0;0,0,1;-22.06,-27,-10]; b=[0;0;40.6]; u=1; for k=1:4 if (k==1) h=0.01; elseif (k==2) h=0.1; elseif (k==3) h=0.223; else (k==4) h=0.447; end x=[0;0;0]; y=x'; t(1)=0; subplot(2,2,k) for i=1:1000 t(i+1)=h*i; k1=a*x+b*u; k2=a*(x+h/2*k1)+b*u; k3=a*(x+h/2*k2)+b*u; k4=a*(x+h*k3)+b*u; x=x+h/6*(k1+2*k2+2*k3+k4); y=[y;x']; end t1=0:0.01:10; y1=1.84-4.95*t1.*exp(-1.88*t1)-1.5*exp(-1.88*t1)-0.34*exp(-6.24*t 1); plot(t,y) hold on plot(t1,y1,'k-.','linewidth',2) xlabel('time(sec)') ylabel('·ùÖµy') legend('y','doty', if (k==1) title('rk h=0.01') 'ddoty','y-shuzhi')
（3）给出仿真结果曲线图。
欧拉法：
改进欧拉法：
四阶龙格库塔法：
四阶亚当姆斯预估校正法：
实验结果得到四阶亚当姆斯预估校正法的临界步长约为 0.21。
5、实验结果及分析：四种方法的精度比较、稳定性比较，步长 稳定区是否符合实验前的理论分析值。其他需要总结的问题。
精度比较： 从保存的数据文件来看， 四阶龙格库塔法和四阶亚当姆斯预估校正法的精度 较高，大约数量级在 104 ，改进欧拉法次之，在 103 数量级上，欧拉法最低，在 102 数量 级上，与理论的结果相符程度较好。 稳定性比较：从实验结果看，四种方法在步长一样且较小的情况下，都能较快的达到稳 定。而且，同一种方法随着步长的增加，达到稳态的时间会增加。但是由于四种方法的临界 步长不相同，因而随着步长的增加，达到稳态的时间有所不同。 步长稳定区：从实验仿真的结果可以看出，实际的临界步长符合实验前理论的分析值。 而且四阶龙格库塔法的临界步长最大，四阶亚当姆斯预估校正法的临界步长最小。 本实验中四阶亚当姆斯预估校正法不能自启动， 需要用四阶龙格库塔法计算其起始的值。
（1）将传递函数化为一阶微分方程组（即状态方程） ； 设： y1 y ， y2 y ， y3 y
1 0 y1 0 y1 0 y 0 0 1 2 y2 0 U 22.06 27 10 40.6 y3 y3
axis([0,10,-1,3]) elseif (k==2) title('rk h=0.1') axis([0,10,-1,3]) elseif (k==3) title('rk h=0.223 axis([0,10,-1,3]) else (k==4) title('rk h=0.447 临界') end end 0.5*临界')
；
（3）理论分析：计算系统特征值。结合系统特征值，对四种方法的稳定性进行分析，确定 每种方法步长的取值范围，即 h 临界。
1.8680 用 matlab 计算系统的特征值为： D 1.8928 ，可见系统的特征值有三个，我们取 6.2392
最大的特征值为： = 6.2392 。由四种方法的稳定区域可得： 方法 欧拉法 改进欧拉法 四阶龙格库塔法 四阶亚当姆斯法-显式 四阶亚当姆斯法-隐式 稳定区域 步长 h 临界步长 0.3206 0.3206 0.4456 0.048 0.2944
(2,0)
(2,0) (2.78,0)
0 h 0.3206
0 h 0.3206 0 h 0.4456
0 h 0.048 0 h 0.2944
( 310 , 0)
( 90 49 , 0)
由于四阶亚当姆斯预估校正法综合受到显式与隐式的影响，其临界步长应在 0.048~0.2944 之间，具体值需用实验验证出来。
此处 U 为单位阶跃信号。
1 0 0 0 0 1 ，b 设： a 0 0 。 22.06 27 10 40.6
（2）分别写出四种方法的计算公式； 欧拉法： yn 1 yn f (tn , yn )h ；
附：matlab 程序 （1）欧拉法：
clear all close all a=[0,1,0;0,0,1;-22.06,-27,-10]; b=[0;0;40.6]; u=1; for k=1:4 if (k==1) h=0.01;
elseif (k==2) h=0.1; elseif (k==3) h=0.16; else (k==4) h=0.321; end x=[0;0;0]; y=x'; t(1)=0; subplot(2,2,k) for i=1:1000 t(i+1)=h*i; x=x+(a*x+b*u)*h; y=[y;x']; end t1=0:0.01:10; y1=1.84-4.95*t1.*exp(-1.88*t1)-1.5*exp(-1.88*t1)-0.34*exp(-6.24*t 1); plot(t,y) hold on plot(t1,y1,'k-.','linewidth',2) xlabel('time(sec)') ylabel('·ùÖµy') legend('y','doty', if (k==1) title('euler h=0.01') axis([0,10,-1.5,3.5]) elseif (k==2) title('euler h=0.1') axis([0,10,-1.5,4.5]) elseif (k==3) title('euler h=0.16 0.5*临界') axis([0,10,-1.5,7]) else (k==4) title('euler h=0.321 临界') axis([0,20,-20,20]) end end 'ddoty','y-shuzhi')