常微分方程组的数值解.ppt

在工程和科学技术的实际问题中，常需求解微分方程，但常微分方程中往往只有少数较简单和典型的常微分方程（例如线性常系数常微分方程等）可求出其解析解,对于变系数常微分方程的解析求解就比较困难，而一般的非线性常微分方程的求解困难就更不用说了。大多数情况下，常微分方程只能用近似方法求解。这种近似解法可分为两大类：一类是近似解析法，如级数解法、逐次逼近法等;另一类是数值解法，它给出方程在一些离散点上的近,第八章 常微分方程数值解法,1、 引 言,似值。,其中 x 是质量，m是离开平衡位o的距离，t为时间，c为弹簧系数。,在具体求解微分方程时，需具备某种定解条件，微分方程和定解条件合在一起组成定解问题。定解条件有两种：一种是给出积分曲线在初始点的状态，称为初始条件，相应的定解问题称为初值问题。另一类是给出积分曲线首尾两端的状态，称为边界条件，相应的定解问题称为边值问题。,我们现在讨论常微分方程的数值解法。先从最简单的一阶常微分方程的初值问题出发开始讨论。,由常微分方程理论可知：只要上式中的函数f(x,y)在区域 G=axb,y内连续，且关于 y 满足Lipschitz条件，即存在与 x, y 无关的常数L，使,至于初值问题（1）的数值解法，常采用差分方法，即把一个连续的初值问题离散化为一个差分方程来求解。具体地，将（1）离散化后，求找其解 y=y(x)在一系列离散点,下面分析均假定满足上述条件。,一、Euler公式,2 Euler方法,因为初值问题中的初始条件 已知，即可利用已知的 来求出下一节点处 的近似值 ；再从 来求 如此继续，直到求出 为止。这种用按节点的排列顺序一步一步地向前推进的方式求解的差分算法称为“步进式”或“递推式”算法，它是初值问题数值解法的各种差分格式的共同特点。因此，只要能写出由前几步已知信息 来计算 的递推公式（即差分格式），即可完全表达这种算法。,若将 和 的近似值分别记为 和 ， 则得： （3）,这就是Euler公式（格式）。 利用它可由初值 出发逐步算出 。 这类形式的方法也称为差分方法。,定义：如果局部截断误差为 ，则这种数值算法的精度为p阶，故Euler格式的精度为一阶。 从几何意义上来看，如图，,当假定 为准确值，即在 的前提下来估计误差 ，这种截断误差称为局部截断误差。,由（2）、（3）知Euler公式在 处的局部截断误差为：,由方程（1）知，其积分曲线 y=f(x)上任一点(x, y)的切线斜率 都等于函数f(x, y)的值。从初始点 （即 点）出发，作积分 曲线 y=y(x) 在 点上的切线 （其斜率为 ）与直线 相交与 点（即 点），得到 作为 的近似值，则有,相比较知，这时用切线 近似代替了曲线段 点近似 代替了 点， 近似代替了 近似代替了 。 递推继续从 点出发，作一斜率为 的直线 与直线 相交于 点（即 点），得到 作为 的近似值 。如此直到 点。这样得出一条折线 近似代替积分曲线 ，当 步数越多时，由于误差的积累，折线 可能会越,解：为便于进行比较，我们后面将用多种数值方法求解上述初值问题。 这里先用Euler 公式，此处具体格式为： 取步长为h=0.1，计算结果略。 由结果可见Euler方法的精度很差。,即为Euler格式（3）。,因为差商是微分的近似，所以Euler格式也可用差商近似代替导数的离散方法来得到。在节点 处有：,二 后退Euler格式,显然Euler格式具有递推性，在计算 时只要用到前一步所得的结果 一个信息就够了，因此是一种单步格式或称一步格式。 若用不同的数值微分计算方法也可导出其它形式的算法。 例如：用向后差商表示的数值微分公式,（6）称为向后Euler公式，又称为隐式Euler公式(后退Euler格式)。后退Euler公式与Euler公式有着本质的区别，后者是关于 的一个直接计算公式，这类公式称作显式的；而前者，即（6）中右端含有未知的 ，它实际上是关于 的一个函数方程。这类公式称作隐式的。,显式与隐式两类方法各有特点，使用显式算法远比隐式算法方便，但考虑数值稳定性等因素，人们常选用隐式算法。,隐试算法（6）常用迭代法来实现，而迭代过程实质上是逐步 显式化。,设用显式Euler格式算出 作迭代初值 ，以此代入（6）右端，使之转化为显示，直接计算得： ，再用 代入（6）右端又有：,从几何上看，梯形公式是取 区间两端点处斜率的 平均斜率。,三 梯形公式,Euler方法是过 点以斜率 引直线交 的点A。 后退Euler方法是以点 处的斜率 为斜率，从 点引直线交 于另一点B。 A、B两点都偏离点 点，梯形方法就是取A、B两点的中点P作为Q2近似，上图表明梯形方法确实改善了精度。,从而也可得二步Euler公式及其 截断误差为：,也可以由导数的中心差分近似式得到：,四 二步欧拉公式,将区间a,bn等分，得子区间 在第i+1 个子区间上，积分,例如：对于初值问题：,对右端利用左矩形公式可得 即 Euler格式,各公式的截断误差可直接利用数值积分截断误差估计而得。从而可知梯形公式（8）的截断误差为：,梯形公式也是隐式的，可用迭代法求解，与后退Euler方法一样，仍用Euler方法提供迭代初值，其迭代格式为：,为分析迭代过程的收敛性，将（12）与（8）相减得：,（12）,k=0，1，2，,L为f(x,y)关于y的Lipschitz常数.如果选取h充分小使得 则 . 时有 这表明迭代过程（12）是收敛于（8）的解的。,五 改进的Euler公式,上面已看到Euler公式计算量小但精度差，梯形方法虽然提高了 计算精度，但算法复杂计算量大，在应用（12）进行迭代时，每次 均要计算函数f的值，而迭代又要反复进行多次，计算量很大，难以 预测。为了控制计算量，通常只迭代一两次就转入 下一步计算， 这就简化了计算。,具体地，先用Euler公式求得一个初步的近似值 称之为预测 值。预测值 的精度可能很差。再用梯形公式（8）将它校正 一次，即按（12）式迭代一次得 ，这个结果称为校正值，这 样建立的 预测校正系统称为改进的Euler公式。,由表可见，与精确解 相比，改进的Euler公式的精度较Euler 公式有明显的提高。 下面再看两步Euler公式（9），除了给出初值 外，还需要借助 其它单步法（如Euler公式，后退Euler公式及梯形公式等）再提供一个,Euler公式,改进的Euler公式,精确解,0,1,1,1,0.1,0.2,0.3,1,0.9000000,0.8100000,0.7290000,0.3486784,0.9050000,0.8190250,0.7412176,0.3685410,0.9048374,0.8187308,0.7408182,0.3678794,启动值 然后才能启动计算公式依次计算,用两步Euler公式与梯形公式相匹配，又可得到下面预测-校正 系统：,校正：,（17）,（18）,两步法优美是由于它调用了两个节点上的信息，从而能以较少的 计算量获得较高的精度。,预测：,与改进的,Euler公式（13）（14）相比较易见（17）（18）的一个,突出特点是它的预测公式与校正公式具有同阶精度。据此可以比较 方便的估计截断误差，并基于这种估计，可以提供一种提高精度的 简易方法。,再由梯形公式截断误差公式（11）知：校正公式（18）的截断误差 为： 比较（19）（20）可见，校正值的误差 大约只是预测值 的误差 的1/4(符号相反).即 ,由此导出误差的事后估计:,若预测公式（17）中的 和 都是准确的。即 则由两步Euler公式的截断误差公式（10）知：,1预测: 2改进: 3计算:,可以期望利用这样估计出的误差作为计算结果的一种 补偿,有 可能使精度得到改善. 以 和 分别表示第I步的预测值和校正值,按估计式(21), 和 分别作为 和 的改进 值,在校正值 尚未求出之前,可用上一步的偏差 替代 来改进预测值 .这样设计的计算方案有如下六个 环节:,3 Runge-Kutta方法,运用上述方案计算 时，要用到先一步的信息 ， ， 和 更前一步的 。因此启动算法之步必须给出启动值 和 。 可用其它单步法（例如改进的Euler方法）来计算， 则一般取为0。计算结果表明，这种简单的处理方法通常 可以获得令人满意的效果。,而具体则有:,式中y（x）的各阶导数可由方程（*）用函数f来表达。引进函数 序列 来描述求导过程：,二 Runge-Kutta公式的导出,例:用Taylor公式求解初值问题:,K可看作是 y=y(x)在区间 上的平均斜率。从而Euler公式 相当于取 点上的斜率 作为平均斜率K的 近似值，这当然十分粗糙，因而精度必然很低。,这个过程启示我们，如果设法在 内多预测几个点的斜率 值，然后将它们加权平均作平均斜率K的近似值，就有可能构造 出更高精度的计算公式。这就是Runge-Kutta方法的基本思想。,根据预测值 再来算出 由此构造出计算格式：,假定 ，分别将 和 作Taylor展开得：,由（2）得：,（7）中三个待定 参数 P ，但只有两个方程，因此还有一个自由度。凡满足条件（7）的一族格式统称为二阶Runge-Kutta格式。,当p=1 ， 时，二阶Runge-Kutta格式（6）即为改进的 Euler格式（15）。 如取p= 1/2 ，则 ，二阶R-K格式（6）成为：,称之为变形的Euler格式 。,由于（8）中的 是由Euler格式预测出来的区 间 中的点 的近似解， 就近似地等于此中点的斜率 ，因此（8）就相当于用中点 的 斜率作为（4） 中平均斜率K的近似值，故格式（8）也称为中点格式。,总之，二阶R-K格式用多算一次函数值f的办法，避开了二 阶Taylor级数法所要求计算的 f的导数值，在这种意义上可以说， R-K方法实质上是Taylor方法的变形。,三 三阶Runge-Kutta方法,斜率值 和 可以利用。我们用 和 线性组合给出区间 上的平均斜率，从而得到 的预测值,于是，再通过计算函数值f得到：,这样设计出的计算格式具有形式：,（9）,我们希望适当选择系数 和p、q、r、 使上述公式具有三阶精 度。为便于数学演算，引进算子：,，,则根据（2）有：,于是三阶展开,可表为：,式中 的下标 i 均表示在 处取值。,（12）,四 四阶Rung-Kutta方法,最常用的一种经典Rung-Kutta格式具体形式如下：,(13),解：三种方法如下：,Euler格式：,改进的Euler格式：,(h=0.05),经典的R-K格式：,(h=0.1),Euler法h=0.025,改进Eulerh=0.05,R-K法h=0.1,精确解y,x=0.1,x=0.2,x=0.3,0.903687890,0.816651803,0.737998345,0.904876562,0.818801593,0.74091437,0.904837500,0.818730901,0.740818,0.90483748,0.81873075,0.74081822,x,这里采用了不同的步长h值，是为了使他们所耗的计算工作量 大致相同，以便于比较。由上表可见，经典的R-K方法的精确度 较改进的Euler方法又有很大的提高。这一结论也可以从理论上 大致的分析出来。,析出来：Euler方法的局部截断误差为：,计算四步后的,而经典R-K方法的局部截断误差则为：,可见，当,为大致相同数量级的常数时有：,但要注意的是：R-K方法的导出利用了Taylor展开，因此要求 所求的解有教好的光滑性，如果解的光滑性差，则采用经典的 R-K方法所的数值解，其精度有可能反而不及改进的Euler方法， 因此在实际计算中应根据问题的具体情况来选择适合的算法。,五 步长的自动选择-变步长的Runge-Kutta方法,在应用数值法求解微分方程中，选择适当的步长是至关重要的。步 长太大则达不到要求，步长太小则步数增多，不但增加计算工作量， 还可能导致舍入误差的严重积累。尤其是当微分方程的解y(x)变化 激烈时，步长的合理取法是在变化激烈处步长取小些，在变化平缓 时取大些，也就是采取自动变步长的方法，即根据精度的要求先估 计出下一步长的合理大小，然后按此计算。,作为近似值，则,的精确度都要高。,当p=4时，可以取,这种修正方法与Romberg的数值积分的思路是一样的。（15） 除以（14）得：,由此可以得出误差,事后估计式：,由上面的分析可见，微分方程数值解的基本思想是，通过,（1） 如果 ，则反复加倍步长进行计算，直到 时为止，并以上一次步长的计算结果作为 （2）若 ， 则反复减半步长进行计算，直到 时为止，并取其最后一次步长的计算作为 这样做时 ，为了选择步长，每一步都要反复判别， 增加 了工作量，但在方程的解 y(x) 变化剧烈的情况时 ，总的计算 工作量可以得到减少，结果还是合算的。,这样就可以从步长减半前后的两次计算结果的偏差 来判断步长选的是否适当，当要求的 数值精度为 时：,4 单步法的收敛性和稳定性,一 单步法的收敛性,本例的Euler公式为 由此式递推可得：,定义 ： 若一种数值方法对任意固定的 当 （同时 ） 时 有 ，则称方法是收敛的。,某种离散化手段，将微分方程转化为差分方程来求解。这种转化是 否合理，还要看差分方程问题的解 当 时是否会收敛到点 对固定的i将趋向 ，这时讨论收敛是没有意义的。,此定理表明，判断单步法（1）的收敛性，归结为验证增量函数 能否满足Lipschitz（3） 对于Euler方法，由于其增量函数 故当f满足 Lipschitz 条件时它是收敛的。,稳定性问题比较复杂，为简化讨论，仅考察下面的模型方程 为保证微分方程本身的稳定性，假设 先考虑Euler方法的稳定性。模型方程 的Euler公式为： 设在节点值 上有一扰动值 , 它的传播使节点值 产生大小为： 的扰动值，假设用 按Euler 公式得出：,定义：若一种数值方法在节点值 上有扰动 ，而对于yi后的各个 节点值 上产生的 偏差均不超过 ,则称该方法是稳定的。,二 单步法的稳定性,0.025,0.050,0.075,0.100,1,2,3,4,5,0,-1,-2,-3,x,y,其方程时间常数为,因此有（10）知，要使Euler法稳定 则步长,如果取步长h=0.025则Euler格式为：,其结果在准确解,上下波动不稳定，再看后退Euler格式,H=0.025时，形式为：,计算结果稳定。具体结果如下：,节点,Eule方法,后退Euler方法,0.025,-1.5,0.2857,0.050.,2.25,0.0816,0.075,-3.375,0.0233,0.100,5.0625,0.0067,前面介绍的几种步进方法，在计算 时大多只用到前一个节点上 的近似值 ，而没有用到前几步的计算所得出的信息，故称为单步 法 。实际上经过多次单步法计算后，已得出一系列近似值 等。为了充分利用这些信息来计算 以减少计算 工作量和获得教高的精度，可采用如下计算公式：,5 线性多步法,一 显式Adams法,（1）,例如 k=3 时有：,（7）,（7）称为Adams四步显式法。它用到了四个节点上的f值，是,一种最常用的多步法，其精度为四阶。,如果利用,共k+1个数据来构造一个Newton内插多次式,，则与上面类似推导可得：,局部截断误差为：,二 Adams隐式公式, 在同一阶数下，隐式的局部截断误差的常数的绝对值 比显式的|Bk|要小。,三 Adzms 显式与隐式的比较,当k=2时，,解：取h=0.1，两种方法具体算式如下：,算法启动值（对四步显式当x1=0.1, x2 =0.2, x3 =0.3时，三步隐式，当x1 =0.1, x2 =0.2时）可用同阶的RungeKutta公式计算。而本例则是由精确解 算出的。计算结果部分如下表,如下列定解问题： 其Ada