第三章：动态系统响应分析的数值法及其程序设计.doc

第三章 动态系统响应分析的数值方法第三章 动态系统响应分析的数值方法 第二章讲述了针对一个动力学系统的数学模型建立的基本方法，接下来的问题是需要针对所建立的方程的求解问题。有些简单问题我们可以通过数学的方法给出其解析解，但是对于更复杂的问题，要得到其解析解是非常困难的事情，然而针对系统的仿真问题，其理论基础是面向数学模型的数值解，因此数值解法奠定了计算机仿真的最重要的计算基础。 本章将介绍连续系统的数学模型，然后介绍几个微分方程的数值解法。在此基础上可以建立面向微分方程的或传递函数的数字仿真程序. 我们知道，仿真模型不是唯一的形式，但他们都是由积分器、加法器和系数器构成的。其中积分器是仿真系统中最重要的环节，N 阶微分方程有N个积分器，为了在数字机器上对它进行仿真，最直观的想法就是构造出N个积分器，也就是进行N次积分运算，因此这里涉及到数值积分和微分的一些方法。3-1 数值积分法 求定积分的近似值的数值方法。即用被积函数的有限个抽样值的离散或加权平均近似值代替定积分的值。求某函数的定积分时，在多数情况下，被积函数的原函数很难用初等函数表达出来， 因此能够借助微积分学的牛顿-莱布尼兹公式计算定积分的机会是不多的。另外，许多实际问题中的被积函数往往是列表函数或其他形式的非连续函数，对这类函数的定积分，也不能用不定积分方法求解。由于以上原因，数值积分的理论与方法一直是计算数学研究的基本课题。对微积分学作出杰出贡献的数学大师，如.牛顿-.欧拉、高斯等人也在数值积分这个领域作出了各自的贡献，并奠定了它的理论基础。 1 欧拉法（折线法） 假设有一阶微分方程为： 选择初时时刻，在此位置做近似做切线有： ，可以解出 ，写成递推形式： 其中是离散点的取值 2 矩形法 假设一阶微分方程： 写成积分的形式： 这里是初值，由于不能得到的值,则积分项是未知的，为了得到近似值，可以将积分项用矩形计算公式来近似计算，即：将其写出递推公式：或： 对比欧拉折线方法可知，两种方法是一样的，3 梯形法为了得到更精确的计算结果，采用梯形方法，将积分表达式写成如下形式 则 问题的提出，在计算时的表达式里要用到的值，实际上，上面的式子的积分项并不能计算结果，通常采用迭代运算方法： 当时，这个式子里面的为初值，实际上，这个表达式是前面所讲的矩形法公式 ，有了这个项，我们就可以计算下一步。 。 在上面这一组式子中，其中的第二项可以认为是初值，第三项是要用到前一次的迭代运算。如果迭代序列是收敛的，则极限肯定存在，当趋向无穷大时，可以得到这个精确结果。 怎样给出迭代的精度呢？可以采用前一次的计算值和后一次的迭代值的差来控制。 即： 为小量 在实际计算中，我们认为，迭代一次或者两次就认为满足精度。所以，可以写出下面的递推公式。4 差分公式基于数值积分的欧拉方法，我们可以得到近似的微分计算方法，用函数的离散数据点来近似表达在节点处的微分，叫做数值微分。 （1） 一阶系统的数值微分公式 设： 应用离散数据点的近似微分，可以表示为： 令： 称为采样时间写成递推公式则有：（差分形式） ： 除此之外还有一阶中心差分形式： （2）二阶系统的数值微分公式 根据一阶差分公式，我们容易得到二阶差分如下： ，根据这个原理，还可以得到高阶数值微分公式。在SIMULINK中，我们经常采用一种叫做单位延迟的离散模块，这样就可以方便的建立基于数值积分和微分的仿真模型了。 例如，设： 分别对的近似积分和近似微分为：近似积分计算公式为 ： 则近似差微商的计算公式为： 在离散模块库中的unit delay (单位延迟模块)，在该模块中的采样周期（sample time）和递推公式中的采样周期相同，在下图中使用了采样周期为0.01秒的正弦波的积分和导数离散仿真模型.近似积分结果 近似差商结果3-2 龙格-库塔方法 龙格-库塔方法由于它具有较高的精度而在计算机仿真中就占有重要的地位，在此做简要介绍1 二阶龙格-库塔方法 将上述梯形法做进一步推广，可以得到龙格-库塔方法，对于： ： 假设我们从跨出一步，在时的解为 ,在附近做泰劳展开，只保留二次项，则有其中： （3-38）我们假设这个解可以写成： （3-39）即有： （3-40）将（3-40）与（3-38）对比可知道， 四个未知数，有三个方程，可以自由的选择一个值令： 于是有： 于是有一组方程为： 由于计算时，只取了和项，所以，上面的公式称为二阶龙格-库塔方法。2 四阶龙格-库塔方法 二阶龙格-库塔法仍然存在精度不高的缺点，通常采用的是四阶龙格-库塔方法。下面直接给出公式：3关于数值积分的精度在实际应用中，为了提高仿真精度，四阶龙格库塔法的计算步距是可以变化的，也就是在计算过程中可根据所要求的精度对步长作适当改变，这对于精确计算有很大的好处。因为一般的系统在阶跃函数作用下的过渡过程一开始变化比较快，而最后变化将趋于平缓（详见第八章），为了达到较高的精度必须按起始段来选择步距，显然这个步距将很小，这样在平缓阶段，就显得十分浪费。因此，在作精确仿真计算时，我们可以采用变步长的办法，具体步骤如下： (1)按系统的过渡过程情况将它分成几段，每一段都预先先给一个步长。 (2)当过渡过程进入该段时，先要做步长的计算。即用作一次计算，然后再用作一次计算，求得二者之差的绝对值，如果它小于某一预定的值，则认为这个步长符合精度的要求，即可用此步长继续算下去，如果二者之差的绝对值大于该予定值，则认为这个步长过大，那么用作一次计算，再用作二次计算，再进行上述比较。当然，这样一来，计算的工作量将会增加。 对于工程设计来讲，计算精度往往并不要求十分高，比如误差不超过0.5也就令人满意了。为此，一般常采用定步长的计算方法，作为一个经验数据，当用四阶龙格库塔法进行计算时，可选为(，。 其中：为系统在阶跃函数作用下的上升时间 为系统在阶跃函数作用下的过渡过程时间最后还要特别指出一点：当采用龙格库塔法对连续系统进行数字仿真时，为了保证计算的稳定性，必须要求所选择的步长应小于系统中最小时间常数之两倍，可保证计算是稳定的。另外，步长的选择会影响系统的稳定性问题，各种方法的稳定性与步距的大小有关，详细情况可以参见有关书籍。33 面向微分方程的仿真程序设计 根据上一节可知，如果能够将一个连续系统表示成一组一阶微分方程，且每个一阶方程都可以用数值积分的方法来求解，则便可以求出整个系统的全部解，包括各阶导数的值都能给出。这种方法就是面向微分方程的仿真程序设计龙格-库塔法： 在数值仿真中，通常采用四阶的龙格-库塔方法。龙格-库塔方法的基本思想是在区间内多预报几个斜率的值，然后将它们加权平均。其表达式如下： 面向微分方程的数字仿真程序结构，由于面向微分方程的数字仿真模型还有状态空间方法和传递函数方法等情况，因此面向微分方程的仿真程序结构也有两种。在此仅介绍关于微分方程（组）与状态空间方法的数字仿真程序设计。 下面我们通过几个例题来说明，怎样建立微分方程的仿真程序。1， 用四阶龙格-库塔方法求解一阶微分方程程序设计在此应用了VB编程方法，对于熟悉其他计算机语言的学者，可以参考其中程序设计的基本思想。计算微分方程： 在 ，取步长为 的数值解。从到各结点上的数值解，请读者将本结果和欧拉方法以及梯形法相比较。注：这个方程有精确解为：新建立一窗体，在窗体上加一个列表框，一文本框数组（0-3），分别存放开始时间、终了时间、采样时间间隔和初始值在窗体的通用事件过程中添加函数过程（给出微分方程模型）在窗体上添加一个命令按钮 Commmand1，将原代码写在命令按钮的单击事件过程中。调试程序，观察运行的结果的正确性Private Function fn1(x, y)fn1 = y - 2 * x / yEnd FunctionPrivate Sub Command1_Click()Dim t0 As LongDim tn As LongDim h As SingleDim y0 As LongDim i As IntegerList1.ClearDim x, y, n As Integer, X1, Y1, k1, k2, k3, k4 As LongList1.AddItem (-龙格库塔方法-)t0 = Val(Text1(0).Text) 数组控件tn = Val(Text1(1).Text) Val转换函数h = Val(Text1(2).Text) 循环步长 - If h = 0 Then MsgBox (请输入步长h) Exit Sub End If -给定初值y0 = Val(Text1(3).Text)n = (tn - t0) / h 循环次数x = Val(t0)y = y0=List1.AddItem (时间T & & Y值)List1.AddItem (-)List1.AddItem (t0 & & y0)For i = 1 To nk1 = fn1(x, y)X1 = x + h / 2Y1 = y + h * k1 / 2k2 = fn1(X1, Y1)Y1 = y + h * k2 / 2k3 = fn1(X1, Y1)Y1 = y + h * k3X1 = x + hk4 = fn1(X1, Y1)x = x + hy = y + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6List1.AddItem (Format(x, 0.0000) & ; & Format(y, 0.0000)Next iEnd SubPrivate Sub Form Load()Text1 (0).Text = 0Text1 (1).Text = 3Text1 (2).Text = 0.1Text1 (3).Text = 1End Sub为了和仿真结果进行比较，下面建立了仿真模型，在该模型中，将积分器的初始条件设置为1，仿真时间设置为0-3秒，计算结果如下。根据VB的数值计算和Sinulink仿真结果相比较，有较高的精度。2 一阶微分方程组的四阶的龙格库塔法程序设计设一阶微分方程组形式如下： 下面我们仅考虑有两个方程的情况，设初值为： 根据上面推导则可以得到四阶龙格库塔法的格式为：其中： 不难将上述情况推广到多维一阶方程组的情况，例如：设有个一微分方程组,下面给出了计算的参考子程序程序代码。初值为： 首先按微分方程的个函数定义四个子函数分别来计算,(在这里使用了,来代替了，),VB程序代码如下：-Private Function fxn1(x1, x2, x3, x4, T) As Singlefxn1 = x3End FunctionPrivate Function fxn2(x1, x2, x3, x4, T) As Singlefxn2 = x4End FunctionPrivate Function fxn3(x1, x2, x3, x4, T) As Singlefxn3 = -1 / x1 / x1 + x1 * x4End FunctionPrivate Function fxn4(x1, x2, x3, x4, T) As Singlefxn4 = -2 * x3 * x4 / x1End Function-Private Sub Command1_Click()Dim x1 As Single, x2 As Single, x3 As Single, x4 As SingleDim T0 As Single, TN As Single, H As Single, T As SingleDim N As Integer, i As IntegerDim K1, K2, K3, K4 As SingleDim S1, S2, S3, S4 As SingleDim T1, T2, T3, T4 As SingleDim U1, U2, U3, U4 As Single-T0 = 0: TN = 3: x1 = 1: x2 = 0.5: x3 = 1: x4 = 1.5H = 0.01N = (TN - T0) / HT = T0-For i = 1 To NK1 = 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 + 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 / 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 + 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) / 6x2 = x2 + (S1 + 2 * S2 + 2 * S3 + S4) / 6x3 = x3 + (T1 + 2 * T2 + 2 * T3 + T4) / 6x4 = x4 + (U1 + 2 * U2 + 2 * U3 + U4) / 6List1.AddItem (Format(T, 0.0000) & & Format(x1, 0.0000) & & Format(x2, 0.0000) & & Format(x3, 0.0000) & & Format(x4, 0.0000)T = T + HNext iEnd Sub3 高阶微分方程的四阶龙格库塔法程序设计 对于高阶微分方程，我们总可以通过引进状态变量(见第五章状态空间)，将高阶微分方程转化为微分方程组来处理。下面是模拟火箭飞行中的一个例子，这个例子的数学模型是一个二阶微分方程，通过这个例子来说明高阶微分方程的仿真程序设计。设有一有阻尼的弹簧质量系统，动力学方程为：假定 ，解：(1) 先将原方程化为一阶方程组：即令： 则一阶微分方程组为： 假定初值为 这是一个由两个方程组成的一阶微分方程组，不难根据前一个例子进行程序设计，在四阶龙格库塔法程序设计中，只需要使用变量，即可，值得注意的是求得解代表系统的速度，代表系统的位移。(2) 程序设计 在VB环境下，建立一个窗体，在窗体上添加一个列表框、一个图片框和一个命令按钮，源程序如下：Option Explicitfunction过程-定义函数 function fxn1 和function fxn2 Private Function fnx1(x1, x2, t)fnx1 = x2End FunctionPrivate Function fnx2(x1, x2, t)fnx2 = (-20 * x1 - 4 * x2 + 0.2* Sin(2 * t)End Function-Private Sub Command1_Click()Dim x1 As Single, x2 As Single, t As Single, H As Single, tn As SingleDim N As Integer, i, j As IntegerDim k1, k2, k3, k4 As SingleDim s1, s2,s3, s4 As SingleDim x0, y0, t0 As Single初始化Me.Picture1.Scale (0, 1.5)-(10, -1.5) 限制图片框的范围 (左上角)-(右下角)For j = -1 To Abs(Me.Picture1.ScaleHeight)Me.Picture1.Line (0, j)-(10, j), RGB(255, 0, 0)Next jFor i = 1 To Me.Picture1.ScaleWidth Step 2Me.Picture1.Line (i, -1.5)-(i, 1.5), RGB(0, 0, 255)Next i设置仿真参数（初始值，开始时间、结束时间、（自动并计算循环步长）x1 = 0: x2 = 0: t = 0H = 0.001: tn = 10: N = (tn - t) / HList1.AddItem ( 时间 & & 位移 & & 速度)List1.AddItem (-) 四阶龙格库塔法-For i = 1 To Nk1 = fnx1(x1, x2, t)s1 = fnx2(x1, x2, t)t = t + H / 2k2 = fnx1(x1 + H * k1 / 2, x2 + H * s1 / 2, t)s2 = fnx2(x1 + H * k1 / 2, x2 + H * s1 / 2, t)k3 = fnx1(x1 + H * k2 / 2, x2 + H * s2 / 2, t)s3 = fnx2(x1 + H * k2 / 2, x2 + H * s2 / 2, t)t = t + H / 2k4 = fnx1(x1 + H * k3, x2 + H * s3, t)s4 = fnx2(x1 + H * k3, x2 + H * s3, t)x1 = x1 + H * (k1 + 2 * k2 + 2 * k3 + k4) / 6x2 = x2 + H * (s1 + 2 * s2 + 2 * s3 + s4) / 6List1.AddItem (Format(t, 00.0000) & & x1 & & x2)-绘图Me.Picture1.Line (t0, x0)-(t, x1)Me.Picture1.Line (t0, y0)-(t, x2)t0 = tx0 = x1y0 = x2Next iEnd Sub=将以上程序拷贝到窗体的通用事件过程中，单击命令按钮，则可以得到一下结果。请读者自己设计一个基于微分方程的SIMULINK仿真模型，并和上面的结果进行比较。对于其他微分方程组，可以先设定函数子程序，然后利用本方法求解。对于高阶微分方程，可以先简化成微分方程组再求解。 34 关于Simulink环境中的求解器Solver在利用Simulink进行仿真时，其中的积分模块的核心是利用了各个不同形式的四阶的龙格库塔法，了解Simulink环境中的求解器Solver，有助于正确建立仿真模型。Simulink所提供的求解器都是当今国际上数值计算研究的最新成果，是一种速度最快，精度最高的计算方法。能够非常理想地求解各类微分方程。不同的系统需要利用不同的求解器，故了解系统的特性是非常重要的，如系统方程是否是刚性方程(stiff equation)等。下面来具体介绍各种计算方法，这对于不同的系统选择不同的方法是至关重要的。对于用数值方法求解常系数微分方程（Ordinary Differential Equation，简写为ODE）或微分方程组，Simulink提供了7种求解函数，他们是： 1 Ode45 这种求解器采用Runge-Kutta方法，这也是利用Simulink求解微分方程时最常用的一种方法。这种算法精度适中，是计算方程中的首选项。 它是利用有限项的Taylor级数取近似解函数，而误差的来源就是Taylor的截断项，误差就是截断误差。 Ode45分别采用4阶、5阶Taylor级数计算每个积分步长终端的状态变量近似值，并利用这个级数的值相减，得到的误差作为计算误差的判断标准。如果误差估计值大于这个系统的设定值，那么就把该积分步长缩短，然后重新计算；如果误差远小于系统的设定值，那么就将积分步长放长。 2 Ode23 这种求解器采用Runge-Kutta方法，为了能够达到ode45同样的精度，ode23的积分步长总要比ode45取的小。因此，ode23处理“中度Stiff”问题的能力优于ode45。ode23是利用有限项的Taylor级数取近似解函数，而误差的来源就是Taylor的截断项，其中，误差就是指截断误差。 ode45分别采用2阶、3阶Taylor级数计算每个积分步长终端的状态变量近似值，并利用这个级数的值相减，得到的误差作和计算误差的判断标准。如果误差估计值大于这个系统的设定值，那么就把该积分步长缩短，然后重新计算。如果误差远小于系统的设定值，那么就将积分步长扩大。ode23和ode45邵是变步K算法。 3 Ode113 Ode113与ode45和ode23不同，它采用的变阶Adams法，是一种多步预报校正算法。 odell3在执行过程中还自动地调整近似多项式的阶数，以平衡其精确性和有效性。 ode45和ode23采用的是Taylor级数方法，而odell3采用的是多项式方法，计算导数的次数也比前面两种方法次数少，所以在计算光滑系统时，odell3的速度更快。 4odel5s是专门刷来求解刚性(Stiff)方程的变阶多步算法，包含一种对系统动态转换进行检测的机理。这种检测使这一算法对非刚性(Stiff)系统计算效率低下，尤其是对那种有快速变化模式的系统情况。 5 ode23s ode23s问odel5s一样都是用来求解刚性方程的，是基于Rosenbrok公式建立起来的定阶单步算法。由于计算阶数不变，所以计算效率要比odel5s效率高。 6 ode23t 用来求解中度刚性方程。7 ode23tb 用来求解中度刚性方程。值得注意的是当系统是刚性方程时，要选择适当的求解器才能得到正确结果， 例如 对刚性微分方程，利用不同的求解器仿真结果和精确值进行比较：并说明求解器使用不当，产生的结果差别非常大。例如我们需要求微分方程 。在初始条件，的解 解：首先求解方程的解析解，然后通过不同求解器所得到的数值方法进行比较，来说明不同算法之间的不同，以及不合适方法所带来的错误结果。【所谓的刚性微分方程的含义是指：在用微分方程描述的一个变化过程中，若往往又包含着多个相互作用但变化速度相差十分悬殊的子过程，这样一类过程就认为具有 “刚性”。描述这类过程的微分方程初值问题称为“刚性问题”。例如，宇航飞行器自动控制系统一般包含两个相互作用但相对速度相差十分悬殊的子系统，一个是控制飞行器质心运动的系统，当飞行器速度较大时，质心运动惯性较大，因而相对来说变化缓慢；另一个是控制飞行器运动姿态的系统，由于惯性小，相对来说变化很快，因而整个系统就是一个刚性系统。】35 关于用MATLAB中命令求解微分方程的解析解在MATLAB中的命令空间中直接可以求解微分方程的函数(dsolve)，其格式为： fs3=dsolve(D2y+999*Dy+0.999*y,y(0)=1,Dy(0)=0,x)% 其中第一项为微分方程的形式，第二第三项是初始条件，最后一项是自变量，回车后显示如下的信息 fs3 = 1/554442780*(-16650+277221390(1/2)*277221390(1/2)*exp(-3/100*(16650+277221390(1/2)*x)+(15/499498*277221390(1/2)+1/2)*exp(3/100*(-16650+277221390(1/2)*x)为了得到速度表达式，我们可以通过微分命令(diff)： df3=diff(fs3,t) df3 = 1/554442780*(16650+277221390(1/2)*277221390(1/2)*(-999/2+3/100*277221390(1/2)*exp(3/100*(-16650+277221390(1/2)*t)+(-15/499498*277221390(1/2)+1/2)*(-999/2-3/100*277221390(1/2)*exp(-3/100*(16650+277221390(1/2)*t) （2）新建立SINMULINK仿真模型，和脚本文件存放在同体格目录中，取名为“modefs3”建立脚本文件，来分析不同的求解器的仿真结果， 仿真模型脚本文件代码如下：fs3=dsolve(D2y+999*Dy+0.999*y,y(0)=1,Dy(0)=0,t); df3=diff(fs3,t) % 对变量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),b); 绘图，其中的因变量取模型的第1个输出序列hold onopts=simset(solver,ode15s); %重新设定ode15s为当前求解器t2,y2,s=sim(modefs3,500,opts);plot(t2,y2(:,1),r-);hold onaxis(249.7 250 0.7785,0.7795) %显示局部区域t=0:0.1:500;y=eval(char(df3);figure(2)plot(t,y,k);hold onopts=simset(solver,ode45);%t1