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

1.5　不需要导数的根求解

除了重根，牛顿方法比二分法和FPI方法的收敛速度更快.它达到了这种更快的速度是因为使用了更多的信息，尤其是通过函数导数得到的函数切线方向的信息.在某些情况下，可能难以计算导数.
在这种情况下，割线方法是牛顿方法的一个非常好的替代.它使用近似值割线替代了切线，并且收敛速度差不多快.割线方法的变体使用抛物线替换了直线，抛物线可能具有垂直轴（Muller方法）或者水平轴（逆二次插值）.本节最后描述Brent方法，这是一种结合了迭代和括号方法的优良特征的混合方法.

1.5.1　割线方法及其变体

割线方法和牛顿方法近似，但是使用差商替换了导数.从几何上来看，切线被通过前面两次估计点的直线替换.“割线”与x轴的交点被看做是新的估计.
在当前估计xi处导数的近似可以写为差商f(xi)-f(xi-1)xi-xi-1使用这种近似直接替换牛顿方法中的f′（xi），就得到割线方法.

和不动点迭代以及牛顿方法不同，割线方法需要两个初始估计.
如果假设割线方法收敛到函数f的根r，且f′(r)≠0，近似误差关系ei+1≈f″(r)2f′(r)eiei-1成立并且ei+1≈f″（r)2f′(r)α-1eαi其中α=(1+5)/2≈1.62（见习题6）.割线方法以超线性的速度收敛到一个单根，意味着它在线性和二次收敛方法之间.61

割线方法有三种推广形式，它们也很重要.试位方法（或称为Regula Falsi）和二分法相似，但是其中的中点被类似割线方法的近似所替换.给定区间［a，b］，该区间包含根（假设f（a）f（b）<0），使用割线方法定义下一个点为c=a-f(a)(a-b)f(a)-f(b)=bf(a)-af(b)f(a)-f(b)但是和割线方法不同，由于点（a，f（a））和（b，f（b））在x轴的两侧，62新的点保证在区间［a，b］内.根据f（a）f（c）<0或者f（c）f（b）<0，分别选择新的区间［a，c］或者［c，b］，新的区间仍然可以括住根.

试位方法开始表现得比二分法和割线方法都要好，具有二者最好的性质.但是，二分法在每一步中可以确保消除1/2的不确定性，试位方法却没有能力做出这样的保证，而且在一些例子中可能收敛很慢.
例1.17　初始区间是［-1，1］，使用试位方法找到如下函数的根r=0：f(x)=x3-2x2+32x给定x0=-1，x1=1作为初始区间，我们计算新的点x2=x1f(x0)-x0f(x1)f(x0)-f(x1)=1(-9/2)-(-1)1/2-9/2-1/2=45由于f（-1）f（4/5）<0，新的括住区间为［x0，x2］=［-1，0.8］.这样就结束了第一步.注意到第一步中解的不确定性减小的因子远小于1/2.如图1.12b所示，以后连续进行的多步趋近到x=0的速度很慢.

Muller方法是割线方法在不同方向的推广.该方法不是计算经过先前两个点的直线和x轴的交点，而是使用三个前面生成的点x0、x1、x2，画出通过它们的63抛物线y=p（x），并计算抛物线和x轴的交点.一般来讲，抛物线会生成0个或者2个交点.如果有两个交点，接近上一步中的x2的点会被选作x3.通过简单的二次公式计算，就可以确定两种可能.如果抛物线和x轴不相交，就会出现复数解.能够处理复数代数的软件就可以计算对应的解.我们不会继续探索这个问题，尽管在这个方向的相关文献中有很多资源.
逆二次插值（IQI）是割线方法到抛物线的一种相近的泛化方法.但是使用形如x=p（y）的抛物线，而不是Muller方法所使用的y=p（x）.我们的问题可以立刻求解:这个抛物线和x轴只有一个交点，所以从上一步中的三个估计xi、xi+1和xi+2寻找xi+3，这个过程中没有混淆.
经过三点（a，A），（b，B），（c，C）的二阶多项式x=P(y)为P(y)=a(y-B)(y-C)(A-B)(A-C)+b(y-A)(y-C)(B-A)(B-C)+c(y-A)(y-B)(C-A)(C-B)(1.35)这是一个拉格朗日插值的例子，在第3章中将进行讲述.现在我们只要知道P（A）=a，P（B）=b，P（C）=c就足够了.用y=0代入，得到抛物线和x轴的交点.经过重新组合与替代，我们得到P(0)=c-r(r-q)(c-b)+(1-r)s(c-a)(q-1)(r-1)(s-1)(1.36)其中q=f（a）/f（b），r=f（c）/f（b），s=f（c）/f（a）.
对于IQI，通过设置a=xi,b=xi+1,c=xi+2,A=f(xi),B=f(xi+1),C=f(xi+2)下一步的估计xi+3=P(0)为xi+3=xi+2-r(r-q)(xi+2-xi+1)+(1-r)s(xi+2-xi)(q-1)(r-1)(s-1)(1.37)其中q=f(xi)/f(xi+1),r=f(xi+2)/f(xi+1),s=f(xi+2)/f(xi).给定三个初始估计，IQI方法对式（1.37）进图1.13　比较Muller方法以及逆向二次迭代方法.前者由插值抛物线y=p(x)确定;后者由插值抛物线x=p(y)确定行迭代，使用最新的估计xi+3替换最旧的估计xi.IQI的另一种方式使用最新估计替换最近的三个估计中的后向误差最大的那个.
图1.13比较了Muller方法与逆二次插值的几何.两种方法比割线方法的收敛速度都快，这是由于更高阶的插值计算.在第3章，我们会更详细地研究插值.割线方法的概念和推广，以及二分法是Brent方法的重要组成部分，该方法将在下节中进行描述.

1.5.2　Brent方法

Brent方法［Brent，1973］是一种混合方法，该方法使用前面介绍的迭代技术，推出一个新的方法，并保留前面方法中有用的性质.如果能把二分法的保证收敛的性质以及更加复杂方法的快速收敛性质结合起来就好了.这个方法最早是由Dekker和Van Wijngaarden在20世纪60年代提出来的.
该方法用于连续函数f，区间的边界是a和b，同时f（a）f（b）<0.Brent方法记录当前点xi，该点具有最优的后向误差，同时有包含根的区间［ai，bi］.简单来讲，64尝试使用逆二次方法，并在下述情况下，使用结果来替代xi、ai、bi中的一个：（1）后向误差得到改进；（2）包含根的区间至少减小一半.否则，尝试使用割线方法以实现相同的目的.如果割线方法也失败了，则使用二分法，保证至少减少一半的不确定性.
MATLAB的命令fzero可以实现Brent方法的一种形式，如果用户没有给定一个初始空间，该命令会提供一个预处理步骤，找到一个好的括住区间.终止条件混合了前向和后向误差的标准.当从xi到一个新点xi+1的变化小于2εmachmax(1,xi)或者后向误差f（xi）变成机器零，算法就会终止.
如果用户提供一个初始括住区间，则不进行预处理步骤.下面使用命令输入函数f(x)=x3+x-1以及初始括住区间［0，1］，并且让MATLAB在每步迭代中显示部分结果:

寻找f（x）在x=1附近的根，在这个过程中，首先确定括住区间，然后使用Brent方法65.
1.5节习题
1.对下列方程使用两步割线方法，初始估计为x0=1，x1=2.
(a)x3=2x+2(b)ex+x=7(c)ex+sinx=4
2.对习题1中的方程使用两步试位方法，初始区间为［1，2］.
3.对习题1中的方程使用两步逆向二次插值方法.初始估计x0=1，x1=2，x2=0，保留三次最近的迭代进行更新.
4.一个商业渔夫想把网撒在温度是10℃的水中.从温度计上的线看到在水深9m时，温度是8℃，水深5m时温度是15℃.使用割线方法找出温度10℃时对应的最合适的水深.
5.通过y=0代入式（1.35），推导公式（1.36）.
6.如果割线方法收敛到r，f′(r)≠0,f″(r)≠0，则近似误差关系ei+1≈f″(r)/(2f′(r))eiei-1成立.证明如果极限limi→∞ei+1/eαi对于α>0成立，并且非0，则α=(1+5)/2，ei+1≈(f″(r)/2f′(r))α-1eαi.
7.考虑下面4种计算21/4，即2开4次方的方法.
（a）根据收敛速度从最快到最慢进行排序.并给出进行这样排序的原因.
（A） 使用二分法f(x)=x4-2
（B） 使用割线方法f(x)=x4-2
（C） 使用不动点迭代g(x)=x2+1x3
（D） 使用不动点迭代g(x)=x3+13x3
（b） 是否还有比以上方法收敛更快的方法？
1.5节编程问题
1.使用割线方法找出习题1中的一个（唯一）根.
2.使用试位方法找出习题1中每个方程的根.
3.使用逆向二次插值找出习题1中每个方程的根.
4.令f(x)=54x6+45x5-102x4-69x3+35x2+16x-4，在区间［-2，2］中画出函数，使用割线方法找出该区间中所有的5个根.哪些根是线性收敛？哪些根是超线性收敛？
5.在习题1.1.6中，曾被问过二分法在区间［-2，1］，函数f（x）=1/x的解是多少.现在比较fzero得到的结果和二分法计算的结果.
6.（a） 如果使用fzero确定函数f(x)=x2在1附近的根会有什么结果（不使用一个括住区间）？解释结果.（b） 该问题对于f（x）=1+cosx在-1附近会得到什么结果？66
事实验证1　Stewart平台运动学
一个Stewart平台包含6个可变长度的支杆，或者棱柱关节，用于支撑负载.棱柱关节通过气动或者水动的方式改变支杆的长度运转.作为一个6自由度的机器人，Stewart平台可以放在任何地方，并趋向它能力范围之内的三维空间.
为了简化问题，本项目考虑Stewart平台的二维版本.该控制器由三个支杆控制的平面上的一个三角形平台构成，如图1.14所示.内部三角形表示　图1.14　平面Stewart平台草图.前向动力学问题中使用p1、p2、p3确定未知的x、y、θ平面Stewart平台，其对应的维数由三个长度L1、L2、L3定义.令γ表示边L1所对的角度.平台位置由三个长度p1、p2、p3控制，这对应三个支杆变化的长度.
给定三个支杆的长度，找到平台位置，被称为控制器的前向，或者方向动力学问题.从名字可以看出来，在给定p1、p2、p3后，计算（x，y）和θ.由于有三个自由度，很自然地使用三个数字表示平台位置.对于运动规划，尽可能快地求解问题非常重要，一般需要实时.但是不幸的是，对于平面Stewart平台的前向运动学问题没有解析解.
当前最好的方法包括如图1.14所示，把几何形式降到一个方程，并使用本章介绍的求解器进行解答.你的工作是计算这个方程的导数并写代码进行求解.
在图1.14中使用简单的三角方法包含下面的三个方程:p21=x2+y2
p22=(x+A2)2+(y+B2)2
p23=(x+A3)2+(y+B3)2(1.38)在方程中，A2=L3cosθ-x1
B2=L3sinθ
A3=L2cos(θ+γ)-x2=L2［cosθcosγ-sinθsinγ］-x2
B3=L2sin(θ+γ)-y2=L2［cosθsinγ+sinθcosγ］-y2注意到式（1.38）求解了平面Stewart平台的逆向运动学问题，通过给定x、y、θ，找出p1、p2、p3.你的目的是求解前向问题，即给定p1、p2、p3，解x、y、θ.67
把式（1.38）中最后两个方程展开，并把第一个方程代入，得到p22=x2+y2+2A2x+2B2y+A22+B22=p21+2A2x+2B2y+A22+B22
p23=x2+y2+2A3x+2B3y+A23+B23=p21+2A3x+2B3y+A23+B23其中只要D=2(A2B3-B2A3)≠0，x和y可以求解为x=N1D=B3(p22-p21-A22-B22)-B2(p23-p21-A23-B23)2(A2B3-B2A3)
y=N2D=-A3(p22-p21-A22-B22)+A2(p23-p21-A23-B23)2(A2B3-B2A3)(1.39)用这些式子替换（1.38）第一个方程中的x和y，乘上D2得到如下的一个方程f=N21+N22-p21D2=0(1.40)其中只有一个未知数θ.（回忆一下p1，p2，p3，L1，L2，L3，γ，x1，x2，y2都是已知的.）如果找到f（θ）的根，对应的x和y可以从（1.39）解出来.
注意到f（θ）是关于sinθ和cosθ的多项式，所以，对于任何给定的根θ，还有其他的根θ+2πk，这些根对于平台也是等价的.由于这个原因，我们把角度θ限制在区间［-π，π］.可以看到f（θ）在这个区间上至多有6个根.
建议活动
1.写出f（θ）的MATLAB函数文件.参数L1，L2，L3，γ，x1，x2，y2为固定的常数，当给定姿态时，可知支杆的长度p1，p2，p3.如果你不了解MATLAB函数文件，可以参考附录B.5.下面为第一行和最后一行:

为了检查代码，由图1.15,将参数设为L1=2,L2=L3=2,γ=π/2,p1=p2=p3=5.然后，用θ=-π/4或者θ=π/4进行替代，这分别对应图1.15a和图1.15b，使得f（θ）=0.
2.画出区间［-π，π］上的f（θ）.你可以使用附录B.5中的@符号为你的函数文件在画图命令中分配一个函数句柄，你可能还需要在算术运算前面加上“.”使得操作矢量化，该计算在附录B.2中描述.作为对你的工作的检查，应该在±π/4位置上有根.
3.重新生成图1.15.MATLAB命令如下

将会画出一个红色三角形，三角形的顶点为（u1,v1）,（u2，v2）,（u3,v3），在支杆锚点（anchor）（0,0）,（0,x1）,（x2,y2）上画上小圆，并且画出支杆.

4.求解平面Stewart平台的前向动力学系统，其中，x1=5,(x2,y2)=(0,6),L1=L3=3,L2=32,γ=π/4,p1=p2=5,p3=3.68首先画出f（θ）.使用方程求解技术找出4个位置，并画出这些位置.通过验证p1、p2、p3是图中支杆的长度，检查你的结果.
5.将支杆长度改变为p2=7并重新求解问题.对于这些参数，有6个姿态.
6.找出支杆长度p2，其他参数和步骤4设置相同，使得其中只有两个姿态.
7.计算p2的区间，其他参数和步骤4设置相同，使得其中分别有0、2、4和6个姿态.
8.推导并确定一个方程,来表示三维、6自由度的Stewart平台前向动力系统.写出MATLAB程序，并验证其用于求解前向动力系统.参考Merlet［2000］的论著可以找出对棱柱机械臂和平台的一个很好的介绍.
软件与进一步阅读
有多种算法可以用于确定非线性方程的解.包括收敛速度慢，但是保证收敛的算法，例如二分法，以及具有更快的收敛速度，但是不能保证收敛的算法，包括牛顿方法及其变种.根据是否需要方程的导数信息，方程求解方法可以分为两类.二分法、割线方法以及逆向二次插值，在计算过程中仅仅需要黑盒子提供给定输入的函数值，而牛顿方法需要方程的导数.Brent方法是个混合方法，结合慢速和快速收敛方法的好的特性，而且不需要导数计算.由于这个原因，Brent方法是最常用的一般方程求解方法，在大量综合的软件包中都包含这个方法.
MATLAB的fzero命令实现了Brent方法，仅仅需要初始区间和一个初始值作为输入.IMSL的ZBREN程序、NAG的程序c05adc，以及netlib FORTRAN的程序fzero.f都依赖这个基本的方法.MATLAB的roots命令可以找出一个多项式的所有根，69
　～
70该命令使用完全不同的方法.该方法计算伴随矩阵的所有特征值，构造得到的矩阵的特征值和多项式的所有根相同.
其他常常提到的算法基于Muller方法和拉格郎日方法，这些方法在正确的条件下，具有三次收敛速度.更多细节可以参考Traub［1964］、Ostrowski［1966］和Householder［1970］关于方程求解的教科书.