常微分方程数值解法模板

第6章常微分方程数值解法《计算方法》第6章常微分方程数值解法§1引言§2欧拉法和改进的欧拉法§3龙格-库塔法§4阿当姆斯方法第6章常微分方程数值解法《计算方法》§1引言在高等数学里我们已经接触过常微分方程，对于一些典型的常微分方程，有求解析解的基本方法，但多数情况下，遇到的问题比较复杂，此时，只能利用近似方法求解，一般有两种近似方法。近似解析方法数值方法第6章常微分方程数值解法《计算方法》实际求解的常微分方程，大多是定解问题┉┉满足指定条件的特解初值问题边值问题本章讨论常微分方程，数值解的最简单问题┉┉一阶方程初值问题，即函数f(x)满足下列微分方程和初值条件：在几何问题是(6-1)表现为一簇曲线，称(6-1)的积分曲线，初值问题(6-1)(6-2)就是要求一条过(x0,y0)的积分曲线00(,)()dyfxydxyxy(6―1)(6－2)第6章常微分方程数值解法《计算方法》方程的精确解y(x)称为积分曲线。方程是否有解，解是否唯一？第6章常微分方程数值解法《计算方法》定理1对初值问题(6-1)(6-2)，若f(x,y)在区域G={a≤x≤b,|y|＜∞}内连续，且关于y满足李普希兹条件，即存在常数L，使|f(x,y1)-f(x,y2)|≤L|y1-y2|（6-3）对G中任意两个y1,y2均成立，其中L是与x,y无关的常数，则初值问题(6-1)(6-2)在(a,b)内存在唯一解，且解是连续可微的。第6章常微分方程数值解法《计算方法》设f(x，y)在带形区域R：｛a≤x≤b，-∞＜y＜+∞｝上为x，y的连续函数，且对任意的y满足李普希茨(Libusize)条件｜f(x，y1)-f(x，y2)｜≤L｜y1-y2｜其中(x，y1)、(x，y2)∈R，L为正常数。在求初值问题(6-1)(6-2)的数值解时，我们通常采用离散化方法，求在自变量x的离散点a=x0＜x1＜x2＜…＜xn=b第6章常微分方程数值解法《计算方法》上的准确解y(x)的近似值y0，y1，y2，…，yn常取离散点x0，x1，x2，…，xn为等距，即xi+1-xi=h，i=0，1，2，…，n-1h称为步长。图6.1表示为初值问题(6―1)(6―2)在n+1个离散点上的准确解y(x)的近似值。第6章常微分方程数值解法《计算方法》图6.1第6章常微分方程数值解法《计算方法》数值解法的重点不在于求准确解（即解析解），而是直接求一系列点上的近似解。求解过程顺着节点排列的顺序一步步向前推进，也即按递推公式由y0，y1…..yi推yi+1，下面各种方法的实质是建立递推公式。初值问题（6.1）(6.2)的数值解法的基本特点是：第6章常微分方程数值解法《计算方法》§2欧拉法和改进的欧拉法一、欧拉方法欧拉法是最早的一种数值方法，由于它不够精度，很少被采用。但是它在某种程度上反映了数值方法的基本思想。从几何上看，初值问题（6.1）（6.2）的解y(x)代表一条过点(x0,y0)的曲线，且曲线上任一点(x,y)的切线斜率为f(x,y)。现在我们从(x0,y0)点出发，作解曲线y(x)的切线，由于这点的切线斜率为f(x0,y0)，故切线方程为y=y0+f(x0,y0)(x-x0)设在x0附近，切线可以作为曲线的近似。令x=x1,我们就得到y(x1)的近似值y1y1=y0+f(x0,y0)(x-x0)=y0+f(x0,y0)h第6章常微分方程数值解法《计算方法》§2欧拉法和改进的欧拉法再从(x1,y1)点出发，以f(x1,y1)为斜率，作直线方程为y=y1+f(x1,y1)(x-x1)令x=x2,我们就得到y(x2)的近似值y2y2=y0+f(x1,y1)h以此类推，一般的y(xn+1)的近似值yn+1为y(xn+1)=yn+hf(xn,yn)n=0,1,……(6.4)称为解初值问题的欧拉方法。其几何意义是用一条折线近似代替积分曲线y=y(x),如图6-1所示，因此欧拉法又称为折线法。第6章常微分方程数值解法《计算方法》Ox0x1x2xn(x0,y0)(x1,y1)(x2,y2)(xn,yn)第6章常微分方程数值解法《计算方法》y(x)过点P0(x0,y0)，从P0出发以f(x0,y0)为斜率做一直线与直线x=x1交于点p1(x1,y1),显然有:y1=y0+hf(x0,y0),再从p1出发,以f(x1,y1)为斜率做一直线推进到x=x2上一点p2(x2,y2),依此类推,这样得到解曲线的一条近似曲线,它就是折线p0p1p2……所以欧拉方法又叫欧拉折线法第6章常微分方程数值解法《计算方法》欧拉法是用yi通过yi+1=yi+hf(xi,yi)i=0,1,……求yi+1,这样利用y0┅y1┅y2┅……计算yi+1用前一步的yi┅┅单步法计算yi+1用前几步的｛yi-n｝┉┉多步法第6章常微分方程数值解法《计算方法》例1：用欧拉法求解方程1)0(2'yxyy0≤x≤1.2h=0.2解:欧拉法的具体形式为:yi+1=yi+hf(xi,yi)=(1-0.4xi)yi所以:y1=y0+hf(x0,y0)y2=y1+hf(x1,y1)=(1-0.4x1)y1=0.920000……=(1-0.4x0)y0=1第6章常微分方程数值解法《计算方法》xiyiy(xi)0110.210.9607890.40.9200000.8521440.60.7728000.6976760.80.5873220.5277921.00.3993830.3678791.20.2396300.236938可见欧拉法的精度是很差的所求值用下表列出,并与精确值对比第6章常微分方程数值解法《计算方法》二、欧拉方法的误差分析定义1(p146)对于初值问题，当假设yi是准确的时，用某种方法求yi+1时所产生的截断误差称为该方法的局部截断误差。第6章常微分方程数值解法《计算方法》我们来看在第i+1步使用欧拉方法所得yi+1的局部截断误差y(xi+1)-yi+1假定yi是准确的,即yi=y(xi)由y(xi+1)=y(xi+h)，应用泰勒展开y(xi+1)=y(xi+h)=y(xi)+hy′(xi)+y″(ξ)/2*h2而由欧拉公式算出yi+1=yi+hf(xi,yi)两式相减得y(xi+1)-yi+1=(h2/2)*y″(ζ)=O(h2)即欧拉方法所得yi+1的局部截断误差为O(h2)注意:只计算了一步，事实上每一步都有可能产生误差，有时误差会原来越大，有时又会得到很好的控制，因此还要考虑整体截断误差。第6章常微分方程数值解法《计算方法》定义2(p147)设yi是用某种方法计算初值问题(6-1)(6-2)在xi点的近似解，而y(xi)是它的精确解，则称为该方法的整体截断误差，也称为该方法的精度。补：若某方法的局部截断误差为O(hp+1)，则该方法的精度为p阶的。iiiyxy)(欧拉方法的整体截断误差是O(h)，一阶的第6章常微分方程数值解法《计算方法》三、改进的欧拉法欧拉法虽然形式简单，计算方便，但比较粗糙，精度也低。特别当y=y(x)的曲线曲率较大时，欧拉法的效果更差。为了构造较高精度的数值解法，对初值问题再做分析00)(),('yxyyxfy第6章常微分方程数值解法《计算方法》对y′=f(x，y)等式两边在（xi，xi+1）上取积分11))(,('iiiixxxxdxxyxfdxy1))(,()()(1iixxiidxxyxfxyxyxixi+1(xi,yi)(xi+1,yi+1)第6章常微分方程数值解法《计算方法》在上图也可利用),(1iiiiyxhfyy矩形公式精度不高，再次说明欧拉精度低后退欧拉公式右矩形公式),(111iiiiyxhfyy这样就利用数值积分公式计算y(xi+1)的近似值。如果用左矩形计算右面的积分：)8.6())(,()()(11iixxiidxxyxfxyxy显式欧拉公式),(1iiiiyxhfyy第6章常微分方程数值解法《计算方法》为了得到更精确的方法我们可使用梯形公式)],(),([2),(111iiiixxyxfyxfhdxyxfii,...2,1,0)],(),([2111iyxfyxfhyyiiiiii此时(6-9)改进的欧拉公式，又称为梯形公式或隐式欧拉公式xixi+1(xi,yi)(xi+1,yi+1)这样得到的点列仍为一折线，只是用平均斜率来代替原来一点处的斜率。第6章常微分方程数值解法《计算方法》第6章常微分方程数值解法《计算方法》不难发现:欧拉公式yi+1=yi+hf(xi,yi)是关于yi+1的显式，只要已知yi，经一次计算可立即得到yi+1的值；而改进欧拉公式,...2,1,0)],(),([2111iyxfyxfhyyiiiiii中的yi+1以隐式给出，且yi+1含在函数f(xi+1，yi+1)中，所以梯形法是隐式单步法，第6章常微分方程数值解法《计算方法》一般来说，这是一个非线性方程（除非f对y是线性的），可用我们前面讲过的非线性方程的各种方法求解，比如用迭代法第6章常微分方程数值解法《计算方法》)],(),([2),()(11)1(1)0(1kiiiiikiiiiiyxfyxfhyyyxhfyy可证明其收敛，局部截断误差O(h3)k=0,1,2,…具体做法是：先用欧拉公式(6―4)求出一个y(0)i+1作为初始近似值，然后再用改进的欧拉公式(6―9)进行迭代，即直到满足(1)()11(1)11kkiikiiyyyy取再转到下一步计算第6章常微分方程数值解法《计算方法》四、预报-校正方法我们看到梯形法虽然提高了精度，但其算法复杂，每算一点，都需进行反复迭代，为了控制计算量，通常只迭代一两次就转入下一步计算，以简化算法。具体地说，我们先用欧拉公式求一个初步的近似值),(1iiiiyxhfyy称为预测值，其精度可能很差，再用梯形公式将其校正一次1iy)],(),([2111iiiiiiyxfyxfhyy第6章常微分方程数值解法《计算方法》预报-校正公式)],(),([2),(1111iiiiiiiiiiyxfyxfhyyyxhfyy(6-11)在实际计算时，还常常将式(6―11)写成下列形式：(6-12),....2,1,0)(2),(),(*2*111*2*1ikkhyyhkyhxfkyxfkiiiiii第6章常微分方程数值解法《计算方法》例3用预报-校正公式求解初值问题2(0)1,01,0.1xyyyyxh解：由预报-校正公式有)]2()2[(2]2[11111iiiiiiiiiiiiiyxyyxyhyyyxyhyyh=0.1,x0=0,x1=0.1,x2=0.2,…..,x9=0.9,x10=1第6章常微分方程数值解法《计算方法》1.1)101(1.01)2(00001yxyhyyi=0095909.1)]2()2[(211100001yxyyxyhyyi=1191818.1)2(11112yxyhyy184097.1)]2()2[(222211112yxyyxyhyyi=2……最后得下表其解析解为21yx第6章常微分方程数值解法《计算方法》表6―2第6章常微分方程数值解法《计算方法》预报-校正公式的截断误差为O(h3)（证明，p150)预报-校正公式的整体截断误差为O(h2)预报-校正公式是单步显式第6章常微分方程数值解法《计算方法》练习1:用Euler法解1)0()6.00(2yxxyydxdy取h=0.2计算到x=0.6第6章常微分方程数值解法《计算方法》练习2:用改进Euler法（梯形公式）解2)1()21(38yxydxdy取h=0.2小数点后至少保留5位第6章常微分方程数值解