1、第九章 常微分方程初值问题的数值解法第一部分 内容提要一、数值解的一般概念常微分方程初值问题 的数值解是指通过一定的近似方法0(),)yxfy得出准确解 在一列离散点 上的近似值()yx12,nx 。数值解的特征是步进式，即 在 点的近似值012,ny ()y1n是由 等若干点处的近似值 的信息给出的递推公式。n1x 1,n若 依赖于前面 步的值 ，则称为 步法； 称为单步1yk1,nky k法。利用 在 的精确解 借助某种()x11,nnkx 11(),()nnkyxyx算法计算出 ，则称 为该方法的局部截断误差。如果一个算1y()y法的局部截断误差是 ，则称该方法是 阶的；而利用数值解1p
2、Ohp得到的 与微分方程的精确解之差 称为整体11,nnky 1ny 1()nyx截断误差，即是该数值方法的误差。对于固定的 ，取 ，用某种算法得到 ，如有0x0xhn=0，则称该方法是收敛的。注意，因 是固定的，随着 ，0lim()nhy x0h数值解的步数 。在实际计算时由于舍入误差不可避免，实际得到数值解是 ，稳定性ny即研究 是否随着计算步骤 的增加而增加。通常所提的稳定性是通ny n过模型方程 来讨论的。若当某一步 有舍入误差时，在以后(0)ny的计算中误差不会逐步扩大，则称这种稳定性为绝对稳定性。二、简单单步法及其收敛性、稳定性Euler 法 的局部截断误差为 ，整体截断误差为1(
3、,)nnyhfxy2()Oh，即一阶收敛。对于模型问题 ，当 时，Euler()Oh (0y法是数值稳定的。 隐式 Euler 法 的误差与 Euler 法相同，但是无条11(,)nnyhfxy件稳定：即对任意步长 ，隐式 Euler 法都是稳定的。0梯形法 的误差比 Euler 法高一阶，1 1(,)(,)2nnnff 也是无条件稳定的。改进 Euler 法 是一种预测1 1(,)(,(,)nnnnhyfxyfyhfx 校正方法： Euler 法预测p梯形法校正1 1(,)(,)2c pnnnyfxyfy 它保持了梯形法的误差阶数，但不是无条件稳定的。三、龙格-库塔方法龙格-库塔类算法采用区
4、间 内若干点的斜率的加权平均来近似整1,nx个区间的平均斜率，一般形式为111 1(,), ),2,snniin iininiijjjyhKKfxyfchchaKis如经典的 4 级 4 阶（局部截断误差为 ）Runge-Kutta 公式为()s5()O112342 13 24 3( )6(,),)2(,)nnnnnnhyKKKfxyhfxyKKh（1）四、线性多步法线性 步法具有一般形式k1011011( )nnknnnknyyyhfff （2）为隐式公式； 为显式公式。构造多步法公式有基于数值积分和00Taylor 展开两种途径。多步法（2）的局部截断误差为 11 1100()()(,)k
5、knnjnjjnjnjj jTyxyxhfxy利用原微分方程后，成为 11 100()()().kknnjnjjnjj jyxyxhyx因此利用 Taylor 公式，分别对 和 作 Taylor 展开，可确定线性多步（2）中的待定参数 和 ，使它达到最高阶精度（或指定精度） 。10kjkj0预测校正格式：不论单步法还是多步法，隐式公式比显式公式的稳定性好，但隐式公式的计算比较困难。预测校正格式是用显式公式进行预测，再用隐式公式进行校正。五、高阶方程和一阶常微分方程组所有应用于一阶常微分方程初值问题的数值方法都可以直接推广到一阶常微分方程组，只需把公式中的未知函数改为向量形式。高阶方程总可以降为
6、一阶方程组，进而用数值方法求解。第二部分 例题精选一、对初值问题 ， ，用梯形公式求解，求数值解 的表达式y 1)0( ny（写成步长 的函数） ，并证明数值解的收敛性。h分析：微分方程用数值方法离散后即变成差分方程，单步法导出的差分方程通常是关于数值解序列 的一个递推公式，因此问题变为已知数列首项和递推ny公式求数列通项公式，收敛性即是数列极限。解：梯形公式 对本问题)(211nnfhnyf于是梯形公式的解 即11nnyy nh2由于 ，易得10y 、 321,2nhyn显然初值问题的准确解为 xe)(对于固定的点 ，准确解在该点的值为nx nxey)(而数值解为 hxn nhy212令 固
7、定，上式对 取极限得：nx0)(2122 nxhxhn yey nn 这证明了数值解的收敛性。二、用隐式单步法（该方法也属于隐式 Runge-Kutta 方法）),(3)(412121hkyxfkynn求解微分方程初值问题 ， 时，试推出其绝对稳定区间。y5)0(分析：将格式应用于所给方程可导出误差传播方程，从而求出绝对稳定区间解：记舍入误差为 。该隐式单步法应用于方程 时n y5hykhyknn5313511 nnnn y10712 nnn hyhy62553041 从而误差传播方程为 nn1021解不等式 得10625h560h故，该方法应用到所给方程的绝对稳定区间为 。,0(三、用 Ta
8、ylor 展开原理构造形如 1101()()nnnyyhf的两步法，试确定系数 使方法具有二阶精度，并推导其局部截断误差0,主项。分析：本题考察构造多步法的方法；二阶精度即局部截断误差为 。3()Oh解：由多步法局部截断误差的定义 11 1011()()()()().(nnnnnnTyxyxhyxyx 将 ， 和 分别在点 展开得：)1x2341()()()()6nnnnnnhyxhyyxOh 2341()()()()()nnnnnnx 231()()()()nnnnnhyxhyxyxOh将此三式代入局部截断误差(3)式，我们有 23 41 0111(2)()()2)()3)(6nnn nnT
9、yxyxyxOh 要使方法具有二阶精度，必需 01112,0.解得 。017,24此时 ，所以局部截断误差主项为 。1383()8nhyx四、取 ，试用 Euler 法求解初值问题 。0.1h240,.2(0)1,()yxyx分析：先将二阶方程写为一阶方程组，再用 Euler 法求解。解：引进 ，并记向量函数 ，则原二阶方程变为一阶方程组,uyvuvy，它满足初值条件 。2(,)4xvxv yf 0(0)1uvyy将 Euler 法 应用到该方程组得：1,nnhy100101(,).2.xfy211 .098(,).204().32h yf即，用 Euler 法求得原二阶方程的数值解为 。12,.y
