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

实验三 利用数值积分算法的仿真实验

学号： 姓名： 学院：

一． 实验目的

1) 熟悉MATLAB 的工作环境；

2) 掌握在MATLAB 命令窗口调试运行程序；

3) 掌握M 文件编写规则及在MATLAB 命令窗口运行程序；

4) 掌握利用欧拉法、梯形法、二阶显式Adams 法及四阶龙格库塔法构建

系统仿真模型的方法,并对仿真结果进行分析。

二．实验内容

电路如图１所示电路进行仿真试验。电路元件参数：V E 1=，Ω=10R ，

H L 01.0=，F C μ1=。电路元件初始值：A i L 0)0(=，V u c 0)0(=。系统输

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

DC

R

C )

(t u c +-

)

(t i L 图１ RLC 串联电路

E

三． 实验步骤

1. 求连续系统传递函数

根据所示电路图，我们利用电路原理建立系统的传递函数模型，根据系统的传递函数是在零初始条件下输出量的拉普拉斯变换与输入量的拉普拉斯变换之比，可得该系统的传递函数：

LC

Ls R s LC

s E s U s G C /1//1)()()(2

++==

2. 离散系统仿真模型

在连续系统的数字仿真算法中,较常用的有欧拉法、梯形法、二阶显式Adams 法及显式四阶Runge-Kutta 法等。欧拉法、梯形法和二阶显式Adams 法是利用离散相似原理构造的仿真算法，而显式四阶Runge-Kutta 法是利用Taylor 级数匹配原理构造的仿真算法。对于线性系统，其状态方程表达式为:

()()()

()()()

t t t t t t ?=+?

=+?x Ax Bu y Cx Du 00)(x x =t 式（3-1）中，[]T

n t x t x t x )()()(21 =x 是系统的n 维状态向

量,[]

T

m t u t u t u t )()()()(21 =u 是系统的m 维输入向量，

[]T

r t y t y t y t )()()

()(21 =y 是系统的r 维输出向量。A 为n n ?阶参数矩阵，

又称动态矩阵，B 为m n ?阶输入矩阵，C 为n r ?阶输出矩阵，D 为m r ?阶交联矩阵。根据图所示电路，系统状态方程模型:

)

()()()(t t E t t Cx y B Ax x

=+=

式中，状态变量T C L T u i x x ],[],[21==x ，输出变量C u t =)(y ，系数矩阵为：

??????--=0/1/1/C L L R A ，??

?

???=0/1L B ，[]10=C 。 (1) 欧拉法

利用前向欧拉法构建线性系统的仿真模型为：

()1111

m m m m m

m m m h h h ++++?=+=++??

=+??x x x 1A x Bu y Cx Du 式中，h 为积分步长，1为单位矩阵。利用后向欧拉法构建线性系统的仿真模型为：

()()111

1111

m m m m m m m m h h h -++++++?=+=-+??

=+??x x x 1A x Bu y Cx Du 对于前向欧拉法，系数矩阵为：A 1A h z +=，B B h z =，C C =z ，D=0。对于后向欧拉法，系数矩阵为：()

1

z h -=-A 1A

，()

1

z h h -=-B 1A

B ，

C C =z 。

(2) 梯形法

利用梯形法构建线性系统的仿真模型为：

()()111

111112222m m

m m m m m m m m m h h h h ++-+++++?=++??

????????

=-+++???

? ????????

?=+???

x x x x x 1A 1A x B u u y Cx Du 对图所示的系统，利用梯形法构造的系统差分方程具有形式：

111

2m t m t m t m E +++=+=x A x B y C x

其系数矩阵为：()()1

/2

/2t h h -=-+A 1A 1A

，，

()

1

/2/2

t h h -=-B 1A B ，t =C C ，D = 0。

（3）二阶显式Adams 法

利用二阶显式Adams 法构建线性系统的仿真模型为：

112111

2316512m m m m m m m m h

+--+++???=+-+????

?=+?x x F F F y

Cx Du

式中： ()()()

111222m m m m m m m m m ------?=+??

=+??=+??F Ax Bu F Ax Bu F Ax Bu

二阶显式Adams 法为多步计算方法，利用多步计算方法对系统进行仿真时，需要与之具有相同计算精度的单步计算方法辅助计算。二阶显式Adams 法的计算精度为二阶，可以采用梯形法或改进的Euler 法等辅助计算。利用改进的Euler 法构建线性系统的仿真模型为：

()1211112(,)(,)()2m m m m m m m m f t x E f t h x h h E h +?

?==+??

=++=++??

?=++??

k Ax B k k A x k B x x k k 其中，0,1m =。由式计算出1x 和2x 后，便可以转入由二阶显式Adams 法构造的离散系统模型计算，即系统差分方程。其计算方程为：

()()()

1122m m m m m m E E E ----?=+??

=+??=+??F Ax B F Ax B F Ax B 11211

2316512m m m m m m m h +--++???=+-+????

?=?x x F F F y

Cx （1m >） （4）显式四阶Runge-Kutta 法

利用显式四阶Runge-Kutta 法构建线性系统的仿真模型为：

()()()121132243311

2341

1

(,)(,)/222(,)/222

(,)(22)6m m m m m m m m m m m m m m m m f t E

h h

f t h E h h

f t h E f t h h h E h

+++?==+??=++=++???=++=++??

?=++=++??=++++???=?k x Ax B k x k A x k B k x k A x k B k x k A x k B x x k k k k y Cx

3. 连续系统解析解

连续系统输出响应)(t u c 的解析解为：

))/sin (cos 1()(ωωωx t t e U t u at s c ?+?-=-

其中，L

R

x 2=，2

21??

? ??-=

L R LC

ω 4. Matlab 仿真

利用Matlab 软件进行编程仿真，分别利用欧拉法、梯形法、二阶显式Adams 法及显式四阶Runge-Kutta 法对电路进行仿真得到图形曲线，比较分析算法特性与积分步长h 的关系。

取不同的积分步长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

四、实验结论

随着步长的增加，各种方法逐渐与曲线分离开来，仿真效果逐渐变差。故仿真步长越小，仿真效果越好。对比分析几种方法，有：

难易性：欧拉法为单步计算法，用到一个过去的值，只要给定初始值就能迭代进行，为自动起步的。梯形法也为单步法，可以自动起步计算。梯形法用了两条折线所谓面积来近似。二阶显示Adams法需要知道k个初始值，不能自起步，二次函数很复杂，常用单步法来得到所需的几个初始值，很为复杂。显式四阶Runge-Kutta法建模最为复杂，仿真时间也较长。

稳定性：我们看到欧拉法最先与解析解分离开来，最后振荡开来，故稳定性最低；显式四阶Runge-Kutta法、二阶显示Adams法和梯形法而则在较长时间可与曲线拟合收敛，故稳定性好。

精度：显式四阶Runge-Kutta法，二阶显示Adams法模型的精度较高，函数曲线与真实曲线较为接近。其次精确的是梯形法，取梯形面积，误差也较小。前向欧拉法和后向欧拉法模型的精度最低。

离散时间间隔：二阶显示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/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,3),plot(t,y,'g',t,y3,'r')

legend('y解析解,','y3梯形法')

title('梯形法')

subplot(2,3,4),plot(t,y,'g',t,y4,'r')

legend('y解析解,','y4Adams法')

title('二阶显式Adams法')

subplot(2,3,5),plot(t,y,'g',t,y5,'r')

legend('y解析解,','y5龙格库塔')

title('显式四阶Runge-Kutta法')