这学期正在学数值分析,看到这个板块有人发了高斯求积,感到非常熟悉,因此我也来发一个(新人首帖,请大家温柔点)。这个帖子主要介绍求解非线性方程。
0. 准备知识:
引理1.
设 g(x) \in C[a,b], 且 g([a,b]) \subset [a,b], 则 \exists x_0 \in [a,b] , 使得 g(x_0)=x_0 。 x_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] 内有唯一的不动点。
注:唯一性可使用反证法。
不动点定理
设
- g(x) \in C^1[a,b], 且 g([a,b]) \subset [a,b] ;
- K>0 是一个常数;
- 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 越小,收敛得越快。 因此 用 g 对 x 求导, 并让它等于 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, 那么 a 是 f 的 M 重根.
加速的 Newton-Raphson 公式可以保证当 a 是 f 的重根时,也是二次收敛的.
想不出例子了,暂时不给例子了,有好的例子可以说一下,有什么说得不清楚的地方也欢迎指出.