计算方法第六章

}第6章常微分方程数值解法6.1欧拉法和改进的欧拉法6.2龙格-库塔法6.3线性多步法引言}一阶常微分方程的初值问题00)(],[),,(yxybaxyxfy当f(x,y)连续,且满足Lipschitz(李普希兹)条件|f(x,y1)–f(x,y2)|≤L|y1–y2|,x∈[a,b],y1,y2∈R则上述初值问题存在连续可微的解函数y=y(x).例如初值问题2)0(1yeyx可求出方程y′=1+ex的通解为y=x+ex+c,将初值条件x=0,y=2代入得2=1+c,故c=1,所以初值问题的解为y=x+ex+1dxeyx)1(引言}又如初值问题0)0(21yxyy可求出它的解为xtxdteey022但要进一步计算指定点的函数值，还需要用数值积分方法。有些微分方程的解是隐函数，例如要求函数值还需要解超越方程。应用中所处理的微分方程往往很复杂且大都得不出一般解，所以一般用数值解法。1xeyyx数值解法：给定节点a=x0x1…xN=b,将初值问题离散化为差分方程,求出解函数y(x)在这些点的近似值y1,y2,…,yn。所求得的近似值称为数值解。}6.1欧拉法和改进的欧拉法6.1.1欧拉法及其截断误差6.1.2改进的欧拉法及预测-校正公式}6.1.1欧拉法及其截断误差1.公式的构造思想：用差商代替导数nxxxx,,,,210设等距，步长为1,,1,0,1nihxxiiyxfhxyxyhxyhxyhxyhxyxy,,)()()(令x=xi,x+h=xi+1,y(xi)≈yi，y(xi+1)≈yi+1,初值问题离散化为00)(),(yxyyxfy初值问题)(,2,1,0,),(001xyyiyxfhyyiiii(欧拉公式)}例6.1取步长h=0.1，用欧拉法求解初值问题)(,2,1,0,),(001xyyiyxfhyyiiii(欧拉公式)解1,0,),(00yxyxyxf1)0(yyxyy1=y0+hf(x0,y0)=1+0.1(0+1)=1.1y2=y1+hf(x1,y1)=1.1+0.1(0.1+1.1)=1.22y3=y2+hf(x2,y2)=1.22+0.1(0.2+1.22)=1.22+0.142=1.362………y10=y9+hf(x9,y9)=y9+0.1(x9+y9)=3.18748}2.几何意义——用折线代替曲线计算解函数的近似值。准确解为xeyx12}3.数值公式的误差来源。(1)局部截断误差（简称截断误差）：假设yi=y(xi)是准确的，计算yi+1所产生的误差y(xi+1)-yi+1若局部截断误差可以表示为O(hk+1),k为正整数，则称公式是k阶公式。（2）由于实际上yi不是准确值，因此它的误差会传播下去。（3）实际计算时，每一步都可能产生舍入误差。}4.欧拉公式的截断误差是O(h2)，公式是1阶的。因为)()(),(1iiiiiixyhxyyxfhyy21)(21)()()(hyhxyxyxyiii（泰勒公式）两式相减，由设yi=y(xi)，有221121hOhyyxyii}6.1.2改进的欧拉法及预测-校正公式对微分方程y′=f(x,y)两边求xi到xi+1的定积分，有1))(,()()(1iixxiidttytfxyxy利用梯形公式计算积分，有))](,())(,([2)()(1111iiiiiiiixyxfxyxfxxxyxy将y(xi)、y(xi+1)分别用yi、yi+1代替，构造相应的数值公式：（改进的欧拉公式）)(,2,1,0,)],(),([200111xyyiyxfyxfhyyiiiiii1.公式的构造32,,,21,hOyxfyxfyxfhyxfhyiiiiyiixiii}2.截断误差假设yi=y(xi),解函数在x=xi处的泰勒公式为32121hOhxyhxyxyxyiiiihyi1)],(),([2111iiiiiiyxfyxfhyy在改进的欧拉公式中，]),(),([2hhxfyxfhyhiiii设则有求出在h=0处的泰勒公式，整理后得00,,021,032hOyxfyxfhyxfhiiyiixiiiy0上式h和h2项的乘数应为零，于是,,0iiixyyxfiiiiiyiixxyyxfyxfyxf,,,0}32,,,2,hOyxfyxfyxfhyxfhyiiiiyiixiii111iiixyhxyy320200hOhh3hO因而改进的欧拉法是二阶的。}3改进的欧拉法的具体使用格式。改进的欧拉法是隐式公式,计算时常用迭代法。一般每一步先由欧拉公式计算出yi+1的初始值yi+1(0),再迭代计算yi+1。,2,1,0)],(),([2),()(11)1(1)0(1kyxfyxfhyyyxfhyykiiiiikiiiii当满足||)(1)1(1kikiyy时，取.)1(11kiiyy可证明当f(x,y)满足一定条件时，迭代是收敛的。}改进的欧拉法的预测－校正公式,2,1,0)],(),([2),()(11)(11)(1)(1iyyyxfyxfhyyyxfhyyciipiiiiiciiiipi可证明预测－校正公式的截断误差也为O(h3)。}例取步长h=0.2,用改进的欧拉法的预测－校正公式求解初值问题的数值解y1,y2.1)0(yyxy解1,0,),(00yxyxyxf)]()[(1.02.12.0)(2.0)(111)(1piiiiiiiiiiipiyxyxyyyxyxyy2.12.12.000)(1yxyp24.1)2.12.010(1.01)]()[(1.0)(110001pyxyxyy528.124.12.12.02.02.12.011)(2yxyp24.11y5768.12y预测-校正公式具体是5768.1)528.14.024.12.0(1.024.1)]()[(1.0)(221112pyxyxyy}}6.2龙格-库塔法6.2.1二阶龙格-库塔（Runge-Kutta)公式6.2.2四阶龙格-库塔公式引言与其他算法的实例比较}思想：从泰勒公式出发，寻找更高阶的数值公式。例如，泰勒公式计算到二阶可得)()(!21)()()(32hOhxyhxyxyhxy))(,())(,())(,())(,()())(,()(xyxfxyxfxyxfdxxyxdfxyxyxfxyyx令nxx则1nxhx，略去余项，得出一个二阶的]),(),(),([2),(21nnnnynnxnnnnyxfyxfyxfhhyxfyy理论上按此方式可以得到更高阶的公式。但需要计算复合函数的高阶导数，使算法复杂而不实用。因数值公式为}龙格和库塔的思想（间接地运用泰勒公式）：利用y(x)在若干个点上的函数值和导数值，作出一个适当的线性组合，使这个线性组合按h展开后的泰勒公式与y(x+h)的泰勒公式有较多的项达到一致，从而得出较高阶的数值公式。}6.2.1二阶龙格-库塔（Runge-Kutta)公式设想一个有二阶精度的数值公式形状为a,b为待定系数。))(,(1ahxyahxbhfyynnnn在h=0处求泰勒公式得仍令x=xn，则x+h=xn+1。如果能找出a,b,使得)()()()(3hOhaxyhbxyhxy那么略去余项就可得到上面所希望的近似计算公式了。因此考虑)()()(haxyhbxyhxyhT)(2)0()0()0(32hOhThTThT由于0)0(T})()()()(2ahxyabahxybhaahxyabhxyhT)()()(ahxyabhahxybhxyhT由T(h)的泰勒公式)()()21()()1(32hOhxybahxybhT为使T(h)=O(h3),令02101abb，解出1,21ba，得32,2hOhxyhxfhxyhxyhT)()()(haxyhbxyhxyhT)()()0(xybxyT)()21()0(xyabT)(2)0()0()0(32hOhThTThT整理得}利用2222,2hOxyhxyhxyhxyhxf)()))(,(2)(,2()()(3hOxyxfhxyhxfhxyhxy可以推出取x=xn并略去O(h3)便得到二阶龙格－库塔公式)),(2,2(1nnnnnnyxfhyhxfhyy或21121)2,2(),,(khyykhyhxfkyxfknnnnnn}6.2.2四阶龙格-库塔公式仿照上述的讨论，可导出四阶龙格－库塔公式:)22(6),(),2,2()2,2(),,(432113423121kkkkhyyhkyhxfkkhyhxfkkhyhxfkyxfknnnnnnnnnn})10(1)0(,2xyyxyy例取步长h=0.2，用四阶龙格－库塔公式求下面初值问题的数值解。解2.0,1,0,2),(00hyxyxyyxf由公式得91818.0)1.1,1.0()2,2(1002fhkyhxfk90864.0)09182.1,1.0()2,2(2003fhkyhxfk1),(001yxfk84324.0)18173.1,2.0(),(3004fhkyhxfk18323.16)22(2.0432101kkkkyy}12xy数值解yi与准确解y(xi)的对照见表准确解是yxi()00.20.40.60.81.011.183231.341671.483281.612511.7321411.183221.341641.483241.612451.73205xiyiy(xi)}yxi()}6.3线性多步法6.3.1四阶阿达姆斯(Adams)外插公式6.3.2四阶阿达姆斯(Adams)内插公式多步法的概念6.3.3初始出发值的计算6.3.4阿达姆斯预测-校正公式与其他算法的实例比较}单步法----计算yn+1时只使用yn的值。多步法----计算yn+1时使用前面的k个yi值，即由yn-k+1,yn-k+2,…,yn-1,yn计算yn+1。(k=1,2,……)线性多步法----计算yn+1的公式由yn-k+1,yn-k+2,…,yn-1,yn的线性组合表达。多步法的概念}6.3.1四阶阿达姆斯(Adams)外插公式设想用yn-3,yn-2,yn-1,yn的值计算yn+1。为方便，讨论由)(),(),2(),3(xyhxyhxyhxy出发计算y(x+h),由初值问题的方程y′=f(x,y(x))两边从x到x+h积分,