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

1实验内容

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

2 算法原理

2.1 不动点迭代法

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

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

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

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

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

迭代公式收敛指迭代序列{xk}收敛,迭代公式发散指迭代序列{xk}不收敛,即发散。迭代公式不一定总是收敛。简单迭代法的几何意义如下:

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

(a)收敛

(b)发散

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

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

2.2 Aitken加速算法

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

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

2.3 牛顿迭代法原理

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

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

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

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

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

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

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


图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

image.png

image.png

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
image.png
image.png
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
image.png
image.png
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-2aggregatepagefirst_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]《数值分析》欧阳洁等。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 204,684评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,143评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,214评论 0 337
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,788评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,796评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,665评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,027评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,679评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 41,346评论 1 299
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,664评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,766评论 1 331
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,412评论 4 321
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,015评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,974评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,203评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,073评论 2 350
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,501评论 2 343

推荐阅读更多精彩内容