数值计算方法 Chapter4. 非线性方程求根本文已参与「新人创作礼」活动，一起开启掘金创作之路。数值计算方法 Ch

本文已参与「新人创作礼」活动，一起开启掘金创作之路。

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

0. 问题描述

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

1. 实根的对分法

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

$f(x)$ 在区间 $[a,b]$ 上连续；
$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^*$ ，且有误差估计式：
$|x^* - x_k| \leq \frac{L^k}{1-L}|x_1 - x_0|$

同样的，我们可以给出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展开公式如下：

f(x) = f(x_0) + f'(x_0) \cdot (x-x_0) + \frac{f''(x_0)}{2!} \cdot (x-x_0)^2 + ... + \frac{f^{(n)}(x_0)}{n!} \cdot (x-x_0)^n + ...

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

x \simeq x_0 + \frac{f(x) - f(x_0)}{f'(x_0)}

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

x \simeq x_0 - \frac{f(x_0)}{f'(x_0)}

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

x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}

其物理含义是说：

每次在曲线的某一个点上做一条切线，取该切线与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迭代法的一个近似版本，具体来说，就是将导数用弦截公式进行替换。

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

x_{k+1} = x_k - f(x_k) \cdot \frac{x_k - x_{k-1}}{f(x_k) - f(x_{k-1})}

同样给出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