《仿真技术与应用》实验报告计算机仿真技术实验报告实验三利用数值积分算法的仿真实验实验三 利用数值积分算法的仿真实验实验目的1) 熟悉MATLAB 勺工作环境；2) 掌握MATLAB 勺.M 文件编写规则，并在命令窗口调试和运行程序;3)掌握利用欧拉法、梯形法、二阶显式Adams 法及四阶龙格库塔法构建系统仿 真模型的方法，并对仿真结果进行分析。

实验内容系统电路如图2.1所示。

电路元件参数：直流电压源E =1V ，电阻R=10门，电感L =0.01H ，电容C 二WF 。

电路元件初始值：电感电流i L (0) =0A ,电容电压u c (0) = 0V 。

系统输出量为电容电压u c (t)。

连续系统输出响应u c (t)的解析解为:u c (t)二 U s (1 — e_at (cos ■ t si nt a/ J)三、要求1) 利用欧拉法、梯形法、二阶显式 Adams 法及显式四阶Runge-Kutta 法构建系统仿真模 型，并求出离散系统的输出量响应曲线；(2-1)其中, a =2L图2.1 RLC 串联电路2) 对比分析利用欧拉法、梯形法、二阶显式 Adams 法及显式四阶Runge-Kutta 法构建系 统仿真模型的仿真精度与模型运行的稳定性问题；3) 分别编写欧拉法、梯形法、二阶显式 Adams 法及显式四阶Runge-Kutta 法的.m 函数文 件，并存入磁盘中。

.m 函数文件要求输入参数为系统状态方程的系数矩阵、仿真时间及仿 真步长。

编写.m 命令文件，在该命令文件中调用已经编写完成的上述 .m 函数文件，完成 仿真实验；4) subplot 和plot 函数将输出结果画在同一个窗口中，每个子图加上对应的标题。

四. 实验原理(1) 连续系统解析解连续系统输出响应u c (t)的解析解为：u c (t)二U s (1-e^t (cos t si nt x/ ，))(2) 原系统的传递函数根据所示电路图，我们利用电路原理建立系统的传递函数模型，根据系统的传递函数 是在零初始条件下输出量的拉普拉斯变换与输入量的拉普拉斯变换之比，可得该系统的传 递函数：(3) 系统的仿真模型在连续系统的数字仿真算法中，较常用的有欧拉法、梯形法、二阶显式Adams 法及显式 四阶Runge-Kutta 法等。

欧拉法、梯形法和二阶显式 Adams 法是利用离散相似原理构造的 仿真算法，而显式四阶 Runge-Kutta 法是利用Taylor 级数匹配原理构造的仿真算法。

对于 线性系统，其状态方程表达式为：；X(t) = Ax (t) + Bu (t)W)二 Cx (t)Du (t)其中: x h X 't)X 2(t)…x n (t) T 是系统的n 维状态向量u (t) = U i (t)比⑴ …U m (t )T 是系统的m 维输入向量其中, G(s)二U c (s) E(s)1/ LC s 2 R/ Ls 1/LCX (t o) = X 0X 页，y(t)二y i(t) y2(t)…y r(t)T是系统的r维输出向量A为n n阶参数矩阵，又称动态矩阵，B为n m阶输入矩阵，C为r n阶输出矩阵, D为r m阶交联矩阵。

根据图所示电路，系统状态方程模型：x(t)二Ax(t) B E y(t)二ex(t)式中，状态变量x = [x1,x2]T二[i L,u C]T，输出变量y(t)=u C，系数矩阵为:A「― R/L -1/U D刁/L l 厂r A=| ，B=| ，C = 101」。

]1/C 0」]0 一(1) 欧拉法利用前向欧拉法构建线性系统的仿真模型为：x m 1 二x m * x m h 二—A h x m - h Bu m,y m 1 二Cx m -1 ' DU m -1式中，h为积分步长，1为单位矩阵。

利用后向欧拉法构建线性系统的仿真模型为：. -1X m ・1 二x m - xm』二1 - A h X g h BU m 1.y m ・1 二Cx m -1 ' Du m -1对于前向欧拉法，系数矩阵为：A z = 1 h A，B z二h B，C z二C，D=0。

对于后向欧拉法，系数矩阵为：£ = 1 -h A ' , B z二1 - h A 一1h B，C z二C。

(2) 梯形法利用梯形法构建线性系统的仿真模型为：对图所示的系统，利用梯形法构造的系统差分方程具有形式:X m 1二 A t X m 2B t Ey m二 C t X m -1二1其系数矩阵为：A = H - h A / 21 h A /2 ,,B t 二 1 - h A /2 J h B 2，C t 二 C ，D = 0。

(3)二阶显式Adams 法利用二阶显式Adams 法构建线性系统的仿真模型为：h --]I x 羽=x + 一 23F —16F q + 5F 21 m 12 - m m _1m -2』m+1 = C x m 41 + DU m+1F m = Ax m - BU m式中：Fm」二 A xmdBUmJFm 上二 Ax m/ BU m/二阶显式Adams 法为多步计算方法，利用多步计算方法对系统进行仿真时， 需要与之 具有相同计算精度的单步计算方法辅助计算。

二阶显式Adams 法的计算精度为二阶，可以 采用梯形法或改进的Euler 法等辅助计算。

利用改进的Euler 法构建线性系统的仿真模型为：k 1 = f (t m ，X m ) = Ax m + B E k 2 二 f(t mh ,x mhkj = A x m k [h B E其中，m = 0,1。

由式计算出 片和x 2后，便可以转入由二阶显式 Adams 法构造 的离散系统模型计算，即系统差分方程。

其计算方程为：xm 单 x m 1y^i匚h /r h jh 口 11 — 一 A < 2丿1 + — A < 2丿x m + ㊁ B (Um + U ^1 )号 k 1 k 2)x m号 X m - Xm .1C X m 1D U m 1x mF m = Ax mB EFm 二 Ax m _1 B EFm^= Ax mJ2'B Ey m 1 二 Cx m 1(4)显式四阶 Runge-Kutta 法利用显式四阶Runge-Kutta 法构建线性系统的仿真模型为:k [=f(t m,xm)=Axm+ BEk 2 = f (t mh + —2xm+ h k 1)=A ( X m + k-i h / 2) + B E k 3= f(t m+h, 2xm+h k 2)=A (X m+ k 2h / 2) + B Ek 4 = f (t m+ h, X m + hk3)= A (X m + k3h ) + B EX m41 =x m h + - 6(k 1 + 2k 2 + 2k 3 + k 4)[_y m 舟=CxmH!五. 实验过程1■实验程序(1) 前向欧拉法fun ctio n []=RLC(R 丄,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]; %前向欧拉法%X m -1 X m—23F m -16F125F m _2for i=1:1: nx1(1: n,1) = 0;endfor k=1:mx1(1: n,k+1) = x1(1: n,k) + (A* x1(1: n,k)+B)*h;endfor k=1:1:my1(k) = D*x1(1: n,k);end%解析解%p = R/(2*L);w=sqrt(1/(L*C)-(R/(2*L)F2);for k=1:1:my(k) = U*(1-exp(-p*(k-1)*h) * ( cos(w*(k-1)*h) + si n(w*(k-1)*h)*p/w)); end %输出曲线%for k=1:1:mt(k) = (k-1)*h;endsubplot(2,3,1),plot(t,y,'g',t,y1,'r')legend('y解析解,','y1前向欧拉')title('前向欧拉法')(2)后向欧拉法fun ctio n []=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: nx2(1: n,1) = 0;endA1 = in v(E-A*h);for k=1:mx2(1: n,k+1) = A1*(x2(1: n,k) + B*h);endy2(k) = D*x2(1: n,k);end%解析解%p = R/(2*L);w=sqrt(1/(L*C)-(R/(2*L)F2);for k=1:1:my(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:mt(k) = (k-1)*h;endsubplot(2,3,2),plot(t,y,'g',t,y2,'r')legend('y解析解,','y2后向欧拉') title('后向欧拉法')(3)梯形法fun ctio n []=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: nx3(1: n,1) = 0;endA2 = in v(E-A*h/2);for k=1:mx3(1: n,k+1) = A2*( x3(1: n,k) + B*h + A*x3(1: n,k)*h/2); end for k=1:1:my3(k) = D*x3(1: n,k);end%解析解%p = R/(2*L);w=sqrt(1/(L*C)-(R/(2*L)F2);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:mt(k) = (k-1)*h;endsubplot(2,3,3),plot(t,y,'g',t,y3,'r')legend('y解析解,','y3梯形法')title('梯形法')(4)二阶显式Adams法fun ctio n []=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];% 二阶显示Adams法%for i=1:1: nx4(1: n,1) = 0;endfor k=1:mx4(1: n, k+1) = A2*(x4(1: n,k) + B*h + A*x4(1: n,k)*h/2);endfor k=3:mfm1 = 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;endfor k=1:1:my4(k) = D*x4(1: n,k);end%解析解%p = R/(2*L);w=sqrt(1/(L*C)-(R/(2*L)F2);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:mt(k) = (k-1)*h;endsubplot(2,3,4),plot(t,y,'g',t,y4,'r')legend('y解析解，','y4Adams 法') title('二阶显式Adams 法')(5)四阶Runge-Kutta 法fun ctio n []=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];% 四阶Runge-Kutta 法%for i=1:1: n % 状态变量初值x5(1: n,1) = 0;endfor k=1:mx5(1: n, k+1) = A2*( x5(1: n,k) + B*h + A*x5(1: n,k)*h/2);endfor k=1:1:mk1= 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;endfor k=1:1:my5(k) = D*x5(1: n,k);end%解析解%《仿真技术与应用》实验报告-11 -p = R/(2*L); w=sqrt(1/(L*C)-(R/(2*L)F2); for k=1:1:my(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:mt(k) = (k-1)*h;end subplot(2,3,5),plot(t,y,'g',t,y5,'r')legend('y 解析解,','y5Runge-Kutta 法') title('显式四阶 Runge-Kutta 法')2■仿真图形取积分步长h=2*10-4s ,可以得到以下几个仿真图形: (1)前向欧拉法1x 10-1-2-3-4-50 0.001 0.002 0.003i* 0.0040.0050.006r0.0070.0080.0090.01(2)后向欧拉法后向欧拉法1.6 1.4 1.2 1 0.8 0.6 0.4 0.2y 解析解，-y2后向欧拉“0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01y 解析解， y1前向欧拉《仿真技术与应用》实验报告-12-显式四阶Runge-Kutta 法1 X 10 16前向欧y%法析解y1前向欧拉后向欧拉法0 -1 22梯形法y 解析解，1.5y2后向欧拉x 10■■二 阶显式Adams 法y 解析解， -21y3梯形法1.5y 解析解， -3 1 0y4Adams法-40.5-11一-50.0050.01-20.5/1 _0 00.005-30.01一11I \0.005 0.01-41 \ F0 0.0020.0040.0060.008 0.01-0.5-10.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01(3)梯形法梯形法21.8 1.6 1.4 1.2 --------- y 解析解,--------- y 3梯形法1 0.8 0.6 0.4 0.2 00.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.01(4)二阶显式Adams 法0.5 (26)X 10二 阶 显 式 Adamsy 解析解，y4Adams 法-0.5-1-1.5-2-2.5-3x 10161前向欧拉法y 解析解， y1前向欧拉后向欧拉法-1-2 -3-4-50.0050.011.50.50.0050.001 0.002 0.003 0.004 0.005 y 解析解， y2后向欧拉梯形法1.50.50.010.0050.006 0.007 0.008 y 解析解， y3梯形法0.010.009 0.01(5)四阶 Runge-Kutta 法 四阶 Runge-Kuttay 解析解，y5龙格库 塔------------- y 解析解，------------- y5龙格库塔-10 0.0051.5《仿真技术与应用》实验报告六.实验结论1. 从仿真的稳定性看，当选取不同的积分步长时，欧拉法稳定性最低，梯形法稳定性其次，而显式四阶Runge-Kutta法、二阶显示Adams法稳定性较好。