计算方法4常微分方程数值解法课件

常微分方程的数值解法NumericalSolutionstoOrdinaryDifferentialEquations

常微分方程的数值解法对象一阶常微分方程初值问题:一阶常微分方程组初值问题:高阶常微分方程初值问题:(4.1)对象一阶常微分方程初值问题:一阶常微分方程组初值问题:高阶常一阶常微分方程初值问题:实际工程技术、生产、科研上会出现大量的微分方程问题很难得到其解析解，有的甚至无法用解析表达式来表示，因此只能依赖于数值方法去获得微分方程的数值解。一阶常微分方程初值问题:实际工程技术、生产、科研上会出现大量ab

x0x1x2...xn-1xn用数值方法，求得y(x)在每个节点xk的值y(xk)的近似值，用yk表示，即yk

≈y(xk)，这样y0,y1,...,yn称为微分方程的数值解求y(x)————>求y0,y1,...,yn

？微分方程的数值解法：不求y=y(x)的精确表达式,而求离散点x0,x1,…xn处的函数值设(4.1)

的解y(x)的存在区间是[a,b]，初始点x0=a，取[a,b]内的一系列节点x0,x1,...,xn。a=x0<x1<…<xn=b，一般采用等距步长abx0x1x2...思路计算过程：方法：采用步进式和递推法将[a,b]n等分,a=x0<x1<…<xn=b,步长h=(b-a)/n,xk=a+kh

怎样建立递推公式？Taylor公式数值积分法思路计算过程：方法：采用步进式和递推法将[a,b]n等分,4.1Euler

公式思想:用向前差商近似代替微商.(4.2)欧拉公式(EulerScheme）4.1Euler公式思想:用向前差商近似代替微商.(几何意义 y(x)过点P0(x0,y0)且在任意点(x,y)的切线斜率为f(x,y)y(x)在点P0(x0,y0)的切线方程为：

y=y0+f(x0,y0)(x-x0)在切线上取点P1(x1,y1)

y1=y0+f(x0,y0)(x1-x0)y1正是Euler公式所求4.类似2，过P1以f(x1,y1)为斜率作y(x)的切线，在其上取点

P2(x2,y2)，依此类推…5.折线P0P1P2…Pn…作为曲线y(x)的近似

——欧拉折线法xp0p1p2p3p4x0x1x2x3x4yy(x)几何意义 y(x)过点P0(x0,y0)且在任意点(x,y)思想:用向后差商近似代替微商.欧拉法（续）用隐式欧拉法，每一步都需解方程（或先解出yn+1的显式表达式），但其稳定性好。隐式欧拉公式(4.3)思想:用向后差商近似代替微商.欧拉法（续）用隐式欧拉法，计算方法4_常微分方程数值解法课件整体误差ek=y(xk)-yk，下面对其加以分析y1=y0+hf(x0,y0)=1+0.1×(1-0/1)=1.1y2=y1+hf(x1,y1)=1.1+0.1×(1.1-2×0.1/1.1)=1.191818y3=y2+hf(x2,y2)=1.277438…其精确解为整体误差ek=y(xk)-yk，下面对其加以分析y1=y0+计算方法4_常微分方程数值解法课件欧拉法（续）思想:用中心差商近似代替微商.注：计算时，先用欧拉法求出y1，以后再用二步欧拉法计算。二步欧拉法(4.4)欧拉法（续）思想:用中心差商近似代替微商.注：计算时，先欧拉法（续）公式单步否显式否单步显式单步隐式二步显式截断误差y(xn+1)-yn+1

欧拉法（续）公式单步否显式否单步显式单步隐式二步显式截断误差截断误差Def4.1设y(xn)

是(4.1)式的精确解，yn是(4.2)式欧拉法得到的近似解，称y(xn)-yn为欧拉法的整体截断误差.

Def4.3若某算法的局部截断误差为O(hp+1),称该算法有p阶精度.Def4.2假设yn=y(xn)

,即第n步计算是精确的前提下,称Rn+1=y(xn+1)-yn+1为欧拉法的局部截断误差.

截断误差Def4.1设y(xn)是(4.1)式的精确解，分析：证明其局部截断误差为O(h2),可通过Taylor展开式分析。证明：Euler

公式为令yn=y(xn)，下证：

y(xn+1)-yn+1=

O(h2)由y’(x)

=f(x,y(x))

定理4.4

欧拉法的精度是一阶。分析：证明其局部截断误差为O(h2),可通过Taylor展开二步欧拉法的局部截断误差Def4.5假设yn=y(xn),yn－1=y(xn－1),称Rn+1=y(xn+1)-yn+1为二步欧拉法的局部截断误差.

定理4.6

隐式欧拉法的精度是一阶,二步欧拉法的精度是二阶。证明：对二步欧拉法进行证明，考虑其局部截断误差，令yn=y(xn),yn－1=y(xn－1),将上两式左右两端同时相减：∴二步欧拉法的局部截断误差为O(h3)，其精度是二阶。二步欧拉法的局部截断误差Def4.5假设yn=y(xn)数值积分法对右端的定积分用数值积分公式求近似值：（1）用左矩形数值积分公式：（2）用梯形公式：——梯形公式梯形公式:将显示欧拉公式,隐式欧拉公式平均可得梯形公式是隐式、单步公式,其精度为二阶数值积分法对右端的定积分用数值积分公式求近似值：（1）用左矩证：令yn=y(xn),由Talor公式有分析：梯形公式是隐式公式，证明其局部截断误差为O(h3)

要用到二元函数的Taylor公式。f(xn+1,yn+1)=f(xn+1,y(xn+1)+(yn+1-y(xn+1))=f(xn+1,y(xn+1))+fy(xn+1,η)(yn+1-y(xn+1))，η∈(xn，xn+1)

=y’(xn+1)+fy(xn+1,η)(yn+1-y(xn+1))=y’(xn)+hy”(xn)+O(h2)+fy(xn+1,η)(yn+1-y(xn+1))=f(xn,yn)+hy”(xn)+fy(xn+1,η)(yn+1-y(xn+1))+O(h2)又y(xn+1)=y(xn+h)=y(xn)+hy’(xn)+h2y”(xn)/2+O(h3)=yn+hf(xn,yn)+h2y”(xn)/2+O(h3)=yn+hf(xn,yn)/2+h[f(xn,yn)+hy”(xn)]/2+O(h3)定理4.7

梯形公式的精度是二阶。证：令yn=y(xn),由Talor公式有分析：梯形公式是隐从而y(xn+1)=yn+1+hfy(xn+1,η)[y(xn+1)-yn+1]/2+O(h3)∴y(xn+1)-yn+1=

hfy(xn+1,η)[y(xn+1)-yn+1]/2+O(h3)∴y(xn+1)-yn+1=O(h3)/[1-hfy(xn+1,η)/2]=O(h3)∴梯形公式的截断误差为O(h3)，其精度是2阶。f(xn+1,yn+1)=f(xn,yn)+hy”(xn)+fy(xn+1,η)(yn+1-y(xn+1))+O(h2)y(xn+1)=yn+hf(xn,yn)/2+h[f(xn,yn)+hy”(xn)]/2+O(h3)=yn+hf(xn,yn)/2+h[f(xn+1,yn+1)-fy(xn+1,η)(yn+1-y(xn+1))+O(h2)]/2+O(h3)=yn+h[f(xn,yn)+f(xn+1,yn+1)]/2+hfy(xn+1,η)(y(xn+1)-yn+1)/2+O(h3)从而y(xn+1)=yn+1+hfy(xn+1,η)[y(解:

取h=0.01x0=0y0=y(0)=1则

y(0.01)≈y1

f(x,y)=y，由梯形公式,基于稳定性考虑解析解解:取h=0.01x0=0y0=y(0)=欧拉公式的比较欧拉法最简单,精度低隐式欧拉法稳定性好二步欧拉法显式,但需要两步初值,且第2个初值只能由其他方法给出,可能对后面的递推精度有影响梯形公式法精度有所提高,但为隐式,需要迭代求解,计算量大欧拉公式的比较欧拉法最简单,精度低隐式欧拉法稳定性好二步欧拉4.2

改进的Euler法

Euler公式计算量小，精度低梯形公式计算量大，精度高综合两个公式，提出预报—校正公式：

预报

校正——改进的Euler法写成嵌套公式：平均化形式：4.2改进的Euler法Euler公式计算量小[例4.4]用改进的Euler法解初值问题在区间[0,0.4]上，步长h=0.1的解，并比较与精确解(y=1/(1-x))的差异。解：Euler法的具体形式为：yn+1=yn+hyn2，改进的Euler法的具体形式为：∵x0=0,h=0.1,则

x1=0.1,x2=0.2,x3=0.3,x4=0.4

计算y1:

yp=y0+0.1y02=1+0.1·12=1.1yc=1+0.1X1.12=1.121y1=(1.1+1.121)/2≈1.1118同样可求y2,、y3、,y4[例4.4]用改进的Euler法解初值问题在区间[0,0.注：（1）令y(xn)=yn，可推导改进的Euler法的局部截断误差

为O(h3),具有二阶精度。（2）

改进的Euler法也可写成如下平均化形式注：（2）改进的Euler法也可写成如下平均化形式4.3龙格—库塔方法如何构造高阶的方法？(4.5)为了构造函数φ使得(4.5)式成为高阶方法，Taylor对于一般的显式单步法若讨论其精度（阶数），即考虑4.3龙格—库塔方法如何构造高阶的方法？(4.5)为了构造令则有考虑到上述方法求高阶微分，实际上不实用令则有考虑到上述方法求高阶微分，实际上不实用改进的Euler公式

:Euler公式：yn+1=yn+hf(xn,yn)写成精度：一阶

精度：二阶

由Lagrange中值定理，称为y(x)在[xn,xn+1]上的平均斜率取—Euler公式

取—改进Euler公式

Euler公式用一个点的值k1作为k*的近似值，而改进的Euler公式用二个点的值k1和k2的平均值作为k*近似值，改进的Euler法比Euler法精度高。改进的Euler公式:Euler公式：yn+1=yn+hf(4.6)Runge-Kutta法的思想：在［xn,xn+1］内多预报几个点的ki值并用其加权平均值作为k*近似值。而构造出具有更高精度的计算公式。多个斜率加权平均可提高精度Runge-Kutta法一般形式：(4.6)Runge-Kutta法的思想：在［xn,xn+1以k1与k2的加权平均来近似k*，设其中c1,c2,α2,b21为待定参数。适当选取参数，使(*)式的精度为二阶，即使其局部截断误差为O(h3)令y(xn)=yn，由泰勒公式：二阶龙格—库塔方法(*)以k1与k2的加权平均来近似k*，设其中c1,c2,α2,b由多元函数的泰勒公式和(*)式则有

比较(a)与(b)要使y(xn+1)-yn+1=O(h3)(a)(b)由多元函数的泰勒公式和(*)式则有比较(a)与(b)要使上述方程组有四个未知量，只有三个方程，有无穷多组解。取任意一组解便得一种二阶龙－库公式。

当c1=c2=1/2,a2=b21=1时二阶Runge-Kutta公式为yn+1=yn+k1/2+k2/2k1=hf(xn,yn)k2=hf(xn+h,yn+k1)此即改进Euler法

取c2=0，c2=1，a2=1/2，b21=1/2yn+1=yn+k2k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)此为中点法或变形的

Euler公式上述方程组有四个未知量，只有三个方程，有无穷多组解。当c1三阶龙格－库塔法是用三个值k1,k2,k3的加权平均来近似k*，即有yn+1=yn+c1k1+c2k2+c3k3k1=hf(xn,yn)k2=hf(xn+a2h,yn+b21k1)k3=hf(xn+a3h,yn+b31k1+b32k2)要使其具有三阶精度，必须使局部截断误差为O(h4)类似二阶龙格－库塔法的推导，c1,c2,c3,a2,a3,b21,b31,b32应满足c1+c2+c3=1a2=b21a3=b31+b32c2a2+c3a3=1/2c2a22+c3a32=1/3c3a32=1/6由该方程组任意解可得三阶龙格-库塔公式例：Kutta公式kn+1=yn+(k1+4k2+k3)/6k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)k3=hf(xn+h,yn-k1+2k2)

三阶龙格－库塔法是用三个值k1,k2,k3的加权平均来近似k类似可推出四阶龙格-库塔公式，常用的有例：经典Runge-Kutta法

yn+1=yn+(k1+2k2+2k3+k4)/6k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)k3=hf(xn+h/2,yn+k2/2)k4=hf(xn+h,yn+k3)局部截断误差O(h5)

还有:Gill公式及m(m>4)阶龙格－库塔法。

m>4时：计算量太大，精确度不一定提高，有时会降低。类似可推出四阶龙格-库塔公式，常用的有yn+1=yn+(k1Gill公式节省存储单元控制舍入误差Gill公式节省存储单元对于经典的四阶Runge-Kutta法给出如下算法：[算法4.2]求解：

dy/dx=f(x,y)a≤x≤by(a)=y0

Step

1:

输入a,b,y0

及NStep

2:

(b-a)/N=>h,a=>x,y0=>yStep

3:

输出

(x,y)Step4:ForI=1T0Nhf(x,y)=>k1hf(x+h/2,y+k1/2)=>k2hf(x+h/2,y+k2/2)=>k3hf(x+h,y＋k3)=>k4y+(k1+2k2+2k3+k4)/6=>yx+h=>x输出(x,y)Step5:停止对于经典的四阶Runge-Kutta法给出如下算法：[算法4[例4.3]用四阶经典Runge－Kutta方法解初值问题：

(1)求，

[例4.3]用四阶经典Runge－Kutta方法解初值问题：(2)求，

(2)求，自适应龙格－库塔法用户提出问题I：问题：①：如何判断|y(xn)-yn|<ε∵精度值y(xn)未知。

②：如何取h=?解①：如用p阶龙格－库塔法计算，局部截断误差为O(hp+1)

y’=f(x,y)y(x0)=y0

要求误差<ε＝10-8求数值解｛｝xnh/2xn+1h如

xn-----------------xn+1

令yn=y(xn)yn+1(h)

则y(xn+1)-yn+1(h)≈chp+1步长折半xnxn+h/2xn+1分两步计算y(xn+1)的近似值yn+1(h/2)。则y(xn+1)-yn+1(h/2)≈2c(h/2)p+1

自适应龙格－库塔法用户提出问题I：问题：①：如何判断|y(定理：对于问题I若用P阶龙格－库塔法计算y(xn+1)在步长折半前后的近似值分别为yn+1(h),yn+1(h/2)则有误差公式

注：10

误差的事后估计法

20

停机准则：△<ε(可保证|y(xn+1)-yn+1(h/2)|<ε)

解②：⑴h取大，局部截断误差chp+1大，不精确⑵h取小，运算量大（步多），舍入误差积累大解决策略：变步长龙格－库塔法

If（△>ε）将步长折半反复计算，直至△<ε为止，取h为最后一次的步长，yn+1为最后一次计算的结果。

Elseif(△<ε)将步长加倍反复计算，直至△>ε为止，最后一次运算的前一次计算结果即为所需。定理：对于问题I若用P阶龙格－库塔法计算y(xn+1)在步长4.5收敛与稳定性对于常微分方程初值问题(4.1)求y=y(x)--------求y(xn)------------

求yn

xn=x0+nh离散化某种公式例：Euler法，改进的Euler法，龙格－库塔法都是单步法。数值解法：单步法：计算yn+1时只用到前一步的结果yn。显式单步法：

yn+1=yn+hφ(xn,yn,h)

φ(x,y,h)为增量函数，它依赖于f,仅是xn,yn,h的函数Def:若某数值方法对于任意固定的xn=x0+nh,当

h->0(n->∞)时，

yn->y(xn)，则称该法是收敛的。即(xn=x0+nh为固定值)4.5收敛与稳定性对于常微分方程初值问题(4.1)求y=改进Euler法的收敛性

判定改进Euler法的收敛性：

yn+1=yn+h[f(xn,yn)+f(xn+h,yn+hf(xn,yn))]/2

则φ(x,y,h)=[f(x,y)+f(x+h,y+hf(x,y))]/2若：①yo=y(x0)

②f(x,y)关于y满足李--条件:则：

限定h≤h0，则即Φ(x,y,h)

满足李普希兹条件,由Th改进的Euler法收敛同样可验证，若f(x,y)关于y满足李普希兹条件,且y0=y(x0)时，Euler法，标准四阶龙格－库法也收敛。改进Euler法的收敛性判定改进Euler法的收敛性：若：4.5.2

稳定性在讨论收敛性时，一般认可：数值方法本身计算过程是准确的。实际上，并非如此：①

初始值y0有误差δ0＝y0-y(x0).②

计算的每一步有舍入误差。这些初始误差在计算过程的传播中，是逐步衰减的，还是恶性增长，这就是稳定性问题。4.5.2稳定性在讨论收敛性时，一般认可：数值方法本身计算Def4.1若一种数值方法在节点xn处的数值解yn的扰动δn≠0，而在以后的各节点值ym(m>n)上有扰动δm.当|δm|≤|δn|，(m=n+1,n+2,……)，则称该数值算法是稳定的。分析某算法的数值稳定性，通常考察模型方程

y’=λy，(λ<0)Def4.1：设在节点xn处用数值法得到的理想数值解为yn，而实际计算得到的近似值为，称为第n步数值解的扰动Def4.1若一种数值方法在节点xn处的数值解yn的扰动δnEuler法的稳定性Euler法：yn+1=yn+hf(xn,yn)考察模型方程

y’=λy，(λ<0)

即yn+1=(1+hλ)yn假设在节点值yn上有扰动δn，在yn+1上有扰动δn+1，且δn+1仅由δn引起(计算过程不再引进新的误差)Euler法稳定Euler法的稳定的条件为0<h≤-2/λEuler法的稳定性Euler法：yn+1=yn+hf(xn隐式Euler法稳定性隐式Euler法：yn+1=yn+hf(xn+1,yn+1)对于模型方程y’=λy(λ<0)则yn+1=yn+hλyn+1=>

yn+1=yn/(1-hλ)yn+1的扰动∵λ<0，∴h>0

上式均成立，隐式Euler法无条件稳定210ReImg隐式Euler法稳定性隐式Euler法：yn+1=yn+hf梯形公式稳定性梯形公式

yn+1=yn+h[f(xn,yn)+f(xn+1,yn+1)]/2，设模型方程为y’=λy(λ<0)则yn+1=yn+h[λyn+λyn+1]/2

λ<0上式对任意h>0时恒成立

梯形公式恒稳定。

yn+1的扰动梯形公式稳定性梯形公式yn+1=yn+h[f(xn,yn3阶Runge－Kutta显式

1~4阶方法的绝对稳定区域为k=1k=2k=3k=4-1-2-3---123ReImg3阶Runge－Kutta显式1~4阶方法的绝对Taylor公式公式单步否显式否精度单步显式一阶单步隐式一阶二步显式二阶数值积分法单步隐式二阶单步显式二阶Runge－Kutta方法局部截断误差整体截断误差收敛与稳定性Taylor公式公式单步否显式否精度单步显式一阶单步隐式一4.6

多步法某一步解和前若干步解的值有关m步线性多步法的一般形式其中a0,…,am-1,b0,…,bm为和h无关的常数初值为y0,y1,…,ym-1若bm≠0，方法是隐式的(implicit)否则是显式的(explicit)(4.7)4.6多步法某一步解和前若干步解的值有关m步线性多步法的一Adams－Bashforth四阶方法Adams－Moulton四阶方法Adams－Bashforth四阶方法Adams－Moult多步法的截断误差多步法的截断误差4.7

常微分方程组和高阶微分方程的数值解法形如欧拉方法的计算公式隐式欧拉方法的计算公式4.7常微分方程组和高阶微分方程的数值解法形如欧拉方法的计梯形公式Adams公式R-K方法梯形公式Adams公式R-K方法对于高阶微分方程，总可以化成方程组的形式，例如可以化成对于高阶微分方程，总可以化成方程组的形式，例如可以化成例令y’(x)=z，则有用4阶RK方法，有例令y’(x)=z，则有用4阶RK方法，有计算方法4_常微分方程数值解法课件4.8

常微分方程边值问题的数值解法形如第一类边界条件(4.8)第二类边界条件第三类边界条件4.8常微分方程边值问题的数值解法形如第一类边界条件(4.定理假定问题（4.8）式中的函数f在集合上连续，且偏导数fy和fy’也在D上连续。若有则边值问题有唯一解。定理假定问题（4.8）式中的函数f在集合上连续，且偏导数f例边值问题有因此边值问题有唯一解。定理若线性边值问题满足则边值问题有唯一解。例边值问题有因此边值问题有唯一解。定理若线性边值问题满足打靶法（ShootingMethod）其中，y1(x)是初值问题的解的解，y2(x)是初值问题打靶法（ShootingMethod）其中，y1(x)是初有限差分法(Finite-DifferenceMethods)以差商代替导数，将微分方程离散成为差分方程用Taylor展开方法有限差分法(Finite-DifferenceMethod两式相加得到由中值定理得线性边值问题可写成两式相加得到由中值定理得线性边值问题可写成然后可得到有限差分公式上式可写成写成矩阵形式然后可得到有限差分公式上式可写成写成矩阵形式其中其中常微分方程的数值解法NumericalSolutionstoOrdinaryDifferentialEquations

常微分方程的数值解法对象一阶常微分方程初值问题:一阶常微分方程组初值问题:高阶常微分方程初值问题:(4.1)对象一阶常微分方程初值问题:一阶常微分方程组初值问题:高阶常一阶常微分方程初值问题:实际工程技术、生产、科研上会出现大量的微分方程问题很难得到其解析解，有的甚至无法用解析表达式来表示，因此只能依赖于数值方法去获得微分方程的数值解。一阶常微分方程初值问题:实际工程技术、生产、科研上会出现大量ab

x0x1x2...xn-1xn用数值方法，求得y(x)在每个节点xk的值y(xk)的近似值，用yk表示，即yk

≈y(xk)，这样y0,y1,...,yn称为微分方程的数值解求y(x)————>求y0,y1,...,yn

？微分方程的数值解法：不求y=y(x)的精确表达式,而求离散点x0,x1,…xn处的函数值设(4.1)

的解y(x)的存在区间是[a,b]，初始点x0=a，取[a,b]内的一系列节点x0,x1,...,xn。a=x0<x1<…<xn=b，一般采用等距步长abx0x1x2...思路计算过程：方法：采用步进式和递推法将[a,b]n等分,a=x0<x1<…<xn=b,步长h=(b-a)/n,xk=a+kh

怎样建立递推公式？Taylor公式数值积分法思路计算过程：方法：采用步进式和递推法将[a,b]n等分,4.1Euler

公式思想:用向前差商近似代替微商.(4.2)欧拉公式(EulerScheme）4.1Euler公式思想:用向前差商近似代替微商.(几何意义 y(x)过点P0(x0,y0)且在任意点(x,y)的切线斜率为f(x,y)y(x)在点P0(x0,y0)的切线方程为：

y=y0+f(x0,y0)(x-x0)在切线上取点P1(x1,y1)

y1=y0+f(x0,y0)(x1-x0)y1正是Euler公式所求4.类似2，过P1以f(x1,y1)为斜率作y(x)的切线，在其上取点

P2(x2,y2)，依此类推…5.折线P0P1P2…Pn…作为曲线y(x)的近似

——欧拉折线法xp0p1p2p3p4x0x1x2x3x4yy(x)几何意义 y(x)过点P0(x0,y0)且在任意点(x,y)思想:用向后差商近似代替微商.欧拉法（续）用隐式欧拉法，每一步都需解方程（或先解出yn+1的显式表达式），但其稳定性好。隐式欧拉公式(4.3)思想:用向后差商近似代替微商.欧拉法（续）用隐式欧拉法，计算方法4_常微分方程数值解法课件整体误差ek=y(xk)-yk，下面对其加以分析y1=y0+hf(x0,y0)=1+0.1×(1-0/1)=1.1y2=y1+hf(x1,y1)=1.1+0.1×(1.1-2×0.1/1.1)=1.191818y3=y2+hf(x2,y2)=1.277438…其精确解为整体误差ek=y(xk)-yk，下面对其加以分析y1=y0+计算方法4_常微分方程数值解法课件欧拉法（续）思想:用中心差商近似代替微商.注：计算时，先用欧拉法求出y1，以后再用二步欧拉法计算。二步欧拉法(4.4)欧拉法（续）思想:用中心差商近似代替微商.注：计算时，先欧拉法（续）公式单步否显式否单步显式单步隐式二步显式截断误差y(xn+1)-yn+1

欧拉法（续）公式单步否显式否单步显式单步隐式二步显式截断误差截断误差Def4.1设y(xn)

是(4.1)式的精确解，yn是(4.2)式欧拉法得到的近似解，称y(xn)-yn为欧拉法的整体截断误差.

Def4.3若某算法的局部截断误差为O(hp+1),称该算法有p阶精度.Def4.2假设yn=y(xn)

,即第n步计算是精确的前提下,称Rn+1=y(xn+1)-yn+1为欧拉法的局部截断误差.

截断误差Def4.1设y(xn)是(4.1)式的精确解，分析：证明其局部截断误差为O(h2),可通过Taylor展开式分析。证明：Euler

公式为令yn=y(xn)，下证：

y(xn+1)-yn+1=

O(h2)由y’(x)

=f(x,y(x))

定理4.4

欧拉法的精度是一阶。分析：证明其局部截断误差为O(h2),可通过Taylor展开二步欧拉法的局部截断误差Def4.5假设yn=y(xn),yn－1=y(xn－1),称Rn+1=y(xn+1)-yn+1为二步欧拉法的局部截断误差.

定理4.6

隐式欧拉法的精度是一阶,二步欧拉法的精度是二阶。证明：对二步欧拉法进行证明，考虑其局部截断误差，令yn=y(xn),yn－1=y(xn－1),将上两式左右两端同时相减：∴二步欧拉法的局部截断误差为O(h3)，其精度是二阶。二步欧拉法的局部截断误差Def4.5假设yn=y(xn)数值积分法对右端的定积分用数值积分公式求近似值：（1）用左矩形数值积分公式：（2）用梯形公式：——梯形公式梯形公式:将显示欧拉公式,隐式欧拉公式平均可得梯形公式是隐式、单步公式,其精度为二阶数值积分法对右端的定积分用数值积分公式求近似值：（1）用左矩证：令yn=y(xn),由Talor公式有分析：梯形公式是隐式公式，证明其局部截断误差为O(h3)

要用到二元函数的Taylor公式。f(xn+1,yn+1)=f(xn+1,y(xn+1)+(yn+1-y(xn+1))=f(xn+1,y(xn+1))+fy(xn+1,η)(yn+1-y(xn+1))，η∈(xn，xn+1)

=y’(xn+1)+fy(xn+1,η)(yn+1-y(xn+1))=y’(xn)+hy”(xn)+O(h2)+fy(xn+1,η)(yn+1-y(xn+1))=f(xn,yn)+hy”(xn)+fy(xn+1,η)(yn+1-y(xn+1))+O(h2)又y(xn+1)=y(xn+h)=y(xn)+hy’(xn)+h2y”(xn)/2+O(h3)=yn+hf(xn,yn)+h2y”(xn)/2+O(h3)=yn+hf(xn,yn)/2+h[f(xn,yn)+hy”(xn)]/2+O(h3)定理4.7

梯形公式的精度是二阶。证：令yn=y(xn),由Talor公式有分析：梯形公式是隐从而y(xn+1)=yn+1+hfy(xn+1,η)[y(xn+1)-yn+1]/2+O(h3)∴y(xn+1)-yn+1=

hfy(xn+1,η)[y(xn+1)-yn+1]/2+O(h3)∴y(xn+1)-yn+1=O(h3)/[1-hfy(xn+1,η)/2]=O(h3)∴梯形公式的截断误差为O(h3)，其精度是2阶。f(xn+1,yn+1)=f(xn,yn)+hy”(xn)+fy(xn+1,η)(yn+1-y(xn+1))+O(h2)y(xn+1)=yn+hf(xn,yn)/2+h[f(xn,yn)+hy”(xn)]/2+O(h3)=yn+hf(xn,yn)/2+h[f(xn+1,yn+1)-fy(xn+1,η)(yn+1-y(xn+1))+O(h2)]/2+O(h3)=yn+h[f(xn,yn)+f(xn+1,yn+1)]/2+hfy(xn+1,η)(y(xn+1)-yn+1)/2+O(h3)从而y(xn+1)=yn+1+hfy(xn+1,η)[y(解:

取h=0.01x0=0y0=y(0)=1则

y(0.01)≈y1

f(x,y)=y，由梯形公式,基于稳定性考虑解析解解:取h=0.01x0=0y0=y(0)=欧拉公式的比较欧拉法最简单,精度低隐式欧拉法稳定性好二步欧拉法显式,但需要两步初值,且第2个初值只能由其他方法给出,可能对后面的递推精度有影响梯形公式法精度有所提高,但为隐式,需要迭代求解,计算量大欧拉公式的比较欧拉法最简单,精度低隐式欧拉法稳定性好二步欧拉4.2

改进的Euler法

Euler公式计算量小，精度低梯形公式计算量大，精度高综合两个公式，提出预报—校正公式：

预报

校正——改进的Euler法写成嵌套公式：平均化形式：4.2改进的Euler法Euler公式计算量小[例4.4]用改进的Euler法解初值问题在区间[0,0.4]上，步长h=0.1的解，并比较与精确解(y=1/(1-x))的差异。解：Euler法的具体形式为：yn+1=yn+hyn2，改进的Euler法的具体形式为：∵x0=0,h=0.1,则

x1=0.1,x2=0.2,x3=0.3,x4=0.4

计算y1:

yp=y0+0.1y02=1+0.1·12=1.1yc=1+0.1X1.12=1.121y1=(1.1+1.121)/2≈1.1118同样可求y2,、y3、,y4[例4.4]用改进的Euler法解初值问题在区间[0,0.注：（1）令y(xn)=yn，可推导改进的Euler法的局部截断误差

为O(h3),具有二阶精度。（2）

改进的Euler法也可写成如下平均化形式注：（2）改进的Euler法也可写成如下平均化形式4.3龙格—库塔方法如何构造高阶的方法？(4.5)为了构造函数φ使得(4.5)式成为高阶方法，Taylor对于一般的显式单步法若讨论其精度（阶数），即考虑4.3龙格—库塔方法如何构造高阶的方法？(4.5)为了构造令则有考虑到上述方法求高阶微分，实际上不实用令则有考虑到上述方法求高阶微分，实际上不实用改进的Euler公式

:Euler公式：yn+1=yn+hf(xn,yn)写成精度：一阶

精度：二阶

由Lagrange中值定理，称为y(x)在[xn,xn+1]上的平均斜率取—Euler公式

取—改进Euler公式

Euler公式用一个点的值k1作为k*的近似值，而改进的Euler公式用二个点的值k1和k2的平均值作为k*近似值，改进的Euler法比Euler法精度高。改进的Euler公式:Euler公式：yn+1=yn+hf(4.6)Runge-Kutta法的思想：在［xn,xn+1］内多预报几个点的ki值并用其加权平均值作为k*近似值。而构造出具有更高精度的计算公式。多个斜率加权平均可提高精度Runge-Kutta法一般形式：(4.6)Runge-Kutta法的思想：在［xn,xn+1以k1与k2的加权平均来近似k*，设其中c1,c2,α2,b21为待定参数。适当选取参数，使(*)式的精度为二阶，即使其局部截断误差为O(h3)令y(xn)=yn，由泰勒公式：二阶龙格—库塔方法(*)以k1与k2的加权平均来近似k*，设其中c1,c2,α2,b由多元函数的泰勒公式和(*)式则有

比较(a)与(b)要使y(xn+1)-yn+1=O(h3)(a)(b)由多元函数的泰勒公式和(*)式则有比较(a)与(b)要使上述方程组有四个未知量，只有三个方程，有无穷多组解。取任意一组解便得一种二阶龙－库公式。

当c1=c2=1/2,a2=b21=1时二阶Runge-Kutta公式为yn+1=yn+k1/2+k2/2k1=hf(xn,yn)k2=hf(xn+h,yn+k1)此即改进Euler法

取c2=0，c2=1，a2=1/2，b21=1/2yn+1=yn+k2k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)此为中点法或变形的

Euler公式上述方程组有四个未知量，只有三个方程，有无穷多组解。当c1三阶龙格－库塔法是用三个值k1,k2,k3的加权平均来近似k*，即有yn+1=yn+c1k1+c2k2+c3k3k1=hf(xn,yn)k2=hf(xn+a2h,yn+b21k1)k3=hf(xn+a3h,yn+b31k1+b32k2)要使其具有三阶精度，必须使局部截断误差为O(h4)类似二阶龙格－库塔法的推导，c1,c2,c3,a2,a3,b21,b31,b32应满足c1+c2+c3=1a2=b21a3=b31+b32c2a2+c3a3=1/2c2a22+c3a32=1/3c3a32=1/6由该方程组任意解可得三阶龙格-库塔公式例：Kutta公式kn+1=yn+(k1+4k2+k3)/6k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)k3=hf(xn+h,yn-k1+2k2)

三阶龙格－库塔法是用三个值k1,k2,k3的加权平均来近似k类似可推出四阶龙格-库塔公式，常用的有例：经典Runge-Kutta法

yn+1=yn+(k1+2k2+2k3+k4)/6k1=hf(xn,yn)k2=hf(xn+h/2,yn+k1/2)k3=hf(xn+h/2,yn+k2/2)k4=hf(xn+h,yn+k3)局部截断误差O(h5)

还有:Gill公式及m(m>4)阶龙格－库塔法。

m>4时：计算量太大，精确度不一定提高，有时会降低。类似可推出四阶龙格-库塔公式，常用的有yn+1=yn+(k1Gill公式节省存储单元控制舍入误差Gill公式节省存储单元对于经典的四阶Runge-Kutta法给出如下算法：[算法4.2]求解：

dy/dx=f(x,y)a≤x≤by(a)=y0

Step

1:

输入a,b,y0

及NStep

2:

(b-a)/N=>h,a=>x,y0=>yStep

3:

输出

(x,y)Step4:ForI=1T0Nhf(x,y)=>k1hf(x+h/2,y+k1/2)=>k2hf(x+h/2,y+k2/2)=>k3hf(x+h,y＋k3)=>k4y+(k1+2k2+2k3+k4)/6=>yx+h=>x输出(x,y)Step5:停止对于经典的四阶Runge-Kutta法给出如下算法：[算法4[例4.3]用四阶经典Runge－Kutta方法解初值问题：

(1)求，

[例4.3]用四阶经典Runge－Kutta方法解初值问题：(2)求，

(2)求，自适应龙格－库塔法用户提出问题I：问题：①：如何判断|y(xn)-yn|<ε∵精度值y(xn)未知。

②：如何取h=?解①：如用p阶龙格－库塔法计算，局部截断误差为O(hp+1)

y’=f(x,y)y(x0)=y0

要求误差<ε＝10-8求数值解｛｝xnh/2xn+1h如

xn-----------------xn+1

令yn=y(xn)yn+1(h)

则y(xn+1)-yn+1(h)≈chp+1步长折半xnxn+h/2xn+1分两步计算y(xn+1)的近似值yn+1(h/2)。则y(xn+1)-yn+1(h/2)≈2c(h/2)p+1

自适应龙格－库塔法用户提出问题I：问题：①：如何判断|y(定理：对于问题I若用P阶龙格－库塔法计算y(xn+1)在步长折半前后的近似值分别为yn+1(h),yn+1(h/2)则有误差公式

注：10

误差的事后估计法

20

停机准则：△<ε(可保证|y(xn+1)-yn+1(h/2)|<ε)

解②：⑴h取大，局部截断误差chp+1大，不精确⑵h取小，运算量大（步多），舍入误差积累大解决策略：变步长龙格－库塔法

If（△>ε）将步长折半反复计算，直至△<ε为止，取h为最后一次的步长，yn+1为最后一次计算的结果。

Elseif(△<ε)将步长加倍反复计算，直至△>ε为止，最后一次运算的前一次计算结果即为所需。定理：对于问题I若用P阶龙格－库塔法计算y(xn+1)在步长4.5收敛与稳定性对于常微分方程初值问题(4.1)求y=y(x)--------求y(xn)------------

求yn

xn=x0+nh离散化某种公式例：Euler法，改进的Euler法，龙格－库塔法都是单步法。数值解法：单步法：计算yn+1时只用到前一步的结果yn。显式单步法：

yn+1=yn+hφ(xn,yn,h)

φ(x,y,h)为增量函数，它依赖于f,仅是xn,yn,h的函数Def:若某数值方法对于任意固定的xn=x0+nh,当

h->0(n->∞)时，

yn->y(xn)，则称该法是收敛的。即(xn=x0+nh为固定值)4.5收敛与稳定性对于常微分方程初值问题(4.1)求y=改进Euler法的收敛性

判定改进Euler法的收敛性：

yn+1=yn+h[f(xn,yn)+f(xn+h,yn+hf(xn,yn))]/2

则φ(x,y,h)=[f(x,y)+f(x+h,y+hf(x,y))]/2若：①yo=y(x0)

②f(x,y)关于y满足李--条件:则：

限定h≤h0，则即Φ(x,y,h)

满足李普希兹条件,由Th改进的Euler法收敛同样可验证，若f(x,y)关于y满足李普希兹条件,且y0=y(x0)时，Euler法，标准四阶龙格－库法也收敛。改进Euler法的收敛性判定改进Euler法的收敛性：若：4.5.2

稳定性在讨论收敛性时，一般认可：数值方法本身计算过程是准确的。实际上，并非如此：①

初始值y0有误差δ0＝