动态系统响应分析的数值法及其程序设计

1、第三章动态系统响应分析的数值方法 第二章讲述了针对一个动力学系统的数学模型建立的基本方法，接下来的问题是需要针对所 建立的方程的求解问题。有些简单问题我们可以通过数学的方法给出其解析解，但是对于更复杂的 问题，要得到其解析解是非常困难的事情，然而针对系统的仿真问题，其理论基础是面向数学模型 的数值解，因此数值解法奠定了计算机仿真的最重要的计算基础。 本章将介绍连续系统的数学模型，然后介绍几个微分方程的数值解法。在此基础上可以建立 面向微分方程的或传递函数的数字仿真程序 我们知道，仿真模型不是唯一的形式，但他们都是由积分器、加法器和系数器构成的。其中 积分器是仿真系统中最重要的环节，N阶微分方程

2、有N个积分器，为了 在数字机器上对它进行仿真，最直观的想法就是构造出N个积分器，也就是进行N次 积分运算，因此这里涉及到数值积分和微分的一些方法。 3-1数值积分法 求定积分的近似值的数值方法。即用被积函数的有限个抽样值的离散或加权平均近似值代替 定积分的值。求某函数的定积分时，在多数情况下，被积函数的原函数很难用初等函数表达出来， 因此能够借助微积分学的牛顿莱布尼兹公式计算定 积分的机会是不多的。另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函 数，对这类函数的定积分，也不能用不定积分方法求解。由于以上原因，数值积分的理论与方法一 直是计算数学研究的基本课题。对微积分学作出杰

3、出贡献的数学大师，如牛顿.欧拉、高斯等人也 在数值积分这个领域作出了各自的贡献，并奠定了它的理论基础。 1欧拉法（折线法） 假设有一阶微分方程为： 业二 f (t,y) dt 选择初时时刻to，在此位置做近似做切线有: 包二3匹二f(t,y。)，可以解出 dt ti -to yA yof(t,y) (t -to)、y (tj,令7才为计算步距(步长)写成递推形式：y(k 1) = y(k) hf (k)其中k是离散点的取值 2矩形法 写成积分的形式: 假设一阶微分方程: y(t)二 y o t ft,y (t) dt 这里y。是初值，由于不能得到y(t)的值，则积分项是未知的，为了得到近似值，

4、可以将积分 项用矩形计算公式来近似计算，即： t ft,y(t)dt: fto,y(to)(t-tO 将其写出递推公式: y(tiJf (ti5y(tj(tiiti) =yi 或:y(k 1) = y(k) hf (k) 对比欧拉折线方法可知，两种方法是一样的， 3梯形法错误！未找到引用源。 为了得到更精确的计算结果，采用梯形方法，将积分表f (y(ti),ti)f(y(tii),ti 达式写成如下形式 ft,y(t)dt： 2fto, y(to) ft-y(t)(t -to) 儿汀丁“)5川 问题的提出，在计算丫 “时的表达式里要用到y J的值，实际上，上面的式子 的积分项并不能计算结果，通

5、常采用迭代运算方法： yh = yi + hf (ti, y.) 当i =0时，这个式子里面的t, y为初值，实际上，这个表达式是前面所讲的矩形法公式， 有了这个项，我们就可以计算下一步。 y(2=yi 十 2h【f(l), yJ + f(i)+ yhy(f =yi+*hf(h, y+ 鮒+ 丫加 在上面这一组式子中，其中的第二项可以认为是初值，第三项是要用到前一次的迭代运算。 如果迭代序列是收敛的，则极限肯定存在，当k趋向无穷大时，可以得 到这个精确结果。 怎样给出迭代的精度呢？可以采用前一次的计算值和后一次的迭代值的差来控制。 即：成6 d为小量 在实际计算中，我们认为，迭代一次或者两次就

6、认为满足精度。所以，可以写出下面的递推 公式。 y(kAy(K-1)-f(k f(k) 4差分公式 基于数值积分的欧拉方法，我们可以得到近似的微分计算方法，用函数 y二f(x)的离散数据点Xkf(Xk) k=1,2,.n来近似表达f(x)在节点处的微分y-叫做数值微 分。 (1)一阶系统的数值微分公式 设：=f (t,y)应用离散数据点的近似微分，可以表示为： dt dy y(k) -y(k-1) dt t(k)-t(k-1) 令：T =t(n) -t(n称为采样时间 写成递推公式则有：(差分形式)：y (k)二啲 除此之外还有一阶中心差分形式：y (k)二”一 ”k丄w 2T (2)二阶系统

7、的数值微分公式 根据一阶差分公式，我们容易得到二阶差分如下： 鱼v(k) v(k1)v(k) v(k1)y(k-1)-y(k-2)_ dt2 %)TTT/T V(k) 2y(k 1) v(k 2) T2 ，根据这个原理，还可以得到高阶数值微分公式。 在SIMULINK中，我们经常采用一种叫做单位延迟的离散模块，这样就可以方便 的建立基于数值积分和微分的仿真模型了。 例如设：f(t)=sint分别对f(t)的近似积分和近似微分为： 则近似差微商的计算公式为： y (k)二丄血來 T 近似积分计算公式为 y(k)z:y(k-1)f(k)T 在离散模块库中的unit delay (单位延迟模块厂在该

8、模块中的采样周期(sample time)和递推公 式中的采样周期T相同，在下图中使用了采样周期为0.01秒的正弦波的积分和导数离散仿真模型 积分模型 f(k) Sine Wave y(k)=f(kH(k-1)/T y(k)=y(k-1)+f(k)*T 差商模型 近似积分结果近 似 差 商 结 果 3-2龙格库塔方法 龙格库塔方法由于它具有较高的精度而在计算机仿真中就占有重要的地位，在此做简要介绍 1二阶龙格库塔方法 将上述梯形法做进一步推广，可以得到龙格 库塔方法，对于: (t, y): y(t = to) = y。 dt 假设我们从to跨出一步，在2 to h时的解为 y, = y(t0

9、h),在to附近做泰 劳展开，只保留二次项，则有 1 yi=y f(t。，yOh ( 2 -t 吋矽 2 )h h由 t y=y (3-38) 我们假设这个解可以写成： %二办f (a,ki azk?- ki 二 f (to, y0) k2 二 f (tobh,y。bzkih) : f(to,y。)bi h Ft t=to b?k, cf |h Gy# (3-39) 即有: Cf yi ny0 ahf (t,y。) a2h f(t。，y。)bi cf h b?ki h 醴 (3-40) 将(340 )与(338 )对比可知道, a a2 =1a2b =1/2 322 =1/2 四个未知数，有三

10、个方程，可以自由的选择一个值令: 0 =b2=1于是有: a a?=1/2 于是有一组方程为： yi =y。- f (ki k2) 2 ki= f (to,yo) kA f (toh,yo kih) 由于计算时，只取了 h和h?项，所以，上面的公式称为二阶龙格库塔方法。 2四阶龙格库塔方法 二阶龙格库塔法仍然存在精度不高的缺点，通常采用的是四阶龙格库塔 方法。下面直接给出公式： yi 二 y0 n (ki 2k2 2k3 k4) 6 k2 二 f (to 2, yohki) h h k3) L Oc，y ck2， 匕二f (toh, yhk3) 3.尖于数值积分的精度 在实际应用中，为了提高仿

11、真精度，四阶龙格库塔法的计算步距是可以变化的，也就是在计 算过程中可根据所要求的精度对步长作适当改变，这对于精确计算有很大的好处。因为一般的系统 在阶跃函数作用下的过渡过程一开始变化比较快，而最后变化将趋于平缓(详见第八章)，为了达到 较高的精度必须按起始段来选择步距h， 显然这个步距将很小，这样在平缓阶段，就显得十分浪费。因此，在作精确仿真计算时，我们可以 采用变步长的办法，具体步骤如下： (1) 按系统的过渡过程情况将它分成几段，每一段都预先先给一个步长h。 (2) 当过渡过程进入该段时，先要做步长的计算。即用h作一次计算，然后再 用h作一次计算，求得二者之差的绝对值，如果它小于某一预定的

12、值，则认为h这 2 个步长符合精度的要求，即可用此步长继续算下去，如果二者之差的绝对值大于该予定值，则认为 h这个步长过大，那么用h/2作一次计算，再用h/4作二次计算，再进行上述比较。当然，这样 一来，计算的工作量将会增加。 对于工程设计来讲，计算精度往往并不要求十分高，比如误差不超过0.5 % 也就令人满意了。为此，一般常采用定步长的计算方法，作为一个经验数据，当用四阶龙格一一库 塔法进行计算时h可选为(5 ,如)。 1040 其中：tc为系统在阶跃函数作用下的上升时间 tn为系统在阶跃函数作用下的过渡过程时间 最后还要特别指出一点：当采用龙格一一库塔法对连续系统进行数字仿真时，为了保证计

13、算 的稳定性，必须要求所选择的步长应小于系统中最小时间常数之两倍，可保证计算是稳定的。 另外，步长h的选择会影响系统的稳定性问题，各种方法的稳定性与步距的大小有尖，详细情 况可以参见有尖书籍。 3-3 面向微分方程的仿真程序设计 根据上一节可知，如果能够将一个连续系统表示成一组一阶微分方程，且每个一阶方程都 可以用数值积分的方法来求解，则便可以求出整个系统的全部解，包 括各阶导数的值都能给出。这种方法就是面向微分方程的仿真程序设计 龙格库塔法： 在数值仿真中，通常采用四阶的龙格库塔方法。龙格-库塔方法的基本思 想是在区间内多预报几个斜率的值，然后将它们加权平均。其表达式如下： ki=f(Xi,

14、 yJ k? = f(X h/2, yi h ki kA f(xi h/2, yi h k2 k4- f (Xih, yih yii* i =0,1,2 h Jki 2k? 2k3 kq 面向微分方程的数字仿真程序结构，由于面向微分方程的数字仿真模型还有状 态空间方法和传递函数方法等情况，因此面向微分方程的仿真程序结构也有两 种。在此仅介绍尖于微分方程(组)与状态空间方法的数字仿真程序设计。下面我们通过几 个 例题来说明，怎样建立微分方程的仿真程序。 1，用四阶龙格库塔方法求解一阶微分方程程序设计 在此应用了 VB编程方法，对于熟悉其他计算机语言的 学者，可以参考其中程序设计的基本思想。 计算

15、微分方程： dy y -2t/y 在 y (0) =1 t =0,3 dt 取步长为h =0.1的数值解。从t =0到t =3各结点上的数 值解，请读者将本结果和欧拉方法以及梯形法相比较。注： 这个方程有精确解为： y = 2t 1 新建立一窗体，在窗体上加一个列表框，一文本框数组 (03 )，分别存放开始时间、终了时间、采样时间间隔和初 始值在窗体的通用事件过程中添加函数过程(给出微分方程 模型) 在窗体上添加一个命令按钮Commma nd1将原代 码写在命令按钮的单击事件过程中。调试程序，观察运行的结 果的正确性 阿间T Y值 1 0 0 1 1 0. 1000: 1.0954 0.200

16、0; L1832 A 0.3000: 1 2B49 0.4000; 1. 3416 11 0. 5000: 1 4142 -I 0.6000; 1.4832 n 7nnn. 1 5432 I i 0.6000; 1.6125 d 0. 9000; 1 5T33 1.0000: L7321 L 1000; 1 7BE9 1.2000： 1-8439 L3000; 1 3974 1.4000; L9494 1.50CO： 2.3000 1.8000; l.TOCO; 2.0976 1 8000; 2.1446 1 9000； z. itn 2.0000; 2.2361 21000： 2.2304

17、2,2QC0; 2,3230 2 3000t 2. 3665 2.4000: 2. 4064 2.5000; 2. 449 2. 5692 2.90CO; 2.6079 3.0000; 2.6460 龙格库塔方法 Private Function fn1 (x, y) fn1 = y 2 * x/y End Function Private Sub Command1_Click() Dim to As Long Dim tn As Long Dim h As Single Dim yO As Long Dim i As Integer Listl .Clear Dim x, y, n As I

18、nteger, X1, Y1, k1, k2, k3, k4 As Long 数组控件 fVal转换函数 循环步长 Listl .Additem (“龙格库塔方法 J tO = Val(Text1 (O).Text) tn = Val(Text1 (1 ).Text) h = Val(Text1(2).Text) If h = 0 Then MsgBox (u请输入步长h”) Exit Sub End If 给定初值 yO = Val(Text1 (3).Text) n = (tn- tO)/ h循环次数 x = Val(tO) y = y Listl.Addltem (n 时间 Tfl Ne

19、xt i End Sub Private Sub Form Load() Textl (O).Text = 0 Textl (1).Text = 3 Textl (2).Text = 0.1 Textl (3).Text = 1 End Sub 为了和仿真结果进行比较，下面建立了仿真模型，在该模型中，将积分器的初 始条件设置为1，仿真时间设置为03秒，计算结果如下。 根据VB的数值计算和Sinulink仿真结果相比较，有较高的精度。 2 一阶微分方程组的四阶的龙格库塔法程序设计设一阶微分方程组形式如下: dxi Xi,X2xji=1,2. n dt 下面我们仅考虑有两个方程的情况，设初值为：

20、y = f (t, y, Z)y(0) = y z=g(t,y,z)z(0)=Zg 根据上面推导则可以得到四阶龙格库塔法的格式为： y“汕 5丫匕2*k2 2*ksk4)/6 z 1 二乙 h* (a 2* P22* p3 Pj /6 其中： ki 二 f (ti, yZ) k2 二 f(ti, Vi h * ki /2, z h * ki /2) 2 k3二 f(ti,yih*k2/2, Zih*k2/2)2 f 魚 h, yih* k3,Z h* k?) Pi 二 g(tyzj pn 2 h*Pi2Zi hp P3=g(ti 卜 h* P2/2,Zh* P2/2) P4=g(tih, yi

21、 h* P3jZ h* Pa) i =0,1,2 不难将上述情况推广到多维一阶方程组的情况, F面给出了计算K(i, j)的参考子程序程序代码。 例如：设有n个一微分方程组 Xi = x3 X2 =X4 277 X3 p X-IX4 Xi x4 二IX3X4 / Xi 初值为：Xi(0)=1 X2(0)=0.5Xa(0) =1 X4(0)=1.5 首先按微分方程的个函数定义四个子函数分别来计算 Kj,(在这里使用了 K,S,T,U来代替了 kii，ki2，K3，ki4),VB程序代码如下: Private Function fXni(xi, x2, x3, x4, T) As Single f

22、Xni = x3 End Function Private Function fxn2(xi, x2, x3, x4, T) As Single fxn2 = x4 End Function Private Fun ction fxn 3(xi, x2, x3, x4, T) As Si ngle fxn3 = -i / xi / xi + xi * x4 End Function Private Function fxn4(x1, x2, x3, x4, T) As Single fxn4 = -2 * x3 * x4 / x1 End Function Private Sub Comman

23、d1_Click() Dim x1 As Single, x2 As Single, x3 As Single, x4 As Single Dim TO As Single, TN As Single, H As Single, T As Single Dim N As Integer, i As Integer Dim K1, K2, K3, K4 As Single Dim S1, S2, S3, S4 As Single Dim T1, T2, T3, T4 As Single Dim U1, U2, U3, U4 As Single T0 = 0:TN = 3:x1 = 1:x2 =

24、0.5: x3 = 1: x4 = 1.5 H = 0.01 N = GN TO)/ H T = TO For i = 1 To N K1 = fxn1(x1, x2, x3, x4, T) S1 = fxn2(x1, x2, x3, x4, T) T1 = fxn3(x1, x2, x3, x4, T) U1 = fxn4(x1, x2, x3, x4, T) K2 = fxn1 (x1 + H * K1 / 2, x2 + H * S1 / 2, x3 + H * T1 / 2, x4 + H * U1 / 2, T + H / 2) S2 = fxn2(x1 +H* K1 / 2,x2

25、+H* S1 / 2,x3 +H* T1 / 2,x4+H * U1 / 2,T +H/ 2)T2 = fxn3(x1 +H* K1 / 2,x2 +H* S1 / 2,x3 +H* T1 / 2,x4+H * U1 / 2,T +H/ 2)U2 = fxn4(x1 +H* K1 / 2,x2 +H* S1 / 2,x3 +H* T1 / 2,x4+H * U1 / 2,T +H/ 2) K3 = fxn1(x1 + H * K2 / 2,x2 + H * S2 / 2,x3 + H * T2/2,x4 + H * U2 / 2, T + H / 2) S3 = fxn2(x1 +H* K2

26、/ 2,x2 +H* S2 / 2,x3 +H* T2 / 2,x4+H * U2 / 2,T +H/ 2)T3 = fxn3(x1 +H* K2 / 2,x2 +H* S2 / 2,x3 +H* T2 / 2,x4+H * U2 / 2,T +H/ 2)U3 = fxn4(x1 +H* K2 / 2,x2 +H* S2 / 2,x3 +H* T2 / 2,x4+H * U2 / 2,T +H/ 2) K4 = fxn1(x1 + H * K3, x2 + H * S3, x3 + H * T3, x4 + H * U3, T + H) S4 = fxn2(x1 + H * K3, x2 +

27、H * S3, x3 + H * T3, x4 + H * U3, T+ H) T4 =fxn3(x1 + H*K3,x2+H* S3, x3 + H * T3, x4 +H * U3, T + H) U4 = fxn4(x1+ H * K3,x2 +H * S3,x3+H *T3,x4 + H * U3, T + H) x1 = x1 + (K1+ 2 * K2 + 2 * K3 + K4) / 6 x2 = x2 + (S1+ 2 * S2 + 2 * S3 + S4) / 6 x3 = x3 + (T1 + 2 * T2 + 2 * T3 + T4) / 6 x4 = x4 + (U1

28、+ 2 * U2 + 2 * U3 + U4) / 6 Listl.Addltem (Format(T, O.OOOO) df3=diff(fs3；f) %对变量t求导数(表示速度表达式) t=0:0.1:500;%设定时间起始、步长和终了时间 y=eval(char(fs3);%求位移的数值解 figure(1)%建立一个绘图 plot(t,y,k)；画图 hold on %保持画好的曲线 opts=simset(solver,ode45,);%设定 ode45 为当前求解器 t1 ,y1 ,s=sim(modefs3,500,opts);调用仿真模型 plot(t1 ,y1 (：,1),V

29、);绘图，其中的因变量取模型的第1个输出序列 hold on opts=simset(solver,ode15s); %重新设定 ode15s 为当前求解器 t2,y2,s=sim(modefs3,500,opts); plot(t2,y2(:51)；r-); hold on axis(249.7 250 0.7785,0.7795) % 显示局部区域 t=0:0.1:500; y=eval(char(df3); figure(2) plot(t,y,k); hold on opts=simset(solver/ode45)；% t1 ,y1 ,s=sim(modefs3,500,opts);

30、 plot(t1 ,y1 (：,2),V);绘图，其中的因变量取模型的第2个输出序列 hold on opts=simset(,solver,ode15s);% t2,y2,s=sim(modefs3,500,opts); plot(t2,y2(:,2),W); hold on axis(249.7 250 -7.793e-4r7.777e-4) % 显示局部区域 位移Lili线 速 度 曲 结果分析：(1)从位移图上可以看到，精确解和ode45、ode15s求解器的结果差另IJ不大。 但是ode15s求解器的运行速度要快，这个结果可以由下面的命令得到： nan ber1=le ngth (y

31、1) nan ber1 = 150508 nan ber2=le ngth (y2) nan ber2 = 86 即ode45的求解点数为150508，而ode15s为86。显然后者的速度比前者的速度要快 (2)从速度曲线上可以看出解析解和ode45的结果差别比较大，(ode45出现明 显的振荡)，可以看出，选择不合适的求解器，会导致仿真结果不理想。 例微分方程组的求解 已知某系统的微分方程组为： 小二月2y y y 2 = - yi 2 _ .2 y y 丿3 i 2 初始条件为：(0)=0, y2(0)=0.5 , y3(0)=0.5求时间区域为t=0,20啲解。 解：建立使用函数过程，然后再通过调用函数来得到仿真结果 (1) 建立子程序 function dy=fa ngc(t,y) % dy(1)=y(2)*y(3); (2) 建立主程序