2收敛性稳定性RK方法

以上我们讨论了求解问题（7-1）,（7-2）的单步法和多步法。应关注三个问题：一、数值方法的局部截断误差和阶二、在离散点tn处的数值解un是否收敛到精确解u(tn)三、数值方法的稳定性具体说，对于上述两类方法求近似解（数值解）还误差估计、收敛性和稳定性。对于第一个问题前面我们已经讨论过，而关于数值方法收敛性问题我们在这里不详细讨论，只给出一些基本结论性的结果，即：对单步法，当方法的阶p≥1时，有整体误差)()(pnnnhOutuE故有0lim0nhE，因此方法是收敛的。对于多步法，若方法是k步p阶法，那么（7-24）是一个k阶差分方程，引入多步法（7-24）的第一特征多项式和第二特征多项式：0(),kjjjkjjj0)(定义7.1若（7-24）的第一特征多项式ρ(λ)的所有根在单位圆内或圆上（︱λ︱≤1），且位于单位圆周上的根都是单根，称多步法（7-24）满足根条件。第二特征多项式第一特征多项式定理7.2若线性多步法(7-24)的阶p≥1，且满足根条件，则方法是收敛的。对于常用的数值方法都是满足收敛性条件的。下面我们着重讨论第三个问题，即数值方法的稳是有误差的，且这些误差将在计算中传递下去。定性问题。110,,,kuuu误差积累无限增长，则会歪曲真解，这样的算法是不如果能用的。用多步法计算时，各种因素如初值精确解为22()(1)utt。考虑二步三阶显式法：21145(42)nnnnnuuuhff例如初值问题124,utu02;t(0)1u取步长h=0.1，初值u0=1，附加值：221(1)(0.1)uhh。精确解数值解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数值结果表在开始几步数值解与精确解符合，但在再往后算，数值解的误差则急剧增长，完全歪曲了真解.通常人们都是通过模型方程来讨论方法的数值稳定性。uu(7-32)而一般形式的一阶微分方程总能化成（7-32）的形式。。因为实际计算时，h是固定的。当某一步un有舍入误差时，若以后的计算中不会逐步扩大，称这种稳定性为绝对稳定性。此后，若不做特殊说明，都是指绝对稳定性。模型方程为：本书中数值方法的稳定性也是如此。前提是求解好条件问题，其中Re(μ)＜0。另外，我们也不考虑h→0时方法的渐近稳定性例如，对最简单的Euler法2,1,0,1nhfuunnn(7-33)用其求解模型方程（7-32）得到1nnnuuhu当un有舍入误差时，其近似解为nu~，从而有nnuhu~)1(~1取nnnuu~，得到误差传播方程,)1(1nnh(1),0,1,2nhun记hh，只要11h都不会恶性发展，此时方法绝对稳定。，则显式Euler方法的解和误差从11,h可得20h。即20h时，(-1,0)为圆心，1为半径的单位圆。又由于实数μ＜0，（7-33）绝对稳定，若μ为复数，在hh的复平面上，则11h表示为以20110,1绝对稳定区域绝对稳定区间定义7.2一个数值方法用于求解模型问题（7-32），若在平面中的某一区域D中方法都是绝对稳定的，而在区域D外，方法是不稳定的，则称D是方法的绝对稳定区域；绝对稳定区间。它与实轴的交称为例如，显式Euler方法的绝对稳定区域、区间。如图现在考察多步法(7-24)，将它用于解模型方程(7-32)得到k阶线性差分方程0kjnjju(7-34)若取hh，则记（7-34）的特征方程为0)()(h(7-35)其中kjjj0)(kjjj0)(0kjnjjhu由k阶线性差分方程的性质我们可以得到如下结论，区域：()1,1,2,,jhhjkD例如，对于k=1时，考虑隐式方法中最简单的后退Euler法111(,)0,1,nnnnuuhftun方程（7-35）的根都在单位圆内（︱λ︱＜1)，则线性多步法（7-4）关于hh绝对稳定,其绝对稳定域是复平面h上的其特征方程为：()()h(1)10h若特征得11,1h当11h时，11,故11h就是隐式Euler法的绝对稳定区域。当μ＜0为实数时，绝对稳定区间为(-∞,0)。h平面上以（1.0）为圆心的单位圆外区域。它是11h11h0Re2,0h0,当Reμ＜0时，它位于平面上y轴左侧区域。h又如，梯形法,1,0)(2111nffhuunnnn其特征方程为：()()h其根1()h当Reμ＜0时，121,12hh故梯形公式h平面的左半平面。绝对稳定区间为(-∞，0)。的绝对稳定域是这样检验绝对稳定性归结为检验特征方程(7-35)的根是否在单位圆内（︱λ︱＜1)。对此有很多判别法，如Schur准则、轨迹法。12,12hh11022hhk=1～4的隐式Adams类方法的绝对稳定区间（μ＜0为实数）。步阶绝对稳定区间12（-∞，0）23（-6.0，0）34（-3.0，0）45（-1.8，0）实系数二次方程λ2-bλ-c=0的根在单位圆内的充要条件为:21cb这里我们给出一种简单的、常用的判别法：例证明求解一阶常微分方程初值问题：),,(utfu0)0(uu的差分格式)85(121212nnnnnfffhuu收敛并求其局部截断误差主项、绝对稳定区间。解：由差分格式可知，2(),()(1)0,则其特征值满足根条件。令得λ1=0,λ2=1。2581(),121212故此为隐式二步三阶法，其局部截断误差主项为：4(4)1()24nhut。0110C1C2C3C4C18521121212013(14)220321185122621212776604311851224!3!1212151612424240注意，0120,1,1,从而012185,,,121212由定理7.2可知，此方法收敛。()()h128125hh而12125hh自然成立。hhhh512412512812得124h128124hh即有332,hh可得其绝对稳定区间：60h。又其特征方程为21280125125hhhh而使得︱λ︱＜1的充要条件为：124125hh12125hh现在再由128h124h2581110121212hhh进一步而自然成立。显式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展开法相等。(,(),)tuthEuler是最简单的单步法。单步法不需要附加初值，所需的存储量小，改变步长灵活，但线性单步法的阶最高为2，Taylor展开法，用在同一点(tn,un)的高阶导数表示，这不便于计算。(,(),)tuth先引进若干记号，首先[t,t+h]取上的m个点:123,mtttttth令,,,2,1mihathattiii11,2,,,iijijbaimmiiicc11,0213132121mmmmbbbbbbRunge-Kutta矩阵B为严格下三角矩阵：满足显式Runge-Kutta公式假设三组系数已给定，则求解（7-1），（7-2）的一般1(,,),nnnnuuhtuh（7-12）其中1(,,)mnniiituhck(7-13)(7-14)显式Runge-Kutta法的计算过程如下：1(,),nnkftu22211(,),nnkfthauhbk33311322(,()),nnkfthauhbkbk31323,bba11(,),mmnmnmjjjkfthauhbk212,ba11mmjmjba0,1,n现在推导一些常用的计算方案，特别地，给出m=3显式首先将u(t+h)在t处展开到h的三次幂,即：3()41()()()()!lllhuthututOhl(7-15)其中(7-16)Runge-Kutta法的推导。()(,(),),uthtuth2311ˆ(,,)()()26utuhfhfhfffOhtuffff2ˆ2tttuuuffffff其次，由二元函数f(t,u(t))在(t,u)点处的Taylor展开式可得：,))(,(1ftutfk),(1222khauhatfk21()tufhafkf)(ˆ21~32222hOfahfhaf)()ˆ21~(~323322233hOfaffbahfhafku于是，将k1,k2,k3代入（7-13）中，即31(,,)iiituhck1cf223323231ˆ2ucfafhabffafh112233ckckck222321112()2tttuuuhafkfkfOh22231ˆ2cfhafafh3()Oh1232233cccfhacacf31(,,)iiituhck2223232322331ˆ2()()2uhabcffacacfOh(7-17)由(7-16)已得(3)211(,,)()()()26tuhututhuth其中,tuffff2ˆ2tttuuuffffff。2311ˆ()26ufhfhfffOh合并f(t,u,h)展开式中的各阶hl(l=0,1,2)的系数，得比较(,,)tuh和(,,)tuh的同次幂系数，231111ˆ()2233ufhfhfffOh可得（一）m=1此时c2=c3=0,f(t,u,h)=c1f,比较h的零次幂，知(,,),tuhf方法（7-21）为一级一阶Runge-Kutta法，实际上为Euler法。（二）m=2,此时c3=0，则2231222221ˆ(,,)()2tuhccfhacfhacfOh与(,,)tuh比较1,h的系数，则2112221cacc它有无穷多组解，从而有无穷多个二级二阶方法。2311ˆ()26ufhfhfffOh(,,)tuh,21,1,0221acc（1）称为中点法。此时12,nnuuhk1(,),nnkftu2111,22nnkfthuhk(7-18)三个常见的方法是：1221,1,2cca（2）称为改进的Euler法。此时112,2nnhuukk1(,),nnkftu21(,)nnkfthuhk(7-