截断误差正比于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 平面的映射关系 由此可以看出，若原微分方程是稳定的，而其特征根在 平面中映射到单位圆