非线性方程(组)的解法

本文关注非线性方程(组)的解法，包括一元情况下非线性方程的bisetion, regula falsi, 不动点迭代, Newton法, 多元情况下非线性方程组的Newton法, Broyden法。与之对比的是线性方程组的求解，如LU分解、QR分解等。

Bisection

求解线性方程组 f(x)=0f(x)=0 相当于求解函数 f(x)f(x) 的零点。

利用零点存在性定理来求解的方法就是二分法。根据闭区间套定理，二分法的收敛性很容易证明。

同时，区间的长度每次减半，精度等于最终区间的长度。

Regula Falsi

Regula falsi 也即 false position bisection，是二分法的一种变种，与二分法的区别是采用线性插值零点作为新的二分点。与二分法一样，它也是能保证收敛性的。

这种方法相当于用线性函数来逼近原函数 f(x)f(x)，那么很显然，当 f(x)f(x) 与线性函数越接近，regula falsi的效果越好；反之，它的效果差于二分法。

Fixed Point Iteration

不动点迭代是一种经典的求解非线性方程的方法。在求解线性方程组的问题中，不动点迭代法相当于 Jacobi 迭代法。

不动点迭代的核心思想是将 f(x)=0f(x)=0 转化为 x=g(x)x=g(x), 然后不断地进行迭代

xk+1=g(xk)x_{k+1}=g(x_k)

对于同一个问题，可以找到多个 g(x)g(x), 不动点迭代的效果取决于 g(x)g(x) 的好坏。

Linear Convergence

假设精确解为 x∗x_*，其满足 x∗=g(x∗)x_*=g(x_*). 那么

xk+1−x∗=g(xk)−g(x∗)=g′(ξk)(xk−x∗)≤∣g′(ξk)∣∣xk−x∗∣x_{k+1}-x_*=g(x_k)-g(x_*)=g'(\xi_k)(x_k-x_*)\le|g'(\xi_k)||x_k-x_*|

当 ∣g′(x)∣<1,∀x|g'(x)|< 1, \forall x 时, 具有全局线性收敛性。

假设仅有 ∣g′(x∗)∣<1|g'(x_*)|<1, 且 g′(x)g'(x) 连续，则有局部线性收敛性。

Quadratic Convergence

当 g′(x∗)=0g'(x_*)=0 时，在 x∗x_* 的邻域进行泰勒展开，有

xk+1−x∗=g(xk)−g(x∗)=g′′(ξ)2!(x−x∗)2≤C(xk−x∗)2，∃Cx_{k+1}-x_*=g(x_k)-g(x_*)=\frac{g''(\xi)}{2!}(x-x_*)^2\le C(x_k-x_*)^2， \exist C

具有局部二次收敛性

Newton’s Method

Basic Newton’s Method

Newton 法的思想是用切线来近似 f(x)f(x), 即 f(x)≈f(x0)+f′(x0)(x−x0)f(x)\approx f(x_0)+f'(x_0)(x-x_0)

Newton 法的迭代公式为

xk+1=xk−f(xk)f′(xk)x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}

实际上，它是一种精心构造的不动点迭代法。令

g(x)=x−f(x)f′(x)g(x)=x-\frac{f(x)}{f'(x)}

g′(x)=1−(f′(x))2−f(x)f′′(x)(f′(x))2=f(x)f′′(x)(f′(x))2g'(x)=1-\frac{(f'(x))^2-f(x)f''(x)}{(f'(x))^2}=\frac{f(x)f''(x)}{(f'(x))^2}

当 f′(x∗)≠0f'(x_*)\not=0 时，根据前面的结论可以知道此时 Newton 法具有局部二次收敛。

可以证明当 x∗x_* 为重根(i.e. f(x∗)=f′(x∗)=⋯=f(m)(x∗)f(x_*)=f'(x_*)=\cdots=f^{(m)}(x_*)) 时，依然具有局部线性收敛性。

对于重根，可以做以下修正

xk+1=xk−(m+1)f(xk)f′(xk)x_{k+1}=x_k-\frac{(m+1)f(x_k)}{f'(x_k)}

这样也有局部二次收敛性，读者可以自行证明。

Secant Method

当求导的操作很贵时，可以考虑用割线来代替割线

f′(xk)≈f(xk)−f(xk−1)xk−xk−1f'(x_k)\approx \frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}

xk+1=xk−f(xk)(xk−xk−1)f(xk)−f(xk−1)x_{k+1}=x_k-\frac{f(x_k)(x_k-x_{k-1})}{f(x_k)-f(x_{k-1})}

可以证明采用割线法，收敛阶为黄金分割比例 1+52≈1.618\frac{1+\sqrt5}{2}\approx1.618

引入 Newton 插值中 nn 阶差商的记号 f[x1,x2,…,xn]=∑i=1k+1f(xi)∏j≠i(xi−xj)−1f[x_1,x_2,\dots,x_n]=\sum\limits_{i=1}^{k+1}f(x_i)\prod\limits_{j\ne i}(x_i-x_j)^{-1}

这样，已知给定 xk,xk−1x_k,x_{k-1}，加入 x∗x_* 进行插值，可以得到

0=f(x∗)=f(xk)+f[xk,xk−1](x∗−xk)+f[xk,xk−1,x∗](x∗−xk)(x∗−xk−1)0=f(x_*)=f(x_k)+f[x_k,x_{k-1}](x_*-x_k)+f[x_k,x_{k-1},x_*](x_*-x_k)(x_*-x_{k-1})

利用差商重写迭代公式

0=f(xk)+f[xk,xk−1](xk+1−xk)0=f(x_k)+f[x_k,x_{k-1}](x_{k+1}-x_k)

两式相减可以得到

x∗−xk+1=−f[xk,xk−1,x∗]f[xk,xk−1](x∗−xk)(x∗−xk−1)x_*-x_{k+1}=-\frac{f[x_k,x_{k-1},x_*]}{f[x_k,x_{k-1}]}(x_*-x_k)(x_*-x_{k-1})

根据 Roole 中值定理，nn 阶差商具有以下性质

f[x1,…,xn+1]=f(n)(ξ)/n!f[x_1,\dots,x_{n+1}]=f^{(n)}(\xi)/n!

所以当 ff 的一阶和二阶导数均有界时，有

∣x∗−xk+1∣≤M∣x∗−xk∣∣x∗−xk−1∣,∃M|x_*-x_{k+1}|\le M |x_*-x_{k}||x_*-x_{k-1}|,\exist M

这与 Fibonacci 数列 FkF_k 十分相关，可以证明

∣x∗−xk∣≤C∣x∗−x0∣Fk|x_*-x_{k}|\le C |x_*-x_{0}|^{F_k}

根据 FkF_k 的通项公式，可以知道收敛阶为黄金分割比例 1+52≈1.618\frac{1+\sqrt5}{2}\approx1.618

Steffensen’s Method

迭代公式如下

xk+1=xk−f(xk)f(xk+f(xk))−f(xk)f(xk)=xk−(f(xk))2f(xk+f(xk))−f(xk)x_{k+1}=x_k-\frac{f(x_k)}{\frac{f(x_k+f(x_k))-f(x_k)}{f(x_k)}}=x_k-\frac{(f(x_k))^2}{f(x_k+f(x_k))-f(x_k)}

这种方法二次收敛

非线性方程组的求解

以上的都是求解一元非线性方程的方法。以下讨论非线性方程组的解法 (特别指 Rn↦Rn\R^n \mapsto \R^n).

Newton’s Method

与一元的类似，将导数换成了Jacobi矩阵

xk+1=xk−Jk−1f(xk)x_{k+1}=x_k-J_k^{-1}f(x_k)

Good Broyden Method

在 Newton 法中，可以看出一个缺点是每一步迭代都要求解一次线性方程组，时间复杂度是 O(n3)O(n^3)，很不划算。Broyden 方法就是为了解决这一问题。

它的目标是近似地求 JkJ_k，从而使得求逆变得方便。很容易联想到 Sherman-Morrison-Woodbury 公式

(A+uvT)−1=A−1−A−1uvTA−11+vTA−1u(A+uv^T)^{-1}=A^{-1}-\frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}

假如每次更新 JkJ_k 都只加入秩一修正，那么就好办了

JK=Jk−1+[rank 1]J_K=J_{k-1}+[rank\ 1]

因为 JkJ_k 为 Jacobi 矩阵，所以它应该满足 JkΔxk=ΔfkJ_k\Delta x_k=\Delta f_k，这样当迭代足够多次后，JkJ_k 能较好地近似真实的 Jacobi 矩阵。

令 rk=Δfk−Jk−1Δxk=ΔJΔxkr_k=\Delta f_k-J_{k-1}\Delta x_k=\Delta J\Delta x_k

因为 rank(ΔJ)=1rank(\Delta J)=1, ∃yk,ΔJ=rkykT,ykTΔxk=1\exist y_k,\Delta J=r_ky_k^T,y_k^T\Delta x_k=1. 同时，要使得 ∣∣ΔJ∣∣F=∣∣rk∣∣2∣∣yk∣∣2||\Delta J||_F=||r_k||_2||y_k||_2 最小.

由 Cauchy 不等式可知 ∣∣yk∣∣2∣∣Δxk∣∣2≥1||y_k||_2||\Delta x_k||_2 \ge 1, 因此 ∣∣yk∣∣2≥1∣∣Δxk∣∣2||y_k||_2\ge \frac1{||\Delta x_k||_2}, 当且仅当 yk=ΔxΔxkTΔxky_k=\frac{\Delta x}{\Delta x_k^T\Delta x_k}时，等号成立

所以可以知道 JkJ_k 的更新公式

Jk=Jk−1+(Δfk−Jk−1Δxk)ΔxkTΔxkTΔxkJ_k=J_{k-1}+\frac{(\Delta f_k-J_{k-1}\Delta x_k)\Delta x_k^T}{\Delta x_k^T\Delta x_k}

代入 Sherman-Morrison-Woodbury 公式可得

Jk−1=Jk−1−1−Jk−1−1(Δfk−Jk−1Δxk)ΔxkTJk−1−1ΔxkTΔxk+ΔxkTΔxk(Δfk−Jk−1Δxk)=Jk−1−1+(Δxk−Jk−1−1Δfk)ΔxkTJk−1−1ΔfkTΔfk\begin{aligned} J_k^{-1}&=J_{k-1}^{-1}-\frac{J_{k-1}^{-1}(\Delta f_k-J_{k-1}\Delta x_k)\Delta x_k^TJ_{k-1}^{-1}}{\Delta x_k^T\Delta x_k+\Delta x_k^T\Delta x_k(\Delta f_k-J_{k-1}\Delta x_k)}\\ &=J^{-1}_{k-1}+\frac{(\Delta x_k-J_{k-1}^{-1}\Delta f_k)\Delta x_k^TJ_{k-1}^{-1}}{\Delta f_k^T\Delta f_k} \end{aligned}

Bad Broyden Method

所谓 Good 和 Bad，其实并不表示实际效果的好坏，只是两种方法名字的区分。Bad Broyden 方法是直接对 Jk−1J_k^{-1} 做秩一修正，也就是要满足

Jk−1=Jk−1−1+ΔJJ_k^{-1}=J_{k-1}^{-1}+\Delta J

Δxk=Jk−1Δfk\Delta x_k=J_k^{-1}\Delta f_k

相同的处理，最终可以得到更新公式

Jk−1=Jk−1−1+(Δxk−Jk−1Δfk)ΔfkTΔfkTΔfkJ_k^{-1}=J_{k-1}^{-1}+\frac{(\Delta x_k-J_k^{-1}\Delta f_k)\Delta f_k^T}{\Delta f_k^T\Delta f_k}

数值算法

本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议 ，转载请注明出处！