第7章-2 收敛性稳定性R-K方法

常微分方程的数值解法的

收敛性、稳定性

第7章--2

以上我们讨论了求解问题（7-1）,（7-2）的单步法和多步法。应关注三个问题：、数值方法的局部截断误差和阶二、在离散点tn处的数值解un是否收敛到精确解u(tn)

三、数值方法的稳定性具体说，对于上述两类方法求近似解（数值解）还误差估计、收敛性和稳定性。对于第一个问题前面我们已经讨论过，而关于数值方法收敛性问题我们在这里不详细讨论，只给出一些基本结论性的结果，即：对单步法，当方法的阶p≥1时，有整体误差故有，因此方法是收敛的。对于多步法，若方法是k步p阶法，那么（7-24）是

一个k阶差分方程，引入多步法（7-24）的第一特征多项式和第二特征多项式：定义7.1

若（7-24）的第一特征多项式ρ(λ)的所有

根在单位圆内或圆上（︱λ︱≤1），且位于单位圆周上的根都是单根，称多步法（7-24）满足根条件。第二特征多项式第一特征多项式定理7.2

若线性多步法(7-24)的阶p≥1，且满足根条件，则方法是收敛的。对于常用的数值方法都是满足收敛性条件的。下面我们着重讨论第三个问题，即数值方法的稳是有误差的，且这些误差将在计算中传递下去。定性问题。误差积累无限增长，则会歪曲真解，这样的算法是不如果能用的。用多步法计算时，各种因素如初值精确解为考虑二步三阶显式法：例如

初值问题取步长h=0.1，初值u0=1，附加值：

精确解数值解01.00000001.00000000.11.02010001.02010000.21.08160001.08120000.31.18810001.18923850.41.34560001.33886600.51.56250001.5929935………1.04.0000000-68.6398041.04.8841000+367.26392………2.025.0000000-6.96×108数值结果表在开始几步数值解与精确解符合，但在再往后算，数值解的误差则急剧增长，完全歪曲了真解.通常人们都是通过模型方程来讨论方法的数值稳定性。(7-32)而一般形式的一阶微分方程总能化成（7-32）的形式。。因为实际计算时，h是固定的。

当某一步un有舍入误差时，

若以后的计算中不会逐步扩大，称这种稳定性为绝对稳定性。此后，若不做特殊说明，都是指绝对稳定性。模型方程为：本书中数值方法的稳定性也是如此。前提是求解好条件问题，其中Re(μ)＜0。另外，我们也不考虑h→0时方法的渐近稳定性

例如，对最简单的Euler法(7-33)用其求解模型方程（7-32）得到

当un有舍入误差时，其近似解为，从而有取，得到误差传播方程记，只要都不会恶性发展，此时方法绝对稳定。，则显式Euler方法的解和误差从可得

即时，(-1,0)为圆心，1为半径的单位圆。又由于实数μ＜0，（7-33）绝对稳定，若μ为复数，在的复平面上，则表示为以绝对稳定区域绝对稳定区间

定义7.2

一个数值方法用于求解模型问题（7-32），若在

平面中的某一区域D中方法都是绝对稳定的，而在区域D外，方法是不稳定的，则称D是方法的绝对稳定区域；绝对稳定区间。它与实轴的交称为例如，显式Euler方法的绝对稳定区域、区间。如图现在考察多步法(7-24)，将它用于解模型方程(7-32)得到k阶线性差分方程(7-34)

若取，则记（7-34）的特征方程为(7-35)

其中由k阶线性差分方程的性质我们可以得到如下结论，区域：

例如，对于k=1时，考虑隐式方法中最简单的后退Euler法方程（7-35）的根都在单位圆内（︱λ︱＜1)，则线性多步法（7-4）关于

绝对稳定,其绝对稳定域是复平面

上的其特征方程为：若特征得当时，故就是隐式Euler法的绝对稳定区域。当μ＜0为实数时，绝对稳定区间为

(-∞,0)。平面上以（1.0）为圆心的单位圆外区域。它是当Reμ＜0时，它位于平面上y轴左侧区域。又如，梯形法其特征方程为：其根当Reμ＜0时，故梯形公式平面的左半平面。绝对稳定区间为(-∞，0)。的绝对稳定域是这样检验绝对稳定性归结为检验特征方程(7-35)的根是否在单位圆内（︱λ︱＜1)。对此有很多判别法，如Schur准则、轨迹法。k=1～4的隐式Adams类方法的绝对稳定区间（μ＜0为实数）。步阶绝对稳定区间12（-∞，0）23（-6.0，0）34（-3.0，0）45（-1.8，0）实系数二次方程λ2-bλ-c=0的根在单位圆内的充要条件为:

这里我们给出一种简单的、常用的判别法：例证明求解一阶常微分方程初值问题：的差分格式收敛并求其局部截断误差主项、绝对稳定区间。解：由差分格式可知，则其特征值满足根条件。令得λ1=0,λ2=1。故此为隐式二步三阶法，其局部截断误差主项为：注意，从而由定理7.2

可知，此方法收敛。而自然成立。得即有可得其绝对稳定区间：又其特征方程为而使得

︱λ︱＜1的充要条件为：现在再由进一步而自然成立。显式Runge-Kutta法

第7章--37.1.4

四阶显式Runge-Kutta法

通过观察我们发现显式Euler法和隐Euler法各用到了u(t)在[t,t+h]上的一个一阶导数值，它们都是一阶方法。改进的Euler法用到了u(t)在[t,t+h]上的两个一阶导数值，它

梯形法和们都是二阶方法。我们要研究的Runge-Kutta方法是一种高阶单步法，它使用u(t)在[t,t+h]上的斜率f在一些点的值非线性表示使得其局部截断误差的阶和Taylor展开法相等。

Euler是最简单的单步法。单步法不需要附加初值，所需的存储量小，改变步长灵活，但线性单步法的阶最高为2，Taylor展开法，用在同一点(tn,un)的高阶导数表示，这不便于计算。先引进若干记号，首先[t,t+h]取上的m个点:令

Runge-Kutta矩阵B为严格下三角矩阵：满足显式Runge-Kutta

公式假设三组系数已给定，则求解（7-1），（7-2）的一般（7-12）其中(7-13)(7-14)显式Runge-Kutta法的计算过程如下：现在推导一些常用的计算方案，特别地，给出m=3显式首先将u(t+h)在t处展开到h的三次幂,即：(7-15)其中(7-16)

Runge-Kutta法的推导。其次，由二元函数f(t,u(t))在(t,u)点处的Taylor展开式可得：

于是，将k1,k2,k3代入（7-13）中，即(7-17)由(7-16)已得其中合并f(t,u,h)展开式中的各阶hl(l=0,1,2)的系数，得比较和的同次幂系数，可得（一）m=1此时c2=c3=0,f(t,u,h)=c1f,

比较h的零次幂，知方法（7-21）为一级一阶Runge-Kutta法，实际上为Euler法。

（二）m=2,此时

c3=0，则与比较1,h的系数，则

它有无穷多组解，从而有无穷多个二级二阶方法。（1）称为中点法。此时(7-18)三个常见的方法是：（2）称为改进的Euler法。此时(7-19)

（3）此时（三）m=3比较（7-16）和（7-17），令

f,h,h2的系数相等，并注意的任意性，得四个方程不能完全确定六个系数，因此这是含两个参数的三级三阶方法类。常见方案有：

Heun三阶方法。

此时取(7-20)（2）Kutta三阶方法，(7-21)此时（四）m=4将（7-16）和（7-17）展开到h3,比较的系数，则含13个待定系数的11个方程，由此得到含两个参数的四级四阶Runge-Kutta方法类，其中最常用的有以下两个方法：经典四阶Runge-Kutta方法：（7-22）Butcher表分别为：

以上讨论的是m级Runge-Kutta法在m=1,2,3,4时，可分别得到最高阶级一、二、三、四阶，但是，通常m级Runge-Kutta

方法最高阶不一定是m阶。若设p(m)是m级Runge-Kutta方法可达到的最高阶，可证：改进的Euler法计算公式为：经典Runge-Kutta法计算公式为：例1

分别用Euler法，改进的Euler法(7-27)和经典Runge-Kutta法(7-30)求解初值问题:解：Euler法计算公式为：

三个方法计算结果比较表n

tn精确解u(tn)Euler法改进Euler法经典Runge-Kutta法数值解un误差数值解un误差数值解un误差00000000010.50.4333330.5000000.0666670.4000000.0333330.4332180.00011521.00.6666670.8000000.1333330.6350000.0316670.6663120.00035531.50.8076920.9000000.0923080.7875960.0200960.8074230.00026942.00.9333530.9856150.0512820.9210250.0123080.9331560.000171作比较

，计算结果见下表：取步长h=0.5,tn=0.5n,n=0,1,2,3。并与精确解：下面考察Runge-Kutta法的绝对稳定性。根据定义，对m级p阶Runge-Kutta法（7-12）取f=mu，则（其中