本节书摘来自华章出版社《数值分析（原书第2版）》一 书中的第2章，第2.7节，作者：（美）Timothy Sauer，更多章节内容可以访问云栖社区“华章计算机”公众号查看。

2.7　非线性方程组

第1章中包含求解一个未知变量的方程，该方程通常是非线性方程.在本章中，我们已经研究了方程组的求解，但是要求方程组是线性的.结合非线性和“多于一个方程”的因素，大大提高了求解问题的难度.本节中我们将描述牛顿方法及其变体，并用于求解非线性方程组.130

2.7.1　多元牛顿方法

单变量的牛顿方法xk+1=xk-f(xk)f′(xk)提供了多元牛顿方法的主要轮廓.两种方法都是根据泰勒展开的线性近似推导得到，例如，令f1(u,v,w)=0
f2(u,v,w)=0(2.49)
f3(u,v,w)=0是三个非线性方程，具有三个未知变量u、v、w.定义向量值函数F(u，v，w)=(f1，f2，f3)，将（2.49）中的问题表示为F（x）=0，其中x=(u，v，w).
单变量情况下函数导数f′现在对应的是雅可比矩阵，定义为DF(x)=f1uf1vf1w
f2uf2vf2w
f3uf3vf3w在点x0附近的向量值函数的泰勒展开是F(x)=F(x0)+DF(x0)·(x-x0)+O(x-x0)2例如，在x0=（0，0）附近函数F(u,v)=(eu+v,sinu)的线性展开是F(x)=1
0+e0e0
cos00u
v+O(x2)=1
0+u+v
u+O(x2)牛顿方法基于线性近似，忽略O（x2）项.正如在一维情况下，令x=r是根，令x0是当前估计.则0=F(r)≈F(x0)+DF(x0)·(r-x0)或者-DF(x0)-1F(x0)≈r-x0(2.50)因而，通过求解（2.50）得到r，可以对根进行更好的近似.
多变量牛顿方法x0=初始向量
xk+1=xk-(DF(xk))-1F(xk),k=0,1,2,…
131　　由于计算矩阵逆的代价很大，我们使用技巧避免该计算.在每一步中，不是按上面描述的过程，而是令xk+1=xk-s，其中s是DF(xk)s=F(xk)的解.现在，仅仅需要在每步中使用高斯消去（n3/3次的乘法），而不是计算矩阵的逆（相对于高斯消去大约3倍的计算次数）.因而，对于多元牛顿方法的迭代步骤如下DF(xk)s=-F(xk)
xk+1=xk+s(2.51)例2.32　使用牛顿方法，初始估计为（1，2），找出如下方程组的解v-u3=0
u2+v2-1=0图2.5显示了一组点，在这些点上f1(u,v)=v-u3和f2(u,v)=u2+v2-1都为0，以及它们的两个交点，交点对应方程组的解.图2.5　对例2.32使用牛顿方法.对应的两个根是圆上的两个点.牛顿方法过程中的点在(0.8260,0.5636)处近似收敛到解雅可比矩阵是DF(u,v)=-3u21
2u2v在第一步中，使用初始点x0=（1，2），我们必须求解矩阵方程（2.51）:-31
24s1
s2=-1
4解是s=（0，-1），因而第一步迭代过程x1=x0+s=(1,1).第二步要求求解-31
22s1
s2=-0
1解是s=(-1/8，-3/8)，x2=x1+s=(7/8，5/8).图2.5显示了两个迭代.后面的步骤得到下表:132
步数uv步数uv01.000000000000002.0000000000000040.826040108170650.5636197735028411.000000000000001.0000000000000050.826031357732410.5636241621316320.875000000000000.6250000000000060.826031357654190.5636241621612630.829036348267120.5643491124260470.826031357654190.56362416216126
小数点后精确数位个数翻倍是熟悉的二次收敛的特征，这在输出的序列中可以明显看到.如图2.5所示，方程的对称性表明，如果（u，v）是一个解，则（-u，-v）也是解.第二个解也可以通过使用附近的初始估计，并应用牛顿方法迭代得到.
例2.33　使用牛顿方法求解系统f1(u,v)=6u3+uv-3v3-4=0
f2(u,v)=u2-18uv2+16v3+1=0注意到（u，v）=（1，1）是一个解.此外还有另外两组解.雅可比矩阵如下DF(u,v)=18u2+vu-9v2
2u-18v2-36uv+48v2正如一维例子中所示，使用牛顿方法能够找到哪个解取决于初始估计.使用初始点(u0，v0)＝（2，2），对前面的公式迭代得到下表：
步数uv步数uv02.000000000000002.0000000000000041.000033678665061.0000224377201011.372580645161291.3403225806451651.000000001119571.0000000005789421.078386812004431.0538012326498461.000000000000001.0000000000000031.005349688965201.0026926187153971.000000000000001.00000000000000其他的初始向量得到另外的两个解，它们近似为（0.865939，0.462168）以及（0.886809，-0.294007）.参见编程问题2.
如果可以计算雅可比矩阵，牛顿方法是一个好的选择.否则，***的替代方法就是Broyden方法，我们将在下一节中讲述.
2.7.2　Broyden方法
牛顿方法求解单变量的方程需要知道导数.在第1章中对于该方法的讨论得到割线方法，割线方法用在没有导数或者导数难以计算的时候.
既然我们有了用于求解非线性方程组F（x）=0的牛顿方法，133我们也面临同样的问题:如果没有雅可比矩阵DF怎么办？尽管没有直接的方法将用于方程的牛顿方法推广到割线方法.但是Broyden［1965］提出一种方法，这通常被认为是次优的方法.
假设Ai是第i步可以得到的雅可比矩阵的最优近似，并被用于生成xi+1=xi-A-1iF(xi)(2.52)为了在下一步中从Ai更新到Ai+1，我们注意到雅可比矩阵DF的导数，满足Ai+1δi+1=Δi+1(2.53)其中δi+1=xi+1-xi，Δi+1=F(xi+1)-F(xi).另一方面，我们没有任何关于δi+1的正交补向量新的信息.因而，我们要求对于所有满足δTi+1w=0的w.Ai+1w=Aiw(2.54)同时满足（2.53）和（2.54）的矩阵如下：Ai+1=Ai+(Δi+1-Aiδi)δTi+1δTi+1δi+1(2.55)Broyden方法使用牛顿方法步骤（2.52）来推进当前的近似估计，使用（2.55）更新雅可比矩阵的近似.总结来说，算法从初始估计x0和初始近似雅可比A0开始.如果没有更好的选择，初始的雅可比矩阵可以使用单位阵.
Broyden方法

其中δi+1=xi+1-xi,Δi+1=F（xi+1）-F（xi）
注意到牛顿类型的步骤通过求解Aiδi+1=F(xi)进行，正如牛顿方法所做的那样.和牛顿方法一样，Broyden方法不保证收敛到根.
Broyden方法的第二种方式避免相对代价较大的矩阵求解步骤Aiδi+1=F(xi).由于我们最多在迭代过程中近似DF导数，我们也可以同样近似DF的逆，这正是牛顿步骤所需要的.
基于Bi=A-1i，重做Broyden导数.我们希望得到δi+1=Bi+1Δi+1(2.56)其中δi+1=xi+1-xi，Δi+1=F(xi+1)-F(xi)，并对于任何满足δTi+1w=0的w，也满足Ai+1w=Aiw，或者Bi+1Aiw=w(2.57)134满足（2.56）和（2.57）的矩阵是Bi+1=Bi+(δi+1-BiΔi+1)δTi+1BiδTi+1BiΔi+1(2.58)不需要矩阵求解的迭代的新版本如下：xi+1=xi-BiF(xi)(2.59)得到的算法称为Broyden方法Ⅱ.
Broyden方法Ⅱ

首先，需要初始向量x0和初始估计B0.如果难以计算导数，可以使用B0=I.Broyden方法Ⅱ的一个可以察觉的缺点是，在一些应用中需要估计雅可比矩阵，但是这个矩阵并不容易得到.矩阵Bi是对雅可比矩阵逆的估计.而Broyden方法Ⅰ正相反，一直记录了Ai用来估计雅可比.正是由于这个原因，在一些圈子里，Broyden方法Ⅰ和方法Ⅱ被分别称为“好的Broyden”和“坏的Broyden”.
两种版本的Broyden方法都超线性收敛（到单根），比牛顿方法的二次收敛要慢一点.如果有雅可比的计算公式，通常使用DF（x0）的逆作为初始矩阵B0，这样通常可以加速收敛.
Broyde方法Ⅱ的MATLAB代码如下:

例如，例2.32的方程组的解可以通过定义如下函数得到

对于Broyden方法Ⅱ的调用如下:
135
Broyden方法的两种实现，在没有雅可比矩阵的时候都非常有用.这种情况的典型例子在事实验证7的管道扭曲问题中得到展示.
2.7节习题
1.找出如下函数的雅可比矩阵.
(a) F(u,v)=(u3,uv3)
(b) F(u,v)=(sinuv,euv)
(c) F(u,v)=(u2+v2-1,(u-1)2+v2-1)
(d) F(u,v,w)=(u2+v-w2,sinuvw,uvw4)
2.使用泰勒展开找出L（x）到F（x）在x0附近的线性近似.
(a) F(u,v)=(1+eu+2v,sin(u+v)),x0=(0,0)
(b) F(u,v)=(u+eu-v,2u+v),x0=(1,1)
3.在uv平面画出两条曲线，通过简单的代数找出所有的解.
(a) u2+v2=1
(u-1)2+v2=1　　(b) u2+4v2=4
4u2+v2=4　　(c) u2-4v2=4
(u-1)2+v2=4
4.对习题3的方程组使用两步牛顿方法，起始点为（1，1）.
5.对习题3的方程组使用两步Broyden方法Ⅰ，起始点为（1，1），使用A0=I.
6.对习题3的方程组使用两步Broyden方法Ⅱ，起始点为（1，1），使用B0=I.
7.证明（2.55）满足（2.53）和（2.54）.
8.证明（2.58）满足（2.56）和（2.57）.2.7节编程问题
1.使用合适的初始点，利用牛顿方法找出所有解.用习题3进行检查确保结果正确.
(a) u2+v2=1
(u-1)2+v2=1(b) u2+4v2=4
4u2+v2=4(c) u2-4v2=4
(u-1)2+v2=4
2.使用牛顿方法找出例2.31的三个解.
3.使用牛顿方法找出如下方程组的两个解u3-v3+u=0，u2+v2=1.
4.使用牛顿方法找出如下三个方程的两个解.2u2-4u+v2+3w2+6w+2=0
u2+v2-2v+2w2-5=0
3u2-12u+v2+3w2+8=05.使用多元牛顿方法找出三维空间中三个给定球的两个交点.（a）每个球的半径都是1，球心在（1，1，0），（1，0，1），以及（0，1，1）.（答案是（1，1，1）和（1/3，1/3，1/3））（b）每个球半径是5，球心在（1，-2，0），（-2，2，-1），以及（4，-2，3）.136
6.尽管三个球在三维空间中的一般交点是两个点，它也可能是唯一点.使用多元牛顿方法找出球心在（1，0，1）半径为8、球心在（0，2，2）半径为2、球心在（0，3，3）半径为2的三个球的唯一交点.这个迭代仍然二次收敛吗？请解释.
7.对习题3的方程组，使用Broyden方法Ⅰ，初始估计x0=（1，1），A0=I.尽可能精确地报告解，以及所需要的步数.
8.对习题3的方程组，使用Broyden方法Ⅱ，初始估计为（1，1），B0=I.报告尽可能精确的解，以及所需要的步数.
9.使用Broyden方法Ⅰ找出编程问题5中的两个交点.
10.使用Broyden方法Ⅰ找出编程问题6中的交点.你所能观测到的收敛率是多少？
11.使用Broyden方法Ⅱ找出编程问题5中的两个交点.
12.使用Broyden方法Ⅱ找出编程问题6中的交点.你所能观测到的收敛率是多少？
软件与进一步阅读
在数值线性代数中有很多优秀的教科书，包含Stewart［1973］和全面的参考材料Golub与Van Loan［1996］.Demmel［1997］和Trefethen与Bau［1997］为两本包含数值线性代数的现代技术的教科书.涵盖迭代技术的书包括Axelsson［1994］、Hackbush［1994］、Kelley［1995］、Saad［1996］、Traub［1964］、Varga［2000］、Young［1971］，以及Dennis与Schnabel［1983］.
LAPACK是一个全面的、公开的软件包，涵盖了高质量的矩阵代数计算程序，包含求解Ax=b、矩阵分解和条件数的方法.代码经过精心设计可以在现在计算机框架下进行移植，包含共享存储向量和并行处理器.这些内容可以参看Anderson等人［1990］的论著.
LAPACK的可移植性依赖于算法在书写方式上最大化使用了基本线性代数子程序（Basic Linear Algebra Subprograms，BLAS），BLAS包括了基本的矩阵/向量计算，这些计算可以通过调整优化在特定的机器和框架上的性能.BLAS可以粗略地分为三个部分:1级，要求O（n）次的操作，例如点积的操作;2级，诸如矩阵/向量乘法的操作，O（n2）次的操作;3级，包含全矩阵/矩阵乘法，具有O（n3）的复杂度.
LAPACK上的DGESV函数，以双精度用PA=LU分解方法求解Ax=b，还有其他用于稀疏和带状矩阵的方法.参见www.netlib.org/lapack可以看到更多细节.LAPACK函数的实现也构成了MATLAB，以及其他IMSL与NAG的矩阵代数计算的基础.137