常微分方程数值解法

数值分析第5章常微分方程数值解法§1引言1.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条件.§1引言1.0基本概念1.常微分方程的初值问题：称为具有初值(1.2)的常微分方程.②(1.1)(1.2)等价于微分方程：

(1.3)注：一般无初等解(解析解)，即使有形式也复杂.§1引言1.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)的计算值.§1引言1.0基本概念3.数值方法①两种转化：由微分出发的数值方法.

由积分出发的数值方法.②计算方法步进法：从初始条件出发，逐步求y1,y2,…,yn.

又有两种：单步法，多步法.注：采用等距节点：§1引言1.1基于数值微分的求解公式.

(1.6)§1引言1.1基于数值微分的求解公式.1.前进欧拉公式

(1.6)的前半部分为：令yi+1=yi

+

hf(xi,yi)

(1.7)其中yi=y(xi),则yi+1

y(xi+1)§1引言1.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

§1引言1.1基于数值微分的求解公式.2.后退欧拉公式

(1.6)的后半部分令yi+1=yi

+hf(xi+1,yi+1)

(1.9)其中yi=y(xi),则yi+1

y(xi+1)

§1引言1.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)§1引言1.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)为关于函数方程称隐式公式.§1引言1.1基于数值微分的求解公式.【例1】取h=0.1求解初值问题：

(1.11).解：，xi=ih=0.1

i,(i=0,1,2,…,10)

①欧拉法：§1引言1.1基于数值微分的求解公式.【例1】取h=0.1求解初值问题：

(1.11).解：，xi=ih=0.1

i,(i=0,1,2,…,10)

②后退欧拉法：

§1引言1.1基于数值微分的求解公式.注：为避免求解函数方程，采用显式与隐式结合的方法：

此方法称为预测——校正系统.求解过程为：§1引言1.1基于数值微分的求解公式.预测——校正系统：【例2】利用预测——校正系统求解例1.§1引言1.1基于数值微分的求解公式.预测——校正系统：注：显式比隐式方便，但有时隐式效果比显式好.(§4介绍).§1引言1.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)做估计.§1引言1.2截断误差定义1.2设ei(h)(i=1,2,…,k)为求解公式第i步的局部截断误差.称为该求解公式在点上的整体截断误差.注：①局部截断误差ek(h)与yk有关.

整体截断误差Ek(h)与y1,y2,…,yk有关.②所有ek(h)都与h有关.§1引言1.2截断误差定义1.3若局部截断误差e(h)=O(hp+1)，则称该求解公式具有p阶精度.注：欧拉法具有一阶精度.(精度越高越好)§1引言作业P2081，2，3.§1引言1.3基于数值积分的求解公式

(1.13)若已知y(xk)=yk,则计算积分可求出y(xk+1)

.

如用矩形公式求积分则有y(xk+1)=y(xk)+hf(xk,yk)令yk+1=y(xk)+hf(xk,yk)即为欧拉公式.故欧拉公式又称矩形法.§1引言1.3基于数值积分的求解公式

(1.13)考虑1.梯形公式记(1.14)§1引言1.3基于数值积分的求解公式1.梯形公式记(1.14)称(1.14)为梯形(求解)公式.简称梯形法.§1引言1.3基于数值积分的求解公式1.梯形公式梯形(求解)公式,简称梯形法:

(1.14)注：①梯形公式的余项：

故是二阶精度.1.3基于数值积分的求解公式1.梯形公式

(1.14)②梯形公式为隐式公式.<显化>预测——校正系统

(1.15)称(1.15)为改进的欧拉公式，也可记为§1引言§1引言1.3基于数值积分的求解公式1.梯形公式

(1.14)③可以证明，改进欧拉公式也具有二阶精度.§1引言1.3基于数值积分的求解公式【例3】用欧拉法，梯形法以及改进欧拉法求解取h=0.1.计算到x=0.5.解：f(x,y)=x–y+1,a=x0=0,b=0.5,y0=1,n=5(Euler法)求解公式：yk=yk–1+h(xk–1–yk–1+1)=hxk–1+(1–h)yk–1+h

=0.1xk–1+0.9yk–1+0.1

§1引言1.3基于数值积分的求解公式【例3】用欧拉法，梯形法以及改进欧拉法求解解：f(x,y)=x–y+1,a=x0=0,b=0.5,y0=1,n=5(梯形法)求解公式：yk=yk–1+h[(xk–1–yk–1+1)+(xk–yk+1)]/2解出yk，得方程§1引言1.3基于数值积分的求解公式【例3】用欧拉法，梯形法以及改进欧拉法求解解：f(x,y)=x–y+1,a=x0=0,b=0.5,y0=1,n=5(改进Euler法)求解公式：yk=yk–1+h[(xk–1–yk–1+1)+xk–(yk+h(xk–yk+1))+1]/2得=0.905yk–1+0.045xk–1+0.05xk+0.095方程§1引言1.3基于数值积分的求解公式2.辛卜生公式

记

(1.17)§1引言1.3基于数值积分的求解公式2.辛卜生公式记

(1.17)其余项§1引言1.3基于数值积分的求解公式2.辛卜生公式记

(1.17)将[xk–1,xk]

对分：调整下标为[xi–2,xi]：xi–2=xk–1,xi–1=xk–1+h1,xi=xk–1+2h1=xk则(1.17)化为

(1.19)称(1.19)为辛卜生求解公式，其中fk–2=f(xk–2,y(xk–2))，fk–1=f(xk–1,y(xk–1))，fk=f(xk,y(xk))§1引言1.3基于数值积分的求解公式2.辛卜生公式记

(1.17)

(1.19)称(1.19)为辛卜生求解公式，其中fi–2=f(xi–2,y(xi–2))，fi–1=f(xi–1,y(xi–1))，fi=f(xi,y(xi))注：①(1.19)的误差：§1引言1.3基于数值积分的求解公式2.辛卜生公式记

(1.17)

(1.19)称(1.19)为辛卜生求解公式，其中fi–2=f(xi–2,y(xi–2))，fi–1=f(xi–1,y(xi–1))，fi=f(xi,y(xi))注：②隐式(需显化)多步——将在§3中讨论.§2Runge-Kutta法2.0原理

其中K

=f(

,y(

))=y'(

)称为y在[xi–1,xi]上的平均斜率.欧拉法：改进欧拉法：(2.1)§2Runge-Kutta法2.0原理

其中K

=f(

,y(

))=y'(

)称为y在[xi–1,xi]上的平均斜率.对(1.17)显化：辛卜生：

(2.4)§2Runge-Kutta法2.0原理其中K

=f(

,y(

))=y'(

)称为y在[xi–1,xi]上的平均斜率.设想：在中多计算(预测)几个点上的值然后可加权取平均值作为的近似值可能构成更高阶的公式.一阶二阶三阶§2Runge-Kutta法2.1Runge-Kutta公式

(*)其中0

j

1，yi–1+

jh是y(xi–1+

jh)的预测值.

称(*)为R-K公式<R-K法>注：①(2.1)(2.4)分别称为二阶，三阶R-K公式.②

j，

j，

j为待定系数.使(*)的阶数尽量高.§2Runge-Kutta法2.1Runge-Kutta公式参数的确定，以m=2为例.

欲求

1，

2，

2

.原则:使ei(h)=y(xi)–yi的阶数尽可能高§2Runge-Kutta法2.1Runge-Kutta公式展开展开

原则:使ei(h)=y(xi)–yi的阶数尽可能高§2Runge-Kutta法2.1Runge-Kutta公式

原则:使ei(h)=y(xi)–yi的阶数尽可能高§2Runge-Kutta法2.1Runge-Kutta公式

欲求截断误差ei(h)=y(xi)–yi关于h的阶数尽可能高，应使无穷多解，从而有许多2阶R-K公式§2Runge-Kutta法2.1Runge-Kutta公式应使注：①取

1=

2=1/2，

2=1，即为改进欧拉公式.§2Runge-Kutta法2.1Runge-Kutta公式应使注：②取

1=0，

2=1，

2=1/2，即为中点公式§2Runge-Kutta法2.1Runge-Kutta公式应使注：③二阶R-K公式的截断误差为故为二阶方法.相仿可得更高阶的R-K公式.§2Runge-Kutta法2.2经典R-K公式在4解R-K公式中最重要的是经典R-K公式<标准R-K公式>.(2.6)注：①(2.6)为4阶方法.§2Runge-Kutta法2.2经典R-K公式在4解R-K公式中最重要的是经典R-K公式<标准R-K公式>.(2.6)注：②R-K法对4阶以上不一定能提高整数阶.§2Runge-Kutta法2.2经典R-K公式【例4】使用三阶，四阶R-K法求解初值问题：

的部分计算值y1，y2，y3，其中h=0.1.解①使用三阶R-K法§2Runge-Kutta法【例4】使用三阶，四阶R-K法求解初值问题：

的部分计算值y1，y2，y3，其中h=0.1.解①使用三阶R-K法§2Runge-Kutta法【例4】使用三阶，四阶R-K法求解初值问题：

的部分计算值y1，y2，y3，其中h=0.1.解②使用四阶R-K法§2Runge-Kutta法【例4】使用三阶，四阶R-K法求解初值问题：

的部分计算值y1，y2，y3，其中h=0.1.解②使用四阶R-K法§2Runge-Kutta法注使用R-K法要求具备较好的光滑性，否则效果不如低阶的.<P189>作业P2098②③9，10.§3线性多步法单步法的优点：简单，计算yk+1只用yk.缺点:没有充分利用前面的信息且计算y(xk+

h)较困难回顾Simpson：(1.19)考虑：(3.1)两种插值求积：①将[xk–1,xk]增加内部节点，改为[xk–2,xk]导出的公式称为闭型求解公式.线性多步§3线性多步法考虑：(3.1)两种插值求积：①将[xk–1,xk]增加内部节点，改为[xk–2,xk]导出的公式称为闭型求解公式.②在[xk–1,xk]外增加插值节点，导出的公式称为开型求解公式.开型有显和隐，闭型也有显和隐.§3线性多步法3.1开型求解公式1.亚当斯显式求解公式取节点xk–3,xk–2,xk–1，在[xk–3,xk]上作F(x)=f(x,y(x))的插值多项式.§3线性多步法3.1开型求解公式1.亚当斯显式求解公式取节点xk–3,xk–2,xk–1，在[xk–3,xk]上记xk–i=xk–ih,x=xk+th，则

§3线性多步法3.1开型求解公式1.亚当斯显式求解公式取节点xk–3,xk–2,xk–1，记xk–i=xk–ih,x=xk+th，则

代入(3.1)得§3线性多步法3.1开型求解公式1.亚当斯显式求解公式取节点xk–3,xk–2,xk–1，记xk–i=xk–ih,x=xk+th，则

令

(3.4)称(3.4)为亚当斯显式求解公式(线性多步).§3线性多步法3.1开型求解公式1.亚当斯显式求解公式取节点xk–3,xk–2,xk–1，记xk–i=xk–ih,x=xk+th，则余项：∵∴从而(3.4)具有3阶精度.称为3阶亚当斯求解公式.§3线性多步法3.1开型求解公式1.亚当斯显式求解公式类似地取xk–4,xk–3,xk–2,xk–1

在[xk–4,xk]上作F(x)=f(x,y(x))的插值多项式，可导出4阶亚当斯显式求解公式：

(3.6)

(3.7)——4阶精度§3线性多步法3.1开型求解公式2.亚当斯隐式求解公式取xk–3,xk–2,xk–1,xk，在[xk–3,xk]上作F(x)=f(x,y(x))的插值多项式用上述方法可导出：

(3.8)

(3.9)称为亚当斯隐式求解公式.§3线性多步法3.1开型求解公式2.亚当斯隐式求解公式

(3.8)

(3.9)称为亚当斯隐式求解公式.注：利用4阶公式(3.6)显化之：

(3.10)称(3.10)为亚当斯预测——校正系统.§3线性多步法3.2闭型求解系统将[xk–1,xk]扩充为[xk–4,xk]，取xk–4，xk–3，xk–2，xk–1为节点，作F(x)=f(x,y(x)