常微分方程数值解法欧拉法

会计学1常微分方程数值解法欧拉法解析解法：（常微分方程理论）只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法。如何求解计算解函数y(x)在一系列节点a=x0<x1<…<xn=b

处的近似值节点间距为步长，通常采用等距节点，即取hi=h

(常数)。数值解法:求解所有的常微分方程步进式：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。因此只需建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。第1页/共32页第2页/共32页第3页/共32页--------Euler’sMethod§1欧拉方法

/*Euler’sMethod*/§1Euler’sMethodTaylor展开法第4页/共32页几何意义亦称为欧拉折线法

/*Euler’spolygonalarcmethod*/

几何直观是帮助我们寻找解决一个问题的思路的好办法哦第5页/共32页定义在假设yn=y(xn)，即第

n

步计算是精确的前提下，考虑公式或方法本身带来的误差:Rn=y(xn+1)

yn+1,称为局部截断误差/*localtruncationerror*/。说明

显然，这种近似有一定误差，而且步长越大，误差越大，如何估计这种误差y(xn+1)

yn+1

？§1Euler’sMethod第6页/共32页截断误差:实际上，y(xn)

yn，yn也有误差，它对yn+1的误差也有影响，见下图。但这里不考虑此误差的影响，仅考虑方法或公式本身带来的误差，因此称为方法误差或截断误差。局部截断误差的分析：由于假设yn=y(xn)

，即yn准确，因此分析局部截断误差时将y(xn+1)和

yn+1都用点xn上的信息来表示，工具：Taylor展开。

欧拉法的局部截断误差：Rn+1

的主项/*leadingterm*/§1Euler’sMethod第7页/共32页定义若某算法的局部截断误差为O(hp+1)，则称该算法有p

阶精度。欧拉法具有1阶精度。在xn点用一阶向前差商近似一阶导数

在第2章讨论牛顿插值公式时介绍了差商的概念和性质，各阶差商可以近似各阶导数，具有不同的精度，且可以用函数值来表示。上一章中数值微分的方法之一就是用差商近似导数Euler’smethod§1Euler’sMethod第8页/共32页§1Euler’sMethod

欧拉公式的改进：隐式欧拉法或后退欧拉法

/*implicitEulermethodorbackwardEulermethod*/xn+1点向后差商近似导数隐式或后退欧拉公式第9页/共32页由于未知数yn+1

同时出现在等式的两边，故称为隐式/*implicit*/

欧拉公式，而前者称为显式/*explicit*/欧拉公式。隐式公式不能直接求解，一般需要用Euler显式公式得到初值，然后用Euler隐式公式迭代求解。因此隐式公式较显式公式计算复杂，但稳定性好收敛性§1Euler’sMethod第10页/共32页

见上图，显然，这种近似也有一定误差，如何估计这种误差y(xn+1)

yn+1

？方法同上，基于Taylor展开估计局部截断误差。但是注意，隐式公式中右边含有f(xn+1

,yn+1),由于yn+1不准确，所以不能直接用y'(xn+1)代替f(xn+1

,yn+1)设已知曲线上一点Pn(xn,yn),过该点作弦线，斜率为(xn+1

,yn+1)点的方向场f(x,y)方向,若步长h充分小，可用弦线和垂线x=xn+1的交点近似曲线与垂线的交点。几何意义xnxn+1PnPn+1xyy(x)§1Euler’sMethod第11页/共32页

隐式欧拉法的局部截断误差：§1Euler’sMethod第12页/共32页§1Euler’sMethod第13页/共32页

隐式欧拉法的局部截断误差：即隐式欧拉公式具有1阶精度。§1Euler’sMethod第14页/共32页比较尤拉显式公式和隐式公式及其局部截断误差显式公式隐式公式§1Euler’sMethod第15页/共32页

若将这两种方法进行算术平均，即可消除误差的主要部分/*leadingterm*/而获得更高的精度,称为梯形法§1Euler’sMethod梯形公式/*trapezoidformula*/—显、隐式两种算法的平均注：的确有局部截断误差，即梯形公式具有2

阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。第16页/共32页梯形法的迭代计算和收敛性收敛性§1Euler’sMethod第17页/共32页梯形法的简化计算

迭代计算量大，且难以预测迭代次数。为了控制计算量，通常只迭代一次就转入下一点的计算。用显式公式作预测，梯形公式作校正，得到如下预测校正系统，也称为改进尤拉法：改进欧拉法/*modifiedEuler’smethod*/Step1:

先用显式欧拉公式作预测，算出),(1nnnnyxfhyy+=+Step2:再将代入隐式梯形公式的右边作校正，得到1+ny)],(),([2111+++++=nnnnnnyxfyxfhyy§1Euler’sMethod第18页/共32页注：此法亦称为预测-校正法/*predictor-correctormethod*/。可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。后面将看到，它的稳定性高于显式欧拉法。§1Euler’sMethod其它形式第19页/共32页几何解释xnxn+1ABPn+1=(A+B)/2尤拉法后退尤拉法梯形法§1Euler’sMethod第20页/共32页令x=x1,得Anotherpointofview对右端积分采用左矩形、右矩形、梯形积分公式，即可得尤拉显式、隐式、梯形公式§1Euler’sMethod第21页/共32页中点欧拉公式/*midpointformula*/中心差商近似导数x0x2x1假设，则可以导出即中点公式也具有2阶精度，且是显式的。需要2个初值y0和y1来启动递推过程，这样的算法称为双步法/*double-stepmethod*/，而前面的三种算法都是单步法/*single-stepmethod*/。§1Euler’sMethod第22页/共32页几何解释xnxn+1尤拉法后退尤拉法中点法xn-1令x=x2,得Anotherpointofview对右端积分采用中矩形公式即得中点公式§1Euler’sMethod第23页/共32页公式局部截断误差精度显隐稳定性步数尤拉显式公式1阶显差单步尤拉隐式公式1阶隐好单步梯形公式2阶隐差单步中点法2阶显好二步summary第24页/共32页第25页/共32页第26页/共32页

算例

分别用显式Euler方法，梯形方法和预估－校正Euler方法解初值问题解：取

h=0.1，梯形方法为：续第27页/共32页

算例

分别用显式Euler方法，梯形方法和预估－校正Euler方法解初值问题解：取

h=0.1，梯形方法为：预估－校正Euler方法：续第28页/共32页

Euler方法

梯形方法

预估－校正方法0.01.0000000.01.0000000.01.0000000.00.11.0000004.8×10－31.0047627.5×10－51.0050001.6×10－40.21.0100008.7×10－31.0185941.4×10－41.0190252.9×10－40.31.0290001.2×10－21.0406331.9×10－41.0412184.0×10－40.41.0561001.4×10－21.0700962.2×10－41.0708004.8×10－40.51.0904901.6×10－21.1062782.5×10－41.1070765.5×10－40.61.1314411.7×10－21.1485372.7×10－41.1494045.9×10－40.71.1782971.8×10－21.1962952.9×10－41.1972106.2×10－40.81.2304671.9×10－21.2490193.0×10－41.2499756.5×10－40.91.2874201.9×10－21.3062643.1×10－41.3072286.6×10－41.01.3486781.9×10－21.3675733.1×