1、1第六章 常微分方程数值方法连续问题的离散处理寻求微分方程的解 在某些离散点上的值)(ty在 处的近似值 ；)(tyi iy记号： 所求函数在 处的（准确）函数值，)(iti计算得到的 处的函数（近似）值，i i节点， 步长； it 1iith等步长： ；以下若不另作说明，一般总记 为等步长。nabhi h 6.1 初值问题的数值方法考察微分方程初值问题： )1.6()(),0yabtatf由微分方程理论可知：若函数 关于 满足 条件，即存在与 无)(,tyfLipschtzy关的常数 ，使L ,),(),( 212121 dybatLytftf 初值问题的解存在且唯一；6.1.1 法及其变形
2、Euler1、 法 由 Taylor 展开式：)(2)(,)(2)()(1 iiiiiiii yhthftyhtytty 作局部化假设： ，并略去 项，便有ii O),()(1iii tft将此式右端作为 的近似 ，便得到公式ity1iyEuler 公式,1iiifh两者的差（即略去的 项）称为局部截断误差，记作 ：)(2 ),(htEi),(iiytEEuler 公式也可由其他方法导出，例如由第四章数值导数公式，可有：；,2)()()( 11 iiiiii thhttty 解出 ，并由 替换，便可得 Euler 公式。1i ,iiityft又如，根据 Newton-leibniz 公式：2d
3、tyttyiitii1)()(1将 以“0 次插值多项式”替换，即以)(t )()(iii tt代入积分，得到数值积分（左矩形）公式，及其误差：)(2)()(1 iit yhtdyii 又得到与 Taylor 展开式相同的表达式，从而又导出 Euler 公式。几何意义 折线法Euler若一个公式的局部截断误差为 ，则称该公式的精度为 阶，或该公)(1phOp式为 阶公式。pEuler 公式是 1 阶公式。注意，以上的截断误差是在局部化假设的前提下得到的，即认定 。)(iity倘若在每一步都按局部化假设，我们有 Euler 公式的总体截断误差： )(2)()(),(12habyhbyhEini
4、2、后退 法uler若取数值导数公式：)(2)()(11 iiii yhttyt 与前相同的推导过程，可以得到)()(,)(11 iiiii tftt 在局部化假设的前提下截去局部截断误差)(2),(iiyhtE便得到后退 法公式：uler),(11iiiitfy注意到此公式中的右端也有 ，需要求解关于 的方程才能得到 。1iy1iy1iy因此将这类公式称为隐式公式，而将可以通过直接计算得到的公式称为显式公式。后退 公式是一阶隐式公式，Euler 公式是 1 阶显式公式Euler6.1.2 多步法1、 梯形公式：在式 右端的积分中，取梯形积分公式，有dtyttyiitii1)()(1 )(1)
5、(231iiiii yhh3由此，并据微分方程，可得：梯形公式 ),(),(211 iiiii ytfytfhy局部截断误差： 1,3iitE这是一个 2 阶隐式公式。2、 Simpson 公式：在式 右端的积分中，取 Simpson 积分公式，有dtytytiitii 1)()(1 )(90)()43)( 511 iiiiii yhtyh由此，并据微分方程，可得：Simpson 公式 ),(),(),( 111 iiiiiii tftftfy局部截断误差： ；90),()5iihtE与以前的公式不同，用 Simpson 公式计算 ，必须有前 2 步的函数值：1iy和 。因此这种方法称为 2
6、步方法，而为启动此算法所需的最初的 2 个函1iyi数值： 称为表头。更一般的，若计算 必须有至少前 2 步函数值，则这10, 1i种方法称为多步法。具体地，若计算 必须有前 k 步的函数值，则这种方法iy称为 k 步方法 ，而为启动该方法所需的最初的 k 个函数值： 称为110,ky该方法的表头。与此相对，以前的方法计算 ，只须前 1 步的函数值 ，便1i i称为单步方法。因此，Simpson 公式是 2 步、4 阶、隐式方法。3、 Adams 方法（线性多步法）在式 右端的积分中，若取具有 k+1 节点的插值多dtyttyiitii1)()(1项式近似替代 作为被积函数，导出初值问题的求解
7、方法称为 Adams 方法。（1）显式 Adams 方法Adams-Bashforth 公式取 处的 构造插值多项式取代 ：kiijt,1, ),(jjytf )(ty)()(tRpty )()!1() ),)2(tykt tfllikikij jiij 4其中 。由于 ，有)()()()( ikiikijj tttt ,0)(1ittdtykytfdl dtykdflyy iiii iiiiii tkkij jtj t iktikij jjtii 11 111 )()!(,() )(!,)2)2(1 由此（以后为方便计，记 ） ，可得显式的多步法 Adams-Bash forth,ii公式及
8、其局部截断误差)(),( )2101iki kiiiii yrhtEfbfbAy这是 步、 阶的显式公式。k下表是 的 的数值（注意： ）4,10 rkjb),1,0(,AbjA02b34r0 1 1 211 2 3 -1 52 12 23 -16 5 833 24 55 -59 37 -9 72014 720 1901 -2774 2616 -1274 251 95以下是 的公式推导过程：2k 2122121121 )()()()( iiiiiiiiiii fttfttftttp作变量代换： ，以 为变量，当 的变化区间为 时， 的变化区shi ,1is间为 ，且 ，有,0dt)51623(
9、 2)()( 210 11 iiit iiiffh dsfsspii51 41(4) ()1204()(,()()23! 3!8iiti iii iihEthyttdtysdshy （2）隐式 Adams 方法Adams-Moulton 公式若取 处的 构造插值多项式取代 ，1,1, kiiijt ),(jjytf )(ty与前一样的方法，可得隐式的多步法 Adams-Moulton 公式及其局部截断误差)(),( ),2 11101iki kiiiiii yhrtEfbftfbAy 这是 步、 阶的隐式公式。k下表是 的 的数值（注意： ）：4,10 rkjb),1,0(, AbjkA012
10、b34r0 1 1 211 2 1 12 12 5 8 -1 43 24 9 19 -5 1 720194 720 251 646 -264 106 -19 63述评：从表可见，对相同的 ， 相同，而 ，特别是：kArbjj,，1Abj，而有若干 ，因而在存在计算误差时，由前步导致的误),0(kj1,bj差显然隐式公式要比显式公式小（显式公式对前步的误差会被放大，而显隐式公式则不会） ，而且局部截断误差也是隐式公式要比显式公式小，结论：隐式公式的稳定性一般比显式公式好。6.1.3 待定系数法6利用 展开式比较有关项的系数，可以直接导出公式待定系数法。Taylor例：求以下数值公式的系数 使公式
11、具有尽可能高的精度：,210)(01 iiiii ffhy解：由于 ，因此由 展开式，)1iitTaylor)(!41!3!2() 5)(1 hOyhty iiiiiii 同时，对所求公式右端各项也作相应的 展开，并乘以相应的系数：lriii yii hfh00 )(!312)( 5)4(312111 hOyhyt iiii 82 2222 hyf iiiiii 由于期望 尽可能准确，比较各对应项的系数，可得方程组：)(11iit，解之可得： ；61212101256310注意到局部截断误差是 ，因此1)(,(iii ythtE)(83)(2!3812!3!4),( 5)4(5)4(4)( h
12、OyhyhtE iiiii 显然，这就是 的显式 Adams-Bashforth 公式。k例：求以下数值公式的系数 使公式具有尽可能高的精度：,210)(1031 iiiii ffhy与上例完全一样，比较系数，可得公式：)2(4131 iiiii ff局部截断误差： (5), 6)(hOyhtEi 这个公式称为 Milne 公式，稍后我们将用到它。综上所述，从公式的构造过程可见，微分方程初值问题的算法公式基本上源于：Taylor 展开式、数值积分或数值微分，而插值方法一般也是通过数值积分或数值微分完成的。7 Taylor 展开式 待定系数法 直接来自于 Taylor 展开式，通过比较系数，形成
13、线性方程组，取得公式与局部截断误差； 数值积分前面的多步法等，基本上都源于此，例如梯形、Simpson方法等，而 Adams 方法用的是 ，当 用不1)()(1iitii dtyty)(ty同的点生成的不同的插值多项式代换时，便得到不同的公式，包括显式和隐式公式； 数值微分由 导出 Euler 公式，而由)(2)()(1iiii yhttty可导出后退 Euler 公式，若取 3 点数2)(11 iiii htty值微分公式可导出更多公式： )(32),(234)()(32)( 111200 yhthfytythty iiiii )(3),(2)(62)(101 yhtfythty iii )
14、(92),(2431)()(21)(31120 yhythfytythty iiiii 其中第二个公式就是“中点公式” ；当然它们也可以由 待定系数法 取得：例如第 3 个公式：用待定系数法求以下隐式公式的系数及局部截断误差： ),(111 iiiii ytChfByAy由 Taylor 展开式：)(21)()(, 6)( 4311 2hOiyCihiyCtyhtyChf BtBiiAiiii 比较： 6iiiii 为使公式具有尽可能高的精度，比较系数，得方程组：8321421CBACBA因此，公式是 ，),(43)( 111 iiiiii ythfyyt与之相应，有局部截断误差： 33343
15、(4)4122(,) ()6 9ii iEthhOOh6.1.4 问题的性态与算法的稳定性1、 问题的性态问题对原始数据的敏感性原始数据 问题的应有的结果)(F原始数据扰动 问题的结果问题的条件数 相对误差之比的上确界： )()()( FCondmF由微积分知识： )()()()( 初值问题：固定步长 ，初值 ； 计算结果：h0y0yFt),( )(,()()( 0020200 tf hOthfOhtyyF因此： 1y ),(1)(),()( 0000 ytfhtthftyCond 由于 ，当 ；)()0tty良 态病 态,yf或： 0)(t扰 动)()(,()( 20200 hOtyhOyt
16、hft ,()0ft又 (),(),(),( 2000 Oytftftf 得： 1hyy9因此： 良态0),(0ytf病态例： 方程的解 病态 （图示）)(2tty tcety2)(例： 方程的解 良态 （图示）t例： ， tftyt y1)(,1)(22当 病态， 当 良态。 （图示）22、 方法的稳定性方法：由初始误差导致的后期误差的可控性。中点法：由数值导数公式 2()()()26yththy 得 可导出3()()ythtt中点法 )(),(,)2 331 hOyhtEyfiii 此为 2 步 2 阶公式，但稳定性差。例： )0(,yy真解： 41 10539.)(, yett取 2 种
17、不同的表头： （即一步 法） ，),(,)2(1)( ythftEuler下表中列出对不同的步长： 计算 的近似值 ， 按后退NhN)0(法Euler计算获得， 按中点法由初值及表头 获得， 按中点法由初值及表头)1(Ny)1(y)2(N)2(1y获得，误差放大则是比值 。)2(1)(yNh后退 法Euler)0()1( )2(Ny误差放大1039765.306.4108.574124117925911909308.2. 6.1112141434210798说明：1、总体误差，由前可知后退 法：Euler)(2)(2).(),(11 yabhyhthbyENi iNii 10因此，对于后退 法
18、，若 ， （ ）可计算Euler210h 1,0,0()(1tet得；若 ，可计算得 显然实21021E4 242E际情况良好；而对于中点法，其总体误差3211(),(.)()()(NNi iiihEybhEtybay 考虑表头 都用准确值（注意在浮点数系中运算，仍是有误差的，这可认为10,是初始误差） ，若 ， （ ）可计算得2 1,0,0()(10tet，而 ，可计算得 ，而实际计算并非如此。13E4h53E2、对于同是中点法：由此表可见，两个中点法的不同结果是由不同的表头导致的，实际的不同仅是 的不同，而它们之间的误差（或不同）导致的)2(1,y计算结果的误差被放大 5000 倍甚至 1
19、 万余倍。说明中点法并不是好的方法，)1(y尽管它是一个 2 阶方法，但实际计算不如后退 法这样的 1 阶方法。Euler定义： 步方法，若存在常数k 00,hK使 ,1,max10 kiyKyjjkjii其中 及 是按此方法分别由表头 及 计算得到，ii 0,ky 110,ky则称此方法是稳定的。此定义因未限制 ，故实际可要求 ，因此本定义又“渐近稳定0h0h性”或“0 稳定性” 。实际计算关心对确定步长 ，一个算法，每步计算误差是否被放大？无法讨论一般方程研究对象：典型方程：其中： 。原因：0)(,)(yabtatyt 0)Re(对 维常微分方程组： ， 为 阶矩阵，它的特n,tYF ,t
20、DYn征值 反映了矩阵的性态，当 时，方程是良态，若 时则是病)Re()(态（工程中，通常称之为“稳定性” ， 时，系统稳定，否则，不稳定）)(。定义：对于典型方程 ，及确定的 ，一个 步方法，若由)(tyt hk表头 及 计算得到 及 ，存在常数110,ky 110,k iyi11使10,C,1,max10 kiyyjjkjii则称此方法对 是绝对稳定的，全体 的集合称为该方法的绝对稳定域。hh事实上，也可将此定义改写成 ,1,ax1 kiyCyjijikjii对于单步法则为：,2,1iiiii对典型方程，单步法 总有iii yhpy)()(1例如： 法： ， Euler iiii )hp1
21、)(后退 法： ， iiiii y,11 hp1)(梯形法： ， ；iiiiii hyyhy )2(,)(2111 2由于单步法 因此，当 则方法绝对稳定：1iiiip 1p法： ，在复平面上 是以 1 为中心，1 为半径的圆内部；Eulerh后退 法： ， 即 ，在复平面上 是以 1 为中心，11h为半径的圆的外部。为稍后介绍方便，对一般的 步方法，写成k )(00101 kjjikjji fhy即： 01101 kiiikiii fffy例如 法： Euler1iiihf后退 法： 1iiiy中点法： 021iiif梯形法： )(iiii法： ；Simpson 0)43111 iiiii
22、ffhy12定义算法 的特征多项式：)()(),h其中 kj kkkkj0 110)( jjk1分别称为算法（*）的第一、第二特征多项式。定理：对 ，若算法 的特征多项式的全体零点 都h)( ),21(ki满足 ，则 使该算法绝对稳定，全体这样的 的集合，称为该算法1i h的绝对稳定域。例如：中点法的特征多项式 ，则 2(,)1h0),(h的根 ，由于22,1h，1:0;1:0 2121 hh可见，中点法总是不稳定的。6.1.5 预估-校正方法原因：显式方法简单；但一般稳定性不好；隐式方法求解复杂；但一般稳定性好；将两者结合问题：如何组合1、误差估计 两个方法获得同一 的不同近似 ，以此估计误
23、差)(ityiy（1） 同阶、同步长的两个不同方法方法 a： )()( 111pii hyt方法 b： 2ii y由于步长 一般较小，因此可认为 ， （实际上就是估计此h )()(211ppy的值，因为一旦获得此估计值，便可得到 的估计)(1py 1,iiythtE值）将两式相减，得)()(11pii yhy即)(1)(iiph因此11)()( iiii yyt13（2） 不同阶、同步长的两个不同方法，设 qp方法 a： )()( 111pii yhyt方法 b： 2qii由于 是 的高阶小量，所以两者相加（减）时，后者可略去；将以上两1qhp式相减，可得)(111pii yhy因此：11)(
24、iiiit（3） 同阶、不同步长的同一方法，设步长 h步长 ：h)()( 11pii yhyt步长 ： 2ii仿（1） ，将两式相减，得)()(111 ppii yhy所以；)()( 111 piiii ht2、预估-校正法 （ Predictor-预估, Evaluation-计算, Correct-校正）（1） Heun 方法改进 Euler 法预估：Euler 法 ),(1iii ythfp校正：梯形法 ),(21iiiii ptfy保持二阶精度： )()(31Otii证明：记 为采用原始的梯形法获得的结果，即： ，由 )()(31hOyti本节初所作的假设：函数 关于变量 满足 条件，
25、),(ytfyLipschz ),(),(2),(,2)()( 11111 tftftfhtyty iiiiiiiii,),( yyytf iiiiiiipLhOi132)()()()( 311hOyttyiii (2) Adams-Bashforth-Moulton 方法取 的 Adams-Bashforth 显式公式作预估，Adams-Moulton 隐式公式校正：3k145937924121 iiiiii ffffhyp),(511iiiiiii pt对应局部截断误差公式：，)(7209)(),(720)( 251151 yhytyhpty iiii 根据前已介绍的误差估计方法，有。)(
26、9)( 11iiii pt由于在计算过程中，右端的值均已得到，故在每一步，均可对误差进行监控。在计算中，对于给定的相对误差界 ，Re若 htneSmypii 2270191tii01此处 Re（Relative error）相对误差界, （Absolute error）绝对误差Ae界， Sm (=Small)是一个与 相当的小值，目的是避免 时出现的估计失真ReA01iy（ 实际上是函数 的数量级，只须注意 ） 。ReA)(ty若相对误差过大，须步长减半（ ） ，新的进程中，将由 等原h21,iy有的值计算 的近似 ，而根据本方法，它是 的组合，)2(htyi21iy 2312,iiiff对此
27、中的原未取得的 与 的值，可采用插值法得到，但应注意，本方法是21if3i4 阶方法，因此，用插值法取得这些数值时应保证误差也应至少是 的，而)(5hO5 点插值多项式的余项： 且 ，因此,)()()(,)( 4321043210 ttttttfR )(tj用此 5 点插值方法可以实现误差是 的。为此，取最近相邻的 5 点构造插5hO值多项式，可得所求点处 的近似值：),(ytf343212(20960)18iiiiiif ff；587435iiiiiif f15(3) Milne-Simpson 方法以 Milne 公式为预估，Simpson 公式做校正。223411 iiiii ffhyp
28、1 1(,)iiiiiitp他们的局部截断误差公式分别为： ),(45)(11yhptyii；902)(ii与前相同，两者相减有：)(2151yhpyii可导出预估公式-Milne 公式的误差估计)(98)(11iiii pt与数值积分的 Romberg 公式的导出思想一样，现在设法将这部分的误差“归还”给预估的 。由于在得到 后，希望对 进行修正，使修正后的值更接近1ip1i 1i再进入校正公式参与计算，这样得到的校正值有望更符合 Simpson 公式)(1ity的实际效果。但由于尚未计算 ，只能用 近似替换 （实际1iyiipy1iipy上，是 ，即： ，因此有)()(1)55iiy1ii
29、iip132111148()94(,)3iiiiiiiiiiiiiiiphffmyyfftm注：这是 4 步方法，必须有表头 才能开始计算。预估步与校正320,y步当 便可执行，但在第一次的修正步计算 中，由于此前无 ，因此只能3i 43p令 ，也即第一次因为无法修正所以不进行修正。4pm（4） 修正 Milne-Hamming 方法解初值问题的 Hamming 公式： ),(283)9(81 112 iiiiiii ytffhyy )40), 5(1ythtEiii 以 Milne 公式为预估，Hamming 公式进行校正。与前相同，两者的局部截断误16差相减有：)(36012151yhpy
30、ii可导出预估公式-Milne 公式的误差估计)(2)( 11iiii pt校正公式-Hamming 公式的误差估计：)(9)( 11iiii yyt与前一样，现在设法将误差“归还”给预估的 ，以及 ；由于通常将最终1ip1iy结果记作 ，因此将 Hamming 校正公式计算结果记作 ；与前 Milne-1i icSimpson 方法一样，在修正 时， 尚未产生，因此只能使用 ，从而形成1ip1ici公式： )(129),(283)(8(12411 1113iiii iiiiiii iiii iiiii pcy mtffhycpmff注：仍与 Milne-Simpson 方法一样，这是一个 4
31、 步方法，必须有表头才能开始计算。预估步与校正步当 便可执行，但在第一次的修3210,y 3i正步计算 中，由于此前无 ，因此只能令 。4m3p4pm注：如果从求解隐式方程的迭代法来看，这校正步可以认为是迭代了一步，而迭代的初值则是预估值，或是预估的一个修正值。因此从迭代收敛的角度看，应对初值提出一定的条件，通过研究，可以得到如下结果，各方法对步长的限制：A-B-M4 方法： ),(75.0ytfhMilne-Simpson 方法： ,4.tfy修正 Milne-Hamming 方法： 。),(69.0tfhy6.1.6 Runge-Kutta 方法将 Heun 方法（改进 Euler 法）改
32、写成如下形式：17),()21112KyhtfKyiiii此为一单步、2 阶、显式公式（比较梯形：2 阶隐式） 。若注意此处的 ，21,K它们是两个不同的增量，Heun 方法则是从 出发加上若干个增量的加权平均iy值后得到 ，将此推广到一般形式， 出发加上 个增量的加权平价取得 ，1iyim1iy便有： -级 Runge-Kutta 公式m ),(),)( 12122211 mmimimii mii KKyhtfKtf 为导出该公式，先引入二元函数 Taylor 展开式： 1,0,),()()!1 ),(!,),),( 2211 kyhtfykthn ytftntfkyhtf n n以二级 R
33、unge-Kutta 公式为例： ),()12121KyhtfKyiiii 先考虑准确值 的展开：)1i )(2() 31 hOyttyiiiii 其中fytyt ffytdy )(,因此 )(21)( 31 hOfhft ytiii 另一方面，考虑近似值 的展开，由二元函数 Taylor 展开式，二级1iRunge-Kutta 公式中，)(212 hfKfhKyti 故应有： )(322211 hOffffy ytiiii 18将准确值 与近似值 展开式作比较，应)(1ity1iy22若令 ，则c2，cc1,1取 ，则 ， Heun 方法（改进 Euler 法） 。12取 ，则 变形 Eul
34、er 法：c 2,01)2,()1212KyhtfKyiiiii（注意此公式与中点法的区别）以上两公式均为二阶二级 Runge-Kutta 公式。常用的四阶四级 Runge-Kutta 公式（称为经典、古典或标准 Runge-Kutta公式Classical Runge-Kutta ）11234122343()6,)(),()iiiiiiiiiyKhftyftyKhK几何意义：若将上述 公式改写为：4R)2(6311 ffyii其中 为分别为原式中 中相应的函数 值。比较)4,32(jf )4,321(j jf在区间 上的数值积分 Simpson 公式,1it )(6),( 11 imiiti
35、i ffhydtfyii其中 分别表示 在区间 的左、中、右端点的值；它们1,imiff ,it分别有： ；因此可简略地得到该 公式4132,ffii 4RK的误差公式是)(80),(5ii yhtE19误差估计：在实际计算中，由于 Runge-Kutta 方法是一个高精度的单步方法。因此，通常作为求解常微分方程的一个独立的方法使用。在计算时，常常要对计算值进行误差估计。以四阶四级的标准 Runge-Kutta 公式为例：若用同一公式的不同步长方法进行误差估计，例如用步长 和 计算获得h2的不同近似进行，为此必须计算 11 个 值（其中： 计算 4)(1ity ),(ytf 1iiy个 值，
36、计算 8 个 值，而 是重复的） 。,f 12iiiyy ),(1itfK这比不估计误差增加 7 个 值的计算量。),(tf一个比较简单的方法是 Runge-Kutta-Fehlberg 方法65431126 4315 324 1321665431 51 20175948360),( )867,2 )0556439,( 197672,1)93,8()4,)( )(,(),20956281356 ,),7402( KKyhtEfKyht KKfyhtKfyt hOtEKKy hiii iiiiiiiiiiiiiii 以此作误差估计，在计算过程中仅需计算 6 个 值的计算量。),(ytf小结：常微
37、分方程求解的理论问题，主要是问题的性态和方法的稳定性。算法包括基本算法的构造（或：形成）及组合。前者就是通过 Taylor 展开式Runge-Kutta 法、待定系数法形成如 Milne 法等，或数值积分（例如：Adams 方法，梯形法，Simpson 法等） ，或数值微分（或称数值导数）形成；当然同一方20法实际上可以通过多种方法导出。后者就是利用前面已得到的基本方法，进行适当组合，一方面利用显式方法简单，作为预估，利用隐式方法比较稳定，作为校正；进一步，同阶的不同方法，或同一方法不同步长的局部截断误差比较，将误差中的估计项 变换成可计算的 ，得到实际误差估)(1py1iipy计值，或直接用
38、来比较误差是否满足要求，或用来对对公式作修正。6.1.7 常微分方程组、高阶常微分方程1、常微分方程组初值问题：常微分方程组初值问题的解法与常微分方程初值问题的解法几乎是一模一样的，唯一区别就在于所有的公式都用向量代替纯量。为此，记向量 )(,)(,)()()(,)()( 2122121 tytytftttfttytt nnnn YFY则可将常微分方程组初值问题 ,)(,)(,)( )(,)(,() 02021012121 nnnnyayaytttft tttft 改写成简单的向量形式：.0)(,()YFYtt对此向量形式的初值问题的解法，只需将以前一维方程的数值方法改为向量形式便可。例如，解
39、常微分方程组的 Euler 方法，写成向量形式：，),YF(Yii1i th标准 Runge-Kutta 公式的向量形式：21)KY,F(tK),(t)K(Ky3ii221ii2i1 4321ii hh62、高阶方程：对于高阶常微分方程初值问题，一般可以将它改写成常微分方程组，然后用前述的向量方法求解。例如：3 阶常微分方程初值问题，000)(,)(,)(,)(,() yayaytytfty 令， )(,)(,)(,)( tvutttvutyftvutY并记 0,)(,)(,)(, ytvutyfttYF则原高阶方程就可以改写成如下形式，从而可按标准方法求解， 0)(,()Y att6.2 边
40、值问题数值方法简介6.2.1 差分方法讨论常微分方程：， btatfytq)()(其中 ；对于一般方程：0)(tttfttp )()()(可以通过以下方法转化为上述形式。对方程乘以 ，得：dtpe)(， dtpdtpdtp feytqey )()()(令 ，使tu，即 dtp)( dtput)(22从而一般方程等价于 dtpdtpefeytqduty)()(在此方程两端同乘以 ，便得t)(dtpdtpefety)(2)(22利用 的反函数（由 知 是关于 的单调函数，故有反函数） ，)(tu0tutt将 转化为 的函数，便可得，)(2RyQd此即为所讨论的形式。第一边值条件 )(,)(ba第二
41、边值条件 yy第三边值条件 ；00,)()(010 1、 差分离散法用数值导数替代导数等分区间 步长： ， 节点： ，nabh)( nihati ,1由 )()(iiii tftyqty及数值微分： )(12)()24yhty得 ()2(11 iiiiiii tftqhtt 略去局部截断误差（Local truncation error）： 注意：)(42iiyhR；)(2hOR得：在内点 处的离散逼近方程：)1,2(nit，1,2,1 nifyqyiiii 或。,)2(211 ihfhiiiii 求方程的数值解， 的近似 尚差 2 个方程；由边界条件得),0ntyiiy到：对第一边值条件，直
42、接代入即可。对第三边值条件（因第二边值为其特殊情况）利用数值导数公式： )(3)(4)(321)( 022100 fhxffxfhxf 23)(3)()(4)(21)( 2210 fhxffxfhxf 略去截断误差部分 项,可得方程:2Oy110)3( hynnn 102 2)3(4至此，已有 个方程，可以求解。2、 差分方程求解：对第一边值条件，方程组： 21211211221 )(1)()()( hffhyyhqhqhqh nnnn 如果将上述系数矩阵按追赶法分解，简记为LUT其中 12112, nnUll 显然： ,43,1,)2(2211 illhqii 回顾线性方程组求解过程，可知，
43、此追赶法与选列主元 消去法是一Gaus致的，因此算法是稳定的。但是，对第三边值条件，按上述介绍的方法离散化，将无法解决算法的稳定性问题。因此需要重新离散化。为此，将节点按下述方法确定：，1,0)2(, nihatnbhi 在内点处（ ） ，仍按以前方法离散，而对边界点，则由ni,1 )(2)()(2)( 62yhtyhtty将公式中的 分别以 和 代入，以 代替 ，分别代入第 3 边值条件，再将tab24项作为截断误差略去，得)(2hO 1011010 2,2 nnyhyyy将此按内点方程形式（非对角元为 1 ）标准化，便形成方程组： 212121)2(1)2(1 0201000 hfhfyyhqhqh nnn 与前面的分析类似，将系数矩阵按追赶法分解，简记为LUT其中 1011, nnUll 注意到 ，因此1200h,1)2(,1201 lhql由此推出 ，最后有 。注意，nilii ,3,0,11nnl的情况仅在 时出现。而此时事实上已出1nii ,00现原方程退化为 形式，可直接由积分得到函数。除此以btatfy)(外，上述解线性方程组，总是稳定的。6.2.2 打靶法（Shooting Method）讨论方程（2.1）)(,)(byaybtattf将
