第5章微分方程01欧拉法02R-K法.

第5章常微分方程初值问题的数值解法00)(],[),,(yxybaxyxfy微分方程常微分方程一阶常微分方程定阶条件：初值问题数值解法：给定点a=x0x1…xn=b,将初值问题离散化为差分方程,求出解函数(积分曲线)y(x)在这些点的近似值y1,y2,…,yn。所求得的近似值y1,y2,…,yn称为微分方程数值解。一阶常微分方程初值问题00(,)()dyfxydxyxy的数值解法。基本思想：常微分方程初值问题的数值解是求微分方程的解()yx（即微分方程初值问题的积分曲线），在区间[,]ab中给定一系列点（节点）1nnnxxh（1,2,n）上的近似值ny。这里nh为1nx到nx的步长，且0nh。第5章常微分方程初值问题的数值解法5.1欧拉法5.2龙格-库塔法5.3线性多步法5.4收敛性与稳定性5.5微分方程组和高阶微分方程5.1欧拉法和改进的欧拉法5.1.1欧拉公式5.1.2局部截断误差和阶5.1.3隐式（后退）欧拉公式和两步欧拉公式5.1.4梯形公式5.1.5改进的欧拉法(预报-校正公式)5.1欧拉法5.1.1欧拉公式一阶常微分方程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。过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欧拉公式的几何意义用一条初始点重合的折线，来近似表示微分方程的解（积分曲线）()yyx。2欧拉法的数学推导泰勒展开法将1()nyx在nx处做泰勒展开21()()()()()2!nnnnnhyxyxhyxhyxy当h充分小时，忽略高次项得22()()2!nhyOh因此，有欧拉公式1(,)nnnnyyhfxy3欧拉法数值微分推导用差商代替导数nxxxx,,,,210设等距，步长1,0,1,nnhxxnyxfhxyxyhxyhxyhxyhxyxy,,)()()(令x=xn,x+h=xn+1,y(xn)≈yn，y(xn+1)≈yn+1,初值问题离散化为00)(),(yxyyxfy初值问题100(,),0,1,2,()nnnnyyhfxynyxy(欧拉公式)4欧拉法的数值积分推导将方程()(,)yxfxy两端从nx到1nx积分，有11()d(,())dnnnnxxxxyxxfxyxx11()()(,())dnnxnnxyxyxfxyxx算出积分项，可得1()nyx。利用左矩形公式1(,())d(,())nnxxfxyxxhfxyx代入，并离散化，有欧拉公式1(,)nnnnyyhfxy用欧拉法解初值问题2'(0)1xyyyy其中0,1x。解欧拉法公式12(,)()nnnnnnnnxyyhfxyyhyy0,1x，取步长h=0.1，有n=001000220()10.1(1)1.11xyyhyy，1n12111220.1()1.10.1(1.1)1.19181.1xyyhyy局部截断误差和阶：数值公式的精度定义局部截断误差：假设第n步是准确的，即y(xn)=yn，将y(xn+1)-yn+1定义为数值方法的局部截断误差。由于实际上yn不是准确值，因此它的误差会传播下去。实际计算时，每一步都可能产生舍入误差。定义若局部截断误差为O(hp+1),p为正整数，则称数值公式是p阶公式。欧拉公式的截断误差是O(h2)，公式是1阶的。1(,)()()nnnnnnyyhfxyyxhyx211()()()()2nnnyxyxyxhyh二阶泰勒公式两式相减，由设yn=y(xn)，有22112nnhyxyyOh欧拉公式的局部截断误差和阶取步长2.0h，用欧拉法解初值问题1)0('2yxyyy其中6.0,0x。解用欧拉法求解公式，得)(),(21nnnnnnnnyxyhyyxhfyy取步长h=0.2时，6.0,0x，有n=0221000010.2(101)0.8yyhyxy1n22211110.80.20.80.20.80.6144yyhyxy2n22322220.61440.20.61440.40.61440.461321yyhyxy隐式（后退）欧拉公式：取1()nyx的向后差商111()[()()]nnnyxyxyxh替代111()(,)nnnyxfxy中的导数项，并离散化，有隐式欧拉公式111(,)nnnnyyhfxy隐式（后退）欧拉公式的局部截断误差：假设()nnyyx，则211()''()2nnnhyxyyx5.1.3两步欧拉公式1隐式（后退）欧拉公式为了提高精度，改用中心差商111[()()]2nnyxyxh替代()(,)nnnyxfxy中的导数项，并离散化，有两步欧拉公式112(,)nnnnyyhfxy两步欧拉公式是两步法，要用前两步的值。两步欧拉公式的局部截断误差231()()()()()()2!3!nnnnnhhyxyxhyxhyxyxy231()()()()()()2!3!nnnnnhhyxyxhyxhyxyxy上二式相减，可得2两步欧拉公式：中点方法311()()2()()3nnnhyxyxhyxy311()()2()()3nnnhyxyxhyxy设()nnyyx，11()nnyyx前两步准确。则上式成为311()2(,)()3nnnnhyxyhfxyy与112(,)nnnnyyhfxy相比较，因此，有311()()3nnhyxyy两步欧拉公式的局部截断误差是3()Oh，是二阶方法。5.1.3梯形法对微分方程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1.公式的构造梯形公式111(,)(,)2nnnnnnhyyfxyfxy梯形公式是将欧拉公式与隐式（后退）欧拉公式的算术平均，也是隐式公式。5.1.4改进的欧拉法(预报-校正公式)欧拉法，显式，计算量小，精度低。梯形公式是隐式公式,计算量大，精度高。实际计算时，将二者综合之，先用欧拉公式计算出yn+1作为初始值,初始值精度不高，取作预报值，代入梯形公式，得到校正值yn+1。写成预报-校正公式1111(,)[(,)(,)]2nnnnnnnnnnyyhfxyhyyfxyfxy1(,)nnnnyyhfxy111[(,)(,)]2nnnnnnhyyfxyfxy预报-校正公式又常常写成一步嵌套显式形式或写成平均化形式11[(,)(,(,))]2nnnnnnnnhyyfxyfxyhfxy11(,)(,)1()2pnnncnnpnpcyyhfxyyyhfxyyyy预报-校正公式的局部截断误差y(xn+1)-yn+1=O(h3)1111(,)[(,)(,)]2nnnnnnnnnnyyhfxyhyyfxyfxy32,,,21,hOyxfyxfyxfhyxfhyiiiiyiixiii预报-校正公式的局部截断误差假设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因而改进的欧拉法是二阶的。预报-校正公式的框图解改进欧拉法公式112()2()1()2npnnnncnppnpcxyyhyyxyyhyyyyy2'(0)1xyyyy0,1x用改进欧拉法解初值问题其中0,1xn01000220()10.1(1)1.11xyyhyy，取步长h=0.1，有=0解由题意知xyyyxfsin,2，改进欧拉法的具体形式为cpnnppncnnnnpyyyxyyhyyxyyhyy21sinsin11221)1(0sin'2yxyyy2.0h)4.1(),2.1(yy用欧拉预报-校正公式求解初值问题要求步长，计算的近似值。由1)1(0yy知10x，取步长2.0h，计算如下0n时22000022011sin10.211sin10.631706sin10.20.6317060.631706sin1.20.799272110.6317060.7992720.71548922pcpppcyyhyyxyyhyyxyyy1n时21111221222sin0.