非线性方程（组）的求解

1. 非线性方程的求解

1.1 简介

非线性方程是指含有指数和余弦函数等非线性函数的方程，例如， e x − c o s ( π x ) = 0 e^x-cos(\pi x)=0 ex−cos(πx)=0
与线性方程相比，无论是解的存在性，还是求解的计算公式，非线性方程问题都比线性问题要复杂的多，对于一般线性方程 f ( x ) = 0 f(x)=0 f(x)=0，既无直接法可用，也无一定章程可寻。常用的方法有：

实根的对分法
迭代法
牛顿跌迭代法
弦截法

通常，非线性方程的根不止一个，因此在求解非线性方程时，要给定初值或求解范围。

1.2 实根对分法

对分法又称二分法，设函数 f ( x ) f(x) f(x)在 [ a , b ] [a,b] [a,b]上连续，且 f ( a ) f ( b ) < 0 f(a)f(b)<0 f(a)f(b)<0，则 f ( x ) f(x) f(x)在 [ a , b ] [a,b] [a,b]上至少有一零点，这是微积分中的介值定理，也是使用对分法的前提条件。算法如下：

　注意：对分法算法简单，然而，若 [ a , b ] [a,b] [a,b]上有几个零点时，只能算出其中一个零点；另一方面，即使 [ a , b ] [a,b] [a,b]上有零点，也未必有 f ( a ) f ( b ) < 0 f(a)f(b)<0 f(a)f(b)<0，这就限制了对分法的使用范围；对分法只能计算方程的实根。

1.3 不动点迭代法

对给定的非线性方程 f ( x ) = 0 f(x)=0 f(x)=0，将它转换成等价形式： x = φ ( x ) . x=\varphi (x). x=φ(x).给定初始值 x 0 x_0 x0，构造迭代序列 x k + 1 = φ ( x k ) , k = 1 , 2 , . . . x_{k+1}=\varphi(x_k), k=1,2,... xk+1=φ(xk),k=1,2,...，如果迭代收敛，
lim ⁡ k → ∞ x k + 1 = lim ⁡ k → ∞ φ ( x k ) = α , \lim_{k\to\infty}x_{k+1}=\lim_{k\to\infty}\varphi(x_k)=\alpha, k→∞limxk+1=k→∞limφ(xk)=α,
有 α = φ ( α ) , \alpha=\varphi(\alpha), α=φ(α),则 α \alpha α就是方程 f ( x ) f(x) f(x)的根。在计算中，当 ∣ x k + 1 − x k ∣ |x_{k+1}-x_k| ∣xk+1−xk∣小于给定精度控制量时，取 x k + 1 x_{k+1} xk+1为方程的根。
例如，代数方程 x 3 − 2 x − 5 = 0 x^3-2x-5=0 x3−2x−5=0的三种等价形式及其迭代格式：
1.等价形式： x 3 = 2 x + 5 , x = 2 x + 5 3 , x^3=2x+5, x=\sqrt[3]{2x+5}, x3=2x+5,x=32x+5 ,
　迭代格式： x k + 1 = 2 x k + 5 3 x_{k+1}=\sqrt[3]{2x_k+5} xk+1=32xk+5
　
2.等价形式： 2 x = x 3 − 5 , x = x 3 − 5 2 2x=x^3-5, x=\frac{x^3-5}{2} 2x=x3−5,x=2x3−5

3.等价形式： x 3 = 2 x + 5 , x = 2 x + 5 x 2 , x^3=2x+5,x=\frac{2x+5}{x^2}, x3=2x+5,x=x22x+5,
迭代形式： x k + 1 = 2 x k + 5 x k 2 x_{k+1}=\frac{2x_k+5}{x_k^2} xk+1=xk22xk+5

对于非线性方程 f ( x ) = 0 f(x)=0 f(x)=0，构造的各种迭代格式 x k + 1 = φ x k x_{k+1}=\varphi{x_k} xk+1=φxk，怎样判断其是否收敛？收敛于迭代的初始值是否有关？

因此可知，对于上面给出的代数方程 x 3 − 2 x − 5 = 0 x^3-2x-5=0 x3−2x−5=0，如果要求在 x 0 = 2 x_0=2 x0=2附近的实根，则只有第一种迭代形式收敛。

1.4 牛顿迭代法

对于非线性方程 f ( x ) = 0 f(x)=0 f(x)=0可构造多种迭代格式 x k + 1 = φ ( x k ) x_{k+1}=\varphi(x_k) xk+1=φ(xk)，牛顿迭代法是借助于对函数 f ( x ) f(x) f(x)在 x 0 x_0 x0出做泰勒展开，取其线性部分构造的一种迭代格式。
将 f ( x ) = 0 f(x)=0 f(x)=0在初始值 x 0 x_0 x0做泰勒展开：
f ( x ) = f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) + f ′ ′ ( x 0 ) 2 ! ( x − x 0 ) 2 + . . . f(x)=f(x_0)+f'(x_0)(x-x_0)+\frac{f''(x_0)}{2!}(x-x_0)^2+... f(x)=f(x0)+f′(x0)(x−x0)+2!f′′(x0)(x−x0)2+...
取展开式的线性部分作为 f ( x ) f(x) f(x)的近似，则有，
f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) = 0 f(x_0)+f'(x_0)(x-x_0)=0 f(x0)+f′(x0)(x−x0)=0,设， f ′ ( x 0 ) ≠ 0 f'(x_0)\neq0 f′(x0)=0，则，
x = x 0 − f ( x 0 ) f ′ ( x 0 ) x=x_0-\frac{f(x_0)}{f'(x_0)} x=x0−f′(x0)f(x0)
因此得到牛顿迭代法的迭代格式：
x k + 1 = x k − f ( x k ) f ′ ( x k ) , k = 1 , 2 , . . . x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}, k=1,2,... xk+1=xk−f′(xk)f(xk),k=1,2,...
　题外话：牛顿法不仅可以用于求非线性方程的根，还可以求最优化问题（等价于求目标函数导数/偏导数等于0的根），而且求解精确，一般用牛顿法求得的解称为ground-truth。

1.5 弦截法

在牛顿迭代法中： x k + 1 = x k − f ( x k ) f ′ ( x k ) , k = 1 , 2 , . . . x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}, k=1,2,... xk+1=xk−f′(xk)f(xk),k=1,2,...需要求解函数 f ( x ) f(x) f(x)的一阶导数，通常比较耗时，因此可用差商代替导数：
f ′ ( x k ) = f ( x k ) − f ( x k − 1 ) x k − x k − 1 f'(x_k)=\frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}} f′(xk)=xk−xk−1f(xk)−f(xk−1)
给定初始值 x 0 x_0 x0和 x 1 x_1 x1，那么弦截法的迭代格式为：
x k + 1 = x k − ( x k − x k − 1 ) f ( x k ) f ( x k ) − f ( x k − 1 , k = 1 , 2 , . . . x_{k+1}=x_k - \frac{(x_k-x_{k-1})f(x_k)}{f(x_k)-f(x_{k-1}}, k=1,2,... xk+1=xk−f(xk)−f(xk−1(xk−xk−1)f(xk),k=1,2,...

2. 非线性方程组的求解

设含有n个未知数的n个方程的非线性方程组为，
F ( x ) = 0 F(x)=0 F(x)=0
其中， x = ( x 1 , x 2 , . . . , x n ) T x=(x_1,x_2,...,x_n)^T x=(x1,x2,...,xn)T为n维列向量，
F ( x ) = ( f 1 ( x ) , f 2 ( x ) , . . . , f n ( x ) ) T F(x)=(f_1(x),f_2(x),...,f_n(x))^T F(x)=(f1(x),f2(x),...,fn(x))T
f i ( x ) ( i = 1 , 2 , . . . , n ) f_i(x)(i=1,2,...,n) fi(x)(i=1,2,...,n)中至少有一个是 x x x的非线性函数，并假设自变量和函数值都是实数。

2.1 不动点迭代法

多元非线性方程组与一元非线性方程 f ( x ) = 0 f(x)=0 f(x)=0具有相同的形式，可以与一元非线性方程并行地讨论它的迭代解法，例如不动点迭代法和Newton型迭代法。
把方程组的每一个方程改写成合适的迭代形式，即可得方程组的不动点迭代形式。
例如，有非线性方程组
{ x 1 2 − 10 x 1 + x 2 2 + 8 = 0 x 1 x 2 2 + x 1 − 10 x 2 + 8 = 0 \begin{cases} x_1^2-10x_1+x_2^2+8=0\\ x_1x_2^2+x_1-10x_2+8=0 \end{cases} {x12−10x1+x22+8=0x1x22+x1−10x2+8=0
写成等价形式：
{ x 1 = φ 1 ( x 1 , x 2 ) = 1 10 ( x 1 2 + x 2 2 + 8 ) x 2 = φ 2 ( x 1 , x 2 ) = 1 10 ( x 1 x 1 2 + x 1 + 8 ) \begin{cases} x_1=\varphi_1(x_1,x_2)=\frac{1}{10}(x_1^2+x_2^2+8)\\ x_2=\varphi_2(x_1,x_2)=\frac{1}{10}(x_1x_1^2+x_1+8) \end{cases} {x1=φ1(x1,x2)=101(x12+x22+8)x2=φ2(x1,x2)=101(x1x12+x1+8)
由此构造不动点迭代法：
{ x 1 ( k + 1 ) = φ 1 ( x 1 ( k ) , x 2 ( k ) ) = 1 10 [ ( x 1 ( k ) ) 2 + ( x 2 ( k ) ) 2 + 8 ] x 2 ( k + 1 ) = φ 2 ( x 1 ( k ) , x 2 ( k ) ) = 1 10 [ x 1 ( k ) ( x 1 ( k ) ) 2 + x 1 ( k ) + 8 ] \begin{cases} x_1^{(k+1)}=\varphi_1(x_1^{(k)},x_2^{(k)})=\frac{1}{10}[(x_1^{(k)})^2+(x_2^{(k)})^2+8]\\ x_2^{(k+1)}=\varphi_2(x_1^{(k)},x_2^{(k)})=\frac{1}{10}[x_1^{(k)}(x_1^{(k)})^2+x_1^{(k)}+8] \end{cases} {x1(k+1)=φ1(x1(k),x2(k))=101[(x1(k))2+(x2(k))2+8]x2(k+1)=φ2(x1(k),x2(k))=101[x1(k)(x1(k))2+x1(k)+8]

2.1 牛顿法

对于非线性方程组，也可以构造类似于一元方程的Newton迭代法。设 x ∗ x^* x∗是方程组的解， x ( k ) x^{(k)} x(k)是方程组的一个近似解。
用点 x ( k ) x^{(k)} x(k)处的一阶 Taylaor 展开式近似每一个分量函数 f i ( x ∗ ) = 0 f_i(x^*)=0 fi(x∗)=0值，有，
f i ( x ∗ ) ≈ f i ( x ( k ) + ∑ j = 1 n ∂ f i ( x ( k ) ) ∂ x j ( x j ∗ − x j ( k ) ) , i = 1 , 2 , . . . , n f_i(x^*)\approx f_i(x^{(k)}+\sum_{j=1}^n\frac{\partial f_i(x^{(k)})}{\partial x_j}(x_j^*-x_j^{(k)}), i=1,2,...,n fi(x∗)≈fi(x(k)+j=1∑n∂xj∂fi(x(k))(xj∗−xj(k)),i=1,2,...,n
写成向量形式有 F ( x ∗ ) ≈ F ( x ( k ) + F ′ ( x ( k ) ) ( x j ∗ − x ( k ) ) F(x^*)\approx F(x^{(k)}+F'(x^{(k)})(x_j^*-x^{(k)}) F(x∗)≈F(x(k)+F′(x(k))(xj∗−x(k))
其中 F ′ ( x ( k ) ) F'(x^{(k)}) F′(x(k))为 F ( x ) F(x) F(x)的Jacobi矩阵 F ′ ( x ) F'(x) F′(x)在 x ( k ) x^{(k)} x(k)的值，若 F ′ ( x ( k ) ) F'(x^{(k)}) F′(x(k))非奇异，可得牛顿迭代法：
x ( k + 1 ) = x ( k ) − ( F ′ ( x ( k ) ) ) − 1 F ( x ( k ) ) , x^{(k+1)}=x^{(k)}-(F'(x^{(k)}))^{-1}F(x^{(k)}), x(k+1)=x(k)−(F′(x(k)))−1F(x(k)),
在实际计算中，由于计算矩阵的逆十分耗时，因此牛顿法的第k步是先解线性方程组
F ′ ( x ( k ) ) △ x ( k ) = − F ( x ( k ) ) F'(x^{(k)})\triangle x^{(k)}=-F(x^{(k)}) F′(x(k))△x(k)=−F(x(k))
解出 △ x ( k ) \triangle x^{(k)} △x(k)后，再令 x ( k + 1 ) = x ( k ) + △ x ( k ) x^{(k+1)}=x^{(k)}+\triangle x^{(k)} x(k+1)=x(k)+△x(k)，其中的计算量包括计算向量 F ( x ( k ) ) F(x^{(k)}) F(x(k))和矩阵 F ′ ( x ( k ) F'(x^{(k)} F′(x(k)。

2.2 阻尼牛顿法

虽然Newton法具有二阶局部收敛性，但它要求 F ′ ( x ∗ ) F'(x^*) F′(x∗)非奇异。如果矩阵 F ′ ( x ∗ ) F'(x^*) F′(x∗)奇异或病态，那么 F ′ ( x ( k ) F'(x^{(k)} F′(x(k)也可能奇异或病态，从而可能导致数值计算失败或产生数值不稳定。这时可使用阻尼牛顿法，即把牛顿法的迭代公式改成，
[ F ′ ( x ( k ) ) + μ k I ] △ x ( k ) = − F ( x ( k ) ) , k = 0 , 1 , . . . [F'(x^{(k)}) + \mu_k I]\triangle x^{(k)}=-F(x^{(k)}), k=0,1,... [F′(x(k))+μkI]△x(k)=−F(x(k)),k=0,1,...
其中的参数 μ k \mu_k μk称为阻尼因子， μ k I \mu_k I μkI称为阻尼项。增加阻尼项的目的，是使线性方程的系数矩阵非奇异并良态。当 μ k \mu_k μk选得很合适时，阻尼Newton法是线性收敛的。

2.3 拟牛顿法

类似于弦截法，导数 F ′ ( x ) F'(x) F′(x)需要近似。
比较流行的有DFP，BFGS，L-BFGS。
to be continued.

3.总结

本文介绍了非线性方程（组）常用的解法。

4.参考资料

中科大精品课程：
http://www.bb.ustc.edu.cn/jpkc/xiaoji/szjsff/jsffkj/chapt4_1.htm
中南大学数学科学与计算技术学院
http://netclass.csu.edu.cn/jpkc2008/China/09WebMaths/Teach_Method/Coursesppt/6/6.4%20%E9%9D%9E%E7%BA%BF%E6%80%A7%E6%96%B9%E7%A8%8B%E7%BB%84%E7%9A%84%E6%95%B0%E5%80%BC%E8%A7%A3%E6%B3%95.pdf