常微分方程数值解法1课件

数数 值值 分分 析析第5章 常微分方程数值解法1111111111112222222222223333333333331 引言v1.0 基本概念基本概念1.常微分方程的初值问题：称为具有初值称为具有初值(1.2)的常微分方程的常微分方程.若若f(x,y)在在a x b,|y|+上上连连续续，且且关关于于y满满足足Lip条件：条件：常数常数L使使|f(x,y1)f(x,y2)|L|y1 y2|则初值问题则初值问题(1.1)(1.2)存在唯一连续可微解存在唯一连续可微解y(x).注：以下总假设注：以下总假设f 满足满足Lip条件条件.1111111111112222222222223333333333331 引言v1.0 基本概念基本概念1.常微分方程的初值问题：称为具有初值称为具有初值(1.2)的常微分方程的常微分方程.(1.1)(1.2)等价于微分方程：等价于微分方程：(1.3)注：一般无初等解注：一般无初等解(解析解解析解)，即使有形式也复杂，即使有形式也复杂.1111111111112222222222223333333333331 引言v1.0 基本概念基本概念2.初值问题的数值解 设设(1.1)(1.2)的解的解y(x)在节点在节点xi处的近似解值为处的近似解值为 yi y(xi),a x1 x2 xn=b则则称称yi (i=1,2,n)为为(1.1)(1.2)的的数数值值解解，又又称称y(xi)的计算值的计算值.1111111111112222222222223333333333331 引言v1.0 基本概念基本概念3.数值方法 两种转化：两种转化：由微分出发的数值方法由微分出发的数值方法.由积分由积分 出发的数值方法出发的数值方法.计算方法计算方法 步进法：从初始条件出发，逐步求步进法：从初始条件出发，逐步求y1,y2,yn.又有两种：单步法，多步法又有两种：单步法，多步法.注：采用等距节点：注：采用等距节点：1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式基于数值微分的求解公式.(1.6)1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式基于数值微分的求解公式.1.前进欧拉公式 (1.6)的前半部分为：的前半部分为：令令 yi+1=yi+hf(xi,yi)(1.7)其中其中yi=y(xi),则则yi+1 y(xi+1)1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式基于数值微分的求解公式.1.前进欧拉公式 令令 yi+1=yi+hf(xi,yi)(1.7)其中其中yi=y(xi),则则yi+1 y(xi+1)记记 (1.8)则则称称(1.7)为为前前进进欧欧拉拉求求解解公公式式.简简称称为为欧欧拉拉公公式式或或欧欧拉拉法法.(1.8)称为欧拉公式的余项：称为欧拉公式的余项：ei+1(h)=y(xi+1)yi+1 1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式基于数值微分的求解公式.2.后退欧拉公式 (1.6)的后半部分的后半部分令令 yi+1=yi+hf(xi+1,yi+1)(1.9)其中其中yi=y(xi),则则yi+1 y(xi+1)1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式基于数值微分的求解公式.2.后退欧拉公式令令 yi+1=yi+hf(xi+1,yi+1)(1.9)其中其中yi=y(xi),则则yi+1 y(xi+1)注：注：(1.9)中中f(xi+1,yi+1)f(xi+1,y(xi+1)余项余项 (1.10)1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式基于数值微分的求解公式.2.后退欧拉公式令令 yi+1=yi+hf(xi+1,yi+1)(1.9)其中其中yi=y(xi),则则yi+1 y(xi+1)注：注：称称(1.9)为后退欧拉公式为后退欧拉公式(后退欧拉法后退欧拉法).称称(1.10)为后退欧拉法的误差近似值为后退欧拉法的误差近似值.欧拉法与后退欧拉公式的区别：欧拉法与后退欧拉公式的区别：(1.7)为直接计算公式称显式公式为直接计算公式称显式公式.(1.9)为关于函数方程称隐式公式为关于函数方程称隐式公式.1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式基于数值微分的求解公式.【例例1】取取h=0.1求解初值问题：求解初值问题：(1.11).解解：，xi=ih=0.1 i,(i=0,1,2,10)欧拉法：欧拉法：1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式基于数值微分的求解公式.【例例1】取取h=0.1求解初值问题：求解初值问题：(1.11).解解：，xi=ih=0.1 i,(i=0,1,2,10)后退欧拉法：后退欧拉法：1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式基于数值微分的求解公式.注注：为为避避免免求求解解函函数数方方程程，采采用用显显式式与与隐隐式式结结合合的的方方法：法：此方法称为此方法称为 预测预测校正系统校正系统.求解过程为：求解过程为：1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式基于数值微分的求解公式.预测预测校正系统：校正系统：【例例2】利用预测利用预测校正系统求解例校正系统求解例1.1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式基于数值微分的求解公式.预测预测校正系统：校正系统：注注：显显式式比比隐隐式式方方便便，但但有有时时隐隐式式效效果果比比显显式式好好.(4介绍介绍).1111111111112222222222223333333333331 引言v1.2 截断误差截断误差定定义义1.1 称称ek(h)=y(xk)yk为为计计算算yk的的公公式式第第k步步的的局局部截断误差部截断误差.注注：“局局部部”是是指指在在计计算算第第k步步时时，假假定定前前面面yi=y(xi)(i k).而而yk y(xk)欧拉法欧拉法.后退欧拉法后退欧拉法.一般根据一般根据y(xk)对对y(k),y(k)做估计做估计.1111111111112222222222223333333333331 引言v1.2 截断误差截断误差定定义义1.2 设设ei(h)(i=1,2,k)为为求求解解公公式式第第i步步的的局局部部截断误差截断误差.称称为该求解公式在点上的整体截断误差为该求解公式在点上的整体截断误差.注：注：局部截断误差局部截断误差ek(h)与与yk有关有关.整体截断误差整体截断误差Ek(h)与与y1,y2,yk有关有关.所有所有ek(h)都与都与h有关有关.1111111111112222222222223333333333331 引言v1.2 截断误差截断误差定定义义1.3 若若局局部部截截断断误误差差e(h)=O(hp+1)，则则称称该该求求解解公公式具有式具有p阶精度阶精度.注：欧拉法具有一阶精度注：欧拉法具有一阶精度.(精度越高越好精度越高越好)1111111111112222222222223333333333331 引言作业作业 P208 1，2，3.1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式基于数值积分的求解公式 (1.13)若已知若已知y(xk)=yk,则计算积分可求出则计算积分可求出y(xk+1).如用矩形公式求积分如用矩形公式求积分则有则有y(xk+1)=y(xk)+hf(xk,yk)令令yk+1=y(xk)+hf(xk,yk)即即为为欧欧拉拉公公式式.故故欧欧拉拉公公式式又又称矩形法称矩形法.1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式基于数值积分的求解公式 (1.13)考虑考虑1.梯形公式记记 (1.14)1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式基于数值积分的求解公式1.梯形公式记记 (1.14)称称(1.14)为梯形为梯形(求解求解)公式公式.简称梯形法简称梯形法.1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式基于数值积分的求解公式1.梯形公式梯形梯形(求解求解)公式公式,简称梯形法简称梯形法:(1.14)注：注：梯形公式的余项：梯形公式的余项：故是二阶精度故是二阶精度.111111111111222222222222333333333333v1.3 基于数值积分的求解公式基于数值积分的求解公式1.梯形公式 (1.14)梯形公式为隐式公式梯形公式为隐式公式.预测预测校正系统校正系统 (1.15)称称(1.15)为改进的欧拉公式，也可记为为改进的欧拉公式，也可记为1 引言1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式基于数值积分的求解公式1.梯形公式 (1.14)可以证明，改进欧拉公式也具有二阶精度可以证明，改进欧拉公式也具有二阶精度.1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式基于数值积分的求解公式【例例3】用欧拉法，梯形法以及改进欧拉法求解用欧拉法，梯形法以及改进欧拉法求解取取h=0.1.计算到计算到x=0.5.解解：f(x,y)=xy+1,a=x0=0,b=0.5,y0=1,n=5(Euler法法)求解公式：求解公式：yk=yk1+h(xk1yk1+1)=hxk1+(1 h)yk1+h=0.1xk1+0.9yk1+0.1 1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式基于数值积分的求解公式【例例3】用欧拉法，梯形法以及改进欧拉法求解用欧拉法，梯形法以及改进欧拉法求解解解：f(x,y)=xy+1,a=x0=0,b=0.5,y0=1,n=5(梯形法梯形法)求解公式：求解公式：yk=yk1+h(xk1yk1+1)+(xkyk+1)/2解出解出yk，得，得方程方程1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式基于数值积分的求解公式【例例3】用欧拉法，梯形法以及改进欧拉法求解用欧拉法，梯形法以及改进欧拉法求解解解：f(x,y)=xy+1,a=x0=0,b=0.5,y0=1,n=5(改进改进Euler法法)求解公式：求解公式：yk=yk1+h(xk1yk1+1)+xk (yk+h(xkyk+1)+1/2得得=0.905yk1+0.045xk1+0.05xk+0.095方程方程1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式基于数值积分的求解公式2.辛卜生公式 记记 (1.17)1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式基于数值积分的求解公式2.辛卜生公式记记 (1.17)其余项其余项1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式基于数值积分的求解公式2.辛卜生公式记记 (1.17)将将xk1,xk 对分：对分：调整下标为调整下标为xi2,xi：xi2=xk1,xi1=xk1+h1,xi=xk1+2h1=xk则则(1.17)化为化为 (1.19)称称(1.19)为为辛辛卜卜生生求求解解公公式式，其其中中fk2=f(xk2,y(xk2)，fk1=f(xk1,y(xk1)，fk=f(xk,y(xk)1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式基于数值积分的求解公式2.辛卜生公式记记 (1.17)(1.19)称称(1.19)为为辛辛卜卜生生求求解解公公式式，其其中中fi2=f(xi2,y(xi2)，fi1=f(xi1,y(xi1)，fi=f(xi,y(xi)注：注：(1.19)的误差：的误差：1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式基于数值积分的求解公式2.辛卜生公式记记 (1.17)(1.19)称称(1.19)为为辛辛卜卜生生求求解解公公式式，其其中中fi2=f(xi2,y(xi2)，fi1=f(xi1,y(xi1)，fi=f(xi,y(xi)注：注：隐式隐式(需显化需显化)多步多步将在将在3中讨论中讨论.1111111111112222222222223333333333332 Runge-Kutta法v2.0 原理原理 其其中中K =f(,y()=y()称称为为y在在xi1,xi上上的的平平均均斜率斜率.欧拉法：欧拉法：改进欧拉法：改进欧拉法：(2.1)1111111111112222222222223333333333332 Runge-Kutta法v2.0 原理原理 其其中中K =f(,y()=y()称称为为y在在xi1,xi上上的的平平均均斜率斜率.对对(1.17)显化：显化：辛卜生：辛卜生：(2.4)1111111111112222222222223333333333332 Runge-Kutta法v2.0 原理原理其其中中K =f(,y()=y()称称为为y在在xi1,xi上上的的平平均均斜率斜率.设设想想：在在中中多多计计算算(预预测测)几几个个点点上上的的值值然然后后可可加加权权取取平均值作为的近似值可能构成更高阶的公式平均值作为的近似值可能构成更高阶的公式.一阶一阶二阶二阶三阶三阶1111111111112222222222223333333333332 Runge-Kutta法v2.1 Runge-Kutta公式公式 (*)其中其中0 j 1，yi1+jh是是y(xi1+jh)的预测值的预测值.称称(*)为为R-K公式公式注：注：(2.1)(2.4)分别称为二阶，三阶分别称为二阶，三阶R-K公式公式.j，j，j为待定系数为待定系数.使使(*)的阶数尽量高的阶数尽量高.1111111111112222222222223333333333332 Runge-Kutta法v2.1 Runge-Kutta公式公式参数的确定，以参数的确定，以m=2为例为例.欲求欲求 1，2，2.原则原则:使使ei(h)=y(xi)yi的阶数尽可能高的阶数尽可能高1111111111112222222222223333333333332 Runge-Kutta法v2.1 Runge-Kutta公式公式展开展开展开展开 原则原则:使使ei(h)=y(xi)yi的阶数尽可能高的阶数尽可能高1111111111112222222222223333333333332 Runge-Kutta法v2.1 Runge-Kutta公式公式 原则原则:使使ei(h)=y(xi)yi的阶数尽可能高的阶数尽可能高1111111111112222222222223333333333332 Runge-Kutta法v2.1 Runge-Kutta公式公式 欲欲求求截截断断误误差差ei(h)=y(xi)yi关关于于h的的阶阶数数尽尽可可能能高高，应使应使无穷多解，从而有许多无穷多解，从而有许多2阶阶R-K公式公式1111111111112222222222223333333333332 Runge-Kutta法v2.1 Runge-Kutta公式公式应使应使注：注：取取 1=2=1/2，2=1，即为改进欧拉公式，即为改进欧拉公式.1111111111112222222222223333333333332 Runge-Kutta法v2.1 Runge-Kutta公式公式应使应使注：注：取取 1=0，2=1，2=1/2，即为中点公式，即为中点公式1111111111112222222222223333333333332 Runge-Kutta法v2.1 Runge-Kutta公式公式应使应使注注：二二阶阶R-K公公式式的的截截断断误误差差为为故故为为二二阶阶方方法法.相相仿仿可得更高阶的可得更高阶的R-K公式公式.1111111111112222222222223333333333332 Runge-Kutta法v2.2 经典经典R-K公式公式 在在4解解R-K公公式式中中最最重重要要的的是是经经典典R-K公公式式.(2.6)注：注：(2.6)为为4阶方法阶方法.1111111111112222222222223333333333332 Runge-Kutta法v2.2 经典经典R-K公式公式 在在4解解R-K公公式式中中最最重重要要的的是是经经典典R-K公公式式.(2.6)注：注：R-K法对法对4阶以上不一定能提高整数阶阶以上不一定能提高整数阶.1111111111112222222222223333333333332 Runge-Kutta法v2.2 经典经典R-K公式公式【例例4】使用三阶，四阶使用三阶，四阶R-K法求解初值问题：法求解初值问题：的部分计算值的部分计算值y1，y2，y3，其中，其中h=0.1.解解 使用三阶使用三阶R-K法法1111111111112222222222223333333333332 Runge-Kutta法【例例4】使用三阶，四阶使用三阶，四阶R-K法求解初值问题：法求解初值问题：的部分计算值的部分计算值y1，y2，y3，其中，其中h=0.1.解解 使用三阶使用三阶R-K法法1111111111112222222222223333333333332 Runge-Kutta法【例例4】使用三阶，四阶使用三阶，四阶R-K法求解初值问题：法求解初值问题：的部分计算值的部分计算值y1，y2，y3，其中，其中h=0.1.解解 使用四阶使用四阶R-K法法1111111111112222222222223333333333332 Runge-Kutta法【例例4】使用三阶，四阶使用三阶，四阶R-K法求解初值问题：法求解初值问题：的部分计算值的部分计算值y1，y2，y3，其中，其中h=0.1.解解 使用四阶使用四阶R-K法法1111111111112222222222223333333333332 Runge-Kutta法注注 使使用用R-K法法要要求求具具备备较较好好的的光光滑滑性性，否否则则效效果果不不如如低阶的低阶的.作业作业P209 8 9，10.1111111111112222222222223333333333333 线性多步法单步法的优点：简单，计算单步法的优点：简单，计算yk+1只用只用yk.缺点缺点:没有充分利用前面的信息且计算没有充分利用前面的信息且计算y(xk+h)较困难较困难回回顾顾Simpson：(1.19)考考虑虑：(3.1)两种插值求积：两种插值求积：将将xk1,xk增增加加内内部部节节点点，改改为为xk2,xk导导出出的的公公式式称为闭型求解公式称为闭型求解公式.线性多步线性多步1111111111112222222222223333333333333 线性多步法考考虑虑：(3.1)两种插值求积：两种插值求积：将将xk1,xk增增加加内内部部节节点点，改改为为xk2,xk导导出出的的公公式式称为闭型求解公式称为闭型求解公式.在在xk1,xk外外增增加加插插值值节节点点，导导出出的的公公式式称称为为开开型型求求解公式解公式.开型有显和隐，闭型也有显和隐开型有显和隐，闭型也有显和隐.1111111111112222222222223333333333333 线性多步法v3.1 开型求解公式开型求解公式1.亚当斯显式求解公式 取取节节点点xk3,xk2,xk1，在在xk3,xk上上作作F(x)=f(x,y(x)的插值多项式的插值多项式.1111111111112222222222223333333333333 线性多步法v3.1 开型求解公式开型求解公式1.亚当斯显式求解公式 取节点取节点xk3,xk2,xk1，在，在xk3,xk上上记记xki=xk ih,x=xk+th，则，则 1111111111112222222222223333333333333 线性多步法v3.1 开型求解公式开型求解公式1.亚当斯显式求解公式 取节点取节点xk3,xk2,xk1，记，记xki=xk ih,x=xk+th，则，则 代入代入(3.1)得得1111111111112222222222223333333333333 线性多步法v3.1 开型求解公式开型求解公式1.亚当斯显式求解公式 取节点取节点xk3,xk2,xk1，记，记xki=xk ih,x=xk+th，则，则 令令 (3.4)称称(3.4)为亚当斯显式求解公式为亚当斯显式求解公式(线性多步线性多步).1111111111112222222222223333333333333 线性多步法v3.1 开型求解公式开型求解公式1.亚当斯显式求解公式 取节点取节点xk3,xk2,xk1，记，记xki=xk ih,x=xk+th，则，则余项：余项：从而从而(3.4)具有具有3阶精度阶精度.称为称为3阶亚当斯求解公式阶亚当斯求解公式.1111111111112222222222223333333333333 线性多步法v3.1 开型求解公式开型求解公式1.亚当斯显式求解公式类似地取类似地取xk4,xk3,xk2,xk1 在在xk4,xk上上作作F(x)=f(x,y(x)的的插插值值多多项项式式，可可导导出出4阶亚当斯显式求解公式：阶亚当斯显式求解公式：(3.6)(3.7)4阶精度阶精度1111111111112222222222223333333333333 线性多步法v3.1 开型求解公式开型求解公式2.亚当斯隐式求解公式 取取xk3,xk2,xk1,xk，在在xk3,xk上上作作F(x)=f(x,y(x)的插值多项式的插值多项式用上述方法可导出：用上述方法可导出：(3.8)(3.9)称为亚当斯隐式求解公式称为亚当斯隐式求解公式.1111111111112222222222223333333333333 线性多步法v3.1 开型求解公式开型求解公式2.亚当斯隐式求解公式 (3.8)(3.9)称为亚当斯隐式求解公式称为亚当斯隐式求解公式.注：利用注：利用4阶公式阶公式(3.6)显化之：显化之：(3.10)称称(3.10)为亚当斯预测为亚当斯预测校正系统校正系统.1111111111112222222222223333333333333 线性多步法v3.2 闭型求解系统闭型求解系统 将将xk1,xk扩扩充充为为xk4,xk，取取xk4，xk3，xk2，xk1为节点，作为节点，作F(x)=f(x,y(x)的牛顿前插多项式的牛顿前插多项式.则则 1111111111112222222222223333333333333 线性多步法v3.2 闭型求解系统闭型求解系统 将将xk1,xk扩扩充充为为xk4,xk，取取xk4，xk3，xk2，xk1为节点，作为节点，作F(x)=f(x,y(x)的牛顿前插多项式的牛顿前插多项式.则则 令令x=xk+(t 4)h 则则 1111111111112222222222223333333333333 线性多步法v3.2 闭型求解系统闭型求解系统令令x=xk+(t 4)h 则则由由1111111111112222222222223333333333333 线性多步法v3.2 闭型求解系统闭型求解系统令令 (3.11)称称(3.11)为米尔恩求解公式为米尔恩求解公式(Miline).余项：余项：1111111111112222222222223333333333333 线性多步法v3.2 闭型求解系统闭型求解系统 (3.11)称称(3.11)为米尔恩求解公式为米尔恩求解公式(Miline).余项：余项：分分别别在在x=xk处处展展开开y(xk 4h)和和y(xk ih)为为台台劳劳级级数并整理得数并整理得 (3.12)1111111111112222222222223333333333333 线性多步法v3.2 闭型求解系统闭型求解系统 (3.11)余项：余项：(3.12)注：注：辛卜生公式求解辛卜生公式求解 (3.13)1111111111112222222222223333333333333 线性多步法v3.2 闭型求解系统闭型求解系统 (3.11)余项：余项：(3.12)注注：多多步步法法的的开开始始须须用用单单步步法法，使使用用高高精精度度多多步步法法时，应用高精度单步法计算开始值时，应用高精度单步法计算开始值.1111111111112222222222223333333333333 线性多步法v3.2 闭型求解系统闭型求解系统【例例5】初值问题初值问题设设h=0.2.已已知知y1=0.181，分分别别用用3阶阶亚亚当当斯斯和和米米尔尔恩恩公公式计算式计算x3=0.6和和x4=0.8的计算值的计算值.解解 用用3阶亚当斯计算阶亚当斯计算y3 y(0.6)3阶亚当斯公式：阶亚当斯公式：1111111111112222222222223333333333333 线性多步法v3.2 闭型求解系统闭型求解系统【例例5】初值问题初值问题设设h=0.2.已已知知y1=0.181，分分别别用用3阶阶亚亚当当斯斯和和米米尔尔恩恩公公式计算式计算x3=0.6和和x4=0.8的计算值的计算值.解解 用用3阶亚当斯计算阶亚当斯计算y3 y(0.6)先用先用3阶阶R-K法求法求y2 y(0.4)，然后求，然后求y31111111111112222222222223333333333333 线性多步法v3.2 闭型求解系统闭型求解系统【例例5】初值问题初值问题设设h=0.2.已已知知y1=0.181，分分别别用用3阶阶亚亚当当斯斯和和米米尔尔恩恩公公式计算式计算x3=0.6和和x4=0.8的计算值的计算值.解解 用米尔恩求用米尔恩求y4 y(0.8)4阶米尔恩公式阶米尔恩公式1111111111112222222222223333333333333 线性多步法v3.2 闭型求解系统闭型求解系统【例例5】初值问题初值问题设设h=0.2.已已知知y1=0.181，分分别别用用3阶阶亚亚当当斯斯和和米米尔尔恩恩公公式计算式计算x3=0.6和和x4=0.8的计算值的计算值.解解 用米尔恩求用米尔恩求y4 y(0.8)4阶米尔恩公式阶米尔恩公式