Newton-Raphson 法简介


#1

这学期正在学数值分析,看到这个板块有人发了高斯求积,感到非常熟悉,因此我也来发一个(新人首帖,请大家温柔点)。这个帖子主要介绍求解非线性方程。
0. 准备知识:
引理1.
g(x) \in C[a,b], 且 g([a,b]) \subset [a,b], 则 \exists x_0 \in [a,b] , 使得 g(x_0)=x_0x_0 称为 g 的不动点。

简单来说就是如果一个连续函数 g[a,b] 映射到 [a,b] 的一个子区间,那么它在 [a,b] 内有不动点。具体的证明略过(使用介值定理)。

引理2
g(x) \in C^1[a,b], 且 g([a,b]) \subset [a,b], 并且 \forall x \in (a,b), 都有 |g'(x) | \leq K <1, 那么 g[a,b] 内有唯一的不动点。
:唯一性可使用反证法。
不动点定理

  1. g(x) \in C^1[a,b], 且 g([a,b]) \subset [a,b]
  2. K>0 是一个常数;
  3. p_0 \in (a,b) 不是 g 的不动点。

那么

  • 如果 \forall x \in (a,b), |g'(x) | \leq K <1, 那么 迭代 p_{n+1}=g(p_n) 将收敛到 g 的不动点.
  • 如果 \forall x \in (a,b), |g'(x) | \geq 1, 那么 迭代 p_{n+1}=g(p_n) 将不会收敛到 g 的不动点.

: 证明结论1 使用 Lagrange 中值定理, 结论2 使用反证法。
下面简单证一下结论1(因为有了小推论后面要用), 设不动点为 x^*
|x_{k+1}-x^*|=|g(x_k)-g(x^*)|=|g'(\eta)||x_k-x^*| \leq K|x_k-x^*|\leq \cdots \leq K^k|x_1-x^*| \rightarrow 0 (k\rightarrow \infty)


下面我们可以开始解方程 f(x) =0 了。
一个简单的想法是,两边同时加上 x, 得到 x=x+f(x), 就可以采用迭代 x_{k+1}=x_k+f(x_k), 即 令 g(x)=x+f(x), 然后用不动点定理进行迭代。 但是聪明的你肯定已经发现了,万一 g 不满足定理条件怎么办?

因此,一种补救的办法是: g(x)=x+\lambda f(x), 我们取合适的 \lambda 使得 g 满足定理的条件就好了。并且要让上面的 K 越小越好,因为 K 越小,收敛得越快。 因此 用 gx 求导, 并让它等于 0. 得到 \lambda =-1/f'(x) .
这样我们就得到了迭代关系式
x_{k+1}=g(x_k)=x_k-f(x_k)/f'(x_k)
这就是 Newton-Raphson 迭代法

先给一个例子吧, 求 e^2 的近似值。 ($f(x)=\ln x-2=0$ 的根)
首先分析下: e^2 大致在4~9 之间, 我们可以取初始值为 x_0=4 来试一下, g(x)=x+x(\ln x-2) ,

x=ones(10)
x.=4

f(x)=log(x)-2
g(x)=x-x*f(x)

for i =1:9
    x[i+1]=g(x[i])
end
x

输出结果为

10-element Array{Float64,1}:
 4.0              
 6.454822555520438
 7.327336834983672
 7.388797614083718
 7.389056094409423
 7.38905609893065 
 7.38905609893065 
 7.38905609893065 
 7.38905609893065 
 7.38905609893065 

可以看到 从第六次结果开始,就已经收敛到了真实值了,收敛速度非常快 (二次收敛, 此处不介绍)。
最后再给出一个加速的 Newton-Raphson 公式

g(x)=x-Mf(x)/f'(x)

其中 M f 的在零点处根的重数。
根的重数定义如下:
f(a)=0, f'(a)=0,\cdots, f^{(M-1)}(a)=0, 但是 f^{(M)}(a)\neq 0, 那么 afM 重根.
加速的 Newton-Raphson 公式可以保证当 af 的重根时,也是二次收敛的.


想不出例子了,暂时不给例子了,有好的例子可以说一下,有什么说得不清楚的地方也欢迎指出.


#2

我就再加一个多变量的版本的例子吧。

julia> using ForwardDiff: jacobian

julia> function f(b)
           x, y = b
          [(x+3)*(y^3-7)+18, sin(y*exp(x)-1)]
       end
f (generic function with 1 method)

julia> x = [1., 1.]
2-element Array{Float64,1}:
 1.0
 1.0

julia> for i in 1:6
           global x
           J = jacobian(f, x)
           x -= J\f(x)
           @info "x: $(x[1])\ty: $(x[2])"
       end
[ Info: x: 2.317487450497444    y: 2.158743725248722
[ Info: x: 2.6221688016462608   y: 1.6851907384212415
[ Info: x: 2.680315902767905    y: 1.5719895457684911
[ Info: x: 2.687043688302479    y: 1.5652758423355826
[ Info: x: 2.6870597848742763   y: 1.565256193249113
[ Info: x: 2.6870597850598585   y: 1.5652561930721185

julia> f(x)
2-element Array{Float64,1}:
 -3.552713678800501e-15
  8.572527594031472e-16

#3

原理是不是这样的? 我们多元的部分讲得比较简单,所以我不是很懂
f: R^n \to R^n ,
0= f(x^*)\approx f(x)+J_f(x)(x^*-x)
那么 J_f(x)(x^*-x)=J_f(x) d=-f(x), 其中 J_f(x)fx 处的 Jacobi 矩阵,
x^*=x+d, d=-J_f(x)\backslash f(x)
个人理解 x_{k+1} 是比 x_k 更接近 x^*, 所以把 x^* 换成了 x_{k+1}.
于是就有迭代过程 x_{k+1}=x_k+d_k .
我感觉基本思路其实还是 g(x)=x+\lambda f(x), 然后让 g 的导数(梯度)等于0.
f 的Jacobi 矩阵在 x^* 处为奇异 (Singular) 矩阵时,可以将采取下面的做法:
(J_f(x_k)+\mu_k I)d_k=-f(x_k), 其中 \mu_k 是一个接近但不是0的数,用于解决奇异性的问题.


#4

Jacobian矩阵。Jacobi矩阵是另一个东西。

牛顿迭代就是迭代解非线性函数的线性逼近。所以有 k(\vec{x}) = f(\vec{x}_n) + f'(\vec{x})(\vec{x}-\vec{x}_n) = 0 ,得到 f'(\vec{x})\Delta x = -f(\vec{x}_n) \implies \Delta x = -f'(\vec{x})\backslash f(\vec{x}_n) 。对于多变量函数,线性逼近就是雅克比矩阵。

写成例子是这样的。

julia> using ForwardDiff: jacobian

julia> using LinearAlgebra

julia> function f(x)
           [10x[1] - 2x[2]^2 + x[2] - 2x[3] - 5,
            8x[2]^2 + 4x[3]^2 - 9,
            8x[2]x[3] + 4]
       end
f (generic function with 1 method)

julia> x = [0., 0., 0.];

julia> rank(jacobian(f, x)) # singular Jacobian
1

julia> report(i, residual) = @info "$i-th iteration: residual = $residual"
report (generic function with 1 method)

julia> for i in 1:80
           global x
           J = jacobian(f, x)
           F = lu(J, check=false)
           if F.info != 0 # singular Jacobian
               F = lu(J+sqrt(eps(maximum(J)))*I)
           end
           x -= F\f(x)
           residual = norm(f(x), Inf)
           i % 5 == 0 && report(i, residual)
           residual <= 100eps() && (report(i, residual); break)
       end
[ Info: 5-th iteration: residual = 1.5657045579530212e15
[ Info: 10-th iteration: residual = 1.5290083573730002e12
[ Info: 15-th iteration: residual = 1.4931722210000002e9
[ Info: 20-th iteration: residual = 1.4581730000037195e6
[ Info: 25-th iteration: residual = 1421.0038062737062
[ Info: 30-th iteration: residual = 0.226533298221808
[ Info: 33-th iteration: residual = 0.0

julia> f(x)
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

#5

要注意的是这个方法只是没有办法的办法,遇到奇异雅克比矩阵那样做收敛很慢。比如说上面的用了33次迭代。

julia> x = [1000,1000,1000.]
3-element Array{Float64,1}:
 1000.0
 1000.0
 1000.0

julia> for i in 1:80
           global x
           J = jacobian(f, x)
           F = lu(J, check=false)
           if F.info != 0 # singular Jacobian
               F = lu(J+sqrt(eps(maximum(J)))*I)
           end
           x -= F\f(x)
           residual = norm(f(x), Inf)
           i % 5 == 0 && report(i, residual)
           residual <= 100eps() && (report(i, residual); break)
       end
[ Info: 5-th iteration: residual = 11715.754804589891
[ Info: 10-th iteration: residual = 10.578288868701229
[ Info: 15-th iteration: residual = 1.1148770795443852e-10
[ Info: 16-th iteration: residual = 0.0

julia> x
3-element Array{Float64,1}:
  0.5
  1.0
 -0.5

用一个很不靠谱,但是不出现奇异雅克比矩阵的最初猜想也就需要16次迭代。