【文档说明】《计算机数值方法教学课件》第二章-常微分方程数值解法.ppt，共(81)页，1.805 MB，由小橙橙上传

转载请保留链接：https://www.ichengzhen.cn/view-77143.html

以下为本文档部分文字说明：

第二章常微分方程数值解法Chapter2NumericalSolutionofOrdinaryDifferentialEquation(s)2§2.1引言Fm假设牵引力F为恒定值dvFdtm为了确定待定常数

，可以给定初始条件：或者给定边界条件：00tvvendtendvv,dvmFtvdt假设牵引力不恒定呢？求速度)(tv3虽然求解微分方程有许多解析方法，但解析方法只能够求解一些特殊类型的方程。还有一类近似方法称为数值方法，它可以给出解在

一些离散点上的近似值。利用计算机解微分方程主要使用数值方法。解析方法与数值方法4主要研究对象：初值问题00)()(),(yxybxayxfdxdy00()(,())xxyxyfxyxdx求数值解：求y(x)在离散数据点xk处的近似值yk

。y=y(x)xyx0=ax1x2x3xk-1xn-1xn=bxk5则称f(x,y)对y满足李普希兹条件，此时初值问题在[a,b]上存在唯一的连续可微的解。定理1：设f(x,y)是定义在区域G={(x,y)|a≤x≤b,y∈R}上的连续函数，若存在正的常数L使：1212|(,)(,)|

||fxyfxyLyy12[,],xabyy使得对任意的及都成立，(Lipschitz)条件有解条件：6则Lipschitz条件成立：在f(x,y)对y可微的情况下，若偏导数有界：(,),(,)fxyLxyGy有解条件的判断：*121212(,)|(,)

(,)|()||fxyfxyfxyyyLyyy7定理2：如果f(x,y)在G={(x,y)|a≤x≤b,y∈R}上满足Lipschitz条件，则初值问题是适定的。适定性：指初值问题中，初始值y0及微分方程的右端函数f(x,y)有微小变化时，只能引起

解的微小变化。适定性条件：8在解的存在区间[a,b]上取n+1个节点bxxxxan210这里把1,...iiihxxi=0,1,,n称为由xi到xi+1的步长一般取成等间距的：nabh求解方法：步进法（分为

单步法和多步法）数值方法的基本思想9本章规定：在处初值问题的理论解用表示，数值解法的近似解用表示。nx()nyxny记，它和是不同的，后者等于。(,)nnnffxy(,())nnfxyx()nyx10§2.2几种简单的数值方法（一）欧拉（Euler）法

00)()(),(yxybxayxfdxdy001()(,)0,1,2,...nnnnyyxyyhfxyn11①、泰勒公式解释②、求导的两点公式解释③、积分公式解释欧拉公式的的分析解释001()

(,)nnnnyxyyyhfxy00(,)()yfxyyxy121)1()(21)!1()(!)(!2)()()()(pppnpnnnnhpyhPxyhxyhxyxyxy

2()(,())()nnnyxhfxyxOh泰勒公式解释其中：1nnhxx可以得到：001()(,)nnnnyxyyyhfxy13求导的两点公式解释1()()()()(,)nnnnn

nyxyxyxhyxfxy可以得到：001()(,)nnnnyxyyyhfxy1400(,)(1.1)()(1.2)yfxyaxbyxy对微分方程(1.1)两端从1nnxx到进行积分11(,)nnn

nxxxxydxfxydx11()()(,)nnxnnxyxyxfxydx积分公式解释1511)(,())((,())nnnnnnnnxxfxyxyyfxyxh右端积分用左矩形数值求积公式：11()()(,)nnxnnxy

xyxfxydx001()(,)nnnnyxyyyhfxy即：(,)fxyfxnx1nx16欧拉公式的的几何描述yxx0x1x2x3x4y=y(x),...2,1,0),()(100nyxhfyyyxynnnn00(,)()yfxyy

xy00,xy11,xy17例题1：（取步长h=0.1）用Euler方法求满足条件的y(t)数值解。22,12(1.0)0.0tdyytetdtty)2(21ntnnnnnetythyy解：27182

81828.0)2(1.00200001tetytyy122111120.1()0.684755578tyyytet18ntnyny(tn)y(tn)-yn01234101.01.11.21.31.43.00.00.271830.684761

.276983.0935515.398240.00.345920.866641.607223.6203618.6830.00.074090.181880.330240.526813.28486数值解列表为)(2eetyt22,12(1.0)0.0tdyytetdtt

y)2(21ntnnnnnetythyy19欧拉方法的误差估计()-nnneyxy通过数值方法进行计算时，考虑每一步产生的误差，从x0开始一步步累积到xn，称为该数值方法在xn点处的整体

截断误差，该误差与xn及之前的各步计算误差都有关系。20欧拉方法的误差估计*1*1111()(,(),)()-()-()(,(),)nnnnnnnnnnnyyxhxyxhRyxyyxyxhxyxh

为了简化分析，着重分析xn点单步计算产生的误差，即把xn点之前的计算当作无误差：称该误差为数值方法在xn+1点处的局部截断误差。12111(,())()ppnnnRHxyxhOh

局部截断误差的第一个非零项为局部截断误差主项。21欧拉方法的误差估计如果求解公式的局部截断误差为R(h)=O(hp+1)，则称该求解公式具有p阶精度，称该方法为p阶方法。定义：欧拉方法：212()()()()2!()(,())()nnnnn

nnyxyxyxyxhhyxhfxyxOh*1()(,())nnnnyyxhfxyx*2111()-()nnnRyxyOh具有1阶精度。22（二）向后欧拉法hxyxyxynnn)()()(1111100(,),0,1,2,(),nnnnyyhf

xynyxy（1）方法其公式为：()(,)yxfxyxfnx1nx23（2）局部截断误差*232111()()()()2nnnnyxRyxyhOhOh231()()()

()()2nnnnyxyxyxyxhhOh*11()()nnnyyxyxh11100(,),0,1,2,(),nnnnyyhfxynyxy21()()()()nnnyxyxyxhOh24例题2：用向后

Euler法解初值问题1)0()10(,2yxyxyy11112()nnnnnxyyhyy向后Euler法的公式为解：x0=0,y0=1,取h=0.125方法比较及推广：Euler方法显式公式向后Euler方法隐式公式解一个非线性方程难求解显

式和隐式相结合隐式的显化01110()(,),nnnnxyyyhfyxy100((),)nnnnxyyyhfyyx26计算公式为：00

1111)(,....2,1,0),,(),(yxynyxhfyyyxhfyynnnnnnnn1ny由显式得到，称为预估值；yn+1由隐式得到，称为校正值。这种求解方法统称为预估－校正方法。其求解过程为：nnyyyyyyy22

11027例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

yxyhyy232222()1.2599xyyhyy332332()1.2547xyyhyy28（三）梯形公式00)()(),(yxybxayxfdxdy11,nnxnnxyxy

xfxyxdx111,,2nnnnnnhyyfxyfxy()(,)yxfxyfxnx1nx29梯形公式局部截断误差*1111(),,()()()22nnnnnnnnnhhyyxfxyfxyyxyxyx

231()()()()()2!3!nnnnnyxyxyxyxyxhhh21()()()()2!nnnnyxyxyxyxhh*111()nnnRyxy

33112nyxhh30111,,2nnnnnnhyyfxyfxy预估-校正方法：称为改进的Euler求解公式或改进Euler法。11,,,2nnnnnn

nnhyyfxyfxyhfxy111100(,)(,)(,),0,1,2,2()nnnnnnnnnnyyhfxyhyyfxyfxynyxy31为了表示方便，可以改写为：112121h()2(,)(,)nnnnnnyyK

KKfxyKfxhyhK11,,,2nnnnnnnnhyyfxyfxyhfxy32（四）欧拉方法的收敛性分析由初值问题的单步法产生的近似解，如果对于任一固定的均有，则称该方法是收敛的。0nxxn

h0lim()nnhnyyxny定义：局部截断误差：若初值问题的一个单步法的局部截断误差为：定理：11(),1,pnROhp且增量函数(,(),)nnxyxh关于y满足Lipschitz条件，则整体截断误差：111()-()pnnneyxyOh

整体截断误差比局部截断误差低1阶*1111()-()-()(,(),)nnnnnnnRyxyyxyxhxyxh33证明：存在常数c，使得*111()Pnnyxych*1111()()PnnnRyxyOh

*11()(,(),)(,,)(1)()nnnnnnnnnnyyyxyhxyxhxyhhLyxy*1()(,(),)nnnnyyxhxyxh1(,,)nnnnyyhxyh

**111111111111211210111()()(1)(1)[(1)][1(1)](1)[1(1)(1)...(1)](1)(1)1(1)(1)1nnnnnnnPnPPnPnPnnnPneyxyyxyyychhLechhLchhLech

hLhLechhLhLhLhLehLchhLehL03411110(1)1(1)(1)1nPnnhLechhLehL

000(),0yxye21011()...2hLhLhLhLe1(1)0(1)nnhLhLe(1)1[1]PnhLnceheL当固定时，1nxx10(1)nnhxxba所以()1

1[1]PbaLpncehechL35（五）欧拉方法的稳定性分析问题：常微分方程初值问题数值解的每步计算都是在前一步计算的结果上进行的，所以必须考虑前面的误差对以后计算结果的影响，误差的积累会不会盖过真解呢？选用代表性试验方程：y'=y(Re()<0)对给定

步长h>0，在计算yn时引入了误差n。若这个误差在计算后面的yn+k(k=1,2,....)中所引的误差n+k按绝对值均不增加，就说这个数值方法对于这个步长h和复数是绝对稳定的。若在区域R内数值方法是绝对稳定的，则R为该数值方法的绝对稳定区域(区间)。36把欧拉方法用于试验方

程：y'=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,(12)(1)3yyxy取步长h=0.2，小数点后至少保

留4位。解：Euler法1831.60.4nnnnyyhyy向后Euler法1180.625113nnnyyhyh1183nnnyyhy40改进的Euler法

13161371nnyy111[(,)(,)]2nnnnnnhyyfxyfxy]3838[22.011nnnnyyyy梯形公式法11,,,2nnnnnnnnhyyfxyfxyhfxy

11624690.581.122nnnnnhyyhyhyy410013xyh=0.23.00003.30083.46593.55653.60623.63351.01.21.41.61.83.0真解改进的欧拉法梯形公式向后EulerEulerxi

3.00003.00003.00003.00003.40003.25003.30773.28003.56003.40633.47343.44243.62403.50393.56263.53663.64963.56493.61063.59123.65983.603

13.63653.622913161371nnyy11.60.4nnyy10.6251nnyy10.581.12nnyy如何得到高精度的求解公式？42§2.3Runge-Kutta

法如何构造更高精度的求解公式？一种思路是：对微分方程右端积分采用高次多项式近似，如用二次多项式近似可得到))(,())(,(4),(62))(,()()(112222nnnnnnxxnnxyxfxyxfyxfhdxxyxfxyxynn

43),(),(4),(311222nnnnnnnnyxfyxfyxfhyy该公式称为Simpson公式。该公式的局部截断误差为：4(4)4()()90khRyOh44rnrnnnnyrhyhyhyy!221

Taylor展开：(,)nnyfxyf(,)(,)xnnynnnxyyfxyfxyyfff222xxxyxyyyyfffffffyff另一种得到高阶方法的想法是直接利用泰勒级数展开。如

果能计算得到y的高阶微商,则可写出r阶的计算公式45例题：取h=0.1，求解初值问题2,(0)1;00.5yyyx解：2324(4)35,22,66,2424yyyyyyyyyyyyyy2231(2)2nnnn

hyyhyy23423451(2)(6)(24)2!3!4!nnnnnnhhhyyhyyyy一阶公式：二阶公式：四阶公式：21(1)nnnnnyyhyyhy461()1yxx真解：p0.10.20.30.40.511.100001.221

001.370081.557791.8004621.110001.246891.421741.662621.9208741.111101.249661.428481.666451.99942y(xn)1.111111.250

021.428571.666672.00000xn问题：1、求微商麻烦；2、计算量大。47R-K方法是通过对不同点上的函数值做线性组合，构造近似公式，把近似公式和泰勒展开相比较，使前面的若干项相吻合，从而使近似公式达到一定的阶数。问题：如何组合函数值？NiiinnKcyy1

1),(1nnyxhfK-11(,)iininijjjKhfxahybK其中，48选择参数ci,ai,bij的原则是，要求的Runge-Kutta法的右端项在(xn,yn)处泰勒展开后按h的幂次重新整理得到的结果231123nnyydhdhdh与微分方程

的解y(xn+1)在xn处的Taylor展开式2311...2!3!nnnnnyxhyxhyxhyxhyx有尽可能多的项重合。49以计算两个函数值为例说明11122122211,,nnnnnn

yychKchKKfxyKfxahybhK选择c1、c2、a2、b21*111R()nnnyxy阶最高1,nnnKfxyxyx222112

2211,,nnnnxyKfxahyxbhKfxyxahfbhKfOh*1ny为了计算50*1122321222()nnnxyyyxcchyxbcahfffOha

于是)(!21321hOxyhxyhxyxynnnn若要求局部截断误差达到)(3hO,则要求有12222121,1/2,/1cccaba122211,12ccab选取51122211,12ccab

1121211122,,nnnnnnyyhKhKKfxyKfxhyhK二阶R-K方法的公式它的局部阶段误差为O(h3)。这是计算两次函数值的情况下所能达到的最高阶。

11122122211,,nnnnnnyychKchKKfxyKfxahybhK521222101,1/2ccab取，12121,11,22nnnnnnyyhKKfxyKfxhyhK

中间点法：12221132,443ccab取，1121211344,22,33nnnnnnyyhKhKKfxyKfxhyhK二阶休恩(Heun)法：53经典

的R-K方法是一个四阶的方法,公式为112341213243226,/2,/2/2,/2,nnnnnnnnnnhyyKKKKKfxyKfxhyhKKfxhyhKKfxhyhK1、一步法，可以自开始；

特点：2、精度较高；3、便于改变步长；4、计算量较小。54例2：用二阶R-K方法和四阶R-K方法(h=0.1)求解解：,101fxyxy000,0,0.1,0.5axyhb1121212101101nnnnnnhyy

KKKxyKxhyhK二阶R-K方法0)0()10()1(10yxyxy11234121324322/610110/21/210/21/2101nn

nnnnnnnnyyhKKKKKxyKxhyhKKxhyhKKxhyhK四阶R-K方法55xn二阶R-K误差四阶R-K误差y=1-e-5x20.10.20.30.40.50.60.70.80.91.

00.05000.00120.04880.00000.04880.18300.00170.18130.00000.18130.36270.00040.36240.00000.36240.5475-0.00310.5507-0.00000.

55070.7059-0.00760.7134-0.00010.71350.8235-0.01120.8346-0.00010.83470.9012-0.01250.9135-0.00020.91370.9476-0.01160.9590

-0.00020.95920.9733-0.00930.9823-0.00020.98260.9866-0.00660.9931-0.00020.9933说明：1、高阶方法比低阶方法精度高，但需要y(x)光滑性好；2、步长h并非越小越好；步长选择从大到小。56数学家Butch

er于1965年证明了运算量与可以达到的最高精度阶数关系如下：每步计算kR的个数2345678910可以达到最高精度阶数23445666≤R-2可见，通常只选用R≤4的公式是合适的。57自动选步长方法设在xn处之前误差不计：)(nnxyy误差：(/2)(/2)()11111()15hhhnnn

nyxyyy()511()hnnnyxyCh(/2)511()2(/2)hnnnyxyCh(/2)(/2)()11111()()15hhhnnnnyxyyy58隐式R-K方法类似于显式R-K公式，稍加改变，就得到隐式R-K

方法。11Lnniiiyyhk1,Lininiijjjkfxchychak显式公式中对系数求和的上限是i-1。而在隐式公式中对系数求和的上限是L，需要用迭代法求出Ki。推导隐式公式的思路和方法与显式R-K法类似。通常，

同级的隐式公式可以获得比显式公式更高的阶。59通常，同级的隐式公式获得比显式公式更高的阶。常用的隐式R-K法有:2,2111nnnnnyyhxhfyy1级2阶中点公式:)],(),([2111nnnnnnyxfyxfhyy2级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阶方法的绝对稳定区域为0

ReImg而显式1~4阶方法的绝对稳定区域为k=1k=2k=3k=4-1-2-3---123ReImg无条件稳定龙格-库塔法稳定区域61§2.4线性多步法单步法的优缺点：简单，可以自开始；提高精度时需要增加中间函数值；没有

充分利用前几步得到的信息。多步法：yn-p,yn-p+1,…,yn-1,ynyn+162yn-p,yn-p+1,…,yn-1,ynyn+1的方法考虑如下形式的求解公式11101kkniniiniiiyyhf

此计算公式为线性的，所以称为线性多步法。(,),(0,1,,1)ijjiffxyik(1,0,,1)iik当时公式含有，这时公式是隐式的，而当时公式是显式的．011nf1063常微分初值问题与积分公式等价11()()(,()

)nnxnnxyxyxfxyxdx基本思想：对f(x,y(x))进行多项式插值，利用插值公式计算右边积分，可以得到常微分初值问题求解公式。64三次多项式插值求积分3121311231213112316121216()()()()()

()()()()()()()()()()()()()nnnnnnnnnnnnnnnnFtLttxtxtxFxhtxtxtxFxhtxtxtxFxhtxtxtxFxh插值余项为(4)311221()()()()()(),4!nn

nnnnFRttxtxtxtxxx令，假设已知，取插值节点，对函数作三次Lagrange多项式插值()(,())Fxfxyx21,,nnnyyy211,,,nnnnxxxx()Fx651133()()()()nnxnnxyxyxLxR

xdx122110,nnxnjnnninixijninjjixxyyfxydxxx四阶Adams内插公式)5199(242111

nnnnnnffffhyy内插公式局部截断误差：15(5)**32119(),720nnxnnxRxdxhyxx661133()()()()nnxnnxyxyxLxRxdx令

，假设已知，取插值节点，对函数作三次Lagrange多项式插值()(,())Fxfxyx321,,,nnnnyyyy321,,,nnnnxxxx()Fx1123(5559379)24nnnn

nnhyyffff外插公式局部截断误差：15(5)331251(),720nnxnnxRxdxhyxx133100,nnxnjnnninixijninjjixxyyfxydxxx

67)5199(242111nnnnnnffffhyy内插公式1123(5559379)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/6nnnnnnnnnnKhftIKhfthIKKhfthIKKhfthIKI

IKKKK691230.1289680.2493910.361332III再利用Adams预测—校正系统1123111112(5559379)24(9(,)195)24nnnnnnnnnnnnnn

hIIffffhIIfxIfffn456…25t0.0080.0100.012…0.050I0.4646910.5592060.644933….1.16486170§2.5常微分方程组、高阶微分方程及

边值问题的数值解考虑如下含二个未知函数的方程组：一、常微分方程组的数值解1112101022122020,,(),,()yfxyyyxyyfxyyyxy00,yfxyyxy则：可依照常微分方程数值方法来求解。令：

1011120202212,,,,,yyfxyyyyfxyyyfxyy710000,,,,,,yfxyzyxyzgxyzzxz则其改进的Euler

格式具有预估形式11,,,,nnnnnnnnnnyyhfxyzzzhgxyz校正公式为11111111,,,,2,,,,6nnnnnnnnnnnnnnnnhyyfxyzfxyzhzzgxyzgxyz

如对于方程组7211121112221211112222112211211321111222231122,,,,,/2,/2,/2/2,/2,/2/2,/2,/2/2,nnnnnnnnnnnnnnnnnnnKfxyyKfxyKfxyyKfx

hyhKyhKKfxhyhKyhKKfxhyhKyhKKfxhyhK24221111222211221,1133433311112,112222432212

4/2,/2,,,,226nnnnnnnnnnnyhKKfxhyhKyhKKfxhyhKyhKyyKKKKhyyKKKK

两方程微分方程组的四阶R-K公式73二、高阶微分方程的数值解对于高阶微分方程，可以把它化成微分方程组后再进行求解。如0000,,yfxyyyxyyxy令：112yyyy则：1210010212200020,,yyyxyy

yfxyyyxyxy使用四阶R-K公式可得求解公式。74四阶R-K公式2112222111222221112222211111121122322334334222,,/2/2,/2,/2/2/2,/2,/2,,nnnnnnnn

nnnnnnnnyKfxyyKyhKKfxhyhKyhKKyhKKfxhyhKyhKKyhKKfxhyhKyhKKy

31,1111112,1222222214314226nnnnyKKKKhyyKKKK75三、边值问题的数值解最简单的边值问题—二阶微分方程及

其定解条件:()(,(),()),(),()yxfxyxyxaxbyayb它的解法可以用数值微分公式代替其导数值，将其变为代数方程，然后求解，这种求解方法通常称为差分法。222()()2()2()()kkkkkkkyx

hyxhyxOhhyxhyxyxhyxOhh761111202(,,)(1,2,,1)2,kkkkkkknyyyyyfxyknhhyayyby其中，h=(b-a)/n，xk=x0+kh，k=1,2,…,n-1，x0=

a如果f是非线性函数，那么差分方程也是非线性的。差分方程为：上述差分方程为n-1个方程组成的n-1元代数方程组。解出y1,y2,…,yn-1即为微分方程的计算解。77特别的，对二阶线性常微分方程()()()ypxyqxyrx其差分方程形式为1111222kkk

kkkkkkyyyyypqyrhh)1(...)1(...2112......

............211222122212121211212222222221212nhnnhnnnnhnnhhprhrhrhprhyyyyqhpphqhqhppqh整理得78边值问题的定解条件分类：第一类边界条件(),()yayb第二类边界条件(

),()yayb第三类边界条件01010000()(),()()0,0,0yayaybyb79例4求解下面的边值问题，取h=0.2。213630,100112xyxyyxxyyy

解：建立差分方程211112110213()63212()kkkkkkkkyyyyyxxyxhhybyhay将k=0代入(a)，联立(b)得：102(2.520.5)2y

yc80将k=1，2，3，4代入(a)，联立(c)得：012342.5220000.5212.07550.9623000.0679012.03330.988900.0200120.91550.016900011

.97671.7724yyyyy解线性方程组得：012341.01321.01671.06931.21

881.5130yyyyy81小结求微分方程（组）数值解的单步法，多步法，但还有许多其它方法。病态问题（刚性问题）求解比较困难。一般采用隐式方法进行求解。实际计算中，方法的选取非常重要。f(x,y)不太复杂时，采用单步法（四阶R-K方法）；但当y(x)的光滑

性不好时，采用低阶单步法（二阶R-K方法）；f(x,y)比较复杂时，采用多步法。