计算机仿真技术实验报告-实验三

1、计算机仿真技术实验报告实验三利用数值积分算法的仿真实验实验三 利用数值积分算法的仿真实验一. 实验目的1) 熟悉MATLAB的工作环境；2) 掌握MATLAB的.M文件编写规则，并在命令窗口调试和运行程序；3) 掌握利用欧拉法、梯形法、二阶显式Adams法及四阶龙格库塔法构建系统仿 真模型的方法，并对仿真结果进行分析。二. 实验容系统电路如图2.1所示。电路元件参数：直流电压源E = W,电阻R = 10G,电感L = 0.01H,电容C = l/zF。电路元件初始值：电感电流z(O) = OA,电容电压c(0) = 0V 0 系统输出量为电容电压i/c(r)o连续系统输出响应的解析解为：(2

2、-1)uc(t) = Ux(cosof + sin 必xa/劲)其中,Ra =2L1R LC69 =9RLM)dcO图2.1 RLC串联电路三、要求1) 利用欧拉法、梯形法、二阶显式Adams法及显式四阶Runge-Kutta法构建系统仿真模 型，并求出离散系统的输出量响应曲线；2) 对比分析利用欧拉法.梯形法.二阶显式Adams法及显式四阶Runge-Kutta法构建系 统仿真模型的仿真精度与模型运行的稳定性问题；3) 分别编写欧拉法.梯形法、二阶显式Adams法及显式四阶Runge-Kutta法的m函数文 件，并存入磁盘中。ni函数文件要求输入参数为系统状态方程的系数矩阵、仿真时间及仿 真

3、步长。编写m命令文件，在该命令文件中调用已经编写完成的上述m函数文件，完成 仿真实验；4) subplot和plot函数将输出结果画在同一个窗口中，每个子图加上对应的标题。四实验原理(1) 连续系统解析解连续系统输出响应/1(0 ：(0儿是系统的厂维输出向量A为nxn阶参数矩阵，又称动态矩阵，为n x m阶输入矩阵，6*为/公阶输出矩阵，D 为rxm阶交联矩阵。根据图所示电路，系统状态方程模型：x(t) = Ax(t) + BEy(t) = Cx(t)式中，状态变量x = xhx2r =iL,ucr.输出变量y(t) = iic,系数矩阵为:-/?/L-1/L 1/LA =,B =1/C00,

4、C = o lo(1) 欧拉法利用前向欧拉法构建线性系统的仿真模型为:卜”+i =兀+ % = ( + Ah)兀” + hB%1九+】=Cxt+i +式中，力为积分步长，2为单位矩阵。利用后向欧拉法构建线性系统的仿真模型为: 1 )h/n+1+ 12(4) 显式四阶Runge-Kutta法利用显式四阶Runge-Kutta法构建线性系统的仿真模型为:K =九,心)=化+ bek2 =+ ，兀” + 人)=A (x)n + kji / 2)+ BEk3 = fg + ，心 + &)= A (xn + k2h / 2)+ BE % = f(tm + h, xm + hk3) = A (xm + k

5、3h) + BE=札+ 匕+汰2 + 2他+心)五.实验过程1. 实验程序(1) 前向欧拉法function =RLC(R,L,C,U,t,h)R=10；L=0.01；C=l. 0e6；U=l；t=0.01；h = 2. Oe-4；m = fix(t/h)；n = 2；A = -R/L -1/L；1/C 0；B = l/L；0；D = 0 1;E = 1 0；0 1;%前向欧拉法%for i=l:1:nxl (l：n,1) = 0；endfor k=l：mxl (l：n,k+l) = xl (l：n,k) + (A* xl (1：n,k)+B)*h；endfor k=l:l：myl(k) =

6、D*xl(l：n,k)；end%解析解%P = R/(2*L)；w=sqrt(1/(L*C)-(R/(2*L)*2);for k=l:l：my (k) = U*(l-exp(-p*(k-l) *h) * ( cos (w*(kl) *h) + sin(w*(kl)*h) *p/w);end%输出曲线%for k=l:l：mt (k) = (k-1) *h；endsubplot (2,3,1) ,plot (t,y,legend(y解析解/ , ryl前向欧拉)titleC前向欧拉法)(2) 后向欧拉法function =RLC(R,L,C,U,t,h)R=10；L二0.01；C=l. Oe-6

7、；U=l；t 二0 01；h = 2.Oe-4；m = fix(t/h)；n = 2；A = -R/L -1/L；1/C 0；B = l/L；0；D = 0 1；E = 1 0；0 1;%后向欧拉法%for i=l：l：nx2(l：n, 1) = 0；endAl = inv(EA*h)；for k=l：mx2(l：n,k+l) = Al*(x2(l：n,k) + B*h);endfor k=l：1：my2(k) = D*x2(l：n,k);end%解析解%P = R/ (2*L);w=sqrt(1/(L*C)-(R/(2*L)2)；for k=l：1：my (k)= LI* (1-exp (-

8、p* (k-1) *h)*( cos (w*(kl) *h)+sin(w*(k-1)*h)*p/w)；end%输出曲线%for k=l:1：mt (k) = (k-1) *h；endsubplot (2,3,2),plot(t,y,* g*,t,y2,r*)legend( y解析解，y2后向欧拉) titleC后向欧拉法)(3) 梯形法function =RLC(R,L,C,U,t,h)R=10；L=0.01；C=l. 0e6；U=l;t 二0. 01；h = 2. Oe-4；m = fix(t/h);n = 2；A = -R/L -1/L；1/C 0；B = l/L；0；D = 0 1；E

9、= 1 0；0 1;%梯形法%for i=l:l：nx3(l：n, 1) = 0；endA2 = inv(E-A*h/2)；for k=l：mx3(l：n,k+l) = A2*( x3(l：n,k) + B*h + A*x3(l：n,k)*h/2)；endfor k=l：1：my3(k) = D*x3(1：n,k)；end%解析解%P = R/ (2*L);w=sqrt (1/ (L*C)-(R/(2*L) 2)；for k=l：l：my (k)= LI* (1-exp (-p* (k-1) *h)*( cos (w*(k-l) *h)+sin(w*(kl)*h)*p/w)；end%输出曲线%

10、for k=l：1：mt (k) = (k-1) *h；endsubplot (2,3,3),plot (t,y,1gf,t,y3,r) legend( y解析解，r, r y3梯形法) title(梯形法)(4) 二阶显式Adams法function =RLC(R,L,C,U,t,h)R=10；0.01；C=l. Oe-6；U=l;t 二0. 01；h = 2. Oe-4；m = fix(t/h)；n = 2；A = -R/L -1/L；1/C 0；B = l/L；0；D = 0 1；E = 1 0；0 1;%二阶显示Adams法%for i=l：1:nx4(l：n, 1) = 0；endf

11、or k=l：mx4(l：n,k+l) = A2*(x4(l：n,k) + B*h + A*x4(l:n,k)*h/2):endfor k=3：mfml = 23*(A*x4(l：n,k) + B)；fm2 = -16*(A*x4(l：n,k-l)+ B)；fm3 = 5*(A*x4(l：n,k-2)+ B)；x4(l：n,k+l) = x4(1：n,k)+(fml+fm2+fm3)*h/12；endfor k=l：1：my4 (k) = D*x4(l：n,k);end%解析解%P = R/ (2*L);w=sqrt (1/ (L*C)-(R/(2*L) 2)；for k=l：l：my (k)

12、= LI* (1-exp (-p* (k-1) *h)*( cos (w*(k-l) *h)+sin(w*(kl)*h)*p/w)；end%输出曲线%for k=l：1：mt (k) = (k-1) *h；endsubplot(2,3,4),plot(t,y,g，t,y4,i)legend ( y 解析解，r,r y4Adams 法 J titleC 二阶显式 Adams 法)(5) 四阶 Runge-Kutta 法function =RLC(R,L,C,U,t,h)R=10；0.01；C=l. Oe-6；U=l;t 二0. 01；h = 2. Oe-4；m = fix(t/h)；n = 2；

13、A = -R/L -1/L；1/C 0；B = l/L；0；D = 0 1；E = 1 0；0 1;% 四阶 Runge-Kutta 法 %for i=l：l：n %状态变量初值x5(l：n, 1) = 0；endfor k=l：mx5(l：n,k+l) = A2*( x5(l：n,k) + B*h + A*x5(l：n,k)*h/2):endfor k=l：1：mkl=A*x5(l：n,k+l)；k2=A*(x5(1：n,k+1)+h*kl/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)

14、+h. *(kl+2*k2+2*k3+k4). /6；endfor k=l：1:my5(k) = D*x5(l：n,k)；end%解析解%P = R/ (2*L);w=sqrt(1/(L*C)-(R/(2*L)2)；for k=l：l：my(k)=U*(1-exp(-p*(k-1)*h)*( cos(w*(k-l)*h)sin (w*(k-1)*h)*p/w)；end%输出曲线%for k=l：1：mt (k) = (k-1) *h；endsubplot (2,3,5),plot(t,y,g*,t,y5,r) legend ( y 解析解，/ y5RungeKutta 法) title(显式四

15、阶 Runge-Kutta 法)2. 仿真图形取积分步长h=2*10 4s,可以得到以下几个仿真图形：(1)前向欧拉法11 1y解侪昶.”仍同吹垃1 1 1 1 1 1r11rrrr11OO.OO10.0020.003O.OCM0.0050.0060.0070.0080.0090.01(2)后向欧拉法y解析解y馮向欧拉0.0010.0020.0030.0040.0050.0060.0070.0080.0090.01(3)梯形法梯形法(4) 二阶显式Adams法x 10二陆 SB re AdamsiAl11ll11lyAdams；/*riirriirrOO.OOl0.0020.0030.0040.0050.0060.0070.0080.009O.Ol(5) 四阶 Runge-Kutta 法六.实验结论1. 从仿真的稳定性看，当选取不同的积分步长时，欧拉法稳定性最低，梯形法稳定性 其次，而显式四阶Runge-Kutta法、二阶显示Adams法稳定性较好。2. 从仿真的难易性看，欧拉法为单步计算法，用到一个过去的值，计算起来比较简单。 而梯形法则是用两条折线所谓面积来近似，与欧拉法相比较为困难。二阶显示Adams法需 要知道k个初始值，不能自起步，二次函数很复杂，因此此方法较复杂。而显式四阶 Runge-Kutta法建模最为复杂，仿真时间也较长。3