计算机数值方法-第五章-常微分方程数值解法-课件

第五章

常微分方程数值解法

(NumericalMethodsforOrdinaryDifferentialEquations)引言1ppt课件引言

考虑一阶常微分方程的初值问题

(Initial-ValueProblem):只要f(x,y)在x∈

[a,b]上连续，且关于y

满足Lipschitz

条件，即存在与x,y无关的常数L

使对任意定义在[a,b]上的y1(x)和y2(x)都成立，则上述常微分方程存在唯一解。2ppt课件引言数值解法就是要计算出解函数y(x)在一系列节点a=x0<x1<…<xn=b处的近似值

y0

y1…yn节点间距称为步长，通常采用等距节点，即取hi=h

(常数)。解决：

数值解法的一个基本特点是“步进式”，即求解时顺着节点排列的次序一步步地向前推进。单步：

yk-1

yk多步：

yk-p…

yk-2

，yk-1

yk注意：与“迭代法”区别3ppt课件第五章常微分方程数值解法欧拉方法（Euler’sMethod

）4ppt课件显式Euler公式由两点公式求导数，在[xj,xj+1]子区间上有：其中j[xj,xj+1]代入方程有显式Euler公式x0x1y(x1)y(x0)h5ppt课件显式Euler公式的误差局部截断误差显式Euler公式定义在假设yj=y(xj)，即第

i

步计算是精确的前提下，考虑的截断误差ej+1=y(xj+1)

yj+1

称为局部截断误差

/*localtruncationerror*/。6ppt课件（P169）例1：取h=0.1，分别用显式Euler法、隐式Euler法、显隐结合的预测校正系统求解初值问题解：用显式Euler法求解，有：…依次下去计算结果见P169

7ppt课件隐式Euler公式由两点公式求导数，在[xj,xj+1]子区间上有：其中j[xj,xj+1]代入方程有局部截断误差隐式Euler公式是一个关于yj+1的方程,要从中解出yj+1

x0x18ppt课件（P169）例1：取h=0.1，分别用显式Euler法、隐式Euler法、显隐结合的预测校正系统求解初值问题计算过程见P170

用隐式Euler法求解，有：解：从中解出yi+1，有：9ppt课件显隐结合的预测校正系统——避免求解方程

（predictor-correctormethod）Step1:

先用显式欧拉公式作预测，算出预测值Step2:再用隐式欧拉公式作校正，得到校正值写成一个公式为：y0

＝y(a)－>y1－>y1－>y2－>y2…－>yn－>yn计算顺序：10ppt课件（P170）例1：取h=0.1，分别用显式Euler法、隐式Euler法、显隐结合的预测校正系统求解初值问题……计算过程见P170

解：用预测校正系统求解，有：y0=111ppt课件p阶精度若某算法的局部截断误差e(h)满足：e(h)=

O(hp+1)，即有：

e(h)/hp+1=c(常数），则称该算法有p

阶精度。定义

欧拉法的局部截断误差：欧拉法具有1阶精度。显式：隐式：12ppt课件梯形公式（trapezoidformula）注：的确有局部截断误差，即梯形公式具有2阶精度，比欧拉方法有了进步。但注意到该公式是隐式公式。从显式、隐式Euler法和Euler局部截断误差来看似乎可以有如下式子——梯形公式：隐式Euler显式Euler13ppt课件基于数值积分的求解思想求出

积分的近似表达，也就求出了y(xk)hxkxf(x,y)f(x,y)xk-1f(xk,yk)f(xk-1,yk-1)14ppt课件基于数值积分的求解思想用矩形求积公式代替hxkxf(x,y)f(x,y)xk-1f(xk,yk)f(xk-1,yk-1)左矩形显式Euler公式右矩形隐式Euler公式15ppt课件基于数值积分的求解思想用梯形求积公式代替hxkxf(x,y)f(x,y)xk-1f(xk,yk)f(xk-1,yk-1)注：梯形公式是隐式公式显式、隐式Euler的平均梯形公式具有2阶精度余项16ppt课件方法显式欧拉隐式欧拉梯形公式计算简单精度低稳定性好精度低,计算不便精度提高计算不便

Can’tyougivemeaformulawithalltheadvantagesyetwithoutanyofthedisadvantages?

Letmetry！Euler公式及梯形公式总结17ppt课件改进的Euler法——预测校正系统

（modifiedEuler’smethod）注：可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单。它的稳定性高于显式欧拉法。误差比较：梯形法≈改进Euler法<显、隐Euler法Step1:

先用显式欧拉公式作预测，算出预测值Step2:再用梯形公式作校正，得到校正值18ppt课件举例

(P173例3)：取h=0.1，分别用Euler法、梯形公式、改进的Euler法求解初值问题解：计算过程参见P173Euler法梯形法解出yk得：19ppt课件举例

(P173例3)：取h=0.1，分别用Euler法、梯形公式、改进的Euler法求解初值问题解：计算过程参见P173改进Euler法20ppt课件Simpson公式（Simpsonformula

）注：Simpson公式是一个隐式公式

Simpson公式是多步法：yk-1,yk-2=>yk应用simpson求积公式：Simpson公式——4阶精度21ppt课件P198习题6用Euler法和改进的Euler法解方程，h＝0.1注意改题目为0<=x<=0.4第五章作业欧拉方法（Euler’sMethod

）22ppt课件第五章常微分方程数值解法龙格－库塔法（Runge－Kutta

Method）23ppt课件二阶Runge－Kutta法——改进的Euler法

考察改进的欧拉法，可以将其改写为：K1K2xi+h步长一定是一个h

吗？斜率一定取K1K2的平均值吗？24ppt课件龙格-库塔法（Runge－Kutta

Method

）首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得

Step1:将K2在(xi,yi)

点作Taylor展开将改进欧拉法推广为：),(),(][12122111phKyphxfKyxfKKKhyyiiiiii++==++=+ll25ppt课件龙格-库塔法（Runge－Kutta

Method

）Step2:将K1

K2代入第1式，得到Step3:将yi+1与y(xi+1)在xi点的泰勒展开作比较要求，则必须有：26ppt课件龙格-库塔法（Runge－Kutta

Method

）这里有

个未知数，

个方程。32存在无穷多个解。所有满足上式的格式统称为2阶龙格-库塔格式。注意:就是改进的欧拉法,也就是2阶R-K法。:

为获得更高的精度，应该如何进一步推广？27ppt课件三阶龙格-库塔法三、四阶R-K法思想：利用各种显化方法，对以上公式显化Simpson公式求积xixi+hxi+h/228ppt课件三阶龙格-库塔法K1K2K3显EulerP177要计算3个函数可以证明(Taylor展开K1，K2，K3）e(h)=O(h4)具有3阶精度29ppt课件

P180例：取h=0.1，用三阶、四阶R-K法求解初值问题举例解：用三阶R-K法30ppt课件取n=1计算K1，K2，K3，=>y1取n=2计算K1，K2，K3，=>y2取n=3计算K1，K2，K3，=>y3n12K1……K2……K3……y1……列表计算31ppt课件四阶龙格-库塔法(经典R-K法)

/*ClassicalRunge-KuttaMethod*/

e(h)=O(h5)4阶精度32ppt课件优缺点：优点是：(1)都是一步法，因此只要给定一个初始值就可以一直计算下去；(2)精度相对较高，如经典R-K法为四阶精度缺点是：(特别是三阶和四阶法)计算量较大。其他问题讨论：R-K法推导基于Taylor展式，因而要求y(x)有较好的光滑性（即有高阶导数）。最常用的是四阶公式，它适用于一般的问题，准确、稳定、易于编程。步长h减小，局部误差O(h5)减小但步数增加，舍入误差积累增加h要适当，总误差才最小四阶龙格-库塔法(经典R-K法)

/*ClassicalRunge-KuttaMethod*/

33ppt课件

P180例：取h=0.1，用三阶、

四阶R-K法求解初值问题解：用三阶R-K法34ppt课件取n=1计算K1，K2，K3，K4

=>y1取n=2计算K1，K2，K3，K4

=>y2取n=3计算K1，K2，K3，K4

=>y3n1K1…K2…K3…K4…y1…列表计算35ppt课件例：取h=0.2，用四阶龙格－库塔法求解初值问题举例解：这里，经典的四阶龙格-库塔公式为36ppt课件举例表中列出了计算结果，同时列出了相应的精确解.比较本章第一个例子的计算结果，显然龙格-库塔方法的精度高.37ppt课件名称RK4（A,B,Y0,H）存储：K[1]~K[4]<=K1~K4Y[N]<=y1,y2…yn算法：(1)

输入a，b，n，y0，设计函数F(x,y)(2)Y[0]<=y0,H=(b-a)/n，H1<=H/2

(3)对i=0，1，2，…,N-1(由Y[i]计算Y[i+1])xi=x0+i*hK[1]<=F(xi,Y[i])K[2]<=F(xi+H1,Y[i]+H1*K[1])K[3]<=F(xi+H1,Y[i]+H1*K[2])K[4]<=F(xi+H,Y[i]+H*K[3])Y[i+1]<=Y[i]+(K[1]+2*K[2]+2*K[3]+K[4])*h/6(4)输出Y[1],Y[2],…,Y[N].四阶龙格-库塔法算法设计（P182)38ppt课件P199习题10（2）第五章作业龙格－库塔法（Runge－Kutta

Method）39ppt课件第五章常微分方程数值解法线性多步法（MultistepMethod）40ppt课件线性多步法（MultistepMethod）单步法：

yk-1

yk对应xk：a=x0<x1<…<xn=b的yk：y0y1…yn优点：不需要其他计算值，计算yk只要用到yk-1缺点：没有充分利用前几步得到的信息处理y(xk-1+ah)的困难多步法：

yk-p…

yk-2

，yk-1

yk用到p步计算的结果41ppt课件开型求解——Adams法思想：计算积分时，插值节点取在积分区间[xk-1,xk]的外部三阶Adams：显式：取xk-1,

xk-2,

xk-3

作节点f(x,y)≈L2(x)误差O(h4)隐式：取xk,xk-1,

xk-2

作节点f(x,y)≈L2(x)误差O(h4)四阶Adams：显式：取xk-1,

xk-2,

xk-3,

xk-4

作节点f(x,y)≈L3(x)误差O(h5)隐式：取xk,

xk-1,

xk-2,

xk-3

作节点f(x,y)≈L3(x)误差O(h5)预测－校正系统：显式Adams＋隐式Adams（三阶、四阶）42ppt课件三阶显式Adams法取xk-1,

xk-2,

xk-3

作节点f(x,y)≈L2(x)误差O(h4)简记为：误差：3阶精度43ppt课件三阶隐式Adams法取xk,xk-1,

xk-2

作节点f(x,y)≈L2(x)误差O(h4)误差：3阶精度44ppt课件四阶显式Adams法取xk-1,

xk-2,

xk-3,

xk-4

作节点f(x,y)≈L3(x)误差O(h5)误差：4阶精度45ppt课件四阶隐式Adams法取xk,

xk-1,

xk-2,

xk-3

作节点f(x,y)≈L3(x)误差O(h5)误差：4阶精度46ppt课件亚当姆斯预测-校正系统

（Adamspredictor-correctorsystem）四阶误差：三阶显式预测隐式校正四阶其中：47ppt课件亚当姆斯预测-校正系统

（Adamspredictor-correctorsystem）Step1:

用Runge-Kutta

法计算前k个初值；Step2:

用Adams显式计算预测值；Step3:

用同阶Adams隐式计算校正值。注意：三步所用公式的精度必须相同。通常用经典Runge-Kutta法（4阶）配合4阶Adams公式。四阶精度48ppt课件闭型求解——Milne公式闭型求解思想：扩大积分区间：[xk-1,xk]－>[xk-4,xk]增加内部节点：xk-1,

xk-2,

xk-3,

xk-4原来的外部节点－>积分区间的内部节点方程形式变化：49ppt课件闭型求解——Milne公式推导：(1)取xk-1,

xk-2,

xk-3,

xk-4为节点作三次Newton前插公式代替f(x,y(x))：f(x,y(x))=N3(x)+R4(x)(2)x＝xk-4+th则x∈[xk-4,xk]时，t∈[0,4]作x->t替换(3)计算积分及误差截断误差：4阶精度Milne公式50ppt课件闭型求解公式Milne公式不能直接使用，一般是首先用四阶R－K法求出前几个节点的值，然后使用Milne公式求解。前面知道Simpson公式（四阶精度）是隐式积分公式，所以利用Milne公式，可以得到如下的精度