计算机仿真实验3数值积分算法

1、仿真技术与应用实验报告实验三 利用数值积分算法的仿真实验学号： 姓名： 学院：一 实验目的1) 熟悉matlab的工作环境；2) 掌握在matlab命令窗口调试运行程序；3) 掌握m文件编写规则及在matlab命令窗口运行程序；4) 掌握利用欧拉法、梯形法、二阶显式adams法及四阶龙格库塔法构建系统仿真模型的方法,并对仿真结果进行分析。二实验内容电路如图所示电路进行仿真试验。电路元件参数：，。电路元件初始值：，。系统输出量为电容电压。三 实验步骤1. 求连续系统传递函数 根据所示电路图，我们利用电路原理建立系统的传递函数模型，根据系统的传递函数是在零初始条件下输出量的拉普拉斯变换与输入量的拉

2、普拉斯变换之比，可得该系统的传递函数：2. 离散系统仿真模型在连续系统的数字仿真算法中,较常用的有欧拉法、梯形法、二阶显式adams法及显式四阶runge-kutta法等。欧拉法、梯形法和二阶显式adams法是利用离散相似原理构造的仿真算法，而显式四阶runge-kutta法是利用taylor级数匹配原理构造的仿真算法。对于线性系统，其状态方程表达式为: 式（3-1）中，是系统的n维状态向量,是系统的m维输入向量，是系统的r维输出向量。a为阶参数矩阵，又称动态矩阵，b为阶输入矩阵，c为阶输出矩阵，d为阶交联矩阵。根据图所示电路，系统状态方程模型: 式中，状态变量，输出变量，系数矩阵为：，。(1

3、) 欧拉法利用前向欧拉法构建线性系统的仿真模型为： 式中，为积分步长，为单位矩阵。利用后向欧拉法构建线性系统的仿真模型为： 对于前向欧拉法，系数矩阵为：，d=0。对于后向欧拉法，系数矩阵为：，。(2) 梯形法利用梯形法构建线性系统的仿真模型为： 对图所示的系统，利用梯形法构造的系统差分方程具有形式： 其系数矩阵为：，d = 0。（3）二阶显式adams法利用二阶显式adams法构建线性系统的仿真模型为： 式中： 二阶显式adams法为多步计算方法，利用多步计算方法对系统进行仿真时，需要与之具有相同计算精度的单步计算方法辅助计算。二阶显式adams法的计算精度为二阶，可以采用梯形法或改进的eul

4、er法等辅助计算。利用改进的euler法构建线性系统的仿真模型为： 其中，。由式计算出和后，便可以转入由二阶显式adams法构造的离散系统模型计算，即系统差分方程。其计算方程为： （）（4）显式四阶runge-kutta法 利用显式四阶runge-kutta法构建线性系统的仿真模型为： 3. 连续系统解析解连续系统输出响应的解析解为： 其中，4. matlab仿真利用matlab软件进行编程仿真，分别利用欧拉法、梯形法、二阶显式adams法及显式四阶runge-kutta法对电路进行仿真得到图形曲线，比较分析算法特性与积分步长的关系。取不同的积分步长h，仿真结果如下：(1) 仿真时间t=0.0

5、1s，积分步长h=10-6s(2) 仿真时间t=0.01s，积分步长h=2*10-6s(3) 仿真时间t=0.01s，积分步长h=10-5s(4) 仿真时间t=0.01s，积分步长h=5*10-5s(5) 仿真时间t=0.01s，积分步长h=10-4s(6) 仿真时间t=0.01s，积分步长h=2*10-4s(7) 仿真时间t=0.01s，积分步长h=4*10-4s4、 实验结论随着步长的增加，各种方法逐渐与曲线分离开来，仿真效果逐渐变差。故仿真步长越小，仿真效果越好。对比分析几种方法，有：难易性：欧拉法为单步计算法，用到一个过去的值，只要给定初始值就能迭代进行，为自动起步的。梯形法也为单步法

6、，可以自动起步计算。梯形法用了两条折线所谓面积来近似。二阶显示adams法需要知道k个初始值，不能自起步，二次函数很复杂，常用单步法来得到所需的几个初始值，很为复杂。显式四阶runge-kutta法建模最为复杂，仿真时间也较长。 稳定性：我们看到欧拉法最先与解析解分离开来，最后振荡开来，故稳定性最低；显式四阶runge-kutta法、二阶显示adams法和梯形法而则在较长时间可与曲线拟合收敛，故稳定性好。 精度：显式四阶runge-kutta法，二阶显示adams法模型的精度较高，函数曲线与真实曲线较为接近。其次精确的是梯形法，取梯形面积，误差也较小。前向欧拉法和后向欧拉法模型的精度最低。离散

7、时间间隔：二阶显示adams法和梯形法对离散时间间隔要求低，前向欧拉法和后向欧拉法由于取的是矩形面积，离散时间间隔要求高。从以上几种方法的仿真，共同规律是在小步长下都收敛，步长越小，误差越小，函数曲线与真实曲线较为接近；步长太大，最后的到的结果不是绝对收敛。同时，同一算法隐式算法要比显式算法稳定性好，计算精度高。matlab仿真程序如下：% 实验三 利用数值积分算法的仿真实验 %function =rlc(r,l,c,u,t,h) r=10; l=0.01; c=1.0e-6; u=1; t=0.01; h = 2.0e-4; m = fix(t/h); n = 2; a = -r/l -1/

8、l;1/c 0; b = 1/l;0; d = 0 1; e = 1 0;0 1;% 前向欧拉法 % for i=1:1:n x1(1:n,1) = 0; end for k=1:m x1(1:n,k+1) = x1(1:n,k) + (a* x1(1:n,k)+b)*h; end for k=1:1:m y1(k) = d*x1(1:n,k); end % 后向欧拉法 % for i=1:1:n x2(1:n,1) = 0; end a1 = inv(e-a*h); for k=1:m x2(1:n,k+1) = a1*(x2(1:n,k) + b*h); end for k=1:1:m y

9、2(k) = d*x2(1:n,k); end % 梯形法 % for i=1:1:n x3(1:n,1) = 0; end a2 = inv(e-a*h/2); for k=1:m x3(1:n,k+1) = a2*( x3(1:n,k) + b*h + a*x3(1:n,k)*h/2); end for k=1:1:m y3(k) = d*x3(1:n,k); end % 二阶显示adams法 % for i=1:1:n x4(1:n,1) = 0; end for k=1:m x4(1:n,k+1) = a2*(x4(1:n,k) + b*h + a*x4(1:n,k)*h/2); en

10、d for k=3:m fm1 = 23*(a*x4(1:n,k)+ b); fm2 = -16*(a*x4(1:n,k-1)+ b); fm3 = 5*(a*x4(1:n,k-2)+ b); x4(1:n,k+1) = x4(1:n,k)+(fm1+fm2+fm3)*h/12; end for k=1:1:m y4(k) = d*x4(1:n,k); end % 四阶runge-kutta法 % for i=1:1:n % 状态变量初值 x5(1:n,1) = 0; end for k=1:m x5(1:n,k+1) = a2*( x5(1:n,k) + b*h + a*x5(1:n,k)*

11、h/2); end for k=1:1:m k1=a*x5(1:n,k+1); k2=a*(x5(1:n,k+1)+h*k1/2); k3=a*(x5(1:n,k+1)+h*k2/2); k4=a*(x5(1:n,k+1)+h*k3); x5(1:n,k+1)=x5(1:n,k+1)+h.*(k1+2*k2+2*k3+k4)./6; end for k=1:1:m y5(k) = d*x5(1:n,k); end % 解析解 % p = r/(2*l); w=sqrt(1/(l*c)-(r/(2*l)2); for k=1:1:m y(k) = u*(1-exp(-p*(k-1)*h) * (

12、 cos(w*(k-1)*h) + sin(w*(k-1)*h)*p/w); end %输出曲线 % for k=1:1:m t(k) = (k-1)*h; end subplot(2,3,1),plot(t,y,'g',t,y1,'r') legend('y解析解,','y1前向欧拉') title('前向欧拉法') subplot(2,3,2),plot(t,y,'g',t,y2,'r') legend('y解析解,','y2后向欧拉') title('后向欧拉法') subplot(2,3,3),plot(t,y,'g',t,y3,'r') legend('y解析解,','y3梯形法') title('梯形法') subplot(2,3,4),plot(t,y,'g',t,y4,'r') legend(