第六章-常微分方程初值问题的数值解法

1计算方法湖南大学电气与信息工程学院第六章常微分方程初值问题的数值解法计算方法课程组200(,)()yfxyyxy微分方程初始条件定理：若f(x,y)在某闭区域R：00||,||(0,0)xxayybab上连续，且在R域内满足李普希兹(Lipschitz)条件，即存在正数L，使得对于R域内的任意两值y1,y2，下列不等式成立：1221|(,)(,)|||fxyfxyLyy则上述初值问题的连续可微的解y(x)存在并且唯一。一阶常微分方程的初值问题:引言3§6.1基本离散方法§6.2Runge-Kutta方法§6.3线性多步法§6.4收敛性与稳定性§6常微分方程数值解法§6.5微分方程组的数值处理§6.6高阶微分方程的数值处理4考虑一阶常微分方程的初值问题:例如：,[,]()0101y=x+yxy其解析解为：12xyxe§6.1基本离散方法0)(],[),(yaybaxyxfdxdy5[,]xdtx01012-tye但是,只有一些特殊类型的微分方程问题能够得到用解析表达式表示的函数解，而大量的微分方程问题很难得到其解析解。因此，只能依赖于数值方法去获得微分方程的数值解。例如：,[,]()01012-xy=exy其解析解为：很难得到其解析解6要计算出解函数y(x)在一系列节点a=x0x1…xn=b处的近似值),...,1()(nixyyii通常取节点间距为步长，通常采用等距节点，即取hi=h(常数)。它适合计算机求解，应用广泛，具有应用价值。)1,...,0(1nixxhiii7Taylor级数复习函数在点作Taylor级数展开：()yyx0xx这里x,x0都可以是任意一点。00000(n)n001()()()()()()2!()()2yx=yx+yxx-x+yxx-x1++yxx-x+n!8Taylor级数复习()yyx0xx()()()()()()()!2i+1iiiii2n!nnhhyxyx+hyxhyxyxyx()()()()()()()!2i-1ii()2nniiih-hyxyx-hyxhyxyxyxn!则：()()()()()!2i+1iiii()2n(n+1)ihhyxyxhyxhyxyxyxn!()()()()()!2i-1iiii()()2n(n+1)ih-hyxyx-hyxhyxyxyxn!(n)n00000001()()()()()()()()2!21yx=yx+yxx-x+yxx-x++yxx-x+n!96.1.1Euler公式方法一：利用等距分割，数值微分来代替导数项，建立差分格式。ibaxim-=1、向前差商公式（单步法或显式欧拉法）1()()''()'2()nnnnyyxyxhyxh1()()'(,())'()2nnnnnyxyxhyxhfyx)(''2))(,()()(21nnnnnyhxyxhfxyxy所以，可以构造差分方程),(1nnnnyxhfyy称为局部截断误差。显然，这个误差在逐步计算过程中会传播，积累。因此还要估计这种由截断误差带来的积累误差10定义:在假设yi=y(xi)，即第i步计算是精确的前提下，考虑的截断误差Ri=y(xi+1)yi+1称为局部截断误差。定义:若某算法的局部截断误差为O(hp+1)，则称该算法有p阶精度。收敛性：)(''2))(,()()(21nnnnnyhxyxhfxyxy),(1nnnnyxhfyy考察局部误差的传播和积累112、向后差商公式（单步法或隐式欧拉法）11()()'()''()2nnnnyxyxhyxyh111()()(,())''()2nnnnnyxyxhfxyxyh2111()()(,())''()2nnnnnhyxyxhfxyxy),(111nnnnyxhfyy是隐格式，要迭代求解)0(1)(11)1(1),(nknnnknyyxhfyy可以由向前差商公式求出(0)1(,)nnnnyyhfxy123、中心差商公式（二步法或中点欧拉公式）311()()'()'''()23!nnnnyxyxhyxyh112(,)nnnnyyhfxy是多步法，具有2次代数精确度，但该迭代格式不稳定13对微分方程],[,),('baxyxfdxdyy两边积分有:11(,)nnnnxxxxdydxfxydxdx)(31hOen类似，可以算出其误差估计式：所以，有111[(,)(,)]2nnnnnnhyyfxyfxy是个隐式的方法，要用迭代法求解具有2次代数精确度11()()(,())nnxnnxyxyxfxyxdx1112()()[(,())(,())]'''()122nnnnnnhyxyxfxyxffxxhy局部截断误差4、梯形公式（单步法或梯形欧拉法）14141111[,()]d()[,(d()()[,()])]()nnnnxxnnnxnnnxnnfxyxxyxyxyxhfxyxxxfxyx左矩公式以近似值yn代替精确值y(xn)可得显式欧拉公式：100(,)()nnnnyyhfxyyyx将微分方程yf(x,y)在区间[xn,xn+1]上积分：方法二：数值积分法1、显式欧拉法15欧拉公式有明显的几何意义依此类推得到一折线由于每步计算只用到，故也称此为单步法，ny或显格式的单步法。000011122(,)()(,)()(,)(,)xyyxxyyxxxxyxy过点的曲线是解在作的切线与直线交于再作切线交于1ny16也称欧拉折线法.就是用这条折线近似地代替曲线欧拉方法()yx()yyxxynx1nxnp1np1npx1717这样便得到了隐式欧拉法：11100(,)()nnnnyyhfxyyyx含有未知的函数值隐式欧拉法没有显式欧拉法方便11111111d()()[,()]d()[,(),()]][nnnnxnnxxnnnnxnnfxyxxxxfxyxyxyxyxhfxyx（右矩公式）2、隐式欧拉法将微分方程yf(x,y)在区间[xn,xn+1]上积分：1811111111[,()]d()d()()2[,(,()])[]nnnnxxxnnxnnnnnnyxyxyxfxyxxxxfhxxyxyfx（中矩公式）以y(x)在xn1,xn上的近似值代替精确值可得：11002(,)()nnnnyyhfxyyyx需要前两步的计算结果将微分方程yf(x,y)在区间[xn-1,xn+1]上积分：3、二步欧拉法（中心欧拉法）19用近似值代替精确值可得梯形公式欧拉法：111[(,)(,)]2nnnnnnhyyfxyfxy注：上式右端出现了未知项，可见梯形法是隐式欧拉法的一种；实际上，梯形公式欧拉法是显式欧拉法与隐式欧拉法的算术平均。11111111[,()]d[,d()[,()][][,()](2(]2),)nnnnxxnnnnnxxnnnnnfxyxxfxyxfyxhfxyxfxyxxxyxx(梯形公式)4、梯形公式欧拉法将微分方程yf(x,y)在区间[xn,xn+1]上积分：201(0)1yyxy取h0.1，计算y(0.5)并与精确解进行比较解：由已知条件可得：h0.1，x00，y01，f(x,y)yx1显式欧拉法：1(1)(1)0.90.10.1nnnnnnnnyyhyxhyhxhyx例题：用显式欧拉法、隐式欧拉法、梯形法求解初值问题：()exyxx其精确解为：21隐式欧拉法：111(1)nnnnyyhyx化简得：11(1)(1)(1)nnnnnhyyhxyhxh11[(1)](0.10.11)1.1(1)nnnnnyyhxhyxh111[(1)(1)]2nnnnnnhyyyxyx梯形公式欧拉法：11[1(2)][1(2)](2)(2)nnnnhyhyhxx11[(2)(2)](2)(1.90.20.21)2.1nnnnnnyhyhxxhyx22xn显式法yn隐式法yn梯形法yn精确解y(xn)0.011110.11.0000001.0090911.0047621.0048370.21.0100001.0264461.0185941.0197310.31.0290001.0513151.0406331.0408180.41.0561001.0830141.0700971.0703200.51.0904901.1209221.1062781.10653110.90.10.1nnnyyx1(0.10.11)1.1nnnyyx1(1.90.20.21)2.1nnnyyx显式法隐式法梯形法计算结果()exyxx其精确解为：23称上述误差Rn1为该常微分方程数值算法的局部截断误差如果某个常微分方程数值算法的局部截断误差可表示为O(hp1)，则称该数值算法的精度是p阶欧拉法的精度为一阶；二步欧拉法的精度为二阶；梯形公式欧拉法的精度为二阶。6.1.2局部截断误差为了简化分析某常微分方程数值算法的误差，现假设yn=y(xn)，即在前一步yn准确的前提下，估计111()nnnRyxy2421()()()()()2!nnnnnhyxyxhyxhyxyx如果只保留线性项，忽略h2及以上各项，则：11()()()nnnnyxyxhyxy()[,()]nnnyxfxyxny1(,)nnnnyyhfxy显式欧拉公式泰勒展开法如果初值问题中的f(x,y)充分可微，则可将y(xn1)在点xn处展开：25将y(xn1)在xn点处用泰勒公式展开：211()()()()()2!(,)nnnnnnyyxyxhyxhyxhxx22111()()()2!nnnyRyxyhOh显式欧拉法的局部截断误差：1()[,()]()()nnnnnnyyxhfxyxyxhyx1(,)nnnnyyhfxy欧拉法()nnyyx令：1阶精度局部截断误差的分析利用泰勒公式展开，比较各算法与展开式的前几项26011100(,)(,)fxyfxy1001010010()[()]()[()]xyxxfxxxyyfyyy0000(,)(,)fxhykfxy00()()xyhfxhkfyk000000(,)(,)()()xyfxhykfxyhfxhkfyk0000(,)(,)fxyhkfxhykxy补充：二元函数微分中值定理271111111()(,)[,](,)[()]nnnynnnnyfxyfxxxfxyy111()[()()(,)]nnnynnyyxhyxOhfxR21()()()()nnnyxyxhyxOhy(xn1)在xn点处展开：211111()(,)()nnnynnRyxyhfxROh111(,)nnnnyyhfxy隐式欧拉法：1阶精度()nyx2111()1(,)nynROhhfx11[(),]nnyxy111()(,)nynn