数值计算方法 Chapter4. 非线性方程求根

• 数值计算方法 Chapter4. 非线性方程求根
  • 0. 问题描述
  • 1. 实根的对分法
  • 2. 迭代法
  • 3. Newton迭代法
  • 4. 弦截法

0. 问题描述

给定一个复杂方程f(x) = 0 ，如果直接求解其解析解非常复杂或者难以求解的话，那么可以通过数值求解的方法得到一定精度条件下的数值解。

1. 实根的对分法

对分法使用的条件需要满足：

1. f(x) 在区间[a,b] 上连续；
2. f(a) \cdot f(b) < 0 ；那么，f(x) 在区间[a,b] 中必然存在至少一个零点，我们可以使用二分法不断地迭代求解。

给出python伪代码如下：

def bisect_solve(fn, a, b, epsilon=1e-9):
    assert(fn(a) * fn(b) < 0)
    while b - a >= epsilon:
        m = (a + b)/2
        if fn(m) * fn(a) > 0:
            a = m
        elif fn(m) * fn(b) > 0:
            b = m
        else:
            return m
    return (a+b)/2

2. 迭代法

迭代法的思路是说，将方程f(x) = 0 改写为方程x = \varphi(x) 。

此时，我们可以构造迭代关系x_{k+1} = \varphi(x_k) ，如果数列x_k 收敛，那么数列x_k 的极限 lim_{k \to \infty} x_k 即为目标方程的解。

而关于数列的收敛条件的判断，我们有定理：

定理 4.1 若\varphi(x) 定义在[a,b] 上，且满足： (1) 当x \in [a, b] 时，有a \leq \varphi(x) \leq b； (2) \varphi(x) 在[a,b] 上可到，且存在正数L < 1 ，使得对于任意x \in [a, b] ，有|\varphi'(x)| \leq L ； 则在[a, b] 上有唯一的点x^* 满足x^* = \varphi(x^*) ，称x^* 为\varphi(x) 的不动点，且迭代格式x_{k+1} = x_{k} 对任意的初值x_0 \in [a,b] 均收敛于\varphi(x) 的不动点x^* ，且有误差估计式：

同样的，我们可以给出python伪代码实现如下：

def iter_solve(fn, x, epsilon=1e-9):
    MAX_ITER_TIME = 10**7
    for _ in range(MAX_ITER_TIME):
        y = fn(x)
        if abs(y-x) <= epsilon:
            return y
        x = y
    return x

3. Newton迭代法

Newton迭代法的思路来源于泰勒展开，给出Talyer展开公式如下：

如果视二阶以上的结果为小量，则有：

当x 为零点时，有f(x) = 0 ，因此，我们近似即有：

由此，我们即可给出Newton迭代公式如下：

其物理含义是说：

• 每次在曲线的某一个点上做一条切线，取该切线与x轴的交点作为下一个迭代点。

给出python伪代码如下：

def newton_solve(fn, dfn, x, epsilon=1e-9):
    MAX_ITER_TIME = 10**7
    for _ in range(MAX_ITER_TIME):
        y = x - fn(x)/dfn(x)
        if abs(y-x) <= epsilon:
            return y
        x = y
    return x

4. 弦截法

弦截法其实就是Newton迭代法的一个近似版本，具体来说，就是将导数用弦截公式进行替换。

给出弦截法的迭代公式如下：

同样给出python伪代码如下：

def secant_solve(fn, x1, x0, epsilon=1e-9):
    MAX_ITER_TIME = 10**7
    for _ in range(MAX_ITER_TIME):
        y = x1 - fn(x1) * (x1 - x0) / (fn(x1) - fn(x0))
        if abs(y-x1) <= epsilon:
            return y
        x1, x0 = y, x1
    return x