单步法

第五章常微分方程的数值方法(NumericalMethodsforOrdinaryDifferentialEquations)§1初值问题的数值方法一阶常微分方程初值问题的一般形式是:称f(x,y)在区域D上对y满足Lipschitz条件是指:0)(),,(yaybxayxfyD={(x,y)|a≤x≤b,c≤y≤d}使得:0L|||),(),(|2121yyyxfyxf],[,],,[21dcyybax一、解的存在性定理若f(x,y)在区域D上连续,且对y满足Lipschitz条件,则初值问题在[a,b]上存在唯一的连续可微解y.利用Picard逼近容易证明:利用Gronwall不等式易证解连续依赖于初值条件:设f(x,y)在区域D上连续,且对y满足Lipschitz条件,若y(x,s)是初值问题saybxayxfy)(),,(的解,则有:|||),(),(|21)(21ssesxysxyaxL定理1定理2定理的意义在于：若初值问题中的初始值有一微小扰动，则解的扰动也是微小的，也就是解连续依赖于初始条件．通常将具有这种特性的初值问题称为是适定的．（稳定的）二、数值解和精确解用数值方法求解初值问题，不是求出它的解析解或其近似解析式，而是给出它的解在某些离散节点上的近似值,称为数值解用y(x)表示问题的准确解y0,y1,…,yN表示数值解，即问题的解y(x)在相应节点处的近似值y(x0),y(x1),…,y(xN)表示解y(x)在节点x0,x1,…,xN处的准确值定义1三、单步法和多步法单步法在计算yi+1时只利用yi在计算yi+1时要用到yi,yi−1,…,yi−k+1k步法显式计算公式可写成:yk+1=yk+hΦf(xk,yk;h）yk+1=yk+hΦf(xk,yk,yk+1;h)它每步求解yk+1需要解一个隐式方程隐式格式定义2多步法yi−1,yi−2,…在计算yi+1时不仅利用yi,还要利用四、基本方法在公式推导或误差估计中主要有三种方法：(1)差商代替导数法：(2)积数值分法：11)()()(ixixiidxxyxyxy1))(,(ixixdxxyxf)()()(11iiiiixyxxxyxy))(,(iixyxf(3)泰勒级数法：)()()()(21hOxyhxyxyiii))(,()(iiixyxhfxy),(yxfybxxxxaNN110§2Euler方法Euler方法是一种最简单的单步法bxxxxaNN110NjNabhjhxxj,,2,1,,0一、Euler公式(1)Euler公式001(),(,),0,1,,1iiiiyyxyyhfxyiN用差商近似导数hxyxyxy)()()(010x0x1),()()()(000001yxfhyxyhxyxy继续这一过程，得到1000(,)yyhfxy记为从而得到求解初值问题的公式)1,...,0(),(1niyxfhyyiiii(2)局部截断误差在假设yi=y(xi)，即第i步计算是精确的前提下，考虑的截断误差Ri=y(xi+1)yi+1称为局部截断误差/*localtruncationerror*/。假定“yi=y(xi)”称为局部化假定对于数值方法),,(1hyxhyyiiii局部截断误差定义为)]),(,()([)(11hxyxhxyxyeiiiii定义3若某算法的局部截断误差为O(hp+1)，则称该算法有p阶精度。估计局部截断误差的主要方法是Taylor展开法Euler方法的局部截断误差1321),(!31)(21)()()(iiiiiiiixxyhxyhxyhxyxy)()(21))(,()(32hOxyhxyxhfxyiiii111)(iiiyxye)],([)]()(21)()([32iiiiiiyxhfyhOxyhxyhxy)()(2132hOxyhi欧拉法具有1阶精度。定义4二、后退的Euler公式以yi代替y(xi)从而得后退Euler公式用差商近似导数hxyxyxyiii)()()(11)()()(11iiixyhxyxy))(,()(11iiixyxfhxy),(111iiiiyxfhyy),()(11100iiiiyxfhyyxyy111)(iiiyxye)()(2132hOxyhi局部截断误差具有1阶精度三、梯形公式两式相加得:称为梯形公式.Euler公式:),(1iiiiyxfhyy向后Euler公式:),(111iiiiyxfhyy)],(),([2111iiiiiiyxfyxfhyy所以)],(),([2)(11100iiiiiiyxfyxfhyyxyy111)(iiiyxye)()(2343hOxyhi局部截断误差具有2阶精度四、改进的Euler公式得到)],(),([2),(1111iiiiiiiiiiyxfyxfhyyyxhfyyEuler公式是显示公式.后退的Euler公式和梯形公式是隐式公式在利用隐式公式进行计算yn+1时,使用迭代法,通常由显示的Euler公式给出迭代初值.称Euler公式给出的迭代初值为预测.称梯形公式迭代一次的值为校正.改进的Euler公式的局部截断误差),(yxfyfffyffyyxyx))(,()()(1iiiixyxhfxyxy)())](,())(,())(,([2132hOxyxfxyxfxyxfhiiiiyiix111)(iiiyxye)))]}(,()(,())(,([2)({)(11iiiiiiiixyxhfxyxfxyxfhxyxy假设yi=y(xi)),(1iiiiyxhfyy)))(,()(,()))(,()(,(1iiiiiiiixyxhfxyhxfxyxhfxyxf)())(,())(,())(,())(,(2hOxyxfxyxhfxyxhfxyxfiiyiiiixii)()(3111hOyxyeiii所以:具有2阶精度欧拉公式和改进欧拉公式分iiipipicicpiyyyyyyyyyyyy9.0)(1.091.0)(1.0905.0211及两种格式的计算结果分别列表如下别为：000,1xy解：取步长h=0.1iiiyyyyi9.0)(1.01ydxdyy1)0(]1,0[x并比较两法所得计算结果的精度。例：试分别用函数公式和改进的欧拉公式求解：由表可见，与精确解相比，改进的Euler公式的精度较Euler公式有明显的提高。下面再看两步Euler公式（9），除了给出初值外，还需要借助其它单步法（如Euler公式，后退Euler公式及梯形公式等）再提供一个xeyoyixEuler公式改进的Euler公式精确解iyiyiiyxyiiyxyixy01110.10.20.310.90000000.81000000.72900000.34867843108374.43107308.821018182.121092010.10.90500000.81902500.74121760.36854100.90483740.81873080.74081820.3678794410626.1410942.2410994.3410616.6§3Runge-Kutta方法考虑一阶常微分方程初值问题0)(),,(yaybxayxfy将区域[a,b]进行分划：NjNabhjhxxaxj,,2,1,,,00h称为步长11))(,()()(ixixiidxxyxfxyxy))(,()(yhfxyi))(,()(hcxyhcxfbhxyjijiji11njjbn个点的加权平均平均高度一、基本思想jjmjmca11))(,(111hcxyhcxfkii))(,(222hcxyhcxfkii))(,(hcxyhcxfkninin……))(,(hcxyhcxfkjijij)))()(,(ijijixyhcxyhcxf))))(,()(,(iijijixyxfhcxyhcxf)),(11jmmjmijikahyhcxf),(),(11111jmmjmijijiinjjjiikahyhcxfkyxfkkbhyyn级显式Runge-Kutta方法jjmjmca11j=2,3,…,n二、二级Runge-Kutta方法),(),()(122122111hkcyhcxfkyxfkkbkbhyyiijiiii特殊地:取n=2记),(iiyxff),(iixxyxff),(iiyyyxff),(iixxxxyxff),(iiyyyyyxff),(iixyxyyxff由此得fk1),(1222hkcyhcxfkii))()(),(212hOfkfhcyxfyxii)()()(31222211hOfkfhbchfbbyyyxii另一方面fyxfyii),(fffyyx)()(21)()()(321hOxyhxyhxyxyiiii)()(2132hOfffhhfyyxi为使局部截断误差为O(h3),应取2112221bcbb特殊地1、改进的Euler方法1,21221cbb1,,1,0Ni),(),()2121(121211hkyhxfkyxfkkkhyyiiiiii2、中点方法21,1,0221cbb21hkyyii),(1iiyxfk)2,2(12khyhxfkii1,,1,0Ni3、二阶Heun方法32,43,41221cbb)32,32(),()4341(121211khyhxfkyxfkkkhyyiiiiii1,,1,0Ni二级Runge-Kutta方法不超过二阶记fffFyx22fffffGyyxyxx则)(21322222hOGhchFcfk)(21)(43222222211hOGhcbFhcbhfbbyyii)()(!3)(!2)()()(4321hOxyhxyhxyhxyxyiiiii)()(62432hOFfGhFhhfyyi因此局部截断误差只能达到O(h3)2112221bcbb当GcbFfGy22221)(61取n=3三、三级Runge-Kutta方法3231323213133122213322111),(),(),()(aachkahkayhcxfkhkcyhcxfkyxfkkbkbkbhyyiiiiiiii记fffFyx22fffffGyyxyxx)(21322222hOGhchFcfk232131kaka)(213232223223hOGhachFacfc),(23213133hkahkayhcxfkii])([2321313yxfkakafchfxyxxfk