牛顿法及其收敛性
- 0. 引论
- 2. 牛顿法
- 2.1 单变量牛顿法
- 2.2 多变量牛顿法
- 3. 拟牛顿法
- 3.1 平行弦法
- 3.2 牛顿下山法
- 3.3 Broyden方法
- 4. 代码及算法对比
- 4.1 原始牛顿法
- 4.2 简单牛顿法
- 4.3 Broyden1算法
- 4.3 Broyden2算法
- 4.4 各种算法迭代速度
0. 引论
在工程中，很多的问题都可以归结为求解一组偏微分方程。描述如下：在某一个区域\(\Omega\)中，物理量\(u={u_1,u_2,\cdots,u_n}\)满足一组控制方程，并且在边界上满足若干条件：
控制方程利用两个微分算子来表示\(A(u),B(u)\)。
\[A(u)=\begin{cases} A_1(u)& =0 \\ A_2(u)& =0 \\ &\vdots \end{cases} \in \Omega \tag{1} \]\[B(u)=\begin{cases} B_1(u)& =0 \\ B_2(u)& =0 \\ &\vdots \end{cases} \in \partial \Omega \tag{2} \]以上的问题，只有在求解域比较简单，控制方程形式比较特殊的时候才能求出解析表达，一般情况下，工程问题中的求解域是比较复杂的，控制方程往往是高阶非线性的，这个时候，只能借助于数值求解。偏微分方程的数值解法，常见的有有限元法，有限差分法，有限体积法等。这些数值解法的基本思想是用有限的自由度插值逼近真实解。不管哪一种数值方式，偏微分方程系统离散之后，形成的是常微分方程组或者代数方程组。
在工程中，对于时间维度往往是特别对待的，含有时间项的问题称为非定常问题，也叫瞬态问题，这类问题在空间离散之后，形成的是关于时间的常微分方程，然后利用一些推进求解的方法，将时间维度离散，转换成纯粹的代数方程。而控制方程不含时间项的问题称为稳态问题，对于稳态问题，将控制方程在空间离散之后，形成代数方程组。
可见，不论是那种问题，最后求解的是代数方程组。因此，如何准确高效的求解代数方程组就显得至关重要。这一部分的内容在数值分析里面有详细的介绍，对于工程师来说，一般不必去了解里面艰深的理论，但是，每一种求解方法的适用范围，求解速度，收敛性等重要特征，最好还是要了解。这样才能在使用的时候不至迷茫。
代数方程组可以简记为如下的形式：
\[F(u)=b\tag{3} \]还可以表示成为矢量的形式：
\[Au=b \tag{4} \]当系数矩阵A与变量u的取值没有关系时，问题是线性的，否在，是非线性问题。众所周知，非线性问题的求解是困难的！！
常用的线性方程组求解方法有：
- 直接法
-
- 高斯消去法、矩阵三角分解法（直接分解法、平方根方法、追赶法）
- 迭代法
-
- 雅可比迭代、高斯-赛德尔迭代、超松弛法、块迭代法、共轭梯度法、最速下降法
常用的非线性方程组的求解方法有：
- 迭代法
-
- 二分法、不动点迭代、牛顿法、简化牛顿法、牛顿下山法、弦截法、抛物线法
- 多变量不动点迭代法、牛顿法、拟牛顿法等。
这里主要介绍牛顿法的一些特性。
2. 牛顿法
2.1 单变量牛顿法
对于单变量的牛顿法，大家已经很熟悉了，这里做一点简要的回顾。
方程
\[f(x)=b \tag{5} \]给定初值\(x^0\)，利用迭代关系式:
\[x^{k+1}=x^k-\frac{f(x)-b}{f'(x)} \tag{6} \]经过一定次数的迭代，就可以收敛到其精确值\(x^*\)。
牛顿法存在两个问题：
- 当函数的二阶导数在精确解的邻域内小于0时，无法收敛到精确解。
- 当初值与精确解相差较大时，收敛困难。
2.2 多变量牛顿法
对于多变量牛顿法，其形式可以利用泰勒公式来推导。有\(F(u)=0\)，假设其精确解为\(u^*\)，则有：
\[F(u^*)=0 \tag{7} \]将F(u)在精确解处展开，如下：
\[F(u^*)=F(u)+\nabla F(u)(u^*-u) \tag{8} \]\[[\nabla F(u)]^{-1}F(u)=u-u^* \tag{9} \]因为精确解是不知道的，所以，利用一个初值代替：
\[u^{k+1}=u^k-[\nabla F(u^k)]^{-1}F(u^k) \tag{10} \]对于方程(10)所示的迭代式，影响其收敛性的一个重要指标是其Hesse矩阵。Hesse矩阵定义如下：
\[G(u)=\nabla_u^2F(u) \tag{11} \]对于多变量牛顿法，\(G(u)\)满足Lipschitz条件：
\[||G(y)-G(z)|| <= \beta||y-z|| \tag{12} \]且\(G(u^*)\)正定，当\(u^0\)充分接近\(u^*\)时，迭代是收敛的，且是二阶收敛。
牛顿法的优点：收敛速度快
缺点：计算量大，收敛条件严格。
3. 拟牛顿法
牛顿法具有收敛速度快的特点，但是其缺点也很明显：一是牛顿法计算量很大，每次都要计算Jacobi矩阵的逆矩阵和\(F(x)\)。二是要求迭代初值在精确解附近，不然很难收敛。针对上面的问题，有一些改进的牛顿法。这里最为常用的是拟牛顿法（又称平行弦法）和牛顿下山法。
3.1 平行弦法
迭代式为：
\[\begin{equation} \begin{cases} x^{k+1} & =\phi(x^k) \\ \phi(x^k)& =x^k-Cf(x^k) \tag{13} \end{cases} \end{equation} \]其中，\(C=\frac{1}{f'(x_0)}\) 。该方法要求迭代式的一阶导数范数小于1才是局部收敛，且是线性收敛。
3.2 牛顿下山法
在迭代过程中，满足单调递减的要求的迭代方法称作下山法。即要求迭代过程中，始终有：
\[||f(x^{k+1})||<||f(x^k)|| \tag{14} \]迭代式：
\[\begin{equation} \begin{cases} x^{k+1} & =\phi(x^k) \\ \phi(x^k)& =x^k-\lambda\frac{f(x^k)}{f'(x^k)} \tag{15} \end{cases} \end{equation} \]其中，\(\lambda(0<\lambda \le 1)\)称为下山因子。选择下山因子的时候，从\(\lambda=1\)开始，逐次减半进行试算，直到下降条件（14）成立为止。
3.3 Broyden方法
Broyden提出的方法如下：
- Broyden秩1迭代公式：
其中：\(s_k=x^{k+1}-x^k\)；\(y_k=F(x^{k+1})-F(x^k)\)
- Broyden秩1第二公式
上面两式中矩阵要求逆，不好看！我们用\(B_k = A_{k}^{-1}\)代替，然后对上面的两个方法分别做改写，就可以得到它们的逆方法！
- 逆Broyden秩1迭代公式：
- 逆Broyden秩1第二公式
一般取\(B_0=(F^{'}(x^0))^{-1}\)。
Broyden方法收敛速度介于线性收敛和二阶收敛之间，算法比较稳定。一般使用过程中，如果方程组的非线性不是很复杂，并且要求收敛速度快，可以选择原始牛顿法。如果方程组非线性比较明显，建议采用更为稳定的拟牛顿法，拟牛顿法中推荐使用Broyden方法。
4. 代码及算法对比
下面给出Broyden算法求解非线性方程组的python代码：
\[\begin{equation} \left\{ \begin{array}{l} x_1^2-10x_1+x_2^2+8=0 \\ x_1x_2^2+x_1-10x_2+8=0 \end{array} \right. \end{equation} \]\[\begin{equation} \frac{\part F}{\part x} = \left[ \begin{array}{ll} 2x_1-10 & 2x_2 \\ x_2^2+1 & 2x_1x_2-10 \end{array} \right] \end{equation} \]主程序如下：
# !/usr/bin/python
# -*- coding:utf-8 -*-
import numpy as np
def normal2(x):
tmp = np.sqrt(np.mean(x ** 2))
return tmp
def F(x):
n = len(x)
tmp = np.zeros(n)
x1, x2 = x
tmp[0] = x1 ** 2 - 10 * x1 + x2 ** 2 + 8
tmp[1] = x1 * x2 ** 2 + x1 - 10 * x2 + 8
return tmp
def dF(x):
n = len(x)
tmp = np.zeros((n, n))
x1, x2 = x
tmp[0, 0] = 2 * x1 - 10
tmp[0, 1] = 2 * x2
tmp[1, 0] = x2 ** 2 + 1
tmp[1, 1] = 2 * x1 * x2 - 10
return tmp
def origin_newton(x0, epsilon1=1e-6, epsilon2=1e-6):
pass
def simple_newton(x0, epsilon1=1e-6, epsilon2=1e-6):
pass
def broyden1(x0, epsilon1=1e-6, epsilon2=1e-6):
pass
def broyden2(x0, epsilon1=1e-6, epsilon2=1e-6):
pass
if __name__ == '__main__':
x0 = np.array([2.0, 1.0])
x = newton(x0)
print(x)
4.1 原始牛顿法
原始牛顿法代码如下：
def origin_newton(x0, epsilon1 = 1e-6, epsilon2=1e-6):
FX = F(x0)
while True:
dFX = dF(x0)
J = np.linalg.inv(dFX)
v = np.matmul(J, FX)
x1 = x0 - v
res1 = normal2(v)
FX = F(x1)
res2 = normal2(FX)
x0 = x1
if res1 < epsilon1 or res2 < epsilon2:
break
return x1
4.2 简单牛顿法
简单牛顿法代码如下：
def simple_newton(x0, epsilon1=1e-6, epsilon2=1e-6):
FX = F(x0)
dFX = dF(x0) # 每次用F'(x0)代替原始牛顿法中的雅可比矩阵的逆
J = np.linalg.inv(dFX)
i = 0
while True:
v = np.matmul(J, FX)
x1 = x0 - v
res1 = normal2(v)
FX = F(x1)
res2 = normal2(FX)
x0 = x1
print(i, res1)
i += 1
if res1 < epsilon1 or res2 < epsilon2:
break
return x1
4.3 Broyden1算法
def broyden2(x0, epsilon1=1e-6, epsilon2=1e-6):
FX = F(x0)
dFX = dF(x0)
B0 = np.linalg.inv(dFX)
i = 0
while True:
s = -np.matmul(B0, FX)
x1 = x0 + s
y = F(x1) - FX
# update B1
t1 = (s - np.matmul(B0, y)).reshape(-1, 1)
t2 = (s - np.matmul(B0, y)).reshape(1, -1)
t3 = np.matmul(t2, y.reshape(-1, 1))
t4 = t2 / t3
B1 = B0 + np.matmul(t1, t4)
FX = F(x1)
# update xk, Bk
x0 = x1
B0 = B0
# residual
res1 = normal2(s)
res2 = normal2(FX)
print(i, res1)
i += 1
if res1 < epsilon1 or res2 < epsilon2:
break
return x1
4.3 Broyden2算法
def broyden2(x0, epsilon1=1e-6, epsilon2=1e-6):
FX = F(x0)
dFX = dF(x0)
B0 = np.linalg.inv(dFX)
i = 0
while True:
s = -np.matmul(B0, FX)
x1 = x0 + s
y = F(x1) - FX
# update B1
t1 = (s - np.matmul(B0, y)).reshape(-1, 1)
t2 = (s - np.matmul(B0, y)).reshape(1, -1)
t3 = np.matmul(t2, y.reshape(-1, 1))
t4 = t2 / t3
B1 = B0 + np.matmul(t1, t4)
FX = F(x1)
# update xk, Bk
x0 = x1
B0 = B0
# residual
res1 = normal2(s)
res2 = normal2(FX)
print(i, res1)
i += 1
if res1 < epsilon1 or res2 < epsilon2:
break
return x1
4.4 各种算法迭代速度
对于上述非线性方程组，各种算法的收敛速度如下：
- 简单牛顿法：
0 0.8408515029421069
1 0.2110161529812992
2 0.10520709965532886
3 0.04958208146353929
4 0.025133245218497215
5 0.012422741268388458
6 0.00624079352688798
7 0.0031120861644459564
8 0.0015580179097346145
9 0.0007785001531163414
10 0.00038937566825012144
11 0.00019465619965536809
12 9.733598330013241e-05
13 4.866601708332767e-05
14 2.4333501790476855e-05
15 1.216662752552937e-05
16 6.083344599211488e-06
17 3.041664589374761e-06
18 1.520834221898411e-06
19 7.604166291096606e-07
[1.00000025 1.00000025]
- 原始牛顿法：
0 0.8408515029421069
1 0.1350278246849652
2 0.004998213941639314
3 9.73120102030655e-06
[1. 1.]
- Broyden1算法
0 0.8408515029421069
1 0.2110161529812992
2 0.10520709965532886
3 0.04958208146353929
4 0.025133245218497215
5 0.012422741268388458
6 0.00624079352688798
7 0.0031120861644459564
8 0.0015580179097346145
9 0.0007785001531163414
10 0.00038937566825012144
11 0.00019465619965536809
12 9.733598330013241e-05
13 4.866601708332767e-05
14 2.4333501790476855e-05
15 1.216662752552937e-05
16 6.083344599211488e-06
17 3.041664589374761e-06
18 1.520834221898411e-06
19 7.604166291096606e-07
[1.00000025 1.00000025]
- Broyden2算法
0 0.8408515029421069
1 0.2110161529812992
2 0.10520709965532886
3 0.04958208146353929
4 0.025133245218497215
5 0.012422741268388458
6 0.00624079352688798
7 0.0031120861644459564
8 0.0015580179097346145
9 0.0007785001531163414
10 0.00038937566825012144
11 0.00019465619965536809
12 9.733598330013241e-05
13 4.866601708332767e-05
14 2.4333501790476855e-05
15 1.216662752552937e-05
16 6.083344599211488e-06
17 3.041664589374761e-06
18 1.520834221898411e-06
19 7.604166291096606e-07
[1.00000025 1.00000025]