本篇详细介绍用于求解一元函数零点的牛顿法及其各种形式的变形,并给出收敛性质。
问题:求函数 f:I⊆R→R f : I ⊆ R → R 的零点 α α ( f(α)=0 f ( α ) = 0 )。 这个问题可以使用如下形式的单点迭代法(one-point iteration),得到迭代函数 g(x) g ( x ) 的不动点(fixed point )即是函数的零点。
xn+1=g(xn),n=0,1,..., x n + 1 = g ( x n ) , n = 0 , 1 , . . . , 其中 x0 x 0 是初始点,如果没有更好的选择(尽量靠近零点 α α )通常设为0。 最常用的单点迭代法是 经典牛顿法: xn+1=xn−f(xn)f′(xn),n=0,1,..., x n + 1 = x n − f ( x n ) f ′ ( x n ) , n = 0 , 1 , . . . , 对于单零点函数( simple zeros),牛顿法二次收敛( converges quadratically);对于多零点函数( multiple zeros),牛顿法线性收敛( converges linearly )( It converges quadratically to simple zeros and linearly to multiple zeros.). 下面介绍收敛性的概念。定义 1: 如果序列 xn|n≥0 x n | n ≥ 0 以下列方式收敛到 α α ,
limxn→αxn+1−α(xn−α)p=C lim x n → α x n + 1 − α ( x n − α ) p = C 其中 C≠0,p≥1 C ≠ 0 , p ≥ 1 ,则称 p p 是序列xnxn收敛的 阶, C C 是渐进误差常数(asymptotic error constant)。当p=1,p=2,p=3p=1,p=2,p=3时,分别称线性收敛( converge linearly)、二次收敛( converge quadratically)、三次收敛( converge cubically)。 p p 称为产生序列xnxn的方法的 收敛阶( order of convergence)。令 en=xn−α. e n = x n − α . 则, en+1=Cepn+O(ep+1n) e n + 1 = C e n p + O ( e n p + 1 ) 称为方法的误差方程( error equation)。 理论上证明一个方法的阶一般都是推导成这种形式。
定义 2: 令 α α 是函数 f f 的零点,xn−1,xn,xn 1xn−1,xn,xn 1为三次连续的迭代生成值,且比较靠近 α α ,则 p p 可用如下公式近似求解: ρ≈ln|(xn 1−α)/(xn−α)|ln|(xn−α)/(xn−1−α)|ρ≈ln|(xn 1−α)/(xn−α)|ln|(xn−α)/(xn−1−α)| 简单证明:反推法。
令 α α 为充分可微函数 f f 的零点。由牛顿-莱布尼茨公式: f(x)=f(xn) ∫xxnf′(t)dt(2)f(x)=f(xn) ∫xnxf′(t)dt(2) 假设我们用常数 f′(xn) f ′ ( x n ) 代替 f′ f ′ 在区间 [xn,x] [ x n , x ] 上的值,则上式右端的定积分值为: (x−xn)f′(xn) ( x − x n ) f ′ ( x n ) 。令, x=α x = α ,则, f(x)=f(α)=0 f ( x ) = f ( α ) = 0 ,则:
0≈f(xn)+(α−xn)f′(xn), 0 ≈ f ( x n ) + ( α − x n ) f ′ ( x n ) , 从而可知对 α α 最新的逼近 xn+1 x n + 1 可有下式决定: xn+1=xn−f(xn)f′(xn) x n + 1 = x n − f ( x n ) f ′ ( x n ) 这就是经典牛顿法( Classical Newton’ Method简称 CN)。如果把(2)式右侧的定积分用梯形公式(trapzoidal rule),则,
xn+1=xn−2f(xn)f′(xn)+f′(xn+1) x n + 1 = x n − 2 f ( x n ) f ′ ( x n ) + f ′ ( x n + 1 ) 显然这是个隐式方程,可先用经典牛顿法求解 xn+1 x n + 1 的值,再带入到上式等式右侧的 f′(xn+1) f ′ ( x n + 1 ) ,即, xn+1=xn−2f(xn)f′(xn)+f′(zn+1),wherezn+1=xn−f(xn)f′(xn)(3) x n + 1 = x n − 2 f ( x n ) f ′ ( x n ) + f ′ ( z n + 1 ) , w h e r e z n + 1 = x n − f ( x n ) f ′ ( x n ) ( 3 ) 上式可改写成:xn+1=xn−f(xn)(f′(xn)+f′(zn+1))/2,n=0,1,...(4) x n + 1 = x n − f ( x n ) ( f ′ ( x n ) + f ′ ( z n + 1 ) ) / 2 , n = 0 , 1 , . . . ( 4 ) 即,使用 f′(xn) f ′ ( x n ) 和 f′(zn+1) f ′ ( z n + 1 ) 的算术平均值( arithmetic mean)而不是 f′(xn). f ′ ( x n ) . 因此这种方法被称为 算术平均值牛顿法(AN)。
在(4)式中,如果使用调和平均值(harmonic mean)取代算术平均值,则,
xn+1=xn−f(xn)(f′(xn)+f′(zn+1))2f′(xn)f′(zn+1) x n + 1 = x n − f ( x n ) ( f ′ ( x n ) + f ′ ( z n + 1 ) ) 2 f ′ ( x n ) f ′ ( z n + 1 )在式(2)中,如果使用中点积分(midpoint integration)代替梯形公式(trapezoidal rule),则,
0≈f(xn)+(α−xn)f′(xn+α2) 0 ≈ f ( x n ) + ( α − x n ) f ′ ( x n + α 2 ) 即, xn+1=xn−f(xn)f′((xn+xn+1)/2). x n + 1 = x n − f ( x n ) f ′ ( ( x n + x n + 1 ) / 2 ) . 显然这也是个隐式方法,可用经典牛顿法近似求解 xn+1 x n + 1 并带入到上式右侧,即, xn+1=xn−f(xn)f′((xn+zn+1)/2),wherezn+1=xn−f(xn)f′(xn),(6) x n + 1 = x n − f ( x n ) f ′ ( ( x n + z n + 1 ) / 2 ) , w h e r e z n + 1 = x n − f ( x n ) f ′ ( x n ) , ( 6 )参考文献[2]提出,如果用Simpson’s formula近似(2)式的定积分,即,
∫xxnf′(t)dt≈x−xn6{f′(x)+4f′(x+xn2)+f′(xn)} ∫ x n x f ′ ( t ) d t ≈ x − x n 6 { f ′ ( x ) + 4 f ′ ( x + x n 2 ) + f ′ ( x n ) } 整理后得, xn+1=xn−6f(xn)f′(xn+1+4f′(xn+1+xn2)+f′(xn) x n + 1 = x n − 6 f ( x n ) f ′ ( x n + 1 + 4 f ′ ( x n + 1 + x n 2 ) + f ′ ( x n ) 同样的,这也是个隐式方程,先用经典牛顿法求得 xn+1 x n + 1 ,再带入上式右侧。最终可得如下公式: xn+1=xn−6f(xn)f′(zn+1+4f′(zn+1+xn2)+f′(xn) x n + 1 = x n − 6 f ( x n ) f ′ ( z n + 1 + 4 f ′ ( z n + 1 + x n 2 ) + f ′ ( x n ) where, zn+1=xn−f(xn)f′(xn). z n + 1 = x n − f ( x n ) f ′ ( x n ) .文献[3]使用反调和平均数(contra harmonic mean)来取代(4)式中的算术平均值,即,
xn+1=xn−f(xn)(f′(xn)+f′(zn+1))f′2(xn)+f′2(zn+1),n=0,1,..., x n + 1 = x n − f ( x n ) ( f ′ ( x n ) + f ′ ( z n + 1 ) ) f ′ 2 ( x n ) + f ′ 2 ( z n + 1 ) , n = 0 , 1 , . . . , where, zn+1=xn−f(xn)f′(xn) z n + 1 = x n − f ( x n ) f ′ ( x n ) 称其为反调和平均值牛顿法( contra harmonic Newton’s method), CHN。 论文声称此方法三阶收敛,适用于求解非线性方程。下面的实验结果来自参考文献[3]。(我很奇怪,结果并没有证明[3]提出的CHN更好,为什么这篇论文还能发表) 1. 测试函数:
停止条件(stopping criterion)为: |xn+1−xn|<ε | x n + 1 − x n | < ε ,取 ε=10−14. ε = 10 − 14 . 2. 满足停止条件(stopping criterion)需要的迭代次数: 3. 满足停止条件需要的函数求值次数(total number of function evaluations, TNFE)(CN迭代次数*2;其它迭代次数*3):
这里只给出收敛性结论,具体证明可见相关参考文献。
一般情况下,都要求如下条件:令 α∈I α ∈ I 是充分可微函数 f:I⊆R f : I ⊆ R 的单零点(simple zero), I I 是开区间,x0x0充分接近 α α 。
经典牛顿法二阶收敛;算术平均值牛顿法三阶收敛;调和平均值牛顿法(HN)三阶收敛(在参考文献[1][3][4]中被证明性能最好);中点牛顿法(MN)三阶收敛(在参考文献[3]中被证明性能最好);Simpson牛顿法三阶收敛;反调和平均数(CHN)牛顿法三阶收敛。注意:如果初始点 x0 x 0 不足够靠近零点 α α ,则所有牛顿法都有可能不收敛。
本节是对参考文献[4]的总结
正实数 a,b a , b 的 α α 次幂平均数 mα m α 定义为:
mα=⎧⎩aα+bα2⎫⎭1α m α = ⟮ a α + b α 2 ⟯ 1 α 显然, - 当 α=−1 α = − 1 时为调和平均数(Harmonic mean) m−1=2aba+b m − 1 = 2 a b a + b - 当 α=12 α = 1 2 时, m12={a√+b√2}2 m 1 2 = { a + b 2 } 2 - 当 α=1 α = 1 时,为算术平均数(Arithmetic mean), m1=a+b2 m 1 = a + b 2 - 当 α→0 α → 0 时, limα→0mα=ab−−√ lim α → 0 m α = a b ,又称为几何平均值(geometric mean)对于给定的整数 a,b a , b ,其它常见的平均数还有, - N(Heronianmean)=a+ab√+b3 N ( H e r o n i a n m e a n ) = a + a b + b 3 - C(Contra−harmonicmean)=a2+b2a+b C ( C o n t r a − h a r m o n i c m e a n ) = a 2 + b 2 a + b - T(Centroidalmean)=2(a2+ab+b2)3(a+b) T ( C e n t r o i d a l m e a n ) = 2 ( a 2 + a b + b 2 ) 3 ( a + b ) - L(Logarithmicmean)=a−blog(a)−log(b) L ( L o g a r i t h m i c m e a n ) = a − b l o g ( a ) − l o g ( b )
用利用上述这些平均数,就可以得到更多的牛顿法变形。例如,使用 α α 次幂平均数,可得一般形式的牛顿法变形,
xn+1=xn−f(xn)sign{f′(x0)}[(f′(xn))α+(f′(zn+1))α2]1α x n + 1 = x n − f ( x n ) s i g n { f ′ ( x 0 ) } [ ( f ′ ( x n ) ) α + ( f ′ ( z n + 1 ) ) α 2 ] 1 α where, zn+1=xn+f(xn)f′(xn) z n + 1 = x n + f ( x n ) f ′ ( x n ) 论文[4]还证明了这种一般形式的变形是 三阶收敛的。可以构造新函数
f^(x)=ep(x−xn)f(x) f ^ ( x ) = e p ( x − x n ) f ( x ) 显然 f^(x) f ^ ( x ) 与 f(x) f ( x ) 有相同的零点,但其导数, f′^(x)=pep(x−xn)f(x)+ep(x−xn)f′(x) f ′ ^ ( x ) = p e p ( x − x n ) f ( x ) + e p ( x − x n ) f ′ ( x ) 在 f′(xn)=0 f ′ ( x n ) = 0 时,不等于0。 得到牛顿法的另一种变形: xn+1=xn−f(xn)f′(xn)+pf(xn) x n + 1 = x n − f ( x n ) f ′ ( x n ) + p f ( x n )本文从求解一元实值函数的零点出发(牛顿法的起源),介绍了牛顿法及其变形。但牛顿法的应用范围远远不局限于此,其还可以用在:
多元函数求零点;非线性方程的求解:等价于一元函数或多元函数求零点;非线性方程组的求解:分别对每一个方程求零点,转换成线性方程组;(多元)函数求极值:通过令导数(偏导)等于0,得到方程(组);上述牛顿法的各种应用本质是一样的。其中函数求极值也称最优化问题,是牛顿法应用最多的地方。另外,
其它和牛顿法相关的概念还包括:
阻尼牛顿法;Levenberg-Marquardt method;高斯牛顿法;拟牛顿法:DFP & BFGS;在接下来的博客中,我会一一介绍。
高阶收敛的牛顿法是否必要? 参考文献[5]给出这样一种观点:
Higher-order methods. The local rate of convergence of Newton’s method is fast enough; however, some researchers construct methods with still higher rate of convergence. This problem looks a bit artificial (actually very few iterations of Newton’s method are required to obtain high precision when we achieve its convergence domain), so there is no need to accelerate the method.
翻译:经典牛顿法的收敛速度已经够快了(通常只需要很少的迭代次数),因此没必要再使用高阶收敛的算法。(orz (:з」∠))