实验三哈工大计算机仿真.docx

实验3 利用数值积分算法的仿真实验(一 实验目的（1）熟悉MATLAB的工作环境；（2）掌握MATLAB的 .M文件编写规则，并在命令窗口调试和运行程序；（3）掌握利用欧拉法、梯形法、二阶显式Adams法及四阶龙格库塔法构建系统仿真模型的方法,并对仿真结果进行分析。二、实验内容系统电路如图3所示。电路元件参数：直流电压源，电阻，电感，电容。电路元件初始值：电感电流，电容电压。系统输出量为电容电压。试利用欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法构建系统仿真模型，并求出离散系统的输出量响应曲线。连续系统输出响应的解析解为： （2-1）其中， ， 。 三、实验要求1）利用欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法构建系统仿真模型，并求出离散系统的输出量响应曲线；2）对比分析利用欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法构建系统仿真模型的仿真精度与模型运行的稳定性问题；3）分别编写欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法的.m函数文件，并存入磁盘中。.m函数文件要求输入参数为系统状态方程的系数矩阵、仿真时间及仿真步长。编写.m命令文件，在该命令文件中调用已经编写完成的上述.m函数文件，完成仿真实验；4）利用subplot和plot函数将输出结果画在同一个窗口中，每个子图加上对应的标题。四、实验原理在连续系统的数字仿真算法中,较常用的有欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法等。欧拉法、梯形法和二阶显式Adams法是利用离散相似原理构造的仿真算法，而显式四阶Runge-Kutta法是利用Taylor级数匹配原理构造的仿真算法。对于线性系统，其状态方程表达式为： (4-1)式（4-1）中，是系统的n维状态向量, 是系统的m维输入向量，是系统的r维输出向量。A为阶参数矩阵，又称动态矩阵，B为阶输入矩阵，C为阶输出矩阵，D为阶交联矩阵。利用前向欧拉法构建线性系统的仿真模型为： （4-2）式中，为积分步长，为单位矩阵。利用后向欧拉法构建线性系统的仿真模型为： （4-3）利用梯形法构建线性系统的仿真模型为： （4-4）利用二阶显式Adams法构建线性系统的仿真模型为： （4-5）式中： （4-6）二阶显式Adams法为多步计算方法，利用多步计算方法对系统进行仿真时，需要与之具有相同计算精度的单步计算方法辅助计算。二阶显式Adams法的计算精度为二阶，可以采用梯形法或改进的Euler法等辅助计算。利用改进的Euler法构建线性系统的仿真模型为： （4-7）利用显式四阶Runge-Kutta法构建线性系统的仿真模型为： （4-8）五、实验结果取不同的积分步长h，仿真结果如下：(1) 仿真时间t=0.01s，积分步长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-4s六、实验结果分析在h=5e-5,h=5e-6,h=5e-7时得到的图像，可以看出，在h=5e-5时，前向欧拉法和后向欧拉法图像明显谁真，而梯形法和二阶显示Adams法图像有轻微失真，步距仍然较大。在步距为时，前向欧拉法和后向欧拉法图像有部分失真，前向欧拉法失真较严重。而梯形法和二阶显示Adams法图像与连续型函数曲线相似度极高，仿真效果非常好。当步距为时，所有方法仿真效果都非常好，与连续型函数曲线相似度都极高。可见，步距时，所有仿真方法都可以应用。从仿真模型实现的难易性、模型的稳定性、模型的精度及离散时间间隔等方面，对比分析上述方法构造的离散系统模型的优缺点。 难易性：通过单个模型的分别仿真，可以得出显式四阶Runge-Kutta法建模最为复杂，仿真时间也较长，对步距要求较低。其次复杂的是二阶显式Adams法和梯形法，仿真时间稍短，梯形法取梯形面积，误差也较小。前向欧拉法和后向欧拉法模型的复杂程度差不多，仿真时间也差不多。模型的稳定性：当步距h=5.0e-5时，前向欧拉法和后向欧拉法明显失真，随着步距的减小，二阶显式Adams法，梯形法和显式四阶Runge-Kutta法的波形变化不大，而前向欧拉法和后向欧拉法的波形得到明显改善。所以显式四阶Runge-Kutta法，二阶显式Adams法和梯形法的稳定性较好，前向欧拉法和后向欧拉法的稳定性较差。模型的精度和离散时间间隔：步距为h=5.0e-6时，显式四阶Runge-Kutta法精度最高，其次是二阶显式Adams法和梯形法。步距为h=5.0e-7时，前向欧拉法和后向欧拉法仿真精度才达到要求。所以，显式四阶Runge-Kutta法，二阶显式Adams法和梯形法模型的精度较高，离散时间间隔要求低，其中，显式四阶Runge-Kutta法模型的精度最高，其次是二阶显式Adams法，由于是二次函数较复杂，函数曲线与真实曲线较为接近；再次精确的是梯形法，取梯形面积，误差也较小；前向欧拉法和后向欧拉法模型的精度较低，由于取的是矩形面积，离散时间间隔要求高。 七、实验结论二阶显示Adams法精度最高，其次是梯形法。在h = 5.0e-6已然如此，当步长在小的时候反而使运行时间延长，效果不一定好。这就说明方法的选取与精度和复杂性，以及可行性等密切相关。八、附录m文件源程序：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/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 y2(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); end 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)*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) * ( 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,