2.机电系统CAD算法基础 2.1.数值积分算法 2.1.1基本数值积分方法 连续系统的动态特性一般可由一个微分方程或一组一阶微分方程加以描述。因 此对系统进行计算机仿真,就需要研究如何利用计算机对微分方程进行求解,目前 采用的方法是应用数值积分法对微分方程求数值解,其中欧拉法、梯形法、龙格 库塔法等数值积分法是几种基本的方法 1.欧拉( Euler)法 欧拉法是一种最简单的数值积分方法,但其精度差,实际应用很少。然而由于 其算法简单且具有明显的几何意义,有利于初学者学习应用数值积分法求解微分方 程,因此在讨论数值积分法时通常先介绍欧拉法。 假设一阶微分方程为 dy dt 初始条件为 y(to)=yo 欧拉法的基本思想就是将(2-1)式中的积分曲线解用直线段所组成的折线加以 近似。根据导数的定义可知,存在 dt △ 故 hny(+△)-y() (2-2) A→0 fIt,y(] 成立,当M足够小时,上式可近似表示为 y(t+△r)-y() ≈f,y( (23) △t 若令M=h,t=lo则得到
24 2．机电系统 CAD 算法基础 2．1．数值积分算法 2．1．1 基本数值积分方法 连续系统的动态特性一般可由一个微分方程或一组一阶微分方程加以描述。因 此对系统进行计算机仿真，就需要研究如何利用计算机对微分方程进行求解，目前 采用的方法是应用数值积分法对微分方程求数值解，其中欧拉法、梯形法、龙格一 库塔法等数值积分法是几种基本的方法。 1．欧拉（Euler）法 欧拉法是一种最简单的数值积分方法，但其精度差，实际应用很少。然而由于 其算法简单且具有明显的几何意义，有利于初学者学习应用数值积分法求解微分方 程，因此在讨论数值积分法时通常先介绍欧拉法。 假设一阶微分方程为 f (t y) dt dy = , （2-1） 初始条件为 ( ) 0 0 y t = y 欧拉法的基本思想就是将（2-1）式中的积分曲线解用直线段所组成的折线加以 近似。根据导数的定义可知，存在 ( ) ( ) t y t t y t dt dy t + − = →0 lim 故 ( ) ( ) t y t t y t t + − →0 lim = f[t, y(t)] （2-2） 成立，当 t 足够小时，上式可近似表示为 ( ) ( ) f [t, y(t)] t y t t y t + − （2-3） 若令 t = h ， 0 t = t 则得到
y(n+h)≈y(a)+h/[n,yn 其中h称为计算步长。 因为y(tn)=y,tn+h=t1,所以有y(t1)≈yo+h/(a 由此可得到所欲求解在t=1处的近似值y1为 y1=y()+h/[n,y( (2-5) 这样就可由已知点(n2y)去确定它在1y平面上的另一点(t1y1)若设P0点是初始 给定点,且微分方程的精确解y(0经过P点,对应1=1+h的y()的精确解为 y(t1),即P点,而由(2-5)式求得的y1值所对应的点为P1点,其坐标为(1,y1) 它是从P点作斜率为f(n,y)的线段与直线【=1的交点。 若从P{(12y1)点出发,可求解在t=12=n+2h上的近似值y2为 y2=y1 (2-6) 即从P1几点作斜率等于f(1,y1)的线段与直线t=t2相交而得到P2点。根据类似的 思想由式(2-6)不难推出 Vk+I=y Kr) P〔! 32 3(r;) 图2-1欧拉数值积分方法
25 ( ) ( ) ( ) 0 0 0 0 y t + h y t + hf t , y t （2-4） 其中 h 称为计算步长。 因为 ( ) 0 0 y t = y ， 0 1 t + h = t ，所以有 ( ) ( ) 1 0 0 0 y t y + hf t , y 由此可得到所欲求解在 1 t = t 处的近似值 1 y 为 ( ) ( ) 1 0 0 0 y = y t + hf t , y t （2-5） 这样就可由已知点 ( ) 0 0 t , y 去确定它在 t-y 平面上的另一点 ( ) 1 1 t , y 。若设 P0 点是初始 给定点，且微分方程的精确解 y(t) 经过 P0 点，对应 t 1 = t 0 + h 的 y(t) 的精确解为 ( ) 1 y t ，即 P 点，而由（2-5）式求得的 1 y 值所对应的点为 P1 点，其坐标为 ( ) 1 1 t , y ， 它是从 P0 点作斜率为 ( ) 0 0 f t , y 的线段与直线 1 t = t 的交点。 若从 ( ) 1 1 1 P t , y 点出发，可求解在 t = t 2 = t 0 + 2h 上的近似值 2 y 为 ( ) 2 1 1 1 y = y + hf t , y （2-6） 即从 P1 几点作斜率等于 ( ) 1 1 f t , y 的线段与直线 2 t = t 相交而得到 P2 点。根据类似的 思想由式（2-6）不难推出 ( ) k k k k y = y + hf t , y +1 （2-7） 图 2-1 欧拉数值积分方法
利用上式不断对微分方程求近似解,直至t=1点P将P,P1,P,…,P诸点连 成一条折线,即得到所求微分方程精确解y()的近似曲线解,如图21所示。(27) 式称为欧拉法递推公式 该方法简单,计算量小,由前一点即可推出后一点值,属于单步法。因其从初 值即可开始进行递推运算,无需其他条件要求,因此能够自启动。该算法精度低 适当减小计算步长h有助于提高计算精度。 2.梯形法 1).欧拉法误差原因的分析 根据欧拉法的分析可知,该方法简单但较粗 ,精度差。考虑对于微分方程 山 存在精确解 ()=y+f(,y 当t=t1时 y)=y+/(y (2-8) 由前述推导可知 y(t1)≈yn+(1-0)/(n,y) (2-9) 对照(2-8)式与(2-9)式得到,在欧拉法中存在 广(y知≈(nx-4) 用∫(ta,yn)代替变量,由矩形近似曲边梯形,因此造成欧拉法误差较大
26 利用上式不断对微分方程求近似解，直至 N t = t 点 PN 。将 P P P PN , , , , 0 1 2 诸点连 成一条折线，即得到所求微分方程精确解 y(t) 的近似曲线解，如图 2-1 所示。（2-7） 式称为欧拉 法递推公式。 该方法简单，计算量小，由前一点即可推出后一点值，属于单步法。因其从初 值即可开始进行递推运算，无需其他条件要求，因此能够自启动。该算法精度低， 适当减小计算步长 h 有助于提高计算精度。 2．梯形法 1）．欧拉法误差原因的分析 根据欧拉法的分析可知，该方法简单但较粗 糙，精度差。考虑对于微分方程 f (t y) dt dy = , ( ) 0 0 y t = y 存在精确解 ( ) ( ) = + t t y t y f t y dt 0 0 , 当 1 t = t 时 ( ) ( ) = + 1 0 1 0 t t y t y f t, y dt （2-8） 由前述推导可知 ( ) ( ) ( ) 1 0 1 0 0 0 y t y + t −t f t , y （2-9） 对照（2-8）式与（2-9）式得到，在欧拉法中存在 ( ) ( )( ) 0 0 1 0 1 0 f t y dt f t y t t t t − , , 用 ( ) 0 0 f t , y 代替变量，由矩形近似曲边梯形，因此造成欧拉法误差较大
2).改进的欧拉法一一梯形法 由前述分析可以看出,进一步改进欧拉法的方法就是利用梯形面积代替矩形面 积,以逼近曲边梯形面积,即 广(,y)≈ h 2 推而广之,得到梯形法公式为 Yl=y +l(,y)+frui (2-10) 由上式可以看出,在计算yx+1时,在等式右边存在有yxn1,因此无法求出。在 这种情况下,一般采用欧拉法进行一次选代,求出yx1的预估值,再将其代人梯形 公式计算出yx+1真值,从而构成预估校正梯形公式为 yo=y +hf(,y) (2-11) 离+1 (2-12) (2-11)式称为预报式,(2-12)式称为校正式。 3.龙格一库塔( Runge-kutt)法 考虑如下一阶微分方程 dy dt f(,y) (a)=y 设1=0+h,这里h为计算步长,在1时刻的解为y1=y{(n)+h,可在1附 近展成泰勒级数,取展开式的前三项得到 VI=yo+ h (2-13) hf(o, yo) h2( 2 at
27 2）．改进的欧拉法——梯形法 由前述分析可以看出，进一步改进欧拉法的方法就是利用梯形面积代替矩形面 积，以逼近曲边梯形面积，即 ( ) ( ) 0 1 2 1 0 f f h f t y dt t t + , 推而广之，得到梯形法公式为 ( ) ( ) 1 1 1 2 + = + + + + f t y f t y h y y , , （2-10） 由上式可以看出，在计算 +1 y 时，在等式右边存在有 +1 y ，因此无法求出。在 这种情况下，一般采用欧拉法进行一次选代，求出 +1 y 的预估值，再将其代人梯形 公式计算出 +1 y 真值，从而构成预估校正梯形公式为 ( ) y = y + hf t , y + 0 1 （2-11） ( ) ( ) 0 1 1 1 2 + = + + + + f t y f t y h y y , , （2-12） （2-11）式称为预报式，（2-12）式称为校正式。 3．龙格一库塔（Runge- kutta）法 考虑如下一阶微分方程 f (t y) dt dy = , ( ) 0 0 y t = y 设 t 1 = t 0 + h ，这里 h 为计算步长，在 1 t 时刻的解为 y1 = y(t 0 )+ h ，可在 0 t 附 近展成泰勒级数，取展开式的前三项得到 ( ) 0 0 0 0 2 2 1 2 0 0 0 2 2 2 1 0 y y t t t t t t t f f t h f y hf t y dt d y h dt dy y y h = = = = + = + + = + + , （2-13）
假设上式的解可表示为 y1=y+(a1K1+a2K2 K,=flo, yo r(,+b, h,y,+b, k,h) (2-16) 将K2在(a,y)附近展开成泰勒级数,并取前三项得到 K2=f(,y)+42+bK,9 将(2-15)式(2-17)式代入(2-14)式中,经整理得 y,=yo +a,hf(o, yo)+arh f(o, yo )+b, h2+b,hi 将(218)式与(213)式相比较得到{2a2b1=1 b2=1 若取h==1,则得到a1=a2=2 将其代入(2-14)式、(2-15)式、(2-16) 式,得 +K K y 将(2-19)式写成递推形式 y +-(K, +K K1=f(t4,y) (2-20) 2=f(t4+h,y4+K1b) 由于(220)式只取了级数展开式的前三项,略去了h2以上的高阶无穷小,其
28 假设上式的解可表示为 y1 = y0 + (a1K1 + a2K2 )h (2-14) ( ) 1 0 0 K = f t , y (2-15) K f (t b h y b K h) 2 = 0 + 1 0 + 2 1 , (2-16) 将 K2 在 ( ) 0 0 t , y 附近展开成泰勒级数，并取前三项得到 ( ) y f b K h t f K f t y b h + 2 0 0 + 1 2 1 , (2-17) 将(2-15)式( 2-17)式代入(2-14)式中，经整理得 ( ) ( ) + = + + + y f b hK t f y1 y0 a1hf t 0 y0 a2 h f t 0 y0 b1h 2 1 , , (2-18) 将（2-18）式与（2-13）式相比较得到 = = + = 2 1 2 1 1 2 2 2 1 1 2 a b a b a a 若取 b1 = b2 = 1 ，则得到 2 1 a1 = a2 = 。将其代入（2-14）式、（2-15）式、（2-16） 式，得 ( ) ( ) ( ) = + + = = + + K f t h y K h K f t y K K h y y 2 0 0 1 1 0 0 1 0 1 2 2 , , （2-19） 将（2-19）式写成递推形式 ( ) ( ) ( ) = + + = + = + + K f t h y K h K f t y K K h y y k k k k k k 2 1 1 1 1 2 2 , , （2-20） 由于（2-20）式只取了级数展开式的前三项，略去了 2 h 以上的高阶无穷小，其
截断误差正比于h2,故这种方法称为二阶龙格库塔法。二阶龙格一库塔法的计算精 度是不高的,为此在实际应用中往往采用精度更高的四阶龙格一库塔法,它的截断 误差正比于h2,由于其推导过程与二阶龙格一库塔法类似,在此就不再推导,而直 接给出四阶龙格一库塔法的递推表达式,即 yk=yk+2(K1+2K2+2K3+K4) K,=f(k,k) h2h +一K K4=f(t4+h,y4+h3) 考察(2-20)式可以清楚地看出,它与前面导出的梯形算法是相同的,因而可以 知道梯形法的精度也是二阶的。当对级数展开式只取前两项时可以得到 y1=y+hf(x,y),这就是前面得到的欧拉法递推式。通过对龙格一库塔法的 分析与推导,可以将欧拉法、梯形法与龙格、库塔法这三种数值积分方法统一起来, 它们都可归于龙格一库塔法,只不过是阶次不同而且。从理论上讲,可以构造任意 阶的计算方法,但截断误差的阶数与所计算函数数值∫的次数并非是等值关系。对 于四阶以上龙格一库塔法,所需计算的∫次数将高于误差阶数,这将大大增加了计 算工作量,因而选择这种高阶龙格-库塔法是不适当的。实际上,对于工程应用四阶 龙格-库塔法已完全能够满足仿真精度要求,所以四阶龙格-库塔法在实际应用中具有 重要地位。 4.四阶龙格一库塔法的向量表示 前面推导过程中假设系统为一阶系统,然而实际上系统往往是高阶的,因而描 述这类系统应是一组一阶微分方程或是状态方程,所以介绍求解状态方程的四阶龙 格一库塔向量式是完全必要的 考虑一n阶系统,其状态方程的一般描述为 Y=F(, Y (n)=Y
29 截断误差正比于 2 h ，故这种方法称为二阶龙格库塔法。二阶龙格一库塔法的计算精 度是不高的，为此在实际应用中往往采用精度更高的四阶龙格一库塔法，它的截断 误差正比于 2 h ，由于其推导过程与二阶龙格一库塔法类似，在此就不再推导，而直 接给出四阶龙格一库塔法的递推表达式，即 ( ) ( ) ( ) = + + = + + = + + = + = + + + + 4 3 3 2 2 1 1 1 1 2 3 4 2 2 2 2 2 2 6 K f t h y hK K h y h K f t K h y h K f t K f t y K K K K h y y k k k k k k k k k k , , , , （2-21） 考察（2-20）式可以清楚地看出,它与前面导出的梯形算法是相同的，因而可以 知道梯 形法 的精 度也 是二 阶的。 当对 级数 展开 式只 取前 两项 时可以 得到 ( ) k k k k y = y + hf t , y +1 ，这就是前面得到的欧拉法递推式。通过对龙格一库塔法的 分析与推导，可以将欧拉法、梯形法与龙格、库塔法这三种数值积分方法统一起来， 它们都可归于龙格一库塔法，只不过是阶次不同而且。从理论上讲，可以构造任意 阶的计算方法，但截断误差的阶数与所计算函数数值 f 的次数并非是等值关系。对 于四阶以上龙格一库塔法，所需计算的 f 次数将高于误差阶数，这将大大增加了计 算工作量，因而选择这种高阶龙格-库塔法是不适当的。实际上，对于工程应用四阶 龙格-库塔法已完全能够满足仿真精度要求，所以四阶龙格-库塔法在实际应用中具有 重要地位。 4．四阶龙格一库塔法的向量表示 前面推导过程中假设系统为一阶系统，然而实际上系统往往是高阶的，因而描 述这类系统应是一组一阶微分方程或是状态方程，所以介绍求解状态方程的四阶龙 格一库塔向量式是完全必要的。 考虑一 n 阶系统，其状态方程的一般描述为 Y = F(t,Y ) ( ) 0 Y0 Y t =
四阶龙格-库塔法的向量表示为 F41=Y4+(k1+2k2+2k3+k k1=F(n,F4) h h-212 其中k=yk2…k人k(=12.…,n;j=1234)是微分方程组中第个方 程的第j个RK系数。有时考虑到计算方便,又可将(2-22)式展为下述形式,即 yAx=yx+2(k1+2k2+2k3+k) k,1 =f(r, yix, y2 ,2 +k1,y2+k1…,ym+km 2,y2 h 4 21x+k13,y2x2 其中i=1,2,…,n:n为系统阶数:下标K为递推下标。 5.亚当斯( Adams)法 前面介绍数值积分算法的特点是,在计算K+1时刻的值yx时,只用到第K 时刻yx和厂x的值。但实际上在逐步递推过程中,计算yx+1:之前已经获得了一系 列的近似值y,y1,…,yx以及f,f1,…,f,如果能够充分利用历史上的
30 四阶龙格-库塔法的向量表示为 ( ) ( ) ( ) = + + = + + = + + = + = + + + + 4 3 3 2 2 1 1 1 1 2 3 4 2 2 2 2 2 2 6 k F t h Y hk k h Y h k F t k h Y h k F t k F t Y k k k k h Y Y m m m m m m m k k k , , , , 其中 ( ) ( 1 2 1 2 3 4) 1 2 k = , , , , i = , , ,n j = , , , j j j nj i j ； 是微分方程组中第 i 个方 程的第 j 个 RK 系数。有时考虑到计算方便，又可将（2-22）式展为下述形式，即 ( ) 1 1 2 2 2 3 4 6 i,k i,k i, i, i, i, k k k k h y + = y + + + + ( ) , , , , , , , , i i n k f t y y y 1 = 1 2 2 = + 1 + 1 1 2 + 2 1 + 1 2 2 2 2 , , , , , , , , , , , i i n n k h k y h k y h y h k f t 3 = + 1 + 1 2 2 + 2 2 + 2 2 2 2 2 , , , , , , , , , , , i n n k h k y h k y h y h k f t 4 = + 1 + 1 3 2 + 2 3 + 3 2 2 2 2 , , , , , , , , , , , i n n k h k y h k y h y h k f t 其中 i=1，2，…，n；n 为系统阶数；下标 为递推下标。 5．亚当斯（Adams）法 前面介绍数值积分算法的特点是，在计算 + 1 时刻的值 +1 y 时，只用到第 时刻 y 和 f 的值。但实际上在逐步递推过程中，计算 +1 y ；之前已经获得了一系 列的近似值 0 y ， 1 y ，…， y 以及 0 f ， 1 f ，…， f ，如果能够充分利用历史上的
一些数据来求解ν’则有望既可加快仿真速度又能获得较高的仿真精度,这就是 构造多步法的基本思路。 Adams法就是利用已经得到的这些信息来实现对yx+1的计 算,由于它充分利用了前面得到的信息,所以在加快仿真速度的同时提高了仿真精 度。在各种多步法中, Adams法最具有代表性,应用也最普遍 本书不再进行公式的推导,直接给出 Adams显式表达式 h .(55f-59/1+372-9/m3) (2-23)式即为 Adams显式表达式,其截断误差为 251 R h5y°(5),5∈(x1,x) 因此(2-23)式为四阶 Adams显式公式。 (2-24)式给出了 Adams隐式表达式,即 h (91+19 f2) (2-24) 其截断误差为 e015),.5∈(x,x) 120 隐式 Adams法的优点主要在于精度较高,然而它无法实现自身递推,需要另外 的显式提供一个估计值,(2-25)式给出了由 Adams显式提供的估计值,隐式加以校 正的四阶 Adams预估校正法公式为 h y (55-59/+37/2-9 9 f1+19x-5Jx1+f8 fP=frl,yp. 显然上式除了由 Adams显式所提供的估计值外,还需要其他方法如四阶龙格-库塔法 计算出∫x、∫x1、fx-2等初值
31 一些数据来求解 +1 y ，则有望既可加快仿真速度又能获得较高的仿真精度，这就是 构造多步法的基本思路。Adams 法就是利用已经得到的这些信息来实现对 +1 y 的计 算，由于它充分利用了前面得到的信息，所以在加快仿真速度的同时提高了仿真精 度。在各种多步法中，Adams 法最具有代表性，应用也最普遍。 本书不再进行公式的推导，直接给出 Adams 显式表达式 ( ) 1 55 59 1 37 2 9 3 24 + = + − − + − − − f f f f h y y （2-23） （2-23）式即为 Adams 显式表达式，其截断误差为 ( ) ( ) 5 5 720 251 R = h y ， ( ) t ,t −1 因此（2-23）式为四阶 Adams 显式公式。 （2-24）式给出了 Adams 隐式表达式，即 ( ) 1 9 1 19 5 1 2 24 + = + + + − − + − f f f f h y y （2-24） 其截断误差为 ( ) ( ) 5 5 120 19 R = h y ， ( ) t ,t −1 隐式 Adams 法的优点主要在于精度较高，然而它无法实现自身递推，需要另外 的显式提供一个估计值，（2-25）式给出了由 Adams 显式提供的估计值，隐式加以校 正的四阶 Adams 预估校正法公式为 ( ) ( ) ( ) = = + + − + = + − + − + + + + + − − + − − − p p c p p f f t y f f f f h y y f f f f h y y 1 1 1 1 1 1 2 1 1 2 3 9 19 5 24 55 59 37 9 24 , （2-25） 显然上式除了由 Adams 显式所提供的估计值外，还需要其他方法如四阶龙格-库塔法 计算出 f 、 −1 f 、 −2 f 等初值
2.1.2数值积分方法的计算稳定性 1.数值积分方法计算稳定性概念 数值积分方法的稳定性与机电系统的稳定性是两个不同的概念。一个稳定的系 统,当用某种数值积分方法进行仿真计算时,若该方法或应用该方法所选参数不适 则会产生意想不到的结果一一仿真结论出错,这就是由于数值积分方法的稳定 性引起的。不同的数值积分方法对应于不同的差分方程,一个数值解是否稳定决定 于该差分方程的特征根是否满足稳定的要求。对于不同的数值积分方法,其稳定性 和约束条件都是仿真过程中的重要参数 考察已知系统,其微分方程为 dy -30y t<1 dt 其精确解为 现在取步长h=0.1,用欧拉法和四阶龙格库塔法计算F=1.5时的y(t)值。 欧拉法y(1.5)=-109225×10 四阶龙格库塔法y(15)=395730×10 精确解y(15)=954173×10 显然此时数值计算结果是错误的。这是因为数值积分法是一种近似积分方法, 它在反复递推计算中将引入误差,若误差积累越来越大,将使计算出现不稳定。所 以原系统稳定与用数值积分法的计算稳定性是不同的概念。对干系统稳定性的讨论, 使用微分方程或传递函数;而后者的稳定性问题需要使用相应的差分方程来讨论 由于对高阶微分方程的数值计算稳定性进行全面的分析十分困难,通常用一简单的 阶微分方程来考察其相应的差分方程的计算稳定性问题。将下述微分方程 dt v(),y(0)=y 称为测试方程。根据稳定性理论,其特征方程的根在S平面的左半部,即根的实部Re λ<0时,原方程稳定。现在以这种条件来讨论数值计算方法的自身稳定性问题
32 2．1．2 数值积分方法的计算稳定性 1．数值积分方法计算稳定性概念 数值积分方法的稳定性与机电系统的稳定性是两个不同的概念。一个稳定的系 统，当用某种数值积分方法进行仿真计算时，若该方法或应用该方法所选参数不适 当，则会产生意想不到的结果——仿真结论出错，这就是由于数值积分方法的稳定 性引起的。不同的数值积分方法对应于不同的差分方程，一个数值解是否稳定决定 于该差分方程的特征根是否满足稳定的要求。对于不同的数值积分方法，其稳定性 和约束条件都是仿真过程中的重要参数。 考察已知系统，其微分方程为 y dt dy = −30 (0 t 1.5) ( ) 3 1 y 0 = 其精确解为 ( ) t y t e 30 3 1 − = 现在取步长 h=0.1，用欧拉法和四阶龙格-库塔法计算 t=1.5 时的 y（t）值。 欧拉法 ( ) 4 y 1.5 = −1.0922510 四阶龙格-库塔法 y(1.5) = 3.9573010 精确解 ( ) 21 1 5 9 54173 10− y . = . 显然此时数值计算结果是错误的。这是因为数值积分法是一种近似积分方法， 它在反复递推计算中将引入误差，若误差积累越来越大，将使计算出现不稳定。所 以原系统稳定与用数值积分法的计算稳定性是不同的概念。对干系统稳定性的讨论， 使用微分方程或传递函数；而后者的稳定性问题需要使用相应的差分方程来讨论。 由于对高阶微分方程的数值计算稳定性进行全面的分析十分困难，通常用一简单的 一阶微分方程来考察其相应的差分方程的计算稳定性问题。将下述微分方程 y(t) dt dy = ， ( ) 0 0 y = y （2-26） 称为测试方程。根据稳定性理论，其特征方程的根在 S 平面的左半部，即根的实部 Re λ＜0 时，原方程稳定。现在以这种条件来讨论数值计算方法的自身稳定性问题
2.欧拉法的计算稳定性分析 对测试方程根据欧拉法得到 Yr=yr +hay (227) 对上式进行Z变换得 (z)=(1+h 显然其特征方程为 (1+h)=0 根据控制理论可知,(2-27)式的稳定域在Z平面上是单位圆。(2-27)式稳定的条 件是它的特征方程的根处于单位圆内,即 =|+b 设原微分方程的特征根元=a+B,并代入上式得到 1+ha+j(<1 经整理得 由上式易见,在a-B坐标系下,即λ平面上的一个圆与Z平面上的稳定域是一 对应的,或者说Z平面上的稳定域映射到原微分方程的参数平面A上也是一个圆 圆心在 (-b0)处,半径为方,如图2所示 图22平面与Z平面的映射关系 由此可以看出,若原微分方程是稳定的,而其特征根在平面中映射到单位圆
33 2．欧拉法的计算稳定性分析 对测试方程根据欧拉法得到 y = y + h y +1 （2-27） 对上式进行 Z 变换得 zy(z) = (1+ h)y(z) 显然其特征方程为 z − (1+ h) = 0 根据控制理论可知，（2-27）式的稳定域在 Z 平面上是单位圆。（2-27）式稳定的条 件是它的特征方程的根处于单位圆内，即 z = 1+ h 1 设原微分方程的特征根 = + j ，并代入上式得到 1+ h + jh 1 经整理得 2 2 2 1 1 h h + + 由上式易见，在 − 坐标系下，即 平面上的一个圆与 平面上的稳定域是—一 对应的，或者说 平面上的稳定域映射到原微分方程的参数平面 上也是一个圆， 圆心在 − 0 1 , h 处，半径为 h 1 ， 如图 2-2 所示。 图 2-2 平面与 Z 平面的映射关系 由此可以看出，若原微分方程是稳定的，而其特征根在 平面中映射到单位圆