第二章非线性方程求根线性方程是方程式中仅包含未知量的一次方项和常数项的方程，除此之外的方程都是非线性方程(nonlinear equation). 例如，大家熟知的“一元二次方程”就是一个非线性方程. 多元线性方程组的求解是数值计算领域的一个重要问题，在后续几章将专门讨论. 本章介绍求解非线性方程的数值方法，主要针对实数域，重点是单个非线性方程的求根问题.2.1引言2.1.1非线性方程的解记要求解的单变量非线性方程为f(x)=0(2.1)其中函数f: ℝ→ℝ. 一般而言，非线性方程的解的存在性和个数是很难确定的，它可能无解，也可能有一个或多个解.例2.1 (非线性方程的解)：分析下列非线性方程的解是否存在和解的个数.(1) e x+1=0. 此方程无解.(2) e−x−x=0. 此方程有一个解.(3) x2−4sinx=0. 此方程有两个解.(4) x3−6x2+5x=0. 此方程有三个解.(5) cosx=0. 此方程有无穷多个解.在实际问题中，往往要求的是自变量在一定范围内的解，比如限定x∈[a,b]. 函数f一般为连续函数，则可记为f(x)∈C[a,b]，C[a,b]表示区间[a,b]上所有连续实函数的集合. 假设在区间[a, b]上方程(2.1)的根为x∗，也称x∗为函数f(x)的零点. 方程的根可能不唯一，而且同一个根x∗也可能是方程(2.1)的多重根.定义2.1：对光滑函数f，若f(x∗)=f′(x∗)=⋯=f(m−1)(x∗)=0，但f(m)(x∗)≠0，则称x∗为方程(2.1)的m重根. 当m=1时，即f(x∗)=0,f′(x∗)≠0时，称x∗为单根.对于多项式函数f(x)，若x∗为m重根，则f(x)可因式分解为f(x)=(x−x∗)m g(x)其中g(x)也是多项式函数，且g(x∗)≠0. 很容易验证，f(x∗)=f′(x∗)=⋯=f(m−1)(x∗)=0，但f(m)(x∗)≠0，即多项式方程重根的概念与定义2.1是一致的. 对一般的函数f，x∗是方程(2.1)的重根的几何含义是，函数曲线在x∗处的斜率为0，且在该点处与x轴相交.非线性方程的一个特例是n次多项式方程（n≥2），根据代数基本定理可知，n次方程在复数域上有n个根（m重根计为m个根）. 当n=1, 2时，方程的求解方法是大家熟知的. 当n=3, 4时，虽然也有求根公式，但已经很复杂，在实际计算时并不一定适用. 当n≥5时，不存在一般的求根公式，只能借助数值求解方法来求根.2.1.2问题的敏感性根据问题敏感性的定义，这里需要考虑输入数据的扰动对方程的根有多大影响. 要分析敏感性首先应假设问题中的数据如何扰动，一种易于分析的情况是将非线性方程写成：f(x)=y的形式，然后讨论y在0值附近的扰动造成的问题敏感性. 此时，求根问题变成了函数求值问题：y =f (x )的反问题. 若函数值f (x )对输入参数x 不敏感（x 在解x ∗附近变化），则求根问题将很敏感；反之，若函数值对参数值敏感的话，求根则不敏感. 这两种情况如图2-1所示.下面分析y 发生扰动Δy 引起的方程的根的扰动Δx . 由于当x =x ∗时，y =0，我们使用绝对（而不是相对）条件数：cond =|Δx |≈1|| 条件数的大小反映方程求根问题(2.1)的敏感程度，若|f ′(x ∗)|很小，则问题很敏感，是一个病态问题；反之，若|f ′(x ∗)|很大，则问题不敏感.一种特殊情况是f ′(x ∗)=0，即x ∗为重根，此时求根问题很敏感，原问题的微小扰动将造成很大的解误差，甚至改变解的存在性和唯一性（如图2-2，问题的扰动可能使解不存在）.对于敏感的非线性方程求根问题，f (x )≈0并不意味着x 很接近x ∗，在后面讨论迭代解法的判停准则时应注意这一点.2.2二分法数值求解非线性方程通常是一个迭代的过程，迭代开始之前要先有个初始的近似解，然后随着迭代步数的增多，近似解越来越接近准确解，当达到一定要求时即停止计算过程. 本节先介绍一种最基本的方法——二分法(interval bisection method).2.2.1 方法原理先介绍有根区间的概念，有根区间就是包含至少一个根的区间，它限定了根存在的范围. 如果能计算出一个非常小的有根区间，那么区间的中点就是一个很好的近似解. 下面的定理给出了有根区间的充分条件.定理2.1：若f (x )∈C[a,b]，且f (a )f (b )<0，则区间(a,b )内至少有一实根.这里省略定理证明过程，只给出图2-3作为一个解释.定理2.1给出了一种获得有根区间的方法，即通过看f (a )、f (b )两个值是否符号相反来判断(a, b)是否为有根区图2-1 方程求根问题的敏感性：(a)不敏感；(b)敏感.(a) (b)图2-2 f . 图2-3 若f (a )f (b )<0，则在区间(a,b )内至少有一实根.间. 在实际操作时，可在一个较大的范围内取多个点计算f(x)函数值，从而得到一个或多个有根区间. 另外应注意，根据定理2.1得到的有根区间内不一定只有一个根，这从图2-3也可以看出.二分法的思想很简单，就是每次将有根区间一分为二，得到长度逐次减半的区间序列{(a k,b k)}，则区间中点x k=(a k+b k)2⁄就是第k步迭代的近似解. 具体算法如下：算法2.1：二分法输入：a, b, 函数f(x) ; 输出：x.While(b−a)> εdox:= a+(b−a)/2;If sign(f(x))=sign(f(a))thena:= x;Elseb:= x;EndEndx:= a+(b−a)/2.在算法2.1中，sign()表示取符号的函数，而二分迭代结束的条件为有根区间(a, b)的长度小于某个阈值ε. 应注意，浮点运算中几乎不可能出现等于0的情况，所以sign()函数的结果只是正号、或负号.假设二分法得到的有根区间序列为{(a k,b k),k=0,1,⋯}，若取解x k=(a k+b k)2⁄，则误差|x k−x∗|<(b k−a k)2⁄=(b0−a0)2k+1⁄,k=0,1,2,⋯.(2.2) 根据公式(2.2)和对解的准确度的要求，也可以事先估算出二分迭代执行的次数，以及相应的计算量. 这里每步迭代的计算量主要是计算一次函数f(x k).例2.2(二分法)：求方程f(x)=x4−x−2=0在区间[1.0, 1.5]上的一个实根，要求准确到小数点后第2位（四舍五入后）.[解] 首先验证(1.0, 1.5)是否是一个有根区间，易知f(1.0)<0,f(1.5)>0. 所以将(1.0, 1.5)作为二分法的初始区间. 利用(2.2)式我们可以估计，若(b−a)2k+1⁄≤0.5×10−2(2.3) 则|x k−x∗|<0.5×10−2，即结果准确到了小数点后第2位. 代入a=1.0, b=1.5,求解(2.3)得，k≥log20.50.5×10−2−1=5.6取最小的整数值k=6. 只需二分6次，可得到满足精度要求的解. 计算过程中的数据和结果列于表2-1. 从中看出，准确到小数点后两位的解为x=1.356(准确解为1.353210).表2-1 采用二分法求解例2.2的过程和结果2.2.2 算法稳定性和结果准确度算法的稳定性考察的是计算过程中的误差对结果的影响. 对于二分法来说，主要的计算步骤是计算函数值，一般采用双精度浮点数计算函数值的误差很小，而其他计算是少量的加减法，因此不至于对有根区间，以及最终结果的准确度造成多大影响. 另外，在计算过程中解的误差限逐次减半，这也说明二分法是稳定的.在实际的浮点算术体系中，二分法运行结果的准确度不可能随迭代过程一直提高. 首先看一个例子.例 2.3(二分法准确度的极限)：编写程序用二分法来求解方程f(x)=x2−2=0，初始区间为[1, 2].下面是用MATLAB语言编写的程序：M =2; a = 1; b= 2; k = 0;while b-a > eps //Matlab中的eps为2倍的机器精度，即2− 2 x = a + (b-a)/2;if x^2 > Mb = x //输出belsea = x //输出aendk = k+1;end这个程序执行了52步就结束了，输出结果如下：b = 1.50000000000000a = 1.25000000000000a = 1.37500000000000b = 1.43750000000000………………a = 1.41421356237309a = 1.41421356237309b = 1.41421356237310b = 1.41421356237310为了看得更清楚，输入MATLAB命令”format hex”，使输出按16进制格式显示，再运行一遍上述程序，最后的四个输出结果为：a = 3ff6a09e667f3bc8a = 3ff6a09e667f3bccb = 3ff6a09e667f3bceb = 3ff6a09e667f3bcd可以看出，最终区间(a, b)的两个端点已经是两个相邻的浮点数，即使让二分过程继续执行下去，区间仍然不会改变（由于a和b平均值的计算结果就是其中的某一个）. 也就是说迭代再多的次数，结果的准确度也无法提高了.上述例子说明了二分法结果准确度的极限情况. 一般地，二分法迭代过程中，有根区间缩小的极限情况是使它的端点a k, b k为两个相邻的机器浮点数. 此时，b k−a k=2⌊log2|x∗|⌋∙2εmacℎ其中εmacℎ为机器精度，⌊∙⌋为下取整符号，而2⌊log2|x∗|⌋为x∗的二进制表示中指数的那部分. 在这种极端情况下，解的误差限就是区间长度，即|e(x k)|=|x k−x∗|≤2⌊log2|x∗|⌋∙2εmacℎ .(2.4) 在IEEE双精度浮点数系统下，εmacℎ=2− 3，则|e(x k)|≤2⌊log2|x∗|⌋∙2εmacℎ≤|x∗|∙2− 2≈|x∗|∙2.22×10−16 .根据(2.4)式，也可得到相对误差的上限：|e r(x k)|=|x k−x∗||x∗|≤2εmacℎ .(2.5)这个相对误差限正好是计算机中用浮点数表示实数的误差限（定理1.5）的两倍.公式(2.4)给出了用二分法求解时绝对误差限可能达到的最小值，它说明在准确解较大的情况下，在执行二分法时无论迭代多少次都无法保证解的误差非常小. 综合上述讨论以及公式(2.2)，得到如下定理.定理2.2：在实际的浮点算术体系下采用二分法解方程f(x)=0，设初始有根区间为(a, b)，则：(1) 结果的误差限最小可达到2⌊log2|x∗|⌋∙2εmacℎ，其中x∗为准确解，相对误差限最小可达2εmacℎ.(2) 若误差阈值ε大于或等于2⌊log2|x∗|⌋∙2εmacℎ，需执行的迭代步数为k=⌈log2(b−a)−1⌉ .(2.6)定理2.2的结论(2)的证明留给读者思考.最后，对二分法说明几点：●二分法是求单变量方程f(x)=0的实根的一种可靠算法，一定能收敛.●二分法解的误差不一定随迭代次数增加一直减小，在实际的有限精度算术体系中，误差限存在最小值.●二分法的缺点是有时不易确定合适的初始有根区间（含两个初始值）、收敛较慢、且无法求解偶数重的根. 因此，实际应用中常将二分法与其他方法结合起来.2.3不动点迭代法二分法的计算效率不够高，在本章后续部分将介绍几种应用广泛、收敛较快的迭代法. 本节介绍不动点迭代法及其收敛性理论，为后续其他方法的讨论建立基础.2.3.1基本原理通过某种等价变换，可将非线性方程(2.1)改写为：x=φ(x)(2.7)其中φ(x)为连续函数. 给定初始值x0后，可构造迭代计算公式：x k+1=φ(x k),(k=0,1,⋯)(2.8) 从而得到近似解序列{x k}. 由于方程(2.1)和(2.7)的等价关系，很容易证明若序列{x k}收敛，其极限必为原方程(2.1)的解x∗. 由于解x∗满足x∗=φ(x∗)，称它为函数φ(x)的不动点(fixed point)，此方法为求解非线性方程(2.1)的不动点迭代法(fixed-point iterative method).不动点迭代法的求解过程如图2-4所示，而算法描述图2-4 采用不动点迭代法，近似∗在下面给出.算法2.2：基于函数φ(x)的不动点迭代法输入：x 0, 函数f (x ), φ(x) ; 输出：x .k:= 0 ;While |f (x k )|> ε1 或|x k −x k−1|> ε2 dox k+1:=φ(x k ) ;k:= k +1 ;Endx:= x k .其中ε1和ε2为用于判断迭代是否应停止的两个阈值. 关于迭代的判停准则，我们在2.4.3节详细讨论.例 2.4(不动点迭代法)：求f (x )=x 4−x −2=0在x 0=1.5附近的根. 以不同的方式得到方程的等价形式，研究相应的不动点迭代法的收敛情况.[解] 将原方程改为等价的(A), (B)两种形式，得到下述两种不动点迭代法.方法(A)：将方程改写为x =x 4−2. 得到的迭代法计算公式为：x 0=1.5, x k+1=x k 4−2,(k =0,1,⋯). 计算出结果如下：x 1=1.54−2=3.0625,x 2=2.3753−1=85.9639,… …从上述计算结果看，序列{x k }有趋于无穷大的趋势，迭代法不收敛，无法求出近似解. 方法(B)：将方程改写为x =√x +24. 得到的迭代法计算公式为：x 0=1.5, x k+1=√x k +24,(k =0,1,⋯). 计算出结果如下：x 1=√1.5+24=1.3678,x 2=√1.3678+24=1.3547,x 3=√1.3547+24=1.3534,x 4=√1.3534+24=1.3532,x =√1.3532+24=1.3532,从上述计算结果看，x 4和x 前5位有效数字均为1.3532，可认为迭代过程是收敛的，要求的根为1.3532.通过例2.4可以看出，用不同的方式改造原方程，可得到多种不动点迭代法计算过程，其收敛性质也是不同的. 因此，判断一个不动点迭代法是否收敛至关重要.2.3.2全局收敛的充分条件下面的定理给出一个函数存在唯一不动点的充分条件.定理2.3：设φ(x)∈C[a,b]，若满足如下两个条件：(1)对任意x ∈[a,b]，有a ≤φ(x)≤b ，(2)存在正常数L ∈(0,1)，使对任意x 1,x 2∈[a,b]，|φ(x 1)−φ(x 2)|≤L |x 1−x 2|则φ(x)在[a,b]上存在不动点，且不动点是唯一的.在证明定理之前，先理解一下定理中两个条件的含义. 首先，采用不动点迭代法的计算公式为x k+1=φ(x k ), k =0, 1, 2, …，因此要使后续迭代步的计算合法，必须要求φ(x k )的值在函数的定义域内，(1)的条件保证了这一点. 其次，(2)中新加的条件表明，φ(x)曲线上任两点连线斜率的绝对值不超过L ，当两点非常靠近时它就是导数，因此φ(x)曲线上任意点的切线斜率的绝对值都小于1，这说明φ(x)曲线变化很平缓，在曲线上任意点处的斜率都比y =x 和y=−x两条直线小. 这个条件也称为L<1的李普希兹（Lipschitz）条件，L为李普希兹系数. [证明] 先证明不动点的存在性，分两种情况：1). 若φ(a)=a，或φ(b)=b，则a或b为不动点;2). 若φ(a)≠a且φ(b)≠b，则φ(a)>a,φ(b)<b. 令f(x)=φ(x)−x，则f(x)为连续函数，且f(a)>0,f(b)<0.根据连续函数性质，必有x∗∈(a,b)，使f(x∗)=0，即φ(x∗)=x∗，x∗为不动点.再证明唯一性，采用反证法. 假设有两个不同的不动点x1∗,x2∗∈[a,b]，它们满足φ(x1∗)=x1∗,φ(x2∗)=x2∗,x1∗≠x2∗ .根据(2)中的条件推出，|x1∗−x2∗|=|φ(x1∗)−φ(x2∗)|≤L|x1∗−x2∗|<|x1∗−x2∗| .产生矛盾！所以假设x1∗≠x2∗不成立，不动点是唯一的.应说明的是，上述证明不动点的存在性只使用了条件(1). 事实上，通过画函数曲线图的方式也可以形象地说明不动点的存在性，这一点留给感兴趣的读者思考.下面的定理给出不动点迭代法收敛的充分条件.定理2.4：设φ(x)∈C[a,b]满足定理2.3的两个条件，则对于任意初值x0∈[a,b]，由不动点迭代法得到的序列{x k}收敛到φ(x)的不动点x∗，并有误差估计：|x k−x∗|≤L k1−L|x1−x0| .[证明] 首先注意到定理条件保证了不动点唯一存在，而且条件(1)保证了不动点迭代法可执行下去，从而得到序列{x k}. 下面证明序列{x k}收敛，其思路是考虑误差序列，证明其极限为0.|x k−x∗|=|φ(x k−1)−φ(x∗)|≤L|x k−1−x∗|≤⋯≤L k|x0−x∗|由于L为小于1的正常数，则lim k→∞L k|x0−x∗|=0,⟹limk→∞|x k−x∗|=0, ⟹limk→∞x k=x∗.这证明了不动点迭代法是收敛的. 剩下|x k−x∗|≤L k1−L|x1−x0|的证明，留给读者思考.定理2.4为判断不动点迭代法的收敛性提供了依据，这种收敛不依赖于初值x0的选取，因此称为全局收敛性. 为了方便应用，也可以将定理2.3, 2.4中的第2个条件替换为：对任意x∈[a,b]，有|φ′(x)|≤L<1，其中L为常数，得到便于使用的定理2.5.定理2.5：设φ(x)∈C[a,b]，且满足如下两个条件：(1)对任意x∈[a,b]，有a≤φ(x)≤b；(2)存在正常数L<1，使对任意x∈[a,b]，有|φ′(x)|≤L<1.则对于任意初值x0∈[a,b]，由不动点迭代法得到的序列{x k}收敛到φ(x)的不动点x∗，并有误差估计：|x k−x∗|≤L k1−L|x1−x0| .此定理可看成定理2.4的推论，其证明留给读者思考.例2.5(不动点迭代法的收敛性)：对于求f(x)=x4−x−2=0在x0=1.5附近的根的问题，使用定理2.5考察例2.4中两种方法的全局收敛性.[解]在区间[1, 2]上考察如下两种不动点迭代法的收敛性：方法(A)：x k+1=x k4−2, (k=0,1,⋯).方法(B)：x k+1=√x k+24, (k=0,1,⋯)很容易看出，方法(B)符合定理中的条件(1)，而φ′(x)=14(x+2)−3/4，也符合条件(2)，因此方法(B)具有全局收敛性. 而对于方法(A)，它不符合定理中的条件(1)，因此无法根据定理2.5说明其具有全局收敛性.关于全局收敛性再说明两点：●定理2.4, 2.5给出的都是不动点迭代法全局收敛的“充分条件”，也就是说，对一些满足条件的方法可以证明其具有全局收敛性，但根据它们并不说明某个方法不具有全局收敛性.●全局收敛性要求初始值x0为定义域内任意值时不动点迭代法都收敛，这常常是很难达到的要求.2.3.3局部收敛性不同于全局收敛性，下面给出重要的局部收敛性的概念.定义2.2：设函数φ(x)存在不动点x∗，若存在x∗的某个邻域D: [x∗−δ,x∗+δ]，对于任意初值x0∈D，迭代法x k+1=φ(x k)产生的解序列{x k}收敛到x∗，则称迭代法局部收敛.这个定义中的邻域是以x∗为中心点的一个对称区间，局部收敛性的定义要求的是存在这样一个邻域，而不关心它的大小. 下面的定理给出迭代法局部收敛的充分条件.定理2.6：设x∗为函数φ(x)的不动点，若φ′(x)在x∗的某个邻域上连续，且|φ′(x∗)|<1，则不动点迭代法x k+1=φ(x k)局部收敛.[证明] 因为φ′(x)在x=x∗附近连续且|φ′(x∗)|<1，则存在x∗的某个邻域D，使得对于任意. 显然L<1, 的xϵD, |φ′(x)|≤L, 其中L是某个介于|φ′(x∗)|和1之间的数, 例如L=|φ′(x∗)|+12即满足定理2.5的条件(2). 另外，对∀x∈D,φ(x)−x∗=φ(x)−φ(x∗)=φ′(ξ)(x−x∗), ξ∈D, ⟹|φ(x)−x∗|≤L|x−x∗|<|x−x∗|, 即φ(x)∈D，满足定理2.5的条件(1).因此，根据定理2.5，此迭代法对区间D内的任意初值都收敛，根据定义2.2知，此迭代法局部收敛.对比定理2.6和定理2.4、2.5可以看出，定理2.6的条件较为宽松，它只需要考察函数φ(x)在x∗这一点上是否满足要求. 因此，不动点迭代法较容易具有局部收敛性，对局部收敛的判断也相对简单.最后说明一点，定理2.5说明李普希兹系数L越小迭代收敛的速度越快，而定理2.6的证明过程说明了L与|φ′(x∗)|的关系. 因此，若|φ′(x∗)|越小，迭代收敛的速度就越快.2.3.4稳定性与收敛阶与二分法类似，不动点迭代法的每步计算都可以通过判停准则（包括考察f(x k)是否接近0）来评估解的准确度，因此解的误差容易被及时发现和纠正. 只要迭代过程是收敛的，误差将随迭代步的增加逐渐趋于零，而不会像某些算法的舍入误差会随迭代过程逐渐累积. 因此收敛的不动点迭代法总是稳定的. 在本章后续算法的讨论中，我们将不再关心稳定性，而将重点放在收敛性的讨论上.对于收敛的迭代法，其收敛速度的快慢也很重要，它关系到达到特定的准确度需要多少步迭代，也就是需要多少计算量. 下面先看一个例子，然后给出收敛阶的概念用于衡量迭代收敛的速度.例2.6(迭代收敛速度)：假设有(1)~(3)三个迭代求解过程，其迭代解的误差|e(x k)|=|x k−x∗|随迭代步变化情况分别为：(1) 10-2, 10-3, 10-4, 10-5, ……。
