第八章节常微分方程组的数值解课件

《第八章节常微分方程组的数值解课件》由会员分享，可在线阅读，更多相关《第八章节常微分方程组的数值解课件（83页珍藏版）》请在金锄头文库上搜索。

1、在工程和科学技术的实际问题中，常需求解微分方程，但常微分方程中往往只有少数较简单和典型的常微分方程（例如线性常系数常微分方程等）可求出其解析解,对于变系数常微分方程的解析求解就比较困难，而一般的非线性常微分方程的求解困难就更不用说了。大多数情况下，常微分方程只能用近似方法求解。这种近似解法可分为两大类：一类是近似解析法，如级数解法、逐次逼近法等;另一类是数值解法，它给出方程在一些离散点上的近,第八章 常微分方程数值解法,1、 引 言,似值。,其中 x 是质量，m是离开平衡位o的距离，t为时间，c为弹簧系数。,在具体求解微分方程时，需具备某种定解条件，微分方程和定解条件合在一起组成定解问题。定解

2、条件有两种：一种是给出积分曲线在初始点的状态，称为初始条件，相应的定解问题称为初值问题。另一类是给出积分曲线首尾两端的状态，称为边界条件，相应的定解问题称为边值问题。,我们现在讨论常微分方程的数值解法。先从最简单的一阶常微分方程的初值问题出发开始讨论。,由常微分方程理论可知：只要上式中的函数f(x,y)在区域 G=axb,y内连续，且关于 y 满足Lipschitz条件，即存在与 x, y 无关的常数L，使,至于初值问题（1）的数值解法，常采用差分方法，即把一个连续的初值问题离散化为一个差分方程来求解。具体地，将（1）离散化后，求找其解 y=y(x)在一系列离散点,下面分析均假定满足上述条件。

3、,一、Euler公式,2 Euler方法,因为初值问题中的初始条件 已知，即可利用已知的 来求出下一节点处 的近似值 ；再从 来求 如此继续，直到求出 为止。这种用按节点的排列顺序一步一步地向前推进的方式求解的差分算法称为“步进式”或“递推式”算法，它是初值问题数值解法的各种差分格式的共同特点。因此，只要能写出由前几步已知信息 来计算 的递推公式（即差分格式），即可完全表达这种算法。,若将 和 的近似值分别记为 和 ， 则得： （3）,这就是Euler公式（格式）。 利用它可由初值 出发逐步算出 。 这类形式的方法也称为差分方法。,定义：如果局部截断误差为 ，则这种数值算法的精度为p阶，故Eu

4、ler格式的精度为一阶。 从几何意义上来看，如图，,当假定 为准确值，即在 的前提下来估计误差 ，这种截断误差称为局部截断误差。,由（2）、（3）知Euler公式在 处的局部截断误差为：,由方程（1）知，其积分曲线 y=f(x)上任一点(x, y)的切线斜率 都等于函数f(x, y)的值。从初始点 （即 点）出发，作积分 曲线 y=y(x) 在 点上的切线 （其斜率为 ）与直线 相交与 点（即 点），得到 作为 的近似值，则有,相比较知，这时用切线 近似代替了曲线段 点近似 代替了 点， 近似代替了 近似代替了 。 递推继续从 点出发，作一斜率为 的直线 与直线 相交于 点（即 点），得到 作

5、为 的近似值 。如此直到 点。这样得出一条折线 近似代替积分曲线 ，当 步数越多时，由于误差的积累，折线 可能会越,解：为便于进行比较，我们后面将用多种数值方法求解上述初值问题。 这里先用Euler 公式，此处具体格式为： 取步长为h=0.1，计算结果略。 由结果可见Euler方法的精度很差。,即为Euler格式（3）。,因为差商是微分的近似，所以Euler格式也可用差商近似代替导数的离散方法来得到。在节点 处有：,二 后退Euler格式,显然Euler格式具有递推性，在计算 时只要用到前一步所得的结果 一个信息就够了，因此是一种单步格式或称一步格式。 若用不同的数值微分计算方法也可导出其它形

6、式的算法。 例如：用向后差商表示的数值微分公式,（6）称为向后Euler公式，又称为隐式Euler公式(后退Euler格式)。后退Euler公式与Euler公式有着本质的区别，后者是关于 的一个直接计算公式，这类公式称作显式的；而前者，即（6）中右端含有未知的 ，它实际上是关于 的一个函数方程。这类公式称作隐式的。,显式与隐式两类方法各有特点，使用显式算法远比隐式算法方便，但考虑数值稳定性等因素，人们常选用隐式算法。,隐试算法（6）常用迭代法来实现，而迭代过程实质上是逐步 显式化。,设用显式Euler格式算出 作迭代初值 ，以此代入（6）右端，使之转化为显示，直接计算得： ，再用 代入（6）右

7、端又有：,从几何上看，梯形公式是取 区间两端点处斜率的 平均斜率。,三 梯形公式,Euler方法是过 点以斜率 引直线交 的点A。 后退Euler方法是以点 处的斜率 为斜率，从 点引直线交 于另一点B。 A、B两点都偏离点 点，梯形方法就是取A、B两点的中点P作为Q2近似，上图表明梯形方法确实改善了精度。,从而也可得二步Euler公式及其 截断误差为：,也可以由导数的中心差分近似式得到：,四 二步欧拉公式,将区间a,bn等分，得子区间 在第i+1 个子区间上，积分,例如：对于初值问题：,对右端利用左矩形公式可得 即 Euler格式,各公式的截断误差可直接利用数值积分截断误差估计而得。从而可知

8、梯形公式（8）的截断误差为：,梯形公式也是隐式的，可用迭代法求解，与后退Euler方法一样，仍用Euler方法提供迭代初值，其迭代格式为：,为分析迭代过程的收敛性，将（12）与（8）相减得：,（12）,k=0，1，2，,L为f(x,y)关于y的Lipschitz常数.如果选取h充分小使得 则 . 时有 这表明迭代过程（12）是收敛于（8）的解的。,五 改进的Euler公式,上面已看到Euler公式计算量小但精度差，梯形方法虽然提高了 计算精度，但算法复杂计算量大，在应用（12）进行迭代时，每次 均要计算函数f的值，而迭代又要反复进行多次，计算量很大，难以 预测。为了控制计算量，通常只迭代一两次

9、就转入 下一步计算， 这就简化了计算。,具体地，先用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

10、.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、基于这种估计，可以提供一种提高精度的 简易方法。,再由梯形公式截断误差公式（11）知：校正公式（18）的截断误差 为： 比较（19）（20）可见，校正值的误差 大约只是预测值 的误差 的1/4(符号相反).即 ,由此导出误差的事后估计:,若预测公式（17）中的 和 都是准确的。即 则由两步Euler公式的截断误差公式（10）知：,1预测: 2改进: 3计算:,可以期望利用这样估计出的误差作为计算结果的一种 补偿,有 可能使精度得到改善. 以 和 分别表示第I步的预测值和校正值,按估计式(21), 和 分别作为 和 的改进 值,在校正值 尚未求出之前,可用上一步的偏差 替代 来改进预测值 .这样

12、设计的计算方案有如下六个 环节:,3 Runge-Kutta方法,运用上述方案计算 时，要用到先一步的信息 ， ， 和 更前一步的 。因此启动算法之步必须给出启动值 和 。 可用其它单步法（例如改进的Euler方法）来计算， 则一般取为0。计算结果表明，这种简单的处理方法通常 可以获得令人满意的效果。,而具体则有:,式中y（x）的各阶导数可由方程（*）用函数f来表达。引进函数 序列 来描述求导过程：,二 Runge-Kutta公式的导出,例:用Taylor公式求解初值问题:,K可看作是 y=y(x)在区间 上的平均斜率。从而Euler公式 相当于取 点上的斜率 作为平均斜率K的 近似值，这当然

13、十分粗糙，因而精度必然很低。,这个过程启示我们，如果设法在 内多预测几个点的斜率 值，然后将它们加权平均作平均斜率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）,我们希望适当选择