MENU

Numerical Analysis: 方程求解的不动点迭代与牛顿拉夫逊法

• September 21, 2019 • Read: 4696 • Knowledge,Numerical Analysis

不动点迭代

不动点迭代法是利用函数的不动点原理来构造迭代关系式的一种方法。函数的不动点的定义为:

对于函数$f(x)$来说,如果存在一个值$x_0$,并且这个值带入函数后,满足这个关系:

$$ f(x_0) = x_0 $$

那么就称$x_0$ 为函数$f(x)$的一个不动点。

比如,对于余弦函数$\cos(x)$来说,0.739085是它的不动点。

对函数来说,不动点原理形成了一个迭代关系式,因此,借助于这个思路,在方程的求解时,努力地构造出这样的形式,把方程变换为一侧为$x$,另一侧为一个函数即可。那么就可以得到迭代关系式:

$$ x_{k+1} = f(x_{k})\quad k=0,1,2,3,....\infty $$

但是,FPI中非常重要的一点就是,这样的迭代关系能否保证收敛?在这里,只要保证迭代式中的函数在不动点附近的斜率绝对值小于1即可,即导数值小于1,这样就能保证收敛。

例如,求方程$x^3 + x - 1 = 0$的根。

如果我们构造这样的形式:

$$ x = 1 - x ^ 3 $$

通过计算,我们会发现这一迭代计算是发散的,不会得到一个有用的估计值。因为函数$1 - x ^ 3$不能完全的保证它的导数的绝对值小于1。

所以,要想得到有效的形式,必须要是一侧函数的斜率绝对值小于1。得到:

$$ x = \sqrt[3]{1 - x} $$

import numpy as np
import matplotlib.pyplot as plt

# 定义一个fpi的函数用于计算方程的根
def myfpi(f, x0, iter_max, eps):
    # x0为初始值
    # f为构造后的函数
    # iter_max最大的迭代次数,避免陷入死循环
    # eps为阀值,当误差小于这个这个值停止计算
    
    number = 0  # 用于记录迭代的次数,先初始化为0
    while np.abs(x0 - f(x0)) > eps:
        print('the error is: ',np.abs(x0 - f(x0)))
        x1 = f(x0)
        x0 = x1
        number += 1
        if number > iter_max:
            print('Warning: The result is not stable!')
            return
    return x0
    
  # 定义迭代式中的函数

def func1(x):  # 这是另外一种,收敛速度更快
    return (1 + 2 * x**3) / (1 + 3 * x**2)

def func2(x):  # (1-x)**(1/3)
    return (1 - x)**(1/3)

x0 = 0.4# 初始值
max_time = 400
eps = 0.0001
res1 = myfpi(func1, x0, max_time, eps)
print(res1)
res2 = myfpi(func2, x0, max_time, eps)
print(res2)

结果的输出为:

the error is:  0.36216216216216224
the error is:  0.07470645955965782
the error is:  0.005105524662041616
0.6823501779404628
the error is:  0.4434326653017492
the error is:  0.30445961169099023
the error is:  0.2335452355366372
the error is:  0.1620700735059284
the error is:  0.11988614320526159
the error is:  0.08427085993241235
the error is:  0.06129859706965046
the error is:  0.04345059424243991
the error is:  0.03135485804934435
the error is:  0.022330818149177367
the error is:  0.01605160099338865
the error is:  0.011460910899801835
the error is:  0.008222159662768846
the error is:  0.00587845620934957
the error is:  0.004213100804023595
the error is:  0.0030142409755113952
the error is:  0.0021592285329744554
the error is:  0.0015453565272561809
the error is:  0.0011067207210083696
the error is:  0.0007922226584874403
the error is:  0.0005672828725480361
the error is:  0.00040611545134028315
the error is:  0.00029078551365002703
the error is:  0.00020818212995765695
the error is:  0.00014905678931631527
the error is:  0.00010671688227037457
0.6822832757973981

误差

不论是之前的二分法,还是上述的FPI法,它们都是一个迭代的方法,只不过形成迭代关系时用到的构造原理不同而已,但是仍然都是迭代法的一种。从计算的角度来看,迭代必须采用循环来让其重复计算,知道获得满意的的结果。我们通过设置一个误差标准,当每一次计算的结果的偏离值满足设定的值时,就停止计算。对于while循环来说,必须有入口条件,而数值结果的误差就是描述、控制循环停止强有力的一个标准。

  • 后向、前向误差

后向误差,指的是依找计算的输入的过程来看的,因为计算最终的步骤就是获得一个解$x_a$,因此,当我们使用$f(x_a)$这一项时,称这个为后向误差。那么只要使用了输入值,就是前向误差,前向误差一般都是指,当前的结果与真实值$r$之间的偏离度,因为真实值往往不知道,所以使用前向误差很有逻辑上的困难。

迭代的收敛

当我们使用后向误差来定义误差时,即$e_{i} = |r - x_{i}|$,那么有

$$ \lim_{i \rightarrow \infty} \frac{e_{i + 1}}{e_{i}} = S $$

当这一结果存在,并且收敛到S,则称为线性收敛,收敛速度为$S$。

因此,可以知道的是,二分法的收敛速度是$\frac{1}{2}$,全程保持这个速度;而FPI由于收敛速度与函数的导数值有关,收敛速度并不确定。当这个速度小于$\frac{1}{2}$时,其速度比二分法要高效。而牛顿拉夫逊法,它的收敛速度$S$为0,计算效率极高,并且使用非常广泛。

牛顿拉夫逊法

牛顿法计算的原理,在于利用当前估计值处的函数导数(斜率)来构造函数的切线,然后利用这个切线方程与$x$轴的交点,来确定一个估计值,一直迭代循环计算,直到满足误差要求。

若我们所要求的方程为$f(x) = 0$,设定一个初始的估计值$x_0$,对应的函数值为$f(x_0)$,由此确定了一个函数曲线上的一个点$(x_0, f(x_0))$,则该点处的斜率(导数值)为$f'(x_0)$,由此可以确定处经过该点的直线方程,然后利用这个直线方程去确定新的估计值。

$$ f'(x_0)(x - x_0) = 0 - f(x_0)\\ x - x_0 = - \frac{f(x_0)}{f'(x_0)}\\ x = x_0 - \frac{f(x_0)}{f'(x_0)} $$

由此得到迭代关系为:

$$ x_{k+1} = x_{k} - \frac{f(x_k)}{f'(x_k)} \quad k=0,1,2,3,\cdots $$

牛顿法的求解方式思路清晰,很容易理解。编程实现时,比较麻烦的是必须要确定出函数的导数形式。

def newton_method(f, fd, x0, iter_max, eps):
    number = 0
    while np.abs(f*(x0)) > eps:
        x1 = x0 - f(x0) / fd(x0)
        x0 = x1
        number += 1
        if number > iter_max:
            print('Warning: The result is not stable!')
            return
    return x0
Last Modified: October 6, 2019
Archives Tip
QR Code for this page
Tipping QR Code