第七章常微分方程数值解法.ppt

第七章 常微分方程数值解法 主 讲：孙 剑 聊城大学计算机学院信息管理系,计算方法 吴筑筑编,本章主要内容：,7.1 欧拉法和改进的欧拉法 7.2 龙格-库塔法 7.3 线性多步法,引言：,可求出方程y=1+ex的通解为 y=x+ex+c,将初值条件 x=0, y=2 代入得 2=1+c, 故 c=1,所以初值问题的解为 y=x+ex+1,求解初值问题,引言：,本章解决的问题：一阶常微分方程的初值问题,引言：,若方程 y=f(x,y)的右端函数f(x,y)在闭矩形域 r：x0ax x0a， y0by y0b上满足： （1） f(x,y)在r上连续, （2）在r上关于y满足lipschitz(李普希兹)条件 即存在常数l，对r上任意点均有以下不等式成立： |f(x,y1 )f(x,y2 )|l|y1y2|, xa,b, y1,y2r 则上述初值问题存在唯一的连续可微的解函数 y = y(x)。,引言：,又如初值问题,可求出它的解为,但要进一步计算指定点的函数值，还需要用数值积分方法。有些微分方程的解是隐函数，例如,要求函数值还需要解超越方程。应用中所处理的微分方程往往很复杂且大都得不出一般解，所以一般用数值解法。,引言：,数值解法： 给定节点a=x0x1xn=b, 将初值问题离散化为差分方程,求出解函数 y(x) 在这些点的近似值y1 ,y2 ,yn 。 所求得的近似值称为数值解。 本章中总假定步长h为定值，节点xi=x0+ih i=1，2，3,7.1.1 欧拉法及其截断误差,7.1.2 改进的欧拉法及预测-校正公式,7.1 欧拉法和改进的欧拉法,7.1.1 欧拉法及其截断误差,初值问题：,1、欧拉公式的构造思想：用差商代替导数,设,等距，步长为,令x=xi , x+h=xi+1 , y(xi )yi ，y(xi+1 ) yi+1 ,初值问题离散化为,(欧拉公式),7.1.1 欧拉法及其截断误差,例 取步长 h=0.1，用欧拉法求解初值问题,解:,y1=y0+h f(x0 ,y0 )=1+0.1(0 + 1 )=1.1,y2=y1+h f(x1 ,y1 )=1.1 + 0.1(0.1 + 1.1 )=1.22 y3=y2+hf(x2 ,y2 )=1.22+0.1 (0.2+1.22)=1.362,y10=y9+h f(x9 ,y9 )=y9+0.1(x9 + y9 )=3.18748,7.1.1 欧拉法及其截断误差,2、欧拉公式几何意义: 用折线代替曲线计算解函数的近似值。,7.1.1 欧拉法及其截断误差,3、数值公式的误差来源。,(1)局部截断误差（简称截断误差）：假设 yi=y(xi )是准确的 ，计算yi+1所产生的误差 y(xi+1 ) - yi+1,若局部截断误差可以表示为o(hk+1), k为正整数，则称公式是k阶公式。,（2）由于实际上yi不是准确值，因此它的误差会 传播下去。,（3）实际计算时，每一步都可能产生舍入误差。,7.1.1 欧拉法及其截断误差,4、欧拉公式的截断误差是o(h2)，公式是1 阶的。,因为,（泰勒公式）,两式相减，由设 yi=y(xi ) ，有,7.1.2 改进的欧拉法及预测-校正公式,对微分方程y=f(x,y) 两边求xi 到xi+1 的定积分，有,利用梯形公式计算积分，有,1、改进的欧拉公式的构造,7.1.2 改进的欧拉法及预测-校正公式,将y(xi ) 、y(xi+1 )分别用yi 、yi+1 代替，构造相应的数值公式：,（改进的欧拉公式）,7.1.2 改进的欧拉法及预测-校正公式,2、改进的欧拉公式的截断误差为o(h3),因而改进的欧拉法是二阶的。,7.1.2 改进的欧拉法及预测-校正公式,3、 改进的欧拉法的具体使用格式。,改进的欧拉法是隐式公式 ,计算时常用迭代法。一般每一步先由欧拉公式计算出yi+1 的初始值yi+1(0),再迭代计算yi+1。,当满足,时，取,可证明当f(x,y)满足一定条件时，迭代是收敛的。,7.1.2 改进的欧拉法及预测-校正公式,改进的欧拉法的预测校正公式,可证明预测校正公式的截断误差也为 o(h3)。,7.1.2 改进的欧拉法及预测-校正公式,例 取步长h=0.2,用改进的欧拉法的预测校正公式求解初值问题的数值解y1 , y2 .,解,预测-校正公式具体是,7.1.2 改进的欧拉法及预测-校正公式,设：改用后差商 替代方程中的导数项,7.1.2 向后 (隐式)欧拉公式,可以得到向后欧拉公式,这是隐式欧拉格式，也是一阶方法，精度与欧拉公式相当。计算yi+1通常用迭代法:,7.1.2 两步欧拉公式,设改用中心差商 替代方程 中的导数项 ，再离散化，即可导出下列格式,无论是显式欧拉公式还是隐式欧拉公式，它们都是单步法，其特点是计算时只用到前一步的信息yi，而该格式却调用了前面两步的信息yi-1,yi ，两步欧拉格式因此而得名。,两步欧拉格式具有更高的精度，它是二阶方法。,引言：（回顾）,本章解决的问题：一阶常微分方程的初值问题,7.1 欧拉法和改进的欧拉法,预测校正公式,改进的欧拉公式,欧拉公式,7.1 欧拉法和改进的欧拉法,两步欧拉公式,向后欧拉公式,7.2 龙格-库塔法(r-k法),7.2.1 二阶龙格-库塔公式,7.2.2 四阶龙格-库塔公式,引言：,公式构造思想：从泰勒公式出发，寻找更高阶的数值公式。,例如，泰勒公式计算到二阶可得,令,则,略去余项，得出一个二阶的数值公式为,因,引言：,理论上按此方式可以得到更高阶的公式。但需要计算复合函数的高阶导数，使算法复杂而不实用。,龙格库塔的思想（间接地运用泰勒公式）： 利用y(x)在若干个点上的函数值和导数值，作出一个适当的线性组合，使这个线性组合按h展开后的泰勒公式与y(x+h)的泰勒公式有较多的项达到一致，从而得出较高阶的数值公式。,引言：,r阶龙格库塔（ runge-kutta)法的一般形式：,7.2.1 二阶龙格-库塔公式,设想一个有二阶精度的数值公式形状为,a, b为待定系数。,仍令x=xi，则x+h=xi+1。如果能找出a,b,使得,略去余项就可得到上面所希望的近似计算公式了。,因此考虑,在h=0处求泰勒公式得,由于,7.2.1 二阶龙格-库塔公式,由t(h)的泰勒公式,7.2.1 二阶龙格-库塔公式,为使t(h)=o(h3),令,解出,，得,整理得,7.2.1 二阶龙格-库塔公式,利用,可以推出,取x=xi并略去o(h3)便得到二阶龙格库塔公式,或,7.2.1 二阶龙格-库塔公式,可以推出,取x=xi并略去o(h3)便得到二阶龙格库塔公式,或,7.2.2 四阶龙格-库塔公式,仿照上述的讨论，可导出四阶龙格库塔公式:,例 取步长h=0.2，用四阶龙格库塔公式求下面初值问题的数值解。,7.2.2 四阶龙格-库塔公式,解,由公式得,7.2.2 四阶龙格-库塔公式,数值解yi与准确解y(xi)的对照见表,准确解是,xi,yi,y(xi),7.3 线性多步法,7.3.1 四阶阿达姆斯(adams)外插公式,7.3.2 四阶阿达姆斯(adams)内插公式,7.3.0 多步法的概念,7.3.3 初始出发值的计算,7.3.4 阿达姆斯预测-校正公式,7.3.0 多步法的概念,单步法 -计算yi+1时只使用yi的值。 多步法 -计算yi+1时使用前面的k个y值，即由 yi-k+1 , yi-k+2 , ,yi-1, yi计算yi+1。(k=1,2,) 线性多步法 -计算yi+1的公式由yi-k+1 , yi-k+2 , ,yi-1, yi 的线性组合表达。,7.3.1 四阶阿达姆斯外插公式,设想用yi-3 ,yi-2 ,yi-1 ,yi 的值计算yi+1。为方便讨论由,出发计算y(x+h), 由初值问题的方程y=f(x,y(x)两边从x到x+h积分, 可得到等价的积分方程,7.3.1 四阶阿达姆斯外插公式,设想运用数值积分方法，取 x-3h, x-2h, x- h, x 为插值基点，做 f(s,y(s) 的三次拉格朗日插值，用它近似计算上式的积分 。 这样得到的数值积分公式是f(s,y(s)在4个插值基点处的函数值的线性组合。,由于 f(x-ih,y(x-ih)= y(x-ih)，所得到的计算 y(x+h) 的近似公式形为：,7.3.1 四阶阿达姆斯外插公式,为达到四阶精度，希望确定参数b0 ,b1 ,b2 ,b3使满足,运用在 h=0 处的泰勒公式得,7.3.1 四阶阿达姆斯外插公式,为达到四阶精度，希望确定参数b0 ,b1 ,b2 ,b3使满足,运用在 h=0 处的泰勒公式得,7.3.1 四阶阿达姆斯外插公式,为使误差等于o(h5)，令h, h2,h3,h4 的系数为0， 得方程组：,求得,代入前面的公式得,7.3.1 四阶阿达姆斯外插公式,令x=xi 并记,7.3.1 四阶阿达姆斯外插公式,略去余项，得到四阶阿达姆斯外插公式 ：,这是显式公式。公式的截断误差为o(h5)。,7.3.2 四阶阿达姆斯内插公式,把四阶阿达姆斯外插公式中使用的yi-3 ,yi-2 ,yi-1 ,yi 改为 ：,yi-2 ,yi-1 ,yi ,yi+1 ,经类似的推导可得近似公式：,确定待定系数的方程组为,7.3.2 四阶阿达姆斯内插公式,解得,得到四阶阿达姆斯内插公式,这是一个隐式公式，截断误差也是o(h5)。,7.3.3 初始出发值的计算,阿达姆斯公式的特点是计算公式简单，只需简单的算术运算，计算量少，结果的精度较高。但需要4个初始出发值。,（1）使用单步法，例如龙格库塔法求出发值。,（2）使用y(x)在x=x0 处的泰勒公式,其中泰勒公式的阶数k按需要选取，各导数值由复合函数 f(x,y(x)的求导得出。,7.3.4 阿达姆斯预测-校正公式,阿达姆斯内插