1、第 3 章 数值算法之一：单步法接下来的章节，我们将详细讲解常微分方程的数值求解方法，特别是编程实现方程的初值问题，这也将成为解决实际问题的必要基础和课程设计的主要内容。本章重点讲解一阶常微分方程的数值算法中最简单的一类方法：单步法。首先介绍 Euler 法、后退 Euler 法与梯形法，并分析单步法的局部截断误差，然后给出了改进的 Euler 法。3.1 简单的单步法及基本概念3.1.1 Euler 法、后退 Euler 法与梯形法求初值问题(1.2.1)的一种最简单方法是将节点 的导数 用差商nx()nyx代替，于是(1.2.1)的方程可近似写成()(nnyxh(3.1.1)从 出发 ，由
2、(3.1.1)求得 再将代入(3.1.1)右端，得到 的近似 ，一般写成(3.1.2)该数值解法称为解微分方程初值问题的 Euler 法.Euler 法的几何意义如图 3-1 所示。初值问题(1.2.1)的解曲线 y=y(x)过点，从 出发，以 为斜率作一段直线，与直线 交点于，显然有 ，再从 出发，以 为斜率作直线推进到 上一点 ，其余类推，这样得到解曲线的一条近似曲线，它就是折线.图 3-1 Euler 法的几何意义显示Euler 法也可利用 的 Taylor 展开式得到，由(3.1.3)略去余项，以 ，就得到近似计算公式(3.1.2).另外，还可对(1.2.1)的方程两端由 到 积分得(
3、3.1.4)若右端积分用左矩形公式，用 ， ，则得(3.1.2).如果在(3.1.4)的积分中用右矩形公式，则得(3.1.5)该算法称为后退(隐式)Euler 法。若在(3.1.4)的积分中用梯形公式，则得(3.1.6)该算法称为梯形方法。上述三个公式(3.1.2)，(3.1.5)及(3.1.6)都是由 计算 ，这种只用前一步即可算出 的公式，我们称之为单步法，其中(3.1.2)可由 逐次求出的值，称为显式方法，而(3.1.5)及(3.1.6)右端含有 当 f对 y 非线性时它不能直接求出 ，此时应把它看作一个方程，求解 ，这类方法称为隐式方法。此时可将(3.1.5)或(3.1.6)写成不动点
4、形式的方程这里对式(3.1.5)有 ，对(3.1.6)则 ，g 与 无关，可构造迭代法(3.1.7)由于 对 y 满足 Lipschitz 条件(1.1.2)，故有当 或 ，迭代法(3.1.7)收敛到 ，因此只要步长 h 足够小，就可保证迭代(3.1.7)收敛。对后退 Euler 法(3.1.5)，当 时迭代收敛，对梯形法(3.1.6)，当 时迭代序列收敛。例 3.1 用 Euler 法、隐式 Euler 法、梯形法解取 h=0.1，计算到 x=0.5，并与精确解比较。解 直接用给出公式计算。由于 ，Euler 法的计算公式为n=0 时， .其余 n=1,2,3,4 的计算结果见表 3-1.对
5、隐式 Euler 法，计算公式为解出当 n=0 时， .其余 n=1,2,3,4 的计算结果见表 3-1.表 3-1 例 3.1 的三种方法及精确解的计算结果对梯形法，计算公式为解得当 n=0 时， .其余 n=1,2,3,4 的计算结果见表7-1.本题的精确解为 ，表 3-1 列出三种方法及精确解的计算结果。附：具体三种算法程序如下：首先定义一阶方程右端函数 f 如下function Y=f(x,y)Y=-y+x+1;（1）Euler 法function y = DEEuler(f, h,a,b,y0,varvec)%这是 Euler法的函数命令%一阶常微分方程的一般表达式的右端函数：f 这
6、里可用 f=inline% 自变量下限 a；上限 b% 函数初值 y0% 积分步长 h% 常微分方程的变量组 varvecformat long;%数据显示方式，不影响计算和存储方式，%是指小数点后 15 位数字表示N = (b-a)/h;y = zeros(N+1,1);x = a:h:b;y(1) = y0;for i=2:N+1y(i) = y(i-1)+h*Funval(f,varvec,x(i-1), y(i-1);%简单 Euler 公式迭代， %本题单步法的格式为：y(i) = y(i-1)+h.*(-y(i-1)+x(i-1)+1); endformat short;本题输入命
7、令为 Y = DEEuler(-v+u+1, 0.1,0,0.5,1,u v)（2）隐式 Euler 法function y = DEimpEuler(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1fx = Funval(f,var(1),x(i);gx = y(i-1)+h*fx - varvec(2);y(i) = NewtonRoot(gx,-10,10,eps);endformat short;（3）梯形法functi
8、on y = DEimpEuler1(f, h,a,b,y0)format long;N = (b-a)/h;y = zeros(N+1,1);x = a:h:b-h;y(1) = y0;var = findsym(f);yd = diff(f,var(4);for i=2:N+1fx = f - yd*var(4);dfy = subs(yd,findsym(yd),x(i-1);cx = subs(fx,findsym(fx),x(i-1);y(i) = (y(i-1)+h*cx)/(1-h*dfy);endformat short;3.1.2 单步法的局部截断误差解初值问题(1.1.1)
9、的单步法可表示为(3.1.8)其中 与 有关，称为增量函数，当 含有 时，是隐式单步法，如(3.1.5)及(3.1.6)均为隐式单步法，而当不含 时，则为显式单步法，它表示为(3.1.9)如 Euler 法(3.1.2)， .为讨论方便，我们只对显式单步法(3.1.9)给出局部截断误差概念.定义 2.1 设 y(x)是初值问题(7.1.1)的精确解，记(3.1.10)称为显式单步法(3.1.9)在 的局部截断误差.之所以称为局部截断误差，可理解为用公式(3.1.9)计算时，前面各步都没有误差，即 ，只考虑由 计算到 这一步的误差，此时由(3.1.10)有局部截断误差(3.1.10)实际上是将精
10、确解 代入(3.1.9)产生的公式误差，利用 Taylor 展开式可得到 .例如对 Euler 法(3.1.2)有，故它表明 Euler 法(3.1.2)的局部截断误差为 ，称为局部截断误差主项.定义 2.2 设 是初值问题(7.1.1)的精确解，若显式单步法(3.1.9)的局部截断误差 ， 是展开式的最大整数，称 为单步法(3.1.9)的阶，含 的项称为局部截断误差主项.根据定义，Euler 法(3.1.2)中的 =1 故此方法为一阶方法.对隐式单步法(3.1.8)也可类似求其局部截断误差和阶，如对后退 Euler 法(3.1.5)有局部截断误差故此方法的局部截断误差主项为 ，也是一阶方法.
11、对梯形法(3.1.6)同样有它的局部误差主项为 ，方法是二阶的.3.1.3 改进 Euler 法上述三种简单的单步法中，梯形法(3.1.6)为二阶方法，且局部截断误差最小，但方法是隐式的，计算要用迭代法。为避免迭代，可先用 Euler 法计算出的近似 ，将(3.1.6)改为(3.1.11)称为改进 Euler 法。改进的 Euler 法是二阶显式方法。即(3.1.12)右端已不含 .可以证明 ， =2,故方法仍为二阶的，与梯形法一样，但用(3.1.11)计算 不用迭代.例 3.2 用改进 Euler 法求例 3.1 的初值问题并与 Euler 法和梯形法比较误差的大小.解 将改进 Euler
12、法用于例 3.1 的计算公式当 n=0 时， .其余结果见表 3-2.表 3-2 改进 Euler 法及三种方法的误差比较从表 3-2 中看到改进 Euler 法的误差数量级与梯形法大致相同，而比 Euler 法小得多，它优于 Euler 法。附：改进的 Euler 法程序function y = DEModifEuler(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1v1 = Funval(f,varvec,x(i-1) y
13、(i-1);t = y(i-1) + h*v1;v2 = Funval(f,varvec,x(i) t);y(i) = y(i-1)+h*(v1+v2)/2;endformat short;讲解： 求初值问题（1.1.1）的数值解就是在假定初值问题解存在唯一的前提下在给定区间 上的一组离散点 上求解析解的一组近似 为此先要建立求数值解 的计算公式，通常称为差分公式，简单的单步法就是由 计算下一步 ，构造差分公式有三种方法，一是用均差（即差商）近似，二是用等价的积分方程（3.1.4）用数值积分方法，三是用函数的 Taylor 展开，其中 Taylor 展开最有普遍性，可以得到任何数值解的计算公式
14、及其局部截断误差。计算公式是微分方程 的一种近似，局部截断误差的概念就是刻划这种逼迫的好坏。当 为微分方程的解，即 ，而用，定义局部截断误差，它表示用精确解代入计算公式（3.1.9）产生的公式误差为 越大表明公式逼近微分方程的精度越高，因此 就定义为公式的阶，通常 的公式才能用于计算初值问题（1.1.1）的数值解。利用 Taylor 展开时，只要将 的表达式在处展开成 Taylor 公式就可得到不同公式的局部截断误差。如 3.1.2 所给出的Euler 法。后退 Euler 法和梯形法，它们只需用一元函数的 Taylor 展开，与后面第 4章的多步法完全一致，而通常单步法（3.1.9）的一般情
15、况则需要用二元函数的Taylor 展开，才能得到公式的具体形式和局部截断误差。例如对改进 Euler 法，其局部截断误差由（3.1.12）可得 要求出它的结果就要用到二元函数 的 Taylor 展开，将在 3.3 节再作介绍。3.2 Runge-Kutta 方法3.2.1 显式 Runge-Kutta 法的一般形式 上节已给出与初值问题(1.1.1)等价的积分形式(3.2.1)只要对右端积分用不同的数值求积公式近似就可得到不同的求解初值问题(1.1.1)的数值方法，若用显式单步法(3.2.2)当 ，即数值求积用左矩形公式，它就是 Euler 法(3.1.2)，方法只有一阶，若取(3.2.3)就
16、是改进 Euler 法，这时数值求积公式是梯形公式的一种近似，计算时要用二个右端函数 f 的值，但方法是二阶的。若要得到更高阶的公式，则求积分时必须用更多的 f 值，根据数值积分公式，可将(3.2.1)右端积分表示为注意，右端 f 中 还不能直接得到，需要像改进 Euler 法(3.1.11)一样，用前面已算得的 f 值表示为(3.2.3)，一般情况可将(3.2.2)的 表示为(3.2.4)其中 这里 均为待定常数，公式(3.2.2)，(3.2.4)称为 r 级的显式 Runge-Kutta 法，简称 R-K 方法。它每步计算 r 个 f 值(即 )，而 ki 由前面(i-1)个已算出的 表示
17、，故公式是显式的。例如当 r=2 时，公式可表示为(3.2.5)其中 .改进 Euler 法(3.1.11)就是一个二级显式 R-K 方法。参数 取不同的值，可得到不同公式。3.2.2 二、三级显式 R-K 方法对 r=2 的显式 R-K 方法(3.2.5)，要求选择参数 ，使公式的阶 p尽量高，由局部截断误差定义 (3.2.6)令 ，对(3.2.6)式在 处按 Taylor 公式展开，由于将上述结果代入(3.2.6)得要使公式(3.2.5)具有的阶 p=2,即 ，必须(3.2.7)即 由此三式求 的解不唯一.因 r=2，故 ，于是有解(3.2.8)它表明使(3.2.5)具有二阶的方法很多，只
18、要 都可得到二阶 R-K 方法.若取 ，则 ，则得改进 Euler 法(3.1.11)，若取 ，则得 ，此时(3.2.5)为(3.2.9)其中 称为中点公式.后退 Euler 法(3.1.11)及中点公式(3.2.9)是两个常用的二级 R-K 方法，注意二级 R-K 方法只能达到二阶，而不可能达到三阶.因为 r=2 只有 4 个参数，要达到 p=3 则在(3.2.6)的展开式中要增加 3 项，即增加三个方程.加上(3.2.7)的三个方程求 4 个待定参数是无解的.当然 r=2,p=2 的 R-K 方法(3.2.5)当 取其他数时，也可得到其他公式，但系数较复杂，一般不再给出.对 r=3 的情形
19、，要计算三个 k 值，即其中 将 按二元函数在 处按 Taylor 公式展开，然后代入局部截断误差表达式，可得可得三阶方法，其系数应满足方程(3.2.10)这是 8 个未知数 6 个方程的方程组，解也是不唯一的，通常 .一种常见的三级三阶 R-K 方法是下面的 Kutta 三阶方法：(3.2.11)3.2.3 四阶 R-K 方法及步长的自动选择利用二元函数 Taylor 展开式可以确定(3.2.4)中 r=4,p=4 的 R-K 方法，经典的四阶 R-K 方法是：(3.2.12)它的局部截断误差 ，故 p=4,这是最常用的四阶 R-K 方法，数学库中都有用此方法求解初值问题的软件.这种方法的优
20、点是精度较高，缺点是每步要算 4 个右端函数值，计算量较大.例 3.3 用经典四阶 R-K 方法解例 7.1 的初值问题 ，仍取 h=0.1，计算到 ，并与改进 Euler 法、梯形法在 处比较其误差大小.解 用四阶 R-K 方法公式(3.2.12)，此处，于是当 n=0 时于是 ，按公式(3.2.12)可算出此方法误差： 改进 Euler 法误差： 梯形法误差： 可见四阶 R-K 方法的精度比二阶方法高得多.用四阶 R-K 方法求解初值问题(1.1.1)精度较高，但要从理论上给出误差的估计式则比较困难.那么应如何判断计算结果的精度以及如何选择合适的步长 h?通常是通过不同步长在计算机上的计算
21、结果近似估计.设 在处的值 ,当 时， 的近似为 ，于是由四阶 R-K方法有若以 为步长，计算两步到 ，则有于是得即 或(3.2.13)它给出了误差的近似估计.如果 ( 为给定精度)，则认为以 为步长的计算结果 满足精度要求，若 ，则还可放大步长.因此(3.2.13)提供了自动选择步长的方法.讲解：求初值问题（1.1.1）的单步法主要是指 Runge-Kutta 法，本节主要讨论显式 RK 方法，建立具体的计算公式使用的是 Taylor 展开，形如（3.2.4）的显式 RK 方法，当 r1 时就是 Euler 法，因此只要讨论 的计算公式，在 r确定后如何推导公式都是一样的，只是 r 越大计算
22、越复杂，为了掌握了解公式来源，只要以 r2 为例推导计算公式即可。因此本节重点就是用 Taylor 展开求出 r2 的显式 R-K 方法的计算公式，由于方法的局部截断误差为（3.2.6），的右端有 的项，要对它做 Taylor 展开，就要用到二元函数 的 Taylor 展开，按照二元函数 Taylor 级数（3.2.14）将它用到（3.2.6）的 的展开式中，即可得到按 升幂整理出的结果，对r2 的公式只能得到 2 阶的公式，即 ，于是 2 级 R-K 方法（3.2.5）的系数 必须满足（3.2.7）给出的方程，它的解由（3.2.8）给出，只要 ，求出的公式都是 r=2 的 2 阶 R-K 方
23、法。而常用的就是 得到的改进 Euler 法（3.1.11）和 得到的中点公式（3.2.9）。3.3 单步法的收敛性与绝对稳定性3.3.1 单步法的收敛性定义 3.1 设 y(x)是初值问题(1.1.1)的精确解， 是单步法(3.2.2)在处产生的近似解，若则称方法(3.2.2)产生的数值解 收敛于 .实际上，定义中 是一固定点，当 h0 时 n，n 不是固定的.因显然方法收敛，则在固定点 处的整体误差，当 p1 时 .下面定理给出方法(3.2.2)收敛的条件.定理 3.1 设初值问题(7.1.1)的单步法(3.2.2)是 p 阶方法(p1)，且函数 对 y 满足 Lipschitz 条件，即
24、存在常数 L0,使对 ,均有 则方法(3.2.2)收敛，且 .定理证明略.可见3 。3.3.2 绝对稳定性用单步法(3.2.2)求数值解 ，由于原始数据及计算过程舍入误差影响，实际得到的不是 而是 ，其中 是误差，再计算下一步得到以 Euler 法为例，若令 ,则(3.3.1)如果 ，则从 计算到 误差不增长，它是稳定的.但如果条件不满足就不稳定.例 3.4 y=-100y,y(0)=1,精确解为 ，用 Euler 法求解得若取 h=0.025，则 ,当 ，而，显然计算是不稳定的.如果用后退 Euler 法(3.1.5)解此例，仍取 h=0.025,则,即显然当 ，计算是稳定的.由此看到稳定性
25、与方法有关，也与 有关，在此例中 .在研究方法的稳定性时，通常不必对一般的 f(x,y)进行讨论，而只针对模型方程(3.3.2)这里 可能为复数.规定 是因为 时微分方程(3.3.2)本身是不稳定的，而讨论数值方法(3.2.2)的稳定性，必须在微分方程本身稳定的前提下进行.另一方面，对初值问题(1.1.1)，若将 f(x,y)在 处线性展开，可得于是方程(1.1.1)可近似表示为它表明用模型方程(3.3.2)是合理的，至于模型方程(3.3.2)中所以用复数 是因为初值问题(7.1.1)如果是方程组，即 ，则 是(mm)阶矩阵，其特征值可能是复数.当然对单个方程， 就是实数，此时只要规定0 即可
26、.用单步法(3.2.2)解模型方程(3.3.2)可得到(3.3.3)其中 依赖所选方法，如用 Euler 法则(3.3.4)此时由(3.3.1)看到误差方程也为 ，与(3.3.4)是一样的.因此对一般单步法(3.2.2)误差方程也与(3.3.3)一致.下面再考虑二阶 R-K 方法有对四阶 R-K 方法，可得定义 3.2 将单步法(3.2.2)用于解模型方程(3.3.2)，若得到(3.3.3)中的则称方法是绝对稳定的.在复平面上复变量 满足 的区域，称为方法(3.2.2)的绝对稳定域，它与实轴的交点称为绝对稳定区间.例如对 Euler 法， 在复平面 上是以(-1，0)为圆心，以 1 为半径的单
27、位圆域内部，当 为实数时，则得绝对稳定区间为 ，因 0，故有 .在例 7.4 中 时方法稳定，而例中取h=0.025 故不稳定.对后退 Euler 法(3.1.5)，因 0，故 ，其绝对稳定域是以(1，0)为圆心的单位圆外部，绝对稳定区间为 ，即对任何 h0 方法都是绝对稳定的.二阶 R-K 方法的绝对稳定区间为 .三阶 R-K 方法的绝对稳定区间为 .四阶 R-K 方法的绝对稳定区间为 .例 3.5 用经典四阶 R-K 方法计算初值问题步长取 h=0.1 及 0.2,给出计算误差并分析其稳定性.解 本题直接按 R-K 方法(3.2.12)的公式计算.因精确解为 ，其计算误差 如表所示.从计算
28、结果看到，h=0.2 时误差很大，这是由于在 =-20,h=0.2 时 h=-4，而四阶 R-K 方法的绝对稳定区间为-2.785,0,故 h=0.2 时计算不稳定，误差很大.而 h=0.1 时 =-2，其值在绝对稳定区间-2.785,0内，计算稳定，故结果是可靠的.讲解：由于微分方程初值问题数值解公式求出的解 是一个逐次递推的过程，因此原始数据误差及计算过程舍入误差对解的影响就是数值方法绝对稳定性研究的问题，如果由 计算 误差不增长，方法就是绝对稳定的。为使问题得到简化通常就是将方法用于解模型方程（3.3.2），对于单步法得到的差分方程为 ，由于模型方程的 ，代入 Euler 法，得 ，对二阶 R-K 方法，例如，用改进 Euler 法于是 对三阶 R-K 方法有对四阶 R-K 方法有 只要 方法，就是绝对稳定的，这时 的值当 n 增大式是减少的，故计算稳定。这时舍入误差影响可忽略不计，而当 ，则 增大，方法不稳定，计算结果是不可靠的。因此用显式单步法必须使 ，也就是步长选择要满足这一要求。对于隐式的梯形公式将模型方程 ，即 代入得于是 注意 ，于是有，对 成立。这就表明对任意步长 h，梯形法都是绝对稳定的。
