1、第二章第二章 常微分方程数值解法常微分方程数值解法Chapter 2 Numerical Solution of Ordinary Differential Equation(s)22.1 引引 言言 Fm假设牵引力F为恒定值dvFdtm为了确定待定常数，可以给定初始条件：或者给定边界条件：00tvvendtendvv,dvmF t vdt假设牵引力不恒定呢？求速度)(tv3虽然求解微分方程有许多解析方法，但解析方法只能够求解一些特殊类型的方程。还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。利用计算机解微分方程主要使用数值方法。解析方法与数值方法解析方法与数值方法4主要研究
2、对象：初值问题 00)()(),(yxybxayxfdxdy00()(,()xxy xyf x y xdx求数值解：求数值解：求求y(x)在离散数据点在离散数据点xk处的近似值处的近似值yk。y=y(x)xyx0=ax1x2x3xk-1xn-1xn=b xk5则称 f(x,y)对y 满足李普希兹条件，此时初值问题在a,b上存在唯一的连续可微的解唯一的连续可微的解。定理1：设 f(x,y)是定义在区域 G=(x,y)|axb,yR上的连续函数，若存在正的常数 L 使：1212|(,)(,)|f x yf x yL yy12,xa byy使得对任意的及都成立，(Lipschitz)条件条件有解条件
3、：有解条件：6则Lipschitz条件成立：在 f(x,y)对y可微的情况下，若偏导数有界：(,),(,)f x yLx yGy有解条件的判断：有解条件的判断：*121212(,)|(,)(,)|()|f x yf x yf x yyyL yyy7定理2：如果f(x,y)在G=(x,y)|axb,yR上满足Lipschitz条件，则初值问题是适定的。适定性：指初值问题中，初始值y0及微分方程的右端函数f(x,y)有微小变化时，只能引起解的微小变化。适定性条件：适定性条件：8在解的存在区间 a,b上取n+1个节点 bxxxxan210这里把 1,.iiihxx i=0,1,n称为由xi到xi+1
4、的步长 一般取成等间距的：nabh求解方法：求解方法：步进法步进法（分为（分为单步法单步法和和多步法多步法）数值方法的基本思想数值方法的基本思想9本章规定：本章规定：在 处初值问题的理论解用 表示，数值解法的近似解用 表示。nx()ny xny记 ，它和 是不同的，后者等于 。(,)nnnff xy(,()nnf xy x()ny x102.2 几种简单的数值方法几种简单的数值方法（一）（一）欧拉（欧拉（Euler）法）法 00)()(),(yxybxayxfdxdy001()(,)0,1,2,.nnnnyy xyyhf xyn 11、泰勒公式解释、泰勒公式解释、求导的两点公式解释、求导的两点
5、公式解释、积分公式解释、积分公式解释欧拉公式的的分析解释欧拉公式的的分析解释001()(,)nnnny xyyyhf xy 00(,)()yf x yy xy 121)1()(21)!1()(!)(!2)()()()(pppnpnnnnhpyhPxyhxyhxyxyxy2()(,()()nnny xhf xy xO h泰勒公式解释泰勒公式解释其中：1nnhxx可以得到：001()(,)nnnny xyyyhf xy 13求导的两点公式解释求导的两点公式解释1()()()()(,)nnnnnny xy xyxhyxf xy可以得到：可以得到：001()(,)nnnny xyyyhf xy 140
6、0(,)(1.1)()(1.2)yfx yaxby xy 对微分方程对微分方程(1.1)(1.1)两端从两端从1nnxx到进行积分进行积分11(,)nnnnxxxxy dxfx y dx11()()(,)nnxnnxy xy xf x y dx积分公式解释积分公式解释1511)(,()(,()nnnnnnnnxxf xy xyyf xy xh右端积分用左矩形数值求积公式：右端积分用左矩形数值求积公式：11()()(,)nnxnnxy xy xf x y dx001()(,)nnnny xyyyhf xy 即：即：(,)f x yfxnx1nx16欧拉公式的的几何描述yxx0 x1 x2 x3
7、x4y=y(x),.2,1,0),()(100nyxhfyyyxynnnn00(,)()yf x yy xy 00,xy 11,xy17例题例题1：（取步长：（取步长h=0.1）用）用Euler方法求满足条方法求满足条件件 的的y(t)数值解。数值解。22,12(1.0)0.0tdyyt etdtty )2(21ntnnnnnetythyy 解：解：2718281828.0)2(1.00200001 tetytyy122111120.1()0.684755578tyyyt et 18ntnyny(tn)y(tn)-yn01234101.01.11.21.31.43.00.00.271830.6
8、84761.276983.0935515.398240.00.345920.866641.607223.6203618.6830.00.074090.181880.330240.526813.28486数值解列表为数值解列表为)(2eetyt 22,12(1.0)0.0tdyyt etdtty )2(21ntnnnnnetythyy 19欧拉方法的误差估计欧拉方法的误差估计()-nnney xy 通过数值方法进行计算时，考虑每一步产生的误差，通过数值方法进行计算时，考虑每一步产生的误差，从从x0开始一步步累积到开始一步步累积到xn，称，称 为该数值方法在为该数值方法在xn点处的点处的整体截断误
9、差整体截断误差，该误差，该误差与与xn 及之前的各步计算误差都有关系。及之前的各步计算误差都有关系。20欧拉方法的误差估计欧拉方法的误差估计*1*1111()(,(),)()-()-()(,(),)nnnnnnnnnnnyy xhxy xhRy xyy xy xhxy xh 为了简化分析，着重分析为了简化分析，着重分析xn点单步计算产生的误差，点单步计算产生的误差，即把即把xn点之前的计算当作无误差：点之前的计算当作无误差：称该误差为数值方法在称该误差为数值方法在xn+1点处的点处的局部截断误差局部截断误差。12111(,()()ppnnnRH xy xhO h局部截断误差的第一个非零项为局部
10、截断误差的第一个非零项为局部截断误差主项局部截断误差主项。21欧拉方法的误差估计欧拉方法的误差估计 如果求解公式的局部截断误差为如果求解公式的局部截断误差为R(h)=O(hp+1)，则则称该求解公式称该求解公式具有具有p阶精度阶精度，称该方法为称该方法为p阶方法阶方法。定义：定义：欧拉方法：欧拉方法：212()()()()2!()(,()()nnnnnnny xy xy xy x hhy xhf xy xO h*1()(,()nnnnyy xhf xy x*2111()-()nnnRy xyO h具有具有1阶精度。阶精度。22（二）向后欧拉法（二）向后欧拉法 hxyxyxynnn)()()(1
11、1 11100(,),0,1,2,(),nnnnyyhf xyny xy （1）方法）方法 其公式为：其公式为：()(,)y xf x yxfnx1nx23（2）局部截断误差）局部截断误差*232111()()()()2nnnnyxRy xyhO hO h 231()()()()()2nnnny xy xy xy x hhO h*11()()nnnyy xy xh11100(,),0,1,2,(),nnnnyyhf xyny xy 21()()()()nnny xy xy x hO h24例题例题2：用向后：用向后Euler法解初值问题法解初值问题 1)0()10(,2yxyxyy11112(
12、)nnnnnxyyh yy 向后向后Euler法的公式为法的公式为 解：解：x0=0,y0=1,取h=0.125方法比较及推广方法比较及推广：nEuler方法方法 显式公式显式公式n向后向后Euler方法方法 隐式公式隐式公式 n解一个非线性方程解一个非线性方程 难求解难求解 n显式和隐式相结合显式和隐式相结合 隐式的显化隐式的显化 01110()(,),nnnnxyyyhfy xy 100(),)nnnnxyyyhfyy x 26计算公式为：计算公式为：001111)(,.2,1,0),(),(yxynyxhfyyyxhfyynnnnnnnn1 ny 由显式得到，称为由显式得到，称为预估值；
13、预估值；yn+1由隐式得到，称为由隐式得到，称为校正值。校正值。这种求解方法统称为这种求解方法统称为预估校正方法预估校正方法。其求解过。其求解过程为：程为：nnyyyyyyy2211027例3 用预估校正方法求解微分方程（取h=0.1）：1)0()10(,2yxyxyy解：1000.1)2(00001 yxyhyy0918.1)1000.11.021000.1(1.01)2(11101 yxyhyy1827.1)0918.11.020918.1(1.00918.1)2(11112 yxyhyy1763.1)1827.12.021827.1(1.00918.1)2(22212 yxyhyy232
14、222()1.2599xyyh yy 332332()1.2547xyyh yy 28（三）梯形公式（三）梯形公式 00)()(),(yxybxayxfdxdy 11,nnxnnxy xy xfx y xdx 111,2nnnnnnhyyfxyfxy()(,)y xf x yfxnx1nx29梯形公式局部截断误差梯形公式局部截断误差 *1111(),()()()22nnnnnnnnnhhyy xfxyfxyy xy xy x231()()()()()2!3!nnnnny xyxy xy xy x hhh21()()()()2!nnnnyxy xy xy x hh*111()nnnRy xy 3
15、3112nyxhh 30 111,2nnnnnnhyyfxyfxy预估-校正方法：称为改进的Euler求解公式或改进Euler法。11,2nnnnnnnnhyyfxyfxyhfxy111100(,)(,)(,),0,1,2,2()nnnnnnnnnnyyhf x yhyyf x yf xyny xy 31为了表示方便，可以改写为：112121h()2(,)(,)nnnnnnyyKKKf xyKf xh yhK 11,2nnnnnnnnhyyfxyfxyhfxy32（四）欧拉方法的收敛性分析（四）欧拉方法的收敛性分析由初值问题的单步法产生的近似解 ，如果对于任一固定的 均有 ，则称该方法是收敛的
16、。0nxxnh0lim()nnhnyy x ny定义：局部截断误差：局部截断误差：若初值问题的一个单步法的局部截断误差为：定理：11(),1,pnRO hp 且增量函数(,(),)nnxy xh 关于y满足Lipschitz条件，则整体截断误差：111()-()pnnney xyO h整体截断误差比局整体截断误差比局部截断误差低部截断误差低1 1阶阶*1111()-()-()(,(),)nnnnnnnRy xyy xy xhxy xh 33证明：存在常数c，使得*111()Pnny xych*1111()()PnnnRy xyO h*11()(,(),)(,)(1)()nnnnnnnnnnyy
17、y xyhxy xhxyhhL y xy*1()(,(),)nnnnyy xhxy xh 1(,)nnnnyyhxyh *111111111111211210111()()(1)(1)(1)1(1)(1)1(1)(1).(1)(1)(1)1(1)(1)1nnnnnnnPnPPnPnPnnnPney xyy xyyychhL echhL chhL echhLhLechhLhLhLhLehLchhLehL 03411110(1)1(1)(1)1nPnnhLechhLehL 000(),0y xye21011().2hLhLhLhLe1(1)0(1)nnhLhLe(1)11PnhLnceheL 当
18、固定时，1nxx 10(1)nnhxxba 所以()111Pb a Lpncehec hL 35（五）欧拉方法的稳定性分析（五）欧拉方法的稳定性分析问题：常微分方程初值问题数值解的每步计算都是在前一步计算的结果上进行的，所以必须考虑前面的误差对以后计算结果的影响，误差的积累会不会盖过真解呢？选用代表性试验方程：y=y (Re()0，在计算yn时引入了误差n。若这个误差在计算后面的yn+k(k=1,2,.)中所引的误差n+k按绝对值均不增加，就说这个数值方法对于这个步长h和复数 是绝对稳定的。若在区域R内数值方法是绝对稳定的，则R为该数值方法的绝对稳定区域(区间)。36把欧拉方法用于试验方程：y
19、=y1nnnyyhy 1nnnh 误差方程：11nnh 11h 要求误差不增加：O-2-1Re(h)Im(h)37把向后欧拉方法用于试验方程：y=y111nnh O21Re(h)Im(h)11h 要求误差不增加：可见隐式的向后Euler方法比显式的Euler方法的绝对稳定域要大得多。同阶精度的数值方法，往往隐式方法比显式方法的绝对稳定域大。38例4：以 y=y为例判断梯形公式的稳定性：11()2nnnnhyyyy解出：11/21/2nnhyyh 1/211/2hh 0 这种稳定性称为无条件稳定！39例5：用Euler法、向后法、向后Euler法、法、改进的欧拉法（梯形公式）解初值问题 83,(
20、12)(1)3yyxy 取步长h=0.2，小数点后至少保留4位。解：Euler法 1831.60.4nnnnyyhyy 向后Euler法 1180.625113nnnyyhyh 1183nnnyyhy40改进的Euler法13161371 nnyy111(,)(,)2nnnnnnhyyf xyf xy 383822.011 nnnnyyyy梯形公式法 11,2nnnnnnnnhyyfxyfxyhfxy 11624690.581.122nnnnnhyyhyhyy 410013xyh=0.23.00003.30083.46593.55653.60623.63351.01.21.41.61.83.0
21、真解真解改进的欧拉法梯形公式梯形公式向后向后EulerEulerxi 3.0000 3.0000 3.0000 3.0000 3.4000 3.2500 3.3077 3.2800 3.5600 3.4063 3.4734 3.4424 3.6240 3.5039 3.5626 3.5366 3.6496 3.5649 3.6106 3.5912 3.6598 3.6031 3.6365 3.622913161371 nnyy11.60.4nnyy 10.6251nnyy 10.581.12nnyy 如何得到高精度的求解公式？422.3 Runge-Kutta法法 如何构造更高精度的求解公式？
22、如何构造更高精度的求解公式？一种思路是：对微分方程右端积分采用高次多项式近似，如用二次多项式近似可得到 )(,()(,(4),(62)(,()()(112222nnnnnnxxnnxyxfxyxfyxfhdxxyxfxyxynn 43 ),(),(4),(311222nnnnnnnnyxfyxfyxfhyy 该公式称为Simpson公式。该公式的局部截断误差为：4(4)4()()90khRyO h 44 rnrnnnnyrhyhyhyy!221 Taylor展开：(,)nnyf xyf (,)(,)xnnynnnxyyfxyfxyyfff 222xxxyxyyyyfffffffyff 另一种得
23、到高阶方法的想法是直接利用泰勒级数展开。如果能计算得到y的高阶微商,则可写出r阶的计算公式45例题：取h=0.1，求解初值问题 2,(0)1;00.5yyyx 解：2324(4)35,22,66,2424yyyy yyyyyyyyyy2231(2)2nnnnhyyhyy 23423451(2)(6)(24)2!3!4!nnnnnnhhhyyhyyyy 一阶公式：二阶公式：四阶公式：21(1)nnnnnyyhyyhy 461()1y xx 真解：p0.10.20.30.40.511.100001.221001.370081.557791.8004621.110001.246891.421741.
24、662621.9208741.111101.249661.428481.666451.99942y(xn)1.111111.250021.428571.666672.00000 xn问题：1、求微商麻烦；2、计算量大。47R-K方法是通过对不同点上的函数值做线性组合，构造近似公式，把近似公式和泰勒展开相比较，使前面的若干项相吻合，从而使近似公式达到一定的阶数。问题：如何组合函数值？NiiinnKcyy11),(1nnyxhfK -11(,)iininijjjKhf xa h yb K其中，48选择参数ci,ai,bij 的原则是，要求的Runge-Kutta法的右端项在(xn,yn)处泰勒展开
25、后按h的幂次重新整理得到的结果231123nnyyd hd hd h 与微分方程的解y(xn+1)在xn处的Taylor展开式 2311.2!3!nnnnny xhy xhyxh yxh yx 有尽可能多的项重合。49以计算两个函数值为例说明 11122122211,nnnnnnyyc hKc hKKfxyKfxa h yb hK 选择c1、c2、a2、b21*111R()nnny xy阶最高 1,nnnKfxy xyx 2221122211,nnnnxyKfxa h y xb hKfxy xa hfb hK fO h*1ny 为了计算50*1122321222()nnnxyyy xcchyx
26、bc a hff fO ha于是 )(!21321hOxyhxyhxyxynnnn 若要求局部截断误差达到)(3hO,则要求有 12222121,1/2,/1ccc aba 122211,12ccab 选取51122211,12ccab 1121211122,nnnnnnyyhKhKKfxyKfxh yhK 二阶R-K方法的公式 它的局部阶段误差为O(h3)。这是计算两次函数值的情况下所能达到的最高阶。11122122211,nnnnnnyyc hKc hKKfxyKfxa h yb hK 521222101,1/2ccab 取取，12121,11,22nnnnnnyyhKKfxyKfxh y
27、hK 中间点法：12221132,443ccab 取取，1121211344,22,33nnnnnnyyhKhKKfxyKfxh yhK 二阶休恩(Heun)法：53经典的R-K方法是一个四阶的方法,公式为 112341213243226,/2,/2/2,/2,nnnnnnnnnnhyyKKKKKfxyKfxhyhKKfxhyhKKfxh yhK 1、一步法，可以自开始；特点：2、精度较高；3、便于改变步长；4、计算量较小。54例2：用二阶R-K方法和四阶R-K方法(h=0.1)求解解：,101fx yxy000,0,0.1,0.5axyhb 1121212101101nnnnnnhyyKKK
28、xyKxhyhK 二阶R-K方法 0)0()10()1(10yxyxy 11234121324322/610110/21/210/21/2101nnnnnnnnnnyyh KKKKKxyKxhyhKKxhyhKKxhyhK 四阶R-K方法55xn二阶二阶R-K误差四阶四阶R-K误差y=1-e-5x20.10.20.30.40.50.60.70.80.91.0 0.0500 0.0012 0.0488 0.0000 0.0488 0.1830 0.0017 0.1813 0.0000 0.1813 0.3627 0.0004 0.3624 0.0000 0.3624 0.5475 -0.0031
29、 0.5507 -0.0000 0.5507 0.7059 -0.0076 0.7134 -0.0001 0.7135 0.8235 -0.0112 0.8346 -0.0001 0.8347 0.9012 -0.0125 0.9135 -0.0002 0.9137 0.9476 -0.0116 0.9590 -0.0002 0.9592 0.9733 -0.0093 0.9823 -0.0002 0.9826 0.9866 -0.0066 0.9931 -0.0002 0.9933说明：1、高阶方法比低阶方法精度高，但需要y(x)光滑性好；2、步长h并非越小越好；步长选择从大到小。56数学家
30、Butcher于1965年证明了运算量与可以达到的最高精度阶数关系如下：每步计算每步计算kR的个数的个数2345678910可以达到最可以达到最高精度阶数高精度阶数23445666R-2可见，通常只选用可见，通常只选用R 4的公式是合适的。的公式是合适的。57自动选步长方法自动选步长方法 设在xn处之前误差不计：)(nnxyy 误差：(/2)(/2)()11111()15hhhnnnny xyyy ()511()hnnny xyC h(/2)511()2(/2)hnnny xyCh(/2)(/2)()11111()()15hhhnnnny xyyy 58隐式隐式R-KR-K方法方法类似于显式类
31、似于显式R-K公式，稍加改变，就得到公式，稍加改变，就得到隐式隐式R-K方法方法。11Lnniiiyyhk1,Lininii jjjkfxc h yc ha k显式公式中对系数求和的上限是显式公式中对系数求和的上限是i-1。而在隐式公式中。而在隐式公式中对系数求和的上限是对系数求和的上限是L，需要用迭代法求出，需要用迭代法求出Ki。推导。推导隐式公式的思路和方法与显式隐式公式的思路和方法与显式R-KR-K法类似。通常，同法类似。通常，同级的隐式公式可以获得比显式公式更高的阶。级的隐式公式可以获得比显式公式更高的阶。59通常，同级的隐式公式获得比显式公式更高的阶。常通常，同级的隐式公式获得比显式
32、公式更高的阶。常用的隐式用的隐式R-K法有法有:2,2111nnnnnyyhxhfyy1级级2阶中点公式阶中点公式:),(),(2111nnnnnnyxfyxfhyy2级级2阶梯形公式：阶梯形公式：n+1n121nn122nn12hy=y+(k+k)23+3h3+23k=fx+h,y+k+hk64123-33-23hk=fx+h,y+hk+k6124 2级级4阶阶R-K公式：公式：60隐式龙格隐式龙格-库塔法库塔法 ),.,1().,(.11111mjhKhKyhxfKKKhyymmjjijijmmii )2,2(1111KhyhxfKhKyyiiii其中其中2阶方法阶方法 的绝对稳定区域为的
33、绝对稳定区域为0ReImg而而显式显式 1 4 阶方法的绝对阶方法的绝对稳定区域为稳定区域为k=1k=2k=3k=4-1-2-3-123ReImg无条件稳定无条件稳定龙格龙格-库塔法稳定区域库塔法稳定区域612.4 线性多步法线性多步法 单步法的优缺点：简单，可以自开始；提高精度时需要增加中间函数值；没有充分利用前几步得到的信息。多步法：多步法：yn-p,yn-p+1,yn-1 ,yn yn+162yn-p,yn-p+1,yn-1,yn yn+1的方法的方法考虑如下形式的求解公式11101kknin iin iiiyyhf此计算公式为线性的，所以称为线性多步法。(,),(0,1,1)ijjif
34、f x yik(1,0,1)iik当 时公式含有 ，这时公式是隐式的，而当 时公式是显式的011nf1063常微分初值问题与积分公式等价11()()(,()nnxnnxy xy xf x y xdx 基本思想：基本思想：对f(x,y(x)进行多项式插值，利用插值公式计算右边积分，可以得到常微分初值问题求解公式。64三次多项式插值求积分3121311231213112316121216()()()()()()()()()()()()()()()()()()nnnnnnnnnnnnnnnnF tLttxtxtxFxhtxtxtxFxhtxtxtxFxhtxtxtxFxh插值余项为插值余项为(4)3
35、11221()()()()()(),4!nnnnnnFRttxtxtxtxxx 令令 ，假设已知，假设已知 ，取插值节，取插值节点点 ，对函数，对函数 作三次作三次Lagrange多项式插值多项式插值()(,()F xf x y x21,nnnyyy211,nnnnxxx x()F x651133()()()()nnxnnxy xy xLxRx dx 122110,nnxnjnnn in ixijn injj ixxyyfxydxxx 四阶Adams内插公式)5199(242111 nnnnnnffffhyy内插公式局部截断误差：15(5)*32119(),720nnxnnxRx dxh yx
36、x 661133()()()()nnxnnxy xy xLxRx dx 令令 ，假设已知，假设已知 ，取插值节，取插值节点点 ，对函数，对函数 作三次作三次Lagrange多项式插值多项式插值()(,()F xf x y x321,nnnnyyyy321,nnnnxxxx()F x1123(5559379)24nnnnnnhyyffff外插公式局部截断误差：15(5)331251(),720nnxnnxRx dxh yxx 133100,nnxnjnnn in ixijn injj ixxyyfxydxxx 67)5199(242111 nnnnnnffffhyy内插公式1123(555937
37、9)24nnnnnnhyyffff外插公式通常计算公式为（预估校正系统）：)519),(9(24)9375955(242111113211nnnnnnnnnnnnnnfffyxfhyyffffhyy特点：精度高，但须与同精度的单步法配合使用。68例3：用Adams预测校正系统求初值问题 3,00.0500dIEIIdtLLLtI 其中，E=200,L=3,=100,=50解：由于Adams预测校正系统为四阶公式，选用四阶R-K公式计算开始值。a=t0=0;b=0.05I0=0;h=0.002f(t,I)=(200-50I3-100I)/3 121324311234,/2,/2/2,/2,22/
38、6nnnnnnnnnnKhf tIKhf thIKKhf thIKKhf th IKIIKKKK 691230.1289680.2493910.361332III再利用Adams预测校正系统1123111112(5559379)24(9(,)195)24nnnnnnnnnnnnnnhIIffffhIIfxIfff n45625t0.0080.0100.0120.050I0.4646910.5592060.644933.1.164861702.5 常微分方程组、高阶微分方程及常微分方程组、高阶微分方程及边值问题的数值解边值问题的数值解 考虑如下含二个未知函数的方程组：一、常微分方程组的数值解一、
39、常微分方程组的数值解 1112101022122020,(),()yfx yyyxyyfx yyyxy 00,yfx yy xy 则：可依照常微分方程数值方法来求解。令：1011120202212,yyfx yyyyfx yyyfx yy710000,yf x y zy xyzg x y zz xz 则其改进的Euler格式具有预估形式11,nnnnnnnnnnyyhf xyzzzhg xyz校正公式为11111111,2,6nnnnnnnnnnnnnnnnhyyf xyzf xyzhzzg xyzg xyz如对于方程组如对于方程组72 1112111222121111222211221121
40、1321111222231122,/2,/2,/2/2,/2,/2/2,/2,/2/2,nnnnnnnnnnnnnnnnnnnKfxyyKfxyKfxyyKfxhyhKyhKKfxhyhKyhKKfxhyhKyhKKfxhyhK 24221111222211221,1133433311112,1122224322124/2,/2,226nnnnnnnnnnnyhKKfxh yhKyhKKfxh yhKyhKyyKKKKhyyKKKK 两方程微分方程组的四阶R-K公式73二、高阶微分方程的数值解二、高阶微分方程的数值解 对于高阶微分方程，可以把它化成微分方程组后再进行求解。如 0000,yfx
41、y yy xyyxy 令：112yyyy 则：1210010212200020 ,yyyxyyyfx yyyxyxy 使用四阶R-K公式可得求解公式。74四阶R-K公式 2112222111222221112222211111121122322334334222,/2/2,/2,/2/2/2,/2,/2,nnnnnnnnnnnnnnnnyKfxyyKyhKKfxhyhKyhKKyhKKfxhyhKyhKKyhKKfxh yhKyhKKy 31,1111112,1222222214314226nnnnyKKKKhyyKKKK 75三、边值问题的数值解三、边值问题的数值解 最简单的边值问题二阶微分
42、方程及其定解条件:()(,(),(),(),()yxf x y xy xaxby ay b 它的解法可以用数值微分公式代替其导数值，将其变为代数方程，然后求解，这种求解方法通常称为差分法。222()()2()2()()kkkkkkky xhy xhyxO hhy xhy xy xhyxO hh 76 1111202(,)(1,2,1)2,kkkkkkknyyyyyf x yknhhy ayy by 其中，其中，h=(b-a)/n，xk=x0+kh，k=1,2,n-1，x0=a如果f是非线性函数，那么差分方程也是非线性的。差分方程为：上述差分方程为n-1个方程组成的n-1元代数方程组。解出y1,
43、y2,yn-1即为微分方程的计算解。77特别的，对二阶线性常微分方程()()()yp x yq x yr x 其差分方程形式为 1111222kkkkkkkkkyyyyypq yrhh )1(.)1(.2112.211222122212121211212222222221212nhnnhnnnnhnnhhprhrhrhprhyyyyqhpphqhqhppqh整理得 78边值问题的定解条件分类：第一类边界条件(),()y ay b第二类边界条件(),()y ay b第三类边界条件01010000()(),()()0,0,0y ay ay by b79例4 求解下面的边值问题，取h=0.2。213
44、630,100112xyxyyxxyyy 解：建立差分方程 211112110213()63212()kkkkkkkkyyyyyxxyxhhybyhay 将k=0代入(a)，联立(b)得：102(2.520.5)2yyc 80将k=1，2，3，4代入(a)，联立(c)得：012342.5220000.5212.07550.9623000.0679 012.03330.988900.0200120.91550.016900011.97671.7724yyyyy 解线性方程组得：012341.01321.01671.06931.21881.5130yyyyy 81小 结n求微分方程（组）数值解的单步法，多步法，但还有许多其它方法。n病态问题（刚性问题）求解比较困难。一般采用隐式方法进行求解。n实际计算中，方法的选取非常重要。f(x,y)不太复杂时，采用单步法（四阶R-K方法）；但当y(x)的光滑性不好时，采用低阶单步法（二阶R-K方法）；f(x,y)比较复杂时，采用多步法。
