7常微分方程数值解法

第七章常微分方程初值问题数值解法数值分析15:24:53NumericalAnalysis2本章内容欧拉法欧拉公式两步欧拉公式梯形法改进欧拉法龙格-库塔法基本思路二阶、三阶龙格-库塔法经典龙格-库塔法隐式龙格-库塔法线性多步法亚当斯法亚当斯预报-校正公式误差修正法收敛性与稳定性微分方程组和高阶微分方程15:24:53NumericalAnalysis3简单的数值方法与基本概念科学技术中常常需要求解常微分方程的定解问题.这类问题最简单的形式，是本章将着重考察的一阶方程的初值问题00(,),().yfxyyxy我们知道，只有f(x,y)适当光滑—譬如关于y满足利普希茨(Lipschitz)条件理论上就可以保证初值问题的解y＝f(x)存在并且唯一..),(),(yyLyxfyxf15:24:53NumericalAnalysis4虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法.所谓数值解法,就是寻求解y(x)在一系列离散节点121nnxxxx上的近似值y1,y2,,yn,yn+1,.相邻两个节点的间距hn=xn+1-xn称为步长.今后如不特别说明，总是假定hi=h(i=1,2,)为定数,这时节点为xn=x0+nh(i=0,1,2,)(等距节点).15:24:53NumericalAnalysis5初值问题的数值解法有个基本特点，他们都采取“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进.描述这类算法，只要给出用已知信息yn,yn-1,yn-2,计算yn+1的递推公式.首先，要对微分方程离散化，建立求解数值解的递推公式.一类是计算yn+1时只用到前一点的值yn，称为单步法.另一类是用到yn+1前面k点的值yn,yn-1,,yn-k+1，称为k步法.其次，要研究公式的局部截断误差和阶，数值解yn与精确解y(xn)的误差估计及收敛性，还有递推公式的计算稳定性等问题.15:24:53NumericalAnalysis6一阶常微分方程00(,)()yfxyyxy的解()yyx是通过点00(,)xy的一条曲线()yyx，称之为微分方程的积分曲线。积分曲线上每一点(,)xy的切线斜率()yx等于函数(,)fxy在这点的值。1几何推导从初始点000(,)Pxy出发，做切线0()yx，与1xx交于111(,)Pxy点，用1y作为曲线()yx上的点11(,())xyx的纵坐标1()yx的近似值。再从1P做切线1()yx，与2xx交于222(,)Pxy点，用2y作为曲线()yx上的点22(,())xyx的纵坐标2()yx的近似值。这样下去便可作出一条折线012PPP。设已作出折线的顶点为nP，再从nP做切线()nyx，推进到111(,)nnnPxy。欧拉法15:24:53NumericalAnalysis7过00(,)xy做以000()(,)yxfxy为切线斜率的方程0000(,)()yyfxyxx当1xx时，得100010(,)()yyfxyxx，取11()yxy。过11(,)xy做以111()(,)yxfxy为切线斜率的方程1111(,)()yyfxyxx当2xx时，得211121(,)()yyfxyxx，取22()yxy。一般地，过(,)nnxy做以()(,)nnnyxfxy为切线斜率的方程(,)()nnnnyyfxyxx当1nxx时，得11(,)()nnnnnnyyfxyxx，取1()nnyxy。从0x出发逐个算出12,,,nxxx，对应的数值解12,,,nyyy。一般取1nnxxh，得欧拉公式1(,)nnnnyyhfxy15:24:53NumericalAnalysis815:24:53NumericalAnalysis92欧拉法的数学推导:泰勒展开法将1()nyx在nx处做泰勒展开21()()()()()2!nnnnnhyxyxhyxhyxy当h充分小时，忽略高次项得22()()2!nhyOh因此，有欧拉公式1(,)nnnnyyhfxy15:24:53NumericalAnalysis102欧拉法数学推导：数值微分（用差商代替导数）nxxxx,,,,210设等距，步长1,0,1,nnhxxnyxfhxyxyhxyhxyhxyhxyxy,,)()()(令x=xn,x+h=xn+1,y(xn)≈yn，y(xn+1)≈yn+1,初值问题离散化为100(,),0,1,2,()nnnnyyhfxynyxy(欧拉公式)15:24:53NumericalAnalysis112欧拉法数学推导:数值积分将方程()(,)yxfxy两端从nx到1nx积分，有11()d(,())dnnnnxxxxyxxfxyxx11()()(,())dnnxnnxyxyxfxyxx算出积分项，可得1()nyx。利用左矩形公式1(,())d(,())nnxxfxyxxhfxyx代入，并离散化，有欧拉公式1(,)nnnnyyhfxy例7-1用欧拉公式求解初值问题2(01),(0)1.xyyxyy解取步长h=0.1，欧拉公式的具体形式为)2(1nnnnnyxyhyy其中xn=nh=0.1n(n=0,1,,10),已知y0=1,由此式可得191818.1)1.12.01.1(1.01.1)2(1.11.01)2(1111200001yxyhyyyxyhyy15:24:53NumericalAnalysis13依次计算下去，部分计算结果见下表.xy21与准确解相比，可看出欧拉公式的计算结果精度很差.xn欧拉公式数值解yn准确解y(xn)误差0.20.40.60.81.01.1918181.3582131.5089661.6497831.7847701.1832161.3416411.4832401.6124521.7320510.0086020.0165720.0257260.0373310.05271915:24:53NumericalAnalysis14局部截断误差和阶：数值公式的精度定义局部截断误差：假设第n步是准确的，即y(xn)=yn，将y(xn+1)-yn+1定义为数值方法的局部截断误差。由于实际上yn不是准确值，因此它的误差会传播下去。实际计算时，每一步都可能产生舍入误差。定义若局部截断误差为O(hp+1),p为正整数，则称数值公式是p阶公式。15:24:53NumericalAnalysis15欧拉公式的截断误差是O(h2)，公式是1阶的。1(,)()()nnnnnnyyhfxyyxhyx211()()()()2nnnyxyxyxhyh二阶泰勒公式两式相减，由设yn=y(xn)，有22112nnhyxyyOh欧拉公式的局部截断误差和阶15:24:53NumericalAnalysis16隐式（后退）欧拉公式：取1()nyx的向后差商111()[()()]nnnyxyxyxh替代111()(,)nnnyxfxy中的导数项，并离散化，有隐式欧拉公式111(,)nnnnyyhfxy隐式（后退）欧拉公式的局部截断误差：假设()nnyyx，则211()''()2nnnhyxyyx隐式（后退）欧拉公式15:24:53NumericalAnalysis17为了提高精度，改用中心差商111[()()]2nnyxyxh替代()(,)nnnyxfxy中的导数项，并离散化，有两步欧拉公式112(,)nnnnyyhfxy两步欧拉公式是两步法，要用前两步的值。两步欧拉公式的局部截断误差231()()()()()()2!3!nnnnnhhyxyxhyxhyxyxy231()()()()()()2!3!nnnnnhhyxyxhyxhyxyxy上二式相减，可得.两步欧拉公式15:24:53NumericalAnalysis18311()()2()()3nnnhyxyxhyxy311()()2()()3nnnhyxyxhyxy设()nnyyx，11()nnyyx前两步准确。则上式成为311()2(,)()3nnnnhyxyhfxyy与112(,)nnnnyyhfxy相比较，因此，有311()()3nnhyxyy两步欧拉公式的局部截断误差是3()Oh，是二阶方法。15:24:53NumericalAnalysis19梯形法对微分方程y′=f(x,y)两边求xn到xn+1的定积分，有11()()(,())dnnxnnxyxyxfxyxx利用梯形公式计算积分，有1111(,())d[(,())(,())]2nnxnnnnnnxxxfxyxxfxyxfxyx将y(xn)、y(xn+1)分别用yn、yn+1代替，构造数值公式11100[(,)(,)],0,1,2,2()nnnnnnhyyfxyfxynyyx15:24:53NumericalAnalysis20改进的欧拉法(预报-校正公式)欧拉法，显式，计算量小，精度低。梯形公式是隐式公式,计算量大，精度高。实际计算时，将二者综合之，先用欧拉公式计算出yn+1作为初始值,初始值精度不高，取作预报值，代入梯形公式，得到校正值yn+1。写成预报-校正公式1111(,)[(,)(,)]2nnnnnnnnnnyyhfxyhyyfxyfxy1(,)nnnnyyhfxy111[(,)(,)]2nnnnnnhyyfxyfxy15:24:53NumericalAnalysis21预报-校正公式又常常写成一步嵌套显式形式或写成平均化形式11[(,)(,(,))]2nnnnnnnnhyyfxyfxyhfxy11(,)(,)1()2pnnncnnpnpcyyhfxyyyhfxyyyy预报-校正公式的局部截断误差y(xn+1)-yn+1=O(h3)1111(,)[(,)(,)]2nnnnnnnnnnyyhfxyhyyfxyfxy15:24:53NumericalAnalysis22开始输出x1,y1结束00,y,,Nh输入x图7-2改进的欧拉法流程图1n10000011(,)y(x,y)y2pcppcxxhyyhfxyyhfyy?nN01011,nnxxyy≠15:24:53NumericalAnalysis23例7-2用改进欧拉法解初值问题2'(0)1xyyyy其中0,1x。解改进欧拉法公式112()2()1()2npnnnncnppnpcxyyhyyxyyhyyyyy0,1x，01000220()10.1(1)1.11xyyhyy取步长h=0.1，有n=015:24:53NumericalAnalysis24龙格-库塔（Runge-Kutta)法基本思想二阶龙格-库塔法经典龙格-库塔法15:24:53NumericalAnalysis25基本思想:一阶常微分方程初值问题00(,)()yfxyyxy欧拉公式1(,)nnnnyyhfxy，写成111(,)nnnnyyhkkfxy，调用1次函数。有1阶精度