数值分析第八章常微分方程数值解法解析课件.ppt

1、数值分析数值分析Numerical Analysis第八章常微分方程数值解法郑州大学研究生课程郑州大学研究生课程 （2014-20152014-2015学年第一学期）学年第一学期）2/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis第八章 常微分方程数值解法 8.1 8.1 引言引言8.2 8.2 欧拉欧拉(Euler)(Euler)法法8.3 8.3 改进欧拉改进欧拉(Euler)(Euler)方法方法8.4 8.4 单步法的稳定性单步法的稳定性3/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.1

2、 引言问题提出问题提出 倒葫芦形状容器壁上的刻度问题倒葫芦形状容器壁上的刻度问题.对于圆柱形状对于圆柱形状容器壁上的容积刻度容器壁上的容积刻度,可以利用圆柱体体积公式可以利用圆柱体体积公式HDV22其中直径D为常数.由于体积V与相对于容器底部的任意高度H的函数关系明确,因此在容器上可以方便地标出容器刻度。对于几何形状不是规则的容器,比如倒葫芦形状容器壁上如何标出刻度呢?4/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.1 引言下表是经过测量得到部分容器高度与直径的关系下表是经过测量得到部分容器高度与直径的关系.H 0 0.2 0.4 0.6 0

3、.8 1.0D 0 0.11 0.26 0.56 1.04 1.17dxDdV241根据上表的数据根据上表的数据,可以拟合出倒可以拟合出倒葫芦形状容器的图葫芦形状容器的图,建立如图所建立如图所示的坐标轴后示的坐标轴后,问题即为问题即为:如何如何根据任意高度根据任意高度x标出容器体积标出容器体积V的刻度的刻度,由由微元思想分析微元思想分析可知可知5/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.1 引言其中其中x表示高度表示高度,直径直径D是高度是高度x的函数的函数,记为记为D(x),),因此得到如下微分方程初值问题因此得到如下微分方程初值问题0

4、)0()(412VxDdxdV只要求解上述方程只要求解上述方程,就可求出体积就可求出体积V与高度与高度x之间的之间的函数关系函数关系,从而可标出容器壁上容积的刻度从而可标出容器壁上容积的刻度,但问题但问题是函数是函数D(x)无解析表达式无解析表达式,我们无法求出其解析解我们无法求出其解析解.6/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.1 引言 包含自变量包含自变量x、未知函数、未知函数y及未知函数的导数或及未知函数的导数或微分的方程称为微分的方程称为微分方程微分方程。在微分方程中。在微分方程中,自变量自变量的个数只有一个的个数只有一个,

5、称为称为常微分方程常微分方程。自变量的个数。自变量的个数为两个或两个以上的微分方程叫为两个或两个以上的微分方程叫偏微分方程偏微分方程。微分。微分方程中出现的未知函数最高阶导数的阶数称为微分方程中出现的未知函数最高阶导数的阶数称为微分方程的方程的阶阶数。如果未知函数数。如果未知函数y及其各阶导数及其各阶导数都是一次的都是一次的,则称它是则称它是线性线性的的,否则称为否则称为非线性非线性的。的。)(,nyyy 7/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis 常微分方程常微分方程(ODEs 未知函数是一元函数未知函数是一元函数)偏微分方程偏微分方程(

6、PDEs 未知函数是多元函数未知函数是多元函数)()dv tcgvdtm22(,)(,)(,)u t xu t xu t xutxx0yx222 8/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis,2)(的的解解xxy cxxy 2)(.,才才能能确确定定方方程程的的解解件件方方程程加加上上适适当当的的定定解解条条边边值值条条件件初初始始值值条条件件定定解解条条件件)2()1(:同一个微分方程,具有不同的初始条件微分方程的定解条件：微分方程的定解条件：9/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis当

7、当x=0时时,y=1,可得可得c=1时特解时特解当当x=0时时,y=1,可得可得c=-1时特解时特解两边积分两边积分通解通解00(,)()yf x yy xy 本章重点讨论一阶常微分方程的初值问题，本章重点讨论一阶常微分方程的初值问题，22001122|1|1(0)1xxdyxydyxdxycxcdxyyyy 例：21yx 220011cossincos(0)1|1|1xxdydyxdxcxcyxyydxyyy例：11 sinyx10/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.1 引言 在高等数学中，对于常微分方程的求解，给出在高等数学中

8、，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如了一些典型方程求解析解的基本方法，如可分离变可分离变量法、常系数齐次线性方程的解法、常系数非齐次量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法线性方程的解法等。但能求解的常微分方程仍然是等。但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。有限的，大多数的常微分方程是不可能给出解析解。11/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.1 引言 待求解的问题待求解的问题：一阶一阶常微分方程的常微分方程的初值问题初值问题/*Initial-Value

9、 Problem*/:00(,),()dyf x yxa bdxy xy 解的存在唯一性解的存在唯一性（“常微分方程常微分方程”理论）：只要理论）：只要 f(x,y)在在a,b R1 上连续，且关于上连续，且关于 y 满足满足 Lipschitz 条件条件，即存在与，即存在与 x,y 无关的常数无关的常数 L 使使则上述则上述IVP存在唯一解存在唯一解。1212|(,)(,)|,f x yf x yL yyxab 12/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis解析解法解析解法：（常微分方程理论）：（常微分方程理论）只能求解极少一类常微分方程；

10、实际中给定的问题不一只能求解极少一类常微分方程；实际中给定的问题不一定是解析表达式，而是函数表，无法用解析解法。定是解析表达式，而是函数表，无法用解析解法。数值解法：数值解法：递推法递推法如何求解如何求解13/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis14/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis记号：()y(,()(,)kkkkkky xf xy xf xy精确值，近似值，是一阶导数精确值,是一阶导数近似值15/66郑州大学研究生2011-2012学年课程 数值分析 Numerical An

11、alysis8.2 欧拉(Euler)法推导推导EulerEuler格式：格式：TaylorTaylor展开法展开法 数值微分数值微分 数值积分法数值积分法对微分方程的离散，可对微分方程的离散，可以有多种思路，但最基以有多种思路，但最基本的想法是本的想法是“以直代曲以直代曲”16/66 郑州大学研究生2014-2015学年课程 数值分析 Numerical Analysis16/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法(1)Taylor展开法展开法)(,()()()()(1nnnnnnxyxhfxyxyhxyxy

12、 100(,)()nnnnyyhf x yyy x方程初值问题方程初值问题Euler公式公式设给定等距剖分设给定等距剖分011,(n0,1,.)nnnxxxhxx，其中步长1()xnny xx将在点=处泰勒展开，当步长h充分小时，略去h2项,得2111()()()(),(,)2nnnnny xy xhyxh yxx17/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法(2)用差商近似导数用差商近似导数)(,()()()()(1nnnnnnxyxhfxyxyhxyxy 100(,)()nnnnyyhf x yyy x差分方

13、程初值问题差分方程初值问题向前向前Euler方法方法1()()(),nnny xy xy xh)(,()()(1nnnnxyxfhxyxy )(,()(nnnxyxfxy 18/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法若若用向后差商近似导数用向后差商近似导数，即，即)(,()()(111 nnnnxyxhfxyxy11100(,)()nnnnyyhf xyyy xhxyxyxynnn)()()(11 向后向后Euler方法方法)(,()()(111 nnnnxyxfhxyxy)(,()(111 nnnxyxfxy1

14、9/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法（3）用数值积分方法）用数值积分方法 1)(,()()(1nnxxnndxxyxfxyxy 11)(,()(nnnnxxxxdxxyxfdxxy1111(,)(,(),(,)(,(),()()nnnnnnnnf x yf xy xf x yf xy xy xyy xy分别用左矩形和右矩形公式，即 代替上式右端的积分,并注意,分别得到1111(,)(,)nnnnnnnnyyh f xyyyh f xy，。向前欧拉公式和向后欧拉公式：20/66 郑州大学研究生2011-201

15、2学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法若对积分用梯形公式，则得若对积分用梯形公式，则得 )(,()(,(2)()(111 nnnnnnxyxfxyxfhxyxy11100(,)(,)2()nnnnnnhyyf x yf xyyy x梯形欧拉公式梯形欧拉公式21/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法 欧拉（欧拉（Euler）方法是解初值问题的最）方法是解初值问题的最简单的数值方法。初值问题简单的数值方法。初值问题的解的解y=y(x)代表通过点代表通过点 的一条称之

16、为的一条称之为微分方程的微分方程的积分曲线积分曲线。积分曲线上每一点。积分曲线上每一点 的切线的斜率的切线的斜率 等于函数等于函数 在在这点的值。这点的值。00)(),(yxyyxfy),(00yx),(yx)(xy),(yxf22/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis Pi+1 Pn y=y(x)Pi Pn P1 Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 Euler法的求解过程是法的求解过程是:从初从初始点始点P0(即点即点(x0 0,y0 0)出发出发,作积分曲线作积分曲线y=y(x)在在P0点上点上切线切线 (

17、其斜率为其斜率为 ),),与与x=x1 1直线直线10PP000()(,()y xf x y x相交于相交于P1点点(即点即点(x1 1,y1 1),),得到得到y1 1作为作为y(x1 1)的近似值的近似值,如上图所示。过点如上图所示。过点(x0 0,y0 0),),以以f(x0 0,y0 0)为为斜率的切线斜率的切线方程为方程为 当当 时时,得得 )(,(0000 xxyxfyy1xx)(,(010001xxyxfyy这样就获得了这样就获得了P1 1点的坐标。点的坐标。23/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis同样同样,过过点点P1

18、(x1 1,y1 1),),作积分曲线作积分曲线y=y(x)的切线的切线交直线交直线x=x2 2于于P2点点,切线切线 的斜率的斜率 直线方程为直线方程为21PP)(1xy11(,)f xy)(,(1111xxyxfyy)(,(121112xxyxfyy当当 时时,得得 2xx Pi+1 Pn y=y(x)P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 24/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis当当 时时,得得由此获得了由此获得了P P2 2的坐标。重复以上过程的坐标。重复以上过程,就可获得一系就可获得一系

19、列的点列的点:P P1 1,P P1 1,P Pn n。对已求得点对已求得点以以 为斜率作直线为斜率作直线 ),(nnnyxP)(nxy(,)nnf xy)(,(nnnnxxyxfyy1nxx)(,(11nxnnnnxxyxfyynnyxy)(取取 Pi+1 Pn y=y(x)P1 Pi Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 25/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis 从图形上看从图形上看,就获得了一条近似于曲线就获得了一条近似于曲线y=y=y(xy(x)的折线的折线 。Pi+1 Pn y=y(x)P1 Pi

20、 Pn Pi+1 P0 x0 x1 xi xi+1 xn Pi P1 这样这样,从从x x0 0逐个算出逐个算出对应的数值解对应的数值解 nxxx,21nyyy,21nPPPP32126/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法00y)y(x );y,x(fydxdy x0 x1x2x30欧拉折线法欧拉折线法27/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法hyynnn:1 单单步步法法一一般般形形式式11100(,)(,)2()nnnnn

21、nhyyf xyf xyyy x1111(,)(,)nnnnnnnnyyh f xyyyh f xy，。28/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis例例8.2.1 用欧拉法解初值问题用欧拉法解初值问题 1)0()6.00(2yxxyyy取步长取步长h=0.2,h=0.2,计算过程保留计算过程保留4 4位小数位小数 解解:h=0.2,:h=0.2,欧拉迭代格式欧拉迭代格式 2),(xyyyxf21(,)nnnnnnnnyyhf xyyhyhx y0.2(4)(0,1,2 3)nnnyx yn，当当 n=0时，已知时，已知x0=0，y0=1，有

22、，有 y(x1)=y(0.2)y1=0.21(401)0.8当当 n=1时，已知时，已知x1=0.2，y1=0.8，有，有 y(0.4)y2=0.20.8(40.20.8)0.6144当当 n=2,时，已知时，已知x2=0.4，y2=0.6144，有，有 y(0.6)y3=0.20.6144(4-0.40.6144)=0.4613 29/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis:,0 x1 y(0)1,dyxydx例2 用欧拉方法求下列方程的数值解 25.1)15.0)(5.0(1)0.1(0.1)1(0.5)(01(0.5)11120001

23、yhxyyyyhxyyy解：解：Euler公式为公式为 ),(1nnnnyxhfyy 当当h=0.5时时,1,0 n nnnyhxy 25.0 5.0 hh分分别别取取30/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis1)1(0.25)(0125.0(0.25)0001 yxyyy),1,0(1 nyhxyynnnn当当h=0.25时时0625.1)125.0)(25.0(1 25.0)5.0(1112 yxyyy191347.1)0625.15.0)(25.0(0625.1 25.0)75.0(2223 yxyyy39600.1)191347.

24、175.0)(25.0(191347.1 )191347.1 ,75.0()75.0(25.0)0.1(3334 hfyyxyyy31/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis1)y(0 ;yxydxdy 00.50.751.010.25h=0.5h=0.25224)/x1(:y解析解为解析解为32/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法欧拉方法的收敛性11212()Taylor,(,),()()()()()2 ()(,()()(1)2nnnnnnnnnnnnn

25、ny xxxxhy xy xhy xhy xyhy xhf xy xy将在点展开1(),(8.2)()(,()(2)nnnnny xyy xh f xy x假定已知准确值 利用欧拉公式，定义假设第n步是准确的33/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法局部截断误差11112()()()(,()().2nnnnnnnnTy xyy xy xhf xy xhy2221 max|()|,|()|().22a x bnnMy xhhTyMO h 令则称为局部截断误差34/66郑州大学研究生2011-2012学年课程 数值

26、分析 Numerical Analysis8.2 欧拉(Euler)法欧拉方法的收敛性定义 若给定方法的局部截断误差满足则称该方法是 P 阶的，或称为具有 P 阶精度。11|(),pnTO h35/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法整体截断误差111().nnney xy记112121112 ,(),(),(),(),.nnnnnnyy yyy xy xy xy xey yy因为计算 时 用到的,是 的近似值 每步产生的误差会累积到计算的误差中因此 与,都有关 称整体截断误差为1111111111|()|()

27、|.(3)nnnnnnnnnney xyy xyyyTyy36/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法欧拉方法的收敛性1111,(8.2)|()(,()(,)|()|(,()(,)|()|()|()(1)nnnnnnnnnnnnnnnnnnnnyyyyy xh f xy xyh f xyy xyh f xy xf xyy xyhL y xyhL对应用欧拉公式得李普希兹条件|.(4)ne37/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis由此知，当 max|,(4)(3

28、)kkTT记 将代入得1121210112|(1)|(1)(1)|(1)(1)|(1)(1)(1)(1)|(1)1(1)1 ()(1)1 (nnnnnnnneThLeThL ThLeThL ThLeThL ThL ThL ThLehLhLTO hhLhLO).h10,|0,nhe欧拉公式是一阶有收敛的。38/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法 注 整体截断误差与局部截断误差的关系：1(1)()0,(1)e,(1)ee.nnznnhLL b azzhLconst111|(|).nneO hT39/66郑州大学研

29、究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法 向后欧拉公式隐式欧拉法或向后欧拉法隐式欧拉法或向后欧拉法/*implicit Euler method or backward Euler method*/11()()()nnny xy xy xhxn+1点向后差商近似导数点向后差商近似导数111111()()()()()(,)nnnnnnnnnny xy xhy xy xyy xyyh f xy代入隐式或后退欧拉公式隐式或后退欧拉公式40/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.

30、2 欧拉(Euler)法 向后欧拉公式由于未知数由于未知数 yn+1 同时出现在等式的两边，故称为同时出现在等式的两边，故称为隐式隐式/*implicit*/欧拉公式，而前者称为欧拉公式，而前者称为显式显式/*explicit*/欧拉公欧拉公式。隐式公式不能直接求解，一般需要用式。隐式公式不能直接求解，一般需要用Euler显式公式显式公式得到初值，然后用得到初值，然后用Euler隐式公式迭代求解隐式公式迭代求解。因此隐式公。因此隐式公式较显式公式计算复杂，但稳定性好（后面分析）。式较显式公式计算复杂，但稳定性好（后面分析）。隐式欧拉公式中的未知数隐式欧拉公式中的未知数 yn+1 可通过以下迭代

31、法可通过以下迭代法求解：求解：0)1(1)()111(,)(,)nnnnkknnnnyyh f xyyyh f xy（41/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis 见上图，见上图，显然，这种近似也有一定误差，显然，这种近似也有一定误差，如何估计这种误差如何估计这种误差y(xn+1)yn+1？方法同上，基于方法同上，基于Taylor展开展开估计估计局部截断误差局部截断误差。但是注意，隐式公式中右边含有但是注意，隐式公式中右边含有f(xn+1,yn+1),由于由于yn+1不准确，所以不能直接用不准确，所以不能直接用y(xn+1)代替代替f(xn

32、+1,yn+1)设已知曲线上一点设已知曲线上一点 Pn(xn,yn),过该过该点作弦线，斜率为点作弦线，斜率为(xn+1,yn+1)点的点的方向场方向场f(x,y)。若步长。若步长h充分小，可充分小，可用弦线和垂线用弦线和垂线x=xn+1的交点近似曲的交点近似曲线与垂线的交点。线与垂线的交点。几何意义xnxn+1PnPn+1xyy(x)42/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis 隐式隐式欧拉法的局部截断误差：欧拉法的局部截断误差：111()nnnTy xy 232()()hny xO h即隐式欧拉公式具有即隐式欧拉公式具有 1 阶精度。阶

33、精度。43/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.2 欧拉(Euler)法 向后欧拉公式1(,)0,1,.nnnnyyh f xyn比较欧拉显式公式和隐式公式及其局部截断误差231112()()()hnnnnTy xyy xO h显式公式111(,)nnnnyyh f xy隐式公式231112()()()hnnnnTy xyy xO h44/66郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis 若将这两种方法进行算术平均，即可消除误差若将这两种方法进行算术平均，即可消除误差的主要部分的主要部分/*

34、leading term*/而获得更高的精度而获得更高的精度,称为梯形法称为梯形法 梯形公式梯形公式/*trapezoid formula*/显、隐式两种算法的显、隐式两种算法的平均平均111(,)(,)2nnnnnnhyyf xyf xy注：注：的确有局部截断误差的确有局部截断误差 ，即梯形公式具有即梯形公式具有2 阶精度，比欧拉方法有了进步。但阶精度，比欧拉方法有了进步。但注意到该公式是注意到该公式是隐式隐式公式，计算时不得不用到迭代公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。法，其迭代收敛性与欧拉公式相似。3111()()nnnRy xyO h45/66郑州大学研究生201

35、1-2012学年课程 数值分析 Numerical Analysis例例8.2.3 对初值问题对初值问题 1)0(0yyy证明用梯形公式求得的近似解为证明用梯形公式求得的近似解为 nnhhy并证明当步长并证明当步长h h0 0时时,y yn n收敛于精确解收敛于精确解证明证明:解初值问题的梯形公式为解初值问题的梯形公式为 xe),(),(nnnnnnyxfyxfhyyyyxf),(211nnnnyyhyy 整理成显式整理成显式 nnyhhy反复迭代反复迭代,得到得到 yhhyhhyhhyhhynnnnn.10y22nxnhyeh 46/66 郑州大学研究生2011-2012学年课程 数值分析

36、Numerical Analysis公式公式局部截断误差局部截断误差精度精度显显隐隐稳定稳定性性步数步数欧拉显欧拉显式公式式公式1 1阶阶显显差差单步单步欧拉隐欧拉隐式公式式公式1 1阶阶隐隐好好单步单步梯形梯形公式公式2 2阶阶隐隐好好单步单步3(3)3nhyx2(2)2nhyx2(2)2nhyx欧拉法小结欧拉法小结47/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.3 改进欧拉(Euler)方法48/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.3 改进欧拉(Euler)方法49/66郑

37、州大学研究生课程 数值分析 Numerical Analysis8.3 改进欧拉(Euler)方法 显式欧拉公式计算工作量小显式欧拉公式计算工作量小,但精度低。梯形公但精度低。梯形公式虽提高了精度式虽提高了精度,但为但为隐式公式隐式公式,需用迭代法求解需用迭代法求解,计算计算工作量大。综合欧拉公式和梯形公式便可得到改进工作量大。综合欧拉公式和梯形公式便可得到改进的欧拉公式。的欧拉公式。结合已有格式的优点，以结合已有格式的优点，以得到计算方便、计算量减得到计算方便、计算量减少且精度保持的数值格式少且精度保持的数值格式50/66 郑州大学研究生2011-2012学年课程 数值分析 Numerica

38、l Analysis8.3 改进欧拉(Euler)方法 先用欧拉公式先用欧拉公式(8.2)求出一个初步的近似值求出一个初步的近似值,称为称为预测值预测值,它的精度不高它的精度不高,再用梯形公式对它校正再用梯形公式对它校正一次一次,即迭代一次即迭代一次,求得求得yn+1,称为称为校正值校正值,这种预测这种预测-校校正方法称为正方法称为改进的欧拉公式改进的欧拉公式：1ny 校校正正预预测测 ),(),(2 ),(1111nnnnnnnnnnyxfyxfhyyyxhfyy称为称为Euler公式与梯形公式的公式与梯形公式的预测预测校正系统校正系统。51/66 郑州大学研究生2011-2012学年课程

39、数值分析 Numerical Analysis8.3 改进欧拉(Euler)方法 )(21),(),(11qpnpnnqnnnpyyyyxhfyyyxhfyy实际计算时，常改写成以下形式实际计算时，常改写成以下形式几何几何解释解释xnxn+1ABPn+1=(A+B)/2欧拉法欧拉法梯形法梯形法改进欧拉法改进欧拉法52/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysispredictorcorrector),(01iiiiyxhfyy ),(),(21011 iiiiyxfyxf ),(),(20111 iiiiiiyxfyxfhyy53/66 郑州大

40、学研究生2011-2012学年课程 数值分析 Numerical Analysis8.3 改进欧拉(Euler)方法 可以证明可以证明,改进的欧拉公式改进的欧拉公式的精度为的精度为二阶二阶。这。这是一种一步是一种一步显式格式显式格式,它可以表示为嵌套形式它可以表示为嵌套形式。),(,(),(211iiiiiiiiyxhfyxfyxfhyy54/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis 校校正正预预测测 ),(),(2 ),(1111nnnnnnnnnnyxfyxfhyyyxhfyy法法计计算算公公式式改改进进的的解解Euler:nnnnnn

41、nnnnnyhhyyhyyyhyyhyyy2)(2 211nnnhhyhhyy)21()21(2021 nhx 令令xhhhhhxhhxhnhehhhhy 22222020022)21(lim)21(limlim例例8.3.18.3.155/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis )(21),(),(11qpnpnnqnnnpyyyyxhfyyyxhfyy.5,2.0,4.12.1 ,1)1(,sin:2位位小小数数点点后后至至少少保保留留取取步步长长处处的的近近似似值值和和的的解解在在方方法法计计算算初初值值问问题题用用改改进进的的例例

42、 hyxyyyEuler方方法法改改进进的的解解Euler:校校正正预预测测 ),(),(2 ),(1111nnnnnnnnnnyxfyxfhyyyxhfyyxyyyxfsin),(2 2211(sin)(sin)1()2pnnnnqnppnnpqyyhyyxyyhyyxyyy2.0,1)0(0 hyy56/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis)(21)sin()sin(1122qpnnppnpnnnnpyyyxyyhyyxyyhyy 1)0(1,00 yyn71549.0)2.1(799272.0 631706.0,01 yyyynp

43、p时时52611.0)4.1(575259.0 476964.0,12 yyyynpp时时8.3 改进欧拉(Euler)方法57/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis 开始 输入 x0,y0,h,N 1 n x0+h x1 y0+hf(x0,y0)yp y0+hf(x1,yp)yq (yp+yq)/2 y1 输出 x1,y1 n+1 n n=N?x1 x0 y1 y0 结束 n y 改进欧拉法的算法58/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.4 单步法的稳定性稳定性稳定性:误差

44、在以后各步的计算中不会无限制扩大误差在以后各步的计算中不会无限制扩大.稳定性在微分方程的数值解法中是一个非常重要的问稳定性在微分方程的数值解法中是一个非常重要的问题。因为微分方程初值问题的数值方法是用差分格式进行题。因为微分方程初值问题的数值方法是用差分格式进行计算的，而在差分方程的求解过程中，存在着各种计算误计算的，而在差分方程的求解过程中，存在着各种计算误差，这些计算误差如差，这些计算误差如舍入误差等引起的扰动，在传播过程舍入误差等引起的扰动，在传播过程中，可能会大量积累，对计算结果的准确性将产生影响中，可能会大量积累，对计算结果的准确性将产生影响。这就涉及到算法稳定性问题。这就涉及到算法

45、稳定性问题。59/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis例：例：考察初值问题考察初值问题 在区间在区间0,0.5上的解。上的解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。分别用欧拉显、隐式格式和改进的欧拉格式计算数值解。1)0()(30)(yxyxy0.00.10.20.30.40.5精确解精确解改进欧拉法改进欧拉法 欧拉隐式欧拉隐式欧拉显式欧拉显式 节点节点 xixey30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 10

46、21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 1013.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 78.4 单步法的稳定性60/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis8.4 单步法的稳定性61/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis定义定义若某算法在计算过程中任一步产生的误差在以后的计若某算法在计算过

47、程中任一步产生的误差在以后的计算中都算中都逐步衰减逐步衰减，则称该算法是，则称该算法是绝对稳定的绝对稳定的/*absolutely stable*/。一般分析某算法的稳定性时，为简单起见，只考虑一般分析某算法的稳定性时，为简单起见，只考虑模型方程模型方程或或试验方程试验方程/*test equation*/yy 8.4 单步法的稳定性62/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis记数值误差为：引进试验方程：*111,nnnyy1*1nnyy为由公式得到的准确值(无舍入误差)，为由公式实际得到的值(有舍入误差).,.yy 为一复常数63/66

48、 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis向前欧拉公式的稳定性试验方程的欧拉公式若每步计算有舍入误差，则 1(,),0,1,2,nnnnyyh f xyn1(,)(1)nnnnnnyyh f xyyh y*1 (2)nnnyyh y11(1)(2):(1),|1|1 nnnnnhhh当 稳定。64/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis当为复数时当为实数时|1|1 Re()11hh 表示 复平面上以为中心，以 为半径的圆,此圆称为向前欧拉方法的绝对稳定区域。|1|1 20.hh 绝对稳定区域为带状。65/66 郑州大学研究生2011-2012学年课程 数值分析 Numerical Analysis本章小结本章小结 本章介绍了常微分方程初值问题的基本数值解法。包括单步法和多步法。单步法主要有欧拉法、改进欧拉法和龙格库塔方法。多步法是亚当姆斯法。它们都是基于把一个连续的定解问题离散化为一个差分方程来求解，是一种步进式的方法。用多步法求常微分方程的数值解可获得较高的精度。实际应用时，选择合适的算法有一定的难度，既要考虑算法的简易性和计算量，又要考虑截断误差和收敛性、稳定性。