数值计算方法 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=\phi (x)$。

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

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

定理 4.1
若$\phi (x)$定义在$[a,b]$上，且满足：
(1) 当$x\in [a,b]$时，有$a\le \phi (x)\le b$；
(2) $\phi (x)$在$[a,b]$上可到，且存在正数$L<1$，使得对于任意$x\in [a,b]$，有$|{\phi }^{\prime }(x)|\le L$；
则在$[a,b]$上有唯一的点$x^{\ast }$满足$x^{\ast }=\phi (x^{\ast })$，称$x^{\ast }$为$\phi (x)$的不动点，且迭代格式$x_{k+1}=x_k$对任意的初值$x_0\in [a,b]$均收敛于$\phi (x)$的不动点$x^{\ast }$，且有误差估计式：
$$|x^{\ast }-x_k|\le \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^{\prime }(x_0)\cdot (x-x_0)+\frac{f^{\prime \prime }(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^{\prime }(x_0)}$$

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

$$x\simeq x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$$

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

$$x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(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