常微分方程数值解法.pptx

《常微分方程数值解法.pptx》由会员分享，可在线阅读，更多相关《常微分方程数值解法.pptx（74页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、会计学1常微分方程数值解法常微分方程数值解法1 引言引言n n1.0 基本概念n n1.1.常微分方程的初值问题：常微分方程的初值问题：n n称为具有初值称为具有初值称为具有初值称为具有初值(1.2)(1.2)的常微分方程的常微分方程的常微分方程的常微分方程.n n 若若若若f f(x x,y y)在在在在 a a x x b b,|,|y y|+|+上连续，且关于上连续，且关于上连续，且关于上连续，且关于y y满满满满足足足足LipLip条件：条件：条件：条件：常数常数常数常数L L使使使使|f f(x x,y y1 1)f f(x x,y y2 2)|)|L L|y y1 1 y y2 2

2、|n n则初值问题则初值问题则初值问题则初值问题(1.1)(1.2)(1.1)(1.2)存在唯一连续可微解存在唯一连续可微解存在唯一连续可微解存在唯一连续可微解y y(x x).).n n注：以下总假设注：以下总假设注：以下总假设注：以下总假设f f 满足满足满足满足LipLip条件条件条件条件.第1页/共73页1 引言引言n n1.0 基本概念n n1.1.常微分方程的初值问题：常微分方程的初值问题：n n称为具有初值称为具有初值称为具有初值称为具有初值(1.2)(1.2)的常微分方程的常微分方程的常微分方程的常微分方程.n n(1.1)(1.2)(1.1)(1.2)等价于微分方程：等价于微

3、分方程：等价于微分方程：等价于微分方程：n n (1.3)(1.3)n n注：一般无初等解注：一般无初等解注：一般无初等解注：一般无初等解(解析解解析解解析解解析解)，即使有形式也复杂，即使有形式也复杂，即使有形式也复杂，即使有形式也复杂.第2页/共73页1 引言引言n n1.0 基本概念n n2.2.初值问题的数值解初值问题的数值解n n n n设设设设(1.1)(1.2)(1.1)(1.2)的解的解的解的解y y(x x)在节点在节点在节点在节点x xi i处的近似解值为处的近似解值为处的近似解值为处的近似解值为n n y yi i y y(x xi i),),a a x x1 1 x x

4、2 2 x xn n=b bn n则称则称则称则称y yi i (i i=1,2,=1,2,n n)为为为为(1.1)(1.2)(1.1)(1.2)的数值解，又称的数值解，又称的数值解，又称的数值解，又称y y(x xi i)的计算值的计算值的计算值的计算值.第3页/共73页1 引言引言n n1.0 基本概念n n3.3.数值方法数值方法n n 两种转化：两种转化：两种转化：两种转化：n n 由微分出发的数值方法由微分出发的数值方法由微分出发的数值方法由微分出发的数值方法.n n 由积分由积分由积分由积分 出发的数值方法出发的数值方法出发的数值方法出发的数值方法.n n 计算方法计算方法计算方

5、法计算方法n n 步进法：从初始条件出发，逐步求步进法：从初始条件出发，逐步求步进法：从初始条件出发，逐步求步进法：从初始条件出发，逐步求y y1 1,y y2 2,y yn n.n n 又有两种：单步法，多步法又有两种：单步法，多步法又有两种：单步法，多步法又有两种：单步法，多步法.n n注：采用等距节点：注：采用等距节点：注：采用等距节点：注：采用等距节点：第4页/共73页1 引言引言n n1.1 基于数值微分的求解公式.n n n n (1.6)(1.6)第5页/共73页1 引言引言n n1.1 基于数值微分的求解公式.n n1.1.前进欧拉公式前进欧拉公式n n (1.6)(1.6)的

6、前半部分为：的前半部分为：的前半部分为：的前半部分为：n n令令令令 y yi i+1+1=y yi i +hfhf(x xi i,y yi i)(1.7)(1.7)n n其中其中其中其中y yi i=y y(x xi i),),则则则则y yi i+1+1 y y(x xi i+1+1)第6页/共73页1 引言引言n n1.1 基于数值微分的求解公式.n n1.1.前进欧拉公式前进欧拉公式n n n n令令令令 y yi i+1+1=y yi i +hfhf(x xi i,y yi i)(1.7)(1.7)n n其中其中其中其中y yi i=y y(x xi i),),则则则则y yi i+

7、1+1 y y(x xi i+1+1)n n记记记记 (1.8)(1.8)n n则则则则n n称称称称(1.7)(1.7)为前进欧拉求解公式为前进欧拉求解公式为前进欧拉求解公式为前进欧拉求解公式.简称为欧拉公式或欧拉简称为欧拉公式或欧拉简称为欧拉公式或欧拉简称为欧拉公式或欧拉法法法法.(1.8).(1.8)称为欧拉公式的余项：称为欧拉公式的余项：称为欧拉公式的余项：称为欧拉公式的余项：e ei i+1+1(h h)=)=y y(x xi i+1+1)y yi i+1+1 第7页/共73页1 引言引言n n1.1 基于数值微分的求解公式.n n2.2.后退欧拉公式后退欧拉公式n n (1.6)(

8、1.6)的后半部分的后半部分的后半部分的后半部分n n令令令令 y yi i+1+1=y yi i +hfhf(x xi i+1+1,y yi i+1+1)(1.9)(1.9)n n其中其中其中其中y yi i=y y(x xi i),),则则则则y yi i+1+1 y y(x xi i+1+1)第8页/共73页1 引言引言n n1.1 基于数值微分的求解公式.n n2.2.后退欧拉公式后退欧拉公式n n令令令令 y yi i+1+1=y yi i +hfhf(x xi i+1+1,y yi i+1+1)(1.9)(1.9)n n其中其中其中其中y yi i=y y(x xi i),),则则

9、则则y yi i+1+1 y y(x xi i+1+1)n n注：注：注：注：(1.9)(1.9)中中中中f f(x xi i+1+1,y yi i+1+1)f f(x xi i+1+1,y y(x xi i+1+1)n n余项余项余项余项n n n n (1.10)(1.10)第9页/共73页1 引言引言n n1.1 基于数值微分的求解公式.n n2.2.后退欧拉公式后退欧拉公式n n令令令令 y yi i+1+1=y yi i +hfhf(x xi i+1+1,y yi i+1+1)(1.9)(1.9)n n其中其中其中其中y yi i=y y(x xi i),),则则则则y yi i+1

10、+1 y y(x xi i+1+1)n n注：注：注：注：称称称称(1.9)(1.9)为后退欧拉公式为后退欧拉公式为后退欧拉公式为后退欧拉公式(后退欧拉法后退欧拉法后退欧拉法后退欧拉法).).n n 称称称称(1.10)(1.10)为后退欧拉法的误差近似值为后退欧拉法的误差近似值为后退欧拉法的误差近似值为后退欧拉法的误差近似值.n n 欧拉法与后退欧拉公式的区别：欧拉法与后退欧拉公式的区别：欧拉法与后退欧拉公式的区别：欧拉法与后退欧拉公式的区别：n n(1.7)(1.7)为直接计算公式称显式公式为直接计算公式称显式公式为直接计算公式称显式公式为直接计算公式称显式公式.n n(1.9)(1.9)

11、为关于函数方程称隐式公式为关于函数方程称隐式公式为关于函数方程称隐式公式为关于函数方程称隐式公式.第10页/共73页1 引言引言n n1.1 基于数值微分的求解公式.n n【例例例例1 1】取取取取h h=0.1=0.1求解初值问题：求解初值问题：求解初值问题：求解初值问题：n n (1.11).(1.11).n n解：解：解：解：，x xi i=ihih=0.1=0.1 i i,(,(i i=0,1,2,=0,1,2,10)10)n n 欧拉法：欧拉法：欧拉法：欧拉法：第11页/共73页1 引言引言n n1.1 基于数值微分的求解公式.n n【例例例例1 1】取取取取h h=0.1=0.1求

12、解初值问题：求解初值问题：求解初值问题：求解初值问题：n n (1.11).(1.11).n n解：解：解：解：，x xi i=ihih=0.1=0.1 i i,(,(i i=0,1,2,=0,1,2,10)10)n n 后退欧拉法：后退欧拉法：后退欧拉法：后退欧拉法：n n 第12页/共73页1 引言引言n n1.1 基于数值微分的求解公式.n n注：为避免求解函数方程，采用显式与隐式结合的注：为避免求解函数方程，采用显式与隐式结合的注：为避免求解函数方程，采用显式与隐式结合的注：为避免求解函数方程，采用显式与隐式结合的方法：方法：方法：方法：n n n n此方法称为此方法称为此方法称为此方

13、法称为 预测预测预测预测校正系统校正系统校正系统校正系统.求解过程为：求解过程为：求解过程为：求解过程为：第13页/共73页1 引言引言n n1.1 基于数值微分的求解公式.n n预测预测预测预测校正系统：校正系统：校正系统：校正系统：n n【例例例例2 2】利用预测利用预测利用预测利用预测校正系统求解例校正系统求解例校正系统求解例校正系统求解例1.1.第14页/共73页1 引言引言n n1.1 基于数值微分的求解公式.n n预测预测预测预测校正系统：校正系统：校正系统：校正系统：n n注：显式比隐式方便，但有时隐式效果比显式好注：显式比隐式方便，但有时隐式效果比显式好注：显式比隐式方便，但有

14、时隐式效果比显式好注：显式比隐式方便，但有时隐式效果比显式好.(4.(4介绍介绍介绍介绍).).第15页/共73页1 引言引言n n1.2 截断误差n n定义定义定义定义1.1 1.1 称称称称e ek k(h h)=)=y y(x xk k)y yk k为计算为计算为计算为计算y yk k的公式第的公式第的公式第的公式第k k步的局部截步的局部截步的局部截步的局部截断误差断误差断误差断误差.n n注：注：注：注：“局部局部局部局部”是指在计算第是指在计算第是指在计算第是指在计算第k k步时，假定前面步时，假定前面步时，假定前面步时，假定前面y yi i=y y(x xi i)(i i k k

15、).).而而而而y yk k y y(x xk k)n n n n 欧拉法欧拉法欧拉法欧拉法.n n 后退欧拉法后退欧拉法后退欧拉法后退欧拉法.n n一般根据一般根据一般根据一般根据yy(x xk k)对对对对yy(k k),),yy(k k)做估计做估计做估计做估计.第16页/共73页1 引言引言n n1.2 截断误差n n定义定义定义定义1.2 1.2 设设设设e ei i(h h)()(i i=1,2,=1,2,k k)为求解公式第为求解公式第为求解公式第为求解公式第i i步的局部截断步的局部截断步的局部截断步的局部截断误差误差误差误差.称称称称n n为该求解公式在点上的整体截断误差为该

16、求解公式在点上的整体截断误差为该求解公式在点上的整体截断误差为该求解公式在点上的整体截断误差.n n注：注：注：注：局部截断误差局部截断误差局部截断误差局部截断误差e ek k(h h)与与与与y yk k有关有关有关有关.n n 整体截断误差整体截断误差整体截断误差整体截断误差E Ek k(h h)与与与与y y1 1,y y2 2,y yk k有关有关有关有关.n n所有所有所有所有e ek k(h h)都与都与都与都与h h有关有关有关有关.第17页/共73页1 引言引言n n1.2 截断误差n n定义定义定义定义1.3 1.3 若局部截断误差若局部截断误差若局部截断误差若局部截断误差e

17、 e(h h)=)=OO(h hp p+1+1)，则称该求解公式具，则称该求解公式具，则称该求解公式具，则称该求解公式具有有有有p p阶精度阶精度阶精度阶精度.n n注：欧拉法具有一阶精度注：欧拉法具有一阶精度注：欧拉法具有一阶精度注：欧拉法具有一阶精度.(.(精度越高越好精度越高越好精度越高越好精度越高越好)第18页/共73页1 引言引言n n作业作业作业作业 P208 1P208 1，2 2，3.3.第19页/共73页1 引言引言n n1.3 基于数值积分的求解公式n n (1.13)(1.13)n n若已知若已知若已知若已知y y(x xk k)=)=y yk k,则计算积分可求出则计算

18、积分可求出则计算积分可求出则计算积分可求出y y(x xk k+1+1).n n 如用矩形公式求积分如用矩形公式求积分如用矩形公式求积分如用矩形公式求积分n n则有则有则有则有y y(x xk k+1+1)=)=y y(x xk k)+)+hfhf(x xk k,y yk k)n n令令令令y yk k+1+1=y y(x xk k)+)+hfhf(x xk k,y yk k)即为欧拉公式即为欧拉公式即为欧拉公式即为欧拉公式.故欧拉公式故欧拉公式故欧拉公式故欧拉公式又称矩形法又称矩形法又称矩形法又称矩形法.第20页/共73页1 引言引言n n1.3 基于数值积分的求解公式n n (1.13)(

19、1.13)n n考虑考虑考虑考虑n n1.1.梯形公式梯形公式n n记记记记 (1.14)(1.14)第21页/共73页1 引言引言n n1.3 基于数值积分的求解公式n n1.1.梯形公式梯形公式n n记记记记 (1.14)(1.14)n n称称称称(1.14)(1.14)为梯形为梯形为梯形为梯形(求解求解求解求解)公式公式公式公式.简称梯形法简称梯形法简称梯形法简称梯形法.第22页/共73页1 引言引言n n1.3 基于数值积分的求解公式n n1.1.梯形公式梯形公式n n梯形梯形梯形梯形(求解求解求解求解)公式公式公式公式,简称梯形法简称梯形法简称梯形法简称梯形法:n n (1.14)(

20、1.14)n n注：注：注：注：梯形公式的余项：梯形公式的余项：梯形公式的余项：梯形公式的余项：n n n n故是二阶精度故是二阶精度故是二阶精度故是二阶精度.第23页/共73页n n1.3 基于数值积分的求解公式n n1.1.梯形公式梯形公式n n (1.14)(1.14)n n 梯形公式为隐式公式梯形公式为隐式公式梯形公式为隐式公式梯形公式为隐式公式.n n 预测预测预测预测校正系统校正系统校正系统校正系统n n (1.15)(1.15)n n称称称称(1.15)(1.15)为改进的欧拉公式，也可记为为改进的欧拉公式，也可记为为改进的欧拉公式，也可记为为改进的欧拉公式，也可记为1 引言引言

21、第24页/共73页1 引言引言n n1.3 基于数值积分的求解公式n n1.1.梯形公式梯形公式n n (1.14)(1.14)n n 可以证明，改进欧拉公式也具有二阶精度可以证明，改进欧拉公式也具有二阶精度可以证明，改进欧拉公式也具有二阶精度可以证明，改进欧拉公式也具有二阶精度.第25页/共73页1 引言引言n n1.3 基于数值积分的求解公式n n【例例例例3 3】用欧拉法，梯形法以及改进欧拉法求解用欧拉法，梯形法以及改进欧拉法求解用欧拉法，梯形法以及改进欧拉法求解用欧拉法，梯形法以及改进欧拉法求解n n取取取取h h=0.1.=0.1.计算到计算到计算到计算到x x=0.5.=0.5.n

22、 n解解解解：f f(x x,y y)=)=x x y y+1,+1,a a=x x0 0=0,=0,b b=0.5,=0.5,y y0 0=1,=1,n n=5=5n n(Euler(Euler法法法法)求解公式：求解公式：求解公式：求解公式：y yk k=y yk k11+h h(x xk k11 y yk k11+1)+1)n n=hxhxk k11+(1 +(1 h h)y yk k11+h h n n=0.1=0.1x xk k11+0.9+0.9y yk k11+0.1+0.1 第26页/共73页1 引言引言n n1.3 基于数值积分的求解公式n n【例例例例3 3】用欧拉法，梯形

23、法以及改进欧拉法求解用欧拉法，梯形法以及改进欧拉法求解用欧拉法，梯形法以及改进欧拉法求解用欧拉法，梯形法以及改进欧拉法求解n n解解解解：f f(x x,y y)=)=x x y y+1,+1,a a=x x0 0=0,=0,b b=0.5,=0.5,y y0 0=1,=1,n n=5=5n n(梯形法梯形法梯形法梯形法)求解公式：求解公式：求解公式：求解公式：n ny yk k=y yk k11+h h(x xk k11 y yk k11+1)+(+1)+(x xk k y yk k+1)/2+1)/2n n解出解出解出解出y yk k，得，得，得，得方程方程第27页/共73页1 引言引言n

24、 n1.3 基于数值积分的求解公式n n【例例例例3 3】用欧拉法，梯形法以及改进欧拉法求解用欧拉法，梯形法以及改进欧拉法求解用欧拉法，梯形法以及改进欧拉法求解用欧拉法，梯形法以及改进欧拉法求解n n解解解解：f f(x x,y y)=)=x x y y+1,+1,a a=x x0 0=0,=0,b b=0.5,=0.5,y y0 0=1,=1,n n=5=5n n(改进改进改进改进EulerEuler法法法法)求解公式：求解公式：求解公式：求解公式：n ny yk k=y yk k11+h h(x xk k11 y yk k11+1)+1)+x xk k (y yk k+h h(x xk k

25、 y yk k+1)+1/2+1)+1/2n n得得得得n n=0.905=0.905y yk k11+0.045+0.045x xk k11+0.05+0.05x xk k+0.095+0.095方程方程第28页/共73页1 引言引言n n1.3 基于数值积分的求解公式n n2.2.辛卜生公式辛卜生公式n n n n记记记记 n n (1.17)(1.17)第29页/共73页1 引言引言n n1.3 基于数值积分的求解公式n n2.2.辛卜生公式辛卜生公式n n记记记记 n n (1.17)(1.17)n n其余项其余项其余项其余项第30页/共73页1 引言引言n n1.3 基于数值积分的求

26、解公式n n2.2.辛卜生公式辛卜生公式n n记记记记 n n (1.17)(1.17)n n将将将将 x xk k11,x xk k 对分：对分：对分：对分：调整下标为调整下标为调整下标为调整下标为 x xi i22,x xi i ：n nx xi i22=x xk k11,x xi i11=x xk k11+h h1 1,x xi i=x xk k11+2+2h h1 1=x xk kn n则则则则(1.17)(1.17)化为化为化为化为n n (1.19)(1.19)n n称称称称(1.19)(1.19)为辛卜生求解公式，其中为辛卜生求解公式，其中为辛卜生求解公式，其中为辛卜生求解公式，

27、其中f fk k22=f f(x xk k22,y y(x xk k22)，f fk k11=f f(x xk k11,y y(x xk k11)，f fk k=f f(x xk k,y y(x xk k)第31页/共73页1 引言引言n n1.3 基于数值积分的求解公式n n2.2.辛卜生公式辛卜生公式n n记记记记 n n (1.17)(1.17)n n (1.19)(1.19)n n称称称称(1.19)(1.19)为辛卜生求解公式，其中为辛卜生求解公式，其中为辛卜生求解公式，其中为辛卜生求解公式，其中f fi i22=f f(x xi i22,y y(x xi i22)，f fi i11

28、=f f(x xi i11,y y(x xi i11)，f fi i=f f(x xi i,y y(x xi i)n n注：注：注：注：(1.19)(1.19)的误差：的误差：的误差：的误差：第32页/共73页1 引言引言n n1.3 基于数值积分的求解公式n n2.2.辛卜生公式辛卜生公式n n记记记记 n n (1.17)(1.17)n n (1.19)(1.19)n n称称称称(1.19)(1.19)为辛卜生求解公式，其中为辛卜生求解公式，其中为辛卜生求解公式，其中为辛卜生求解公式，其中f fi i22=f f(x xi i22,y y(x xi i22)，f fi i11=f f(x

29、xi i11,y y(x xi i11)，f fi i=f f(x xi i,y y(x xi i)n n注：注：注：注：隐式隐式隐式隐式(需显化需显化需显化需显化)n n多步多步多步多步将在将在将在将在33中讨论中讨论中讨论中讨论.第33页/共73页2 Runge-Kutta法法n n2.0 原理n n n n其中其中其中其中K K =f f(,y y()=)=y y()称为称为称为称为y y在在在在 x xi i11,x xi i 上的平上的平上的平上的平均斜率均斜率均斜率均斜率.n n欧拉法：欧拉法：欧拉法：欧拉法：n n改进欧拉法：改进欧拉法：改进欧拉法：改进欧拉法：n n(2.1)(

30、2.1)第34页/共73页2 Runge-Kutta法法n n2.0 原理n n n n其中其中其中其中K K =f f(,y y()=)=y y()称为称为称为称为y y在在在在 x xi i11,x xi i 上的平上的平上的平上的平均斜率均斜率均斜率均斜率.n n对对对对(1.17)(1.17)显化：显化：显化：显化：n n辛卜生：辛卜生：辛卜生：辛卜生：n n (2.4)(2.4)第35页/共73页2 Runge-Kutta法法n n2.0 原理n n其中其中其中其中K K =f f(,y y()=)=y y()称为称为称为称为y y在在在在 x xi i11,x xi i 上的平上的

31、平上的平上的平均斜率均斜率均斜率均斜率.n n设想：在中多计算设想：在中多计算设想：在中多计算设想：在中多计算(预测预测预测预测)几个点上的值然后可加权几个点上的值然后可加权几个点上的值然后可加权几个点上的值然后可加权取平均值作为的近似值可能构成更高阶的公式取平均值作为的近似值可能构成更高阶的公式取平均值作为的近似值可能构成更高阶的公式取平均值作为的近似值可能构成更高阶的公式.一阶一阶二阶二阶三阶三阶第36页/共73页2 Runge-Kutta法法n n2.1 Runge-Kutta公式n n (*)(*)n n其中其中其中其中0 0 j j 1 1，y yi i11+j jh h是是是是y

32、y(x xi i11+j jh h)的预测值的预测值的预测值的预测值.n n 称称称称(*)(*)为为为为R-KR-K公式公式公式公式R-K n n注：注：注：注：(2.1)(2.4)(2.1)(2.4)分别称为二阶，三阶分别称为二阶，三阶分别称为二阶，三阶分别称为二阶，三阶R-KR-K公式公式公式公式.n n j j，j j，j j为待定系数为待定系数为待定系数为待定系数.使使使使(*)(*)的阶数尽量高的阶数尽量高的阶数尽量高的阶数尽量高.第37页/共73页2 Runge-Kutta法法n n2.1 Runge-Kutta公式n n参数的确定，以参数的确定，以参数的确定，以参数的确定，以m

33、 m=2=2为例为例为例为例.n n n n欲求欲求欲求欲求 1 1，2 2，2 2 .原则原则:使使ei(h)=y(xi)yi的阶数尽可能高的阶数尽可能高第38页/共73页2 Runge-Kutta法法n n2.1 Runge-Kutta公式n n展开展开展开展开n n展开展开展开展开n n 原则原则:使使ei(h)=y(xi)yi的阶数尽可能高的阶数尽可能高第39页/共73页2 Runge-Kutta法法n n2.1 Runge-Kutta公式n n 原则原则:使使ei(h)=y(xi)yi的阶数尽可能高的阶数尽可能高第40页/共73页2 Runge-Kutta法法n n2.1 Runge

34、-Kutta公式n n n n欲求截断误差欲求截断误差欲求截断误差欲求截断误差e ei i(h h)=)=y y(x xi i)y yi i关于关于关于关于h h的阶数尽可能高，的阶数尽可能高，的阶数尽可能高，的阶数尽可能高，应使应使应使应使无穷多解，从而有许多无穷多解，从而有许多2阶阶R-K公式公式第41页/共73页2 Runge-Kutta法法n n2.1 Runge-Kutta公式n n应使应使应使应使n n注：注：注：注：取取取取 1 1=2 2=1/2=1/2，2 2=1=1，即为改进欧拉公式，即为改进欧拉公式，即为改进欧拉公式，即为改进欧拉公式.第42页/共73页2 Runge-K

35、utta法法n n2.1 Runge-Kutta公式n n应使应使应使应使n n注：注：注：注：取取取取 1 1=0=0，2 2=1=1，2 2=1/2=1/2，即为中点公式，即为中点公式，即为中点公式，即为中点公式第43页/共73页2 Runge-Kutta法法n n2.1 Runge-Kutta公式n n应使应使应使应使n n注：注：注：注：二阶二阶二阶二阶R-KR-K公式的截断误差为故为二阶方法公式的截断误差为故为二阶方法公式的截断误差为故为二阶方法公式的截断误差为故为二阶方法.相相相相仿可得更高阶的仿可得更高阶的仿可得更高阶的仿可得更高阶的R-KR-K公式公式公式公式.第44页/共73

36、页2 Runge-Kutta法法n n2.2 经典R-K公式n n 在在在在4 4解解解解R-KR-K公式中最重要的是经典公式中最重要的是经典公式中最重要的是经典公式中最重要的是经典R-KR-K公式公式公式公式.n n (2.6)(2.6)n n注：注：注：注：(2.6)(2.6)为为为为4 4阶方法阶方法阶方法阶方法.第45页/共73页2 Runge-Kutta法法n n2.2 经典R-K公式n n 在在在在4 4解解解解R-KR-K公式中最重要的是经典公式中最重要的是经典公式中最重要的是经典公式中最重要的是经典R-KR-K公式公式公式公式.n n (2.6)(2.6)n n注：注：注：注：

37、R-KR-K法对法对法对法对4 4阶以上不一定能提高整数阶阶以上不一定能提高整数阶阶以上不一定能提高整数阶阶以上不一定能提高整数阶.第46页/共73页2 Runge-Kutta法法n n2.2 经典R-K公式n n【例例例例4 4】使用三阶，四阶使用三阶，四阶使用三阶，四阶使用三阶，四阶R-KR-K法求解初值问题：法求解初值问题：法求解初值问题：法求解初值问题：n n n n的部分计算值的部分计算值的部分计算值的部分计算值y y1 1，y y2 2，y y3 3，其中，其中，其中，其中h h=0.1.=0.1.n n解解解解 使用三阶使用三阶使用三阶使用三阶R-KR-K法法法法第47页/共73

38、页2 Runge-Kutta法法n n【例例例例4 4】使用三阶，四阶使用三阶，四阶使用三阶，四阶使用三阶，四阶R-KR-K法求解初值问题：法求解初值问题：法求解初值问题：法求解初值问题：n n n n的部分计算值的部分计算值的部分计算值的部分计算值y y1 1，y y2 2，y y3 3，其中，其中，其中，其中h h=0.1.=0.1.n n解解解解 使用三阶使用三阶使用三阶使用三阶R-KR-K法法法法第48页/共73页2 Runge-Kutta法法n n【例例例例4 4】使用三阶，四阶使用三阶，四阶使用三阶，四阶使用三阶，四阶R-KR-K法求解初值问题：法求解初值问题：法求解初值问题：法求

39、解初值问题：n n n n的部分计算值的部分计算值的部分计算值的部分计算值y y1 1，y y2 2，y y3 3，其中，其中，其中，其中h h=0.1.=0.1.n n解解解解 使用四阶使用四阶使用四阶使用四阶R-KR-K法法法法第49页/共73页2 Runge-Kutta法法n n【例例例例4 4】使用三阶，四阶使用三阶，四阶使用三阶，四阶使用三阶，四阶R-KR-K法求解初值问题：法求解初值问题：法求解初值问题：法求解初值问题：n n n n的部分计算值的部分计算值的部分计算值的部分计算值y y1 1，y y2 2，y y3 3，其中，其中，其中，其中h h=0.1.=0.1.n n解解解

40、解 使用四阶使用四阶使用四阶使用四阶R-KR-K法法法法第50页/共73页2 Runge-Kutta法法n n注注注注 使用使用使用使用R-KR-K法要求具备较好的光滑性，否则效果不法要求具备较好的光滑性，否则效果不法要求具备较好的光滑性，否则效果不法要求具备较好的光滑性，否则效果不如低阶的如低阶的如低阶的如低阶的.n n作业作业作业作业P209 8 P209 8 9 9，10.10.第51页/共73页3 线性多步法线性多步法n n单步法的优点：简单，计算单步法的优点：简单，计算单步法的优点：简单，计算单步法的优点：简单，计算y yk k+1+1只用只用只用只用y yk k.n n缺点缺点缺点

41、缺点:没有充分利用前面的信息且计算没有充分利用前面的信息且计算没有充分利用前面的信息且计算没有充分利用前面的信息且计算y y(x xk k+h h)较困较困较困较困难难难难n n回顾回顾回顾回顾SimpsonSimpson：(1.19)(1.19)n n考虑：考虑：考虑：考虑：(3.1)(3.1)n n两种插值求积：两种插值求积：两种插值求积：两种插值求积：n n 将将将将 x xk k11,x xk k 增加内部节点，改为增加内部节点，改为增加内部节点，改为增加内部节点，改为 x xk k22,x xk k 导出的公导出的公导出的公导出的公式称为闭型求解公式式称为闭型求解公式式称为闭型求解公

42、式式称为闭型求解公式.线性多步线性多步第52页/共73页3 线性多步法线性多步法n n考虑：考虑：考虑：考虑：(3.1)(3.1)n n两种插值求积：两种插值求积：两种插值求积：两种插值求积：n n 将将将将 x xk k11,x xk k 增加内部节点，改为增加内部节点，改为增加内部节点，改为增加内部节点，改为 x xk k22,x xk k 导出的公导出的公导出的公导出的公式称为闭型求解公式式称为闭型求解公式式称为闭型求解公式式称为闭型求解公式.n n在在在在 x xk k11,x xk k 外增加插值节点，导出的公式称为开型外增加插值节点，导出的公式称为开型外增加插值节点，导出的公式称为

43、开型外增加插值节点，导出的公式称为开型求解公式求解公式求解公式求解公式.n n开型有显和隐，闭型也有显和隐开型有显和隐，闭型也有显和隐开型有显和隐，闭型也有显和隐开型有显和隐，闭型也有显和隐.第53页/共73页3 线性多步法线性多步法n n3.1 开型求解公式n n1.1.亚当斯显式求解公式亚当斯显式求解公式n n 取节点取节点取节点取节点x xk k33,x xk k22,x xk k11，在，在，在，在 x xk k33,x xk k 上作上作上作上作F F(x x)=)=f f(x x,y y(x x)的插值多项式的插值多项式的插值多项式的插值多项式.第54页/共73页3 线性多步法线性

44、多步法n n3.1 开型求解公式n n1.1.亚当斯显式求解公式亚当斯显式求解公式n n 取节点取节点取节点取节点x xk k33,x xk k22,x xk k11，在，在，在，在 x xk k33,x xk k 上上上上n n记记记记x xk k i i=x xk k ihih,x x=x xk k+thth，则，则，则，则n n 第55页/共73页3 线性多步法线性多步法n n3.1 开型求解公式n n1.1.亚当斯显式求解公式亚当斯显式求解公式n n 取节点取节点取节点取节点x xk k33,x xk k22,x xk k11，记，记，记，记x xk k i i=x xk k ihih

45、,x x=x xk k+thth，则，则，则，则n n n n代入代入代入代入(3.1)(3.1)得得得得第56页/共73页3 线性多步法线性多步法n n3.1 开型求解公式n n1.1.亚当斯显式求解公式亚当斯显式求解公式n n 取节点取节点取节点取节点x xk k33,x xk k22,x xk k11，记，记，记，记x xk k i i=x xk k ihih,x x=x xk k+thth，则，则，则，则n n n n令令令令n n (3.4)(3.4)n n称称称称(3.4)(3.4)为亚当斯显式求解公式为亚当斯显式求解公式为亚当斯显式求解公式为亚当斯显式求解公式(线性多步线性多步线

46、性多步线性多步).).第57页/共73页3 线性多步法线性多步法n n3.1 开型求解公式n n1.1.亚当斯显式求解公式亚当斯显式求解公式n n 取节点取节点取节点取节点x xk k33,x xk k22,x xk k11，记，记，记，记x xk k i i=x xk k ihih,x x=x xk k+thth，则，则，则，则n n余项：余项：余项：余项：n nn n 从而从而从而从而(3.4)(3.4)具有具有具有具有3 3阶精度阶精度阶精度阶精度.称为称为称为称为3 3阶亚当斯求解公式阶亚当斯求解公式阶亚当斯求解公式阶亚当斯求解公式.第58页/共73页3 线性多步法线性多步法n n3.

47、1 开型求解公式n n1.1.亚当斯显式求解公式亚当斯显式求解公式n n类似地取类似地取类似地取类似地取x xk k44,x xk k33,x xk k22,x xk k11 n n在在在在 x xk k44,x xk k 上作上作上作上作F F(x x)=)=f f(x x,y y(x x)的插值多项式，可导的插值多项式，可导的插值多项式，可导的插值多项式，可导出出出出4 4阶亚当斯显式求解公式：阶亚当斯显式求解公式：阶亚当斯显式求解公式：阶亚当斯显式求解公式：n n n n (3.6)(3.6)n n (3.7)(3.7)n n44阶精度阶精度阶精度阶精度第59页/共73页3 线性多步法线

48、性多步法n n3.1 开型求解公式n n2.2.亚当斯隐式求解公式亚当斯隐式求解公式n n 取取取取x xk k33,x xk k22,x xk k11,x xk k，在，在，在，在 x xk k33,x xk k 上作上作上作上作F F(x x)=)=f f(x x,y y(x x)的插值多项式的插值多项式的插值多项式的插值多项式n n用上述方法可导出：用上述方法可导出：用上述方法可导出：用上述方法可导出：n n (3.8)(3.8)n n (3.9)(3.9)n n称为亚当斯隐式求解公式称为亚当斯隐式求解公式称为亚当斯隐式求解公式称为亚当斯隐式求解公式.第60页/共73页3 线性多步法线性

49、多步法n n3.1 开型求解公式n n2.2.亚当斯隐式求解公式亚当斯隐式求解公式n n (3.8)(3.8)n n (3.9)(3.9)n n称为亚当斯隐式求解公式称为亚当斯隐式求解公式称为亚当斯隐式求解公式称为亚当斯隐式求解公式.n n注：利用注：利用注：利用注：利用4 4阶公式阶公式阶公式阶公式(3.6)(3.6)显化之：显化之：显化之：显化之：n n (3.10)(3.10)n n称称称称(3.10)(3.10)为亚当斯预测为亚当斯预测为亚当斯预测为亚当斯预测校正系统校正系统校正系统校正系统.第61页/共73页3 线性多步法线性多步法n n3.2 闭型求解系统n n 将将将将 x xk

50、 k11,x xk k 扩充为扩充为扩充为扩充为 x xk k44,x xk k，取，取，取，取x xk k44，x xk k33，x xk k22，x xk k11为节点，作为节点，作为节点，作为节点，作F F(x x)=)=f f(x x,y y(x x)的牛顿前插多项式的牛顿前插多项式的牛顿前插多项式的牛顿前插多项式.n n n n则则则则 第62页/共73页3 线性多步法线性多步法n n3.2 闭型求解系统n n 将将将将 x xk k11,x xk k 扩充为扩充为扩充为扩充为 x xk k44,x xk k，取，取，取，取x xk k44，x xk k33，x xk k22，x x