python数值分析——非线性方程迭代求解与Aitken加速

作者: 遇见飖雪 | 来源:发表于2022-04-06 19:11 被阅读0次

1实验内容

用迭代法求下列方程的根
f(x)=x^3+2x2+10x-20=0 (x*=1.368808107821372...)，要求误差< 10-15 。
(1)把方程改写成等价的迭代格式:
x(k+1) =20/(xk^2+2xk+10)，取初值x0=1, 求方程的根。
(2)把方程改写成等价的迭代格式:
xk+1 = 2-(xk^3+2xk2)/10，取初值x0=1, 求方程的根。
(3)对迭代格式(1)和(2)分别用Aitken法进行加速，并对结果进行分析。

2 算法原理

2.1 不动点迭代法

迭代法是一种逐次逼近的方法，用某种固定格式反复校正根的近似值，使之逐步精确，最后达到满足精度要求的结果。

方程f(x)=0化为等价形式的方程x=φ(x) ，构造迭代公式x_k+1=φ(x_k), k=0,1,2,……取初始近似根(迭代初值)x₀，进行迭代计算x1=φ(x₀),x₂=φ(x₁),...,x_k+1=φdf(x_k),……

则有x1，x2 ,……,xk,……，得到迭代序列{x_k}.如果这个序列有极限，则迭代公式是收敛的。

这时则x(x)，x称为函数φ(x)的不动点，等价地有f(x)=0,x^*即为方程的根。连续函数φ(x)称为迭代函数。

实际计算到|xk-xk-1|<ε(ε是预定的精度)，取x^*≈x^k。

迭代公式收敛指迭代序列{x^k}收敛，迭代公式发散指迭代序列{x^k}不收敛，即发散。迭代公式不一定总是收敛。简单迭代法的几何意义如下：

图1 不动点迭代的几何表示

（a）收敛

（b）发散
图2 收敛与发散的几何对比

𝑥=𝑔(𝑥) 的解称为函数 𝑔(𝑥) 的不动点，几何图像为函数图像与正比例函数 𝑦=𝑥的交点。迭代形式为：𝑥𝑘+1=𝑔(𝑥𝑘) 。设函数的不动点为 𝑥∗，当|𝑔′(𝑥∗)|<1 时，在不动点附近的存在一个领域，使得从领域内的某个值开始的不动点迭代收敛，收敛到不动点。通过 𝛿𝑥𝑘+1≈𝑔′(𝑥∗)𝛿𝑥𝑘式，当导数绝对值小于1时，迭代后误差的线性主部较迭代前误差小，以导数绝对值为收敛常数线性收敛。唯一不同的是若导数值为零，此时该迭代至少二次收敛。用不动点迭代解非线性方程的核心在于将非线性方程转化为一个收敛的不动点迭代。

2.2 Aitken加速算法

Aitken加速收敛算法基本原理：

对于收敛的迭代过程，只要迭代足够多次，就可以使结果达到任意的精度。但有时迭代过程收敛缓慢，从而使计算量变得很大，因此，迭代过程的加速是个重要的过程。

2.3 牛顿迭代法原理

牛顿迭代法是用于求解等式方程的一种方法。类似于求解F(x)=0的根，牛顿迭代法求解的是近似根，这个想法准确来说来源于泰勒展开式，需要求解的表达式可能非常复杂，通过一般的方法，我们很难求出它的解。

所以采用了一种近似求解的方法，取泰勒展开式的前几项，队原来的求解函数做一个取代，然后，求解这个取代原方程的方程的解，作为近似解。当然只对原方程做一次近似求解不行，因为第一次近似肯定不会太准确，所以还需要不断地迭代。

首先就要取一个值作为初始的近似值，然后去求解该点的泰勒展开近似项，然后求解根，之后再以此根对原方程进行近似，然后再求解结果不断重复迭代，最终就能求得近似解。

牛顿迭代法迭代公式如下：

设x是f(x) = 0的根，选取x₀作为x初始近似值，过点（x₀, f(x)）做曲线y = f(x)的切线L，则L的方程为y = f(x₀)+f '(x₀)(x-x₀)，求出L与x轴交点的横坐标x₁= x₀-f(x₀) / f '(x₀），称x₁为x的一次近似值。过点（x₁, f(x₁)）做曲线y = f(x)的切线，并求该切线与x轴交点的横坐标x₂= x1-f(x₁)/ f '(x₁)，称x₂为x的二次近似值。重复以上过程，得r的近似值序列，其中x_n+1=x(nf(x_n)/f'(x_n)，称为x的n+1次近似值，上式称为牛顿迭代公式。

解非线性方程f(x)=0的牛顿法是把非线性方程线性化的一种近似方法。把f(x）在x₀点附近展开成泰勒级数f(x) = f(x₀)+(x－x₀)f'(x₀)+(x－x₀)²*f''(x₀)/2! +… 取其线性部分（一次），作为非线性方程f(x) = 0的近似方程，即泰勒展开的前两项，则有f(x₀)+f'(x₀)(x－x₀)=0 设f'(x₀）≠0则其解为x₁=x₀－f(x₀)/f'(x₀) 这样，得到牛顿法的一个迭代序列：x(n+1）= x_n－f(x_n) / f '(x_n)。

牛顿迭代法，取得是泰勒展开式的前两项，也就是线性近似，所以迭代比较快，容易计算。

图3 牛顿迭代法几何表示

3 程序框图

图4 Aitken加速法算法框图

图5 牛顿迭代法算法框图

4 编程实现

利用所编程序，求解下列方程的根:
以下将用Python进行编程实验，利用不动点迭代法与牛顿迭代法求解题目中方程的根，并且利用Aitken加速方法对其进行加速，最后进行结果分析与讨论。
利用不动点迭代法的Python代码详解及运行结果：

（1）把方程改写成等价的迭代格式: xk+1+20/(xk2+2xk+10)，取初值x0=1，求方程的根。
第（1）种迭代格式代码详解：

from sympy import *


def fun():
    """
    求解方程的符号定义
    :return:方程
    error:这里不能用math库中的函数，pow()需要调用符号库中的函数，即sympy.pow()
    或者用from sympy import * 来导入sympy中的所有函数，之后直接写pow()函数即
    是默认的sympy中的函数。
    """
    x = symbols('x')  # 符号变量的定义
    # return 2 * exp(-x) * sin(x) + 2 * cos(x) - 0.25
    return pow(x, 3) + 2.0 * pow(x, 2) + 10.0 * x - 20.0  # 返回函数的值


def fixedPoint(x0, eps, maxIter):
    """
    不动点迭代法的作用是，求解非线性方程的根，采用逐步逼近的方法进行计算
    :param x0:迭代初值
    :param eps:误差精度要求
    :param maxIter:最大迭代次数
    :return:返回值为None
    """
    x = symbols('x')
    fh = fun()
    x_n = x0  # 定义初值
    k = 0  # 初始化迭代次数
    errval = 0  # 初始化误差
    print('%3s %10s %18s' % ('迭代次数', '方程近似值', '迭代误差'))
    for k in range(maxIter):
        x_b = x_n  # 代表x_n
        x_n = 20 / (pow(x_b, 2) + 2 * x_b + 10)  # 写出第一种迭代函数格式
        errval = abs(fh.evalf(subs={x: x_n}))  # 第k次迭代误差的大小
        print('%3d %22.15f %22.15f' % (k + 1, x_n, errval))  # 分别输出迭代次数，方程的近似根以及迭代误差
        if errval <= eps:
            break
    if k + 1 <= maxIter - 1:
        print('方程在满足精度' + str(eps) + '的条件下，近似解为：'
              + str(x_n) + ',误差为：' + str(errval))
    else:
        print('不动点迭代法求解方程的根，已经达到最大迭代次数，也可能不收敛或精度过高...')
    return None


if __name__ == '__main__':
    fh = fun()
    plot(fh)
    x0 = float(input('请输入迭代初值：'))  # input函数总是以字符串的形式返回
    eps = float(input('请输入误差精度要求：'))  # 方程解的精度要求是近似解与真值之间的误差
    maxIter = int(input('请输入最大迭代次数：'))  # 方程一般会迭代无数次，必须定义其迭代的次数，以求收敛
    fixedPoint(x0, eps, maxIter)
    print('方程为：', '%30s' % (str(fun())))
    print('迭代函数为：', '%20s' % ('20 / (x**2 + 2 * x + 10)'))

image.png

（2）把方程改写成等价的迭代格式: xk+1 = 2-(xk3+2xk2)/10，取初值x0=1, 求方程的根。
第二种迭代格式代码详解：

from sympy import *


def fun():
    """
    求解方程的符号定义
    :return:方程
    error:这里不能用math库中的函数，pow()需要调用符号库中的函数，即sympy.pow()
    或者用from sympy import * 来导入sympy中的所有函数，之后直接写pow()函数即
    是默认的sympy中的函数。
    """
    x = symbols('x')  # 符号变量的定义
    # return 2 * exp(-x) * sin(x) + 2 * cos(x) - 0.25
    return pow(x, 3) + 2.0 * pow(x, 2) + 10.0 * x - 20.0  # 返回函数的值


def fixedPoint(x0, eps, maxIter):
    """
    不动点迭代法的作用是，求解非线性方程的根，采用逐步逼近的方法进行计算
    :param x0:迭代初值
    :param eps:误差精度要求
    :param maxIter:最大迭代次数
    :return:返回值为None
    """
    x = symbols('x')
    fh = fun()
    x_n = x0  # 定义初值
    k = 0  # 初始化迭代次数
    errval = 0  # 初始化误差
    print('%3s %10s %18s' % ('迭代次数', '方程近似值', '迭代误差'))
    for k in range(maxIter):
        x_b = x_n  # 代表x_n
        x_n = 2 - (pow(x_b, 3) + 2 * pow(x_b, 2))/10  # 写出第二种迭代函数格式
        errval = abs(fh.evalf(subs={x: x_n}))  # 第k次迭代误差的大小
        print('%3d %22.15f %22.15f' % (k + 1, x_n, errval))  # 分别输出迭代次数，方程的近似根以及迭代误差
        if errval <= eps:
            break
    if k + 1 <= maxIter - 1:
        print('方程在满足精度' + str(eps) + '的条件下，近似解为：'
              + str(x_n) + ',误差为：' + str(errval))
    else:
        print('不动点迭代法求解方程的根，已经达到最大迭代次数，也可能不收敛或精度过高...')
    return None


if __name__ == '__main__':
    fh = fun()
    plot(fh)
    x0 = float(input('请输入迭代初值：'))  # input函数总是以字符串的形式返回
    eps = float(input('请输入误差精度要求：'))  # 方程解的精度要求是近似解与真值之间的误差
    maxIter = int(input('请输入最大迭代次数：'))  # 方程一般会迭代无数次，必须定义其迭代的次数，以求收敛
    fixedPoint(x0, eps, maxIter)
    print('方程为：', '%30s' % (str(fun())))
    print('迭代函数为：', '%20s' % ('2 - (x**3 + 2 * x**2)/10'))

image.png

（3）对迭代格式(1)和(2)分别用Aitken法进行加速，并对结果进行分析。
迭代格式
第（1）种迭代格式Aitken加速法代码详解：

from sympy import *


def Aitken(x0, eps, iterNum, Phi):
    xk_1 = x0
    errval = 0
    print('%3s %10s %18s' % ('迭代次数', '方程近似值', '迭代误差'))
    for k in range(iterNum):
        y = Phi(xk_1)
        z = Phi(y)
        if (z - 2 * y + xk_1) != 0:
            xk = xk_1 - (y - xk_1) ** 2 / (z - 2 * y + xk_1)
            errval = abs(xk - xk_1)
            print('%3d %22.15f %22.15f' % (k + 1, xk, errval))
            if errval < eps:
                return xk
            else:
                xk_1 = xk
        else:
            return xk
    print('方法失败')
    return 0


def Phi(x):
    return 20 / (pow(x, 2) + 2 * x + 10)


if __name__ == '__main__':
    x = symbols('x')
    x0 = float(input('请输入迭代初值：'))  # input函数总是以字符串形式返回
    eps = float(input('请输入迭代误差精度要求：'))  # 方程误差精度要求为迭代近似值与真值之间的差值
    iterNum = int(input('请输入最大迭代次数：'))  # 最大迭代次数限制了方程的收敛
    Phi(x)
    Aitken(x0, eps, iterNum, Phi)

image.png

第（2）种迭代格式Aitken加速法代码详解：

from sympy import *


def Aitken(x0, eps, iterNum, Phi):
    xk_1 = x0
    errval = 0
    print('%3s %10s %18s' % ('迭代次数', '方程近似值', '迭代误差'))
    for k in range(iterNum):
        y = Phi(xk_1)
        z = Phi(y)
        if (z - 2 * y + xk_1) != 0:
            xk = xk_1 - (y - xk_1) ** 2 / (z - 2 * y + xk_1)
            errval = abs(xk - xk_1)
            print('%3d %22.15f %22.15f' % (k + 1, xk, errval))
            if errval < eps:
                return xk
            else:
                xk_1 = xk
        else:
            return xk
    print('方法失败')
    return 0


def Phi(x):
    return 2 - (pow(x, 3) + 2 * pow(x, 2))/10


if __name__ == '__main__':
    x = symbols('x')
    x0 = float(input('请输入迭代初值：'))  # input函数总是以字符串形式返回
    eps = float(input('请输入迭代误差精度要求：'))  # 方程误差精度要求为迭代近似值与真值之间的差值
    iterNum = int(input('请输入最大迭代次数：'))  # 最大迭代次数限制了方程的收敛
    Phi(x)
    Aitken(x0, eps, iterNum, Phi)

image.png

牛顿迭代法代码详解：

from sympy import *


def fun():
    """
    求解方程的符号定义
    :return:方程
    error:这里不能用math库中的函数，pow()需要调用符号库中的函数，即sympy.pow()
    或者用from sympy import * 来导入sympy中的所有函数，之后直接写pow()函数即
    是默认的sympy中的函数。
    """
    x = symbols('x')  # 符号变量的定义
    # return 2 * exp(-x) * sin(x) + 2 * cos(x) - 0.25
    return pow(x, 3) + 2.0 * pow(x, 2) + 10.0 * x - 20.0  # 返回函数的值


def diffFun():
    """
    求解方程的一阶导数
    :return: 一阶导数
    """
    return fun().diff()


def newtonIterative(x0, eps, maxIter):
    """
    牛顿迭代函数的作用，求解非线性方程的根，采用逼近法的思想
    :param x0:迭代初值
    :param eps:误差精度要求
    :param maxIter:最大迭代次数
    :return:返回值为None
    """
    # x1 = x0 - f(x0)/f'(x0)==>x1 = x2 - f(x1)/f'(x1)...
    x = symbols('x')
    fh = fun()  # 引用方程
    dfh = diffFun()  # 引用方程的一阶导数
    x_n = x0
    k = 0
    errval = 0
    print('%3s %10s %18s' % ('迭代次数', '方程近似值', '迭代误差'))
    # 利用牛顿迭代法的思想逐步逼近精确解
    for k in range(maxIter):
        x_b = x_n  # 代表x(n)
        fx = fh.evalf(subs={x: x_b})  # 方程在x_n处的数值
        dfx = dfh.evalf(subs={x: x_b})  # 一阶导数方程在x_n处的方程
        x_n = x_b - fx / dfx  # 牛顿迭代公式
        errval = abs(fh.evalf(subs={x: x_n}))  # 第k次迭代误差的大小
        print('%3d %22.15f %22.15f' % (k + 1, x_n, errval))
        if errval <= eps:
            break
    if k + 1 <= maxIter - 1:
        print('方程在满足精度' + str(eps) + '的条件下，近似解为：'
              + str(x_n) + ',误差为：' + str(errval))
    else:
        print('牛顿迭代法求解数值逼近，已经达到最大迭代次数，也可能不收敛或精度过高...')
    # print(x0, eps, maxIter)
    return None


if __name__ == '__main__':
    fh = fun()
    dfh = diffFun()
    plot(fh)
    plot(dfh)
    x0 = float(input('请输入迭代初值：'))  # input函数总是以字符串的形式返回
    eps = float(input('请输入误差精度要求：'))  # 方程解的精度要求是近似解与真值之间的误差
    maxIter = int(input('请输入最大迭代次数：'))  # 方程一般会迭代无数次，必须定义其迭代的次数，以求收敛
    newtonIterative(x0, eps, maxIter)
    print('方程为：', '%30s' % (str(fun())))
    print('方程的导数为：' '%22s' % (str(diffFun())))

image.png

5 结果分析与讨论

图6 方程的曲线

由方程的曲线可以看出，此方程的根在1附近只有一个根，因此迭代初值取1可以降低迭代次数，节省迭代时间。
（1）从不动点迭代法与牛顿迭代法的结果对比可以看出，在精度1e-15高精度的情况下，方程很难在100次迭代中得到近似解。但是牛顿迭代法对于方程的收敛速度高于不动点迭代法。
（2）在两种不同的迭代格式下，从运行结果来看：第（1）种迭代格式收敛，但由于精度过高，达到最大迭代次数；第（2）种迭代格式是发散的，不论精度，都不会收敛。
（3）分别对两种迭代格式用Aitken法进行加速，结果可以看出两种迭代格式在1e-15的高精度下，均处于收敛的情况。不同的是第（1）种迭代格式在4此迭代后收敛，而第（2）种迭代格式却需要5次迭代才能收敛。但通过Aitken法加速之后，出现了以下两种情况：
1）由于精度过高引起的达到最大迭代次数的问题得以解决，函数收敛；
2）由于迭代格式的发散问题导致无法求解方程，在加速的情况下，也可以很快收敛，进一步说明了Aitken加速法的巨大作用。

6 参考：

[1]https://baike.baidu.com/item/%E7%89%9B%E9%A1%BF%E8%BF%AD%E4%BB%A3%E6%B3%95/10887580?fr=aladdin

[2]https://blog.csdn.net/weixin_45720246/article/details/115916645

[3]https://www.freesion.com/article/6284628111/

[4]https://blog.csdn.net/ReDreamme/article/details/111561395?utm_medium=distribute.pc_aggpage_search_result.none-task-blog-2_{aggregatepage}first_rank_ecpm_v1~rank_v31_ecpm-10-111561395.pc_agg_new_rank&utm_term=%E4%B8%8D%E5%8A%A8%E7%82%B9%E8%BF%AD%E4%BB%A3%E6%B3%95%E5%92%8CAitken%E6%B3%95&spm=1000.2123.3001.4430

[5]https://www.bilibili.com/video/BV1Lq4y1U7Hj?spm_id_from=333.999.0.0

[6]《数值分析》欧阳洁等。

python数值分析——非线性方程迭代求解与Aitken加速

1实验内容

2 算法原理

2.1 不动点迭代法

2.2 Aitken加速算法

2.3 牛顿迭代法原理

3 程序框图

4 编程实现

5 结果分析与讨论

6 参考：

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读