python數(shù)值分析——非線性方程迭代求解與Aitken加速

1實驗內(nèi)容

用迭代法求下列方程的根
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法進(jìn)行加速，并對結(jié)果進(jìn)行分析。

2 算法原理

2.1 不動點迭代法

迭代法是一種逐次逼近的方法，用某種固定格式反復(fù)校正根的近似值，使之逐步精確，最后達(dá)到滿足精度要求的結(jié)果。

方程f(x)=0化為等價形式的方程x=φ(x) ，構(gòu)造迭代公式xk+1=φ(xk), k=0,1,2,……取初始近似根(迭代初值)x0，進(jìn)行迭代計算x1=φ(x0),x2=φ(x1),...,xk+1=φdf(xk),……

則有x1，x2 ,……,xk,……，得到迭代序列{xk}.如果這個序列有極限，則迭代公式是收斂的。

這時則x(x)，x稱為函數(shù)φ(x)的不動點，等價地有f(x)=0,x*即為方程的根。連續(xù)函數(shù)φ(x)稱為迭代函數(shù)。

實際計算到|xk-xk-1|<ε(ε是預(yù)定的精度)，取x*≈xk。

迭代公式收斂指迭代序列{xk}收斂，迭代公式發(fā)散指迭代序列{xk}不收斂，即發(fā)散。迭代公式不一定總是收斂。簡單迭代法的幾何意義如下：

圖1 不動點迭代的幾何表示

（a）收斂

（b）發(fā)散

圖2 收斂與發(fā)散的幾何對比

??=??(??) 的解稱為函數(shù) ??(??) 的不動點，幾何圖像為函數(shù)圖像與正比例函數(shù) ??=??的交點。迭代形式為：????+1=??(????) 。設(shè)函數(shù)的不動點為 ???，當(dāng)|??′(???)|<1 時，在不動點附近的存在一個領(lǐng)域，使得從領(lǐng)域內(nèi)的某個值開始的不動點迭代收斂，收斂到不動點。通過 ??????+1≈??′(???)??????式，當(dāng)導(dǎo)數(shù)絕對值小于1時，迭代后誤差的線性主部較迭代前誤差小，以導(dǎo)數(shù)絕對值為收斂常數(shù)線性收斂。唯一不同的是若導(dǎo)數(shù)值為零，此時該迭代至少二次收斂。用不動點迭代解非線性方程的核心在于將非線性方程轉(zhuǎn)化為一個收斂的不動點迭代。

2.2 Aitken加速算法

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

對于收斂的迭代過程，只要迭代足夠多次，就可以使結(jié)果達(dá)到任意的精度。但有時迭代過程收斂緩慢，從而使計算量變得很大，因此，迭代過程的加速是個重要的過程。

2.3 牛頓迭代法原理

牛頓迭代法是用于求解等式方程的一種方法。類似于求解F(x)=0的根，牛頓迭代法求解的是近似根，這個想法準(zhǔn)確來說來源于泰勒展開式，需要求解的表達(dá)式可能非常復(fù)雜，通過一般的方法，我們很難求出它的解。

所以采用了一種近似求解的方法，取泰勒展開式的前幾項，隊原來的求解函數(shù)做一個取代，然后，求解這個取代原方程的方程的解，作為近似解。當(dāng)然只對原方程做一次近似求解不行，因為第一次近似肯定不會太準(zhǔn)確，所以還需要不斷地迭代。

首先就要取一個值作為初始的近似值，然后去求解該點的泰勒展開近似項，然后求解根，之后再以此根對原方程進(jìn)行近似，然后再求解結(jié)果不斷重復(fù)迭代，最終就能求得近似解。

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

設(shè)x是f(x) = 0的根，選取x0作為x初始近似值，過點（x0, f(x)）做曲線y = f(x)的切線L，則L的方程為y = f(x0)+f '(x0)(x-x0)，求出L與x軸交點的橫坐標(biāo)x1= x0-f(x0) / f '(x0），稱x1為x的一次近似值。過點（x1, f(x1)）做曲線y = f(x)的切線，并求該切線與x軸交點的橫坐標(biāo)x2= x1-f(x1)/ f '(x1)，稱x2為x的二次近似值。重復(fù)以上過程，得r的近似值序列，其中xn+1=x(nf(xn)/f'(xn)，稱為x的n+1次近似值，上式稱為牛頓迭代公式。

解非線性方程f(x)=0的牛頓法是把非線性方程線性化的一種近似方法。把f(x）在x0點附近展開成泰勒級數(shù)f(x) = f(x0)+(x－x0)f'(x0)+(x－x0)2*f''(x0)/2! +… 取其線性部分（一次），作為非線性方程f(x) = 0的近似方程，即泰勒展開的前兩項，則有f(x0)+f'(x0)(x－x0)=0 設(shè)f'(x0）≠0則其解為x1=x0－f(x0)/f'(x0) 這樣，得到牛頓法的一個迭代序列：x(n+1）= xn－f(xn) / f '(xn)。

牛頓迭代法，取得是泰勒展開式的前兩項，也就是線性近似，所以迭代比較快，容易計算。


圖3 牛頓迭代法幾何表示

3 程序框圖

圖4 Aitken加速法算法框圖

圖5 牛頓迭代法算法框圖

4 編程實現(xiàn)

利用所編程序，求解下列方程的根:
以下將用Python進(jìn)行編程實驗，利用不動點迭代法與牛頓迭代法求解題目中方程的根，并且利用Aitken加速方法對其進(jìn)行加速，最后進(jìn)行結(jié)果分析與討論。
利用不動點迭代法的Python代碼詳解及運行結(jié)果：

（1） 把方程改寫成等價的迭代格式: xk+1+20/(xk2+2xk+10)，取初值x0=1，求方程的根。
第（1）種迭代格式代碼詳解：

from sympy import *


def fun():
    """
    求解方程的符號定義
    :return:方程
    error:這里不能用math庫中的函數(shù)，pow()需要調(diào)用符號庫中的函數(shù)，即sympy.pow()
    或者用from sympy import * 來導(dǎo)入sympy中的所有函數(shù)，之后直接寫pow()函數(shù)即
    是默認(rèn)的sympy中的函數(shù)。
    """
    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  # 返回函數(shù)的值


def fixedPoint(x0, eps, maxIter):
    """
    不動點迭代法的作用是，求解非線性方程的根，采用逐步逼近的方法進(jìn)行計算
    :param x0:迭代初值
    :param eps:誤差精度要求
    :param maxIter:最大迭代次數(shù)
    :return:返回值為None
    """
    x = symbols('x')
    fh = fun()
    x_n = x0  # 定義初值
    k = 0  # 初始化迭代次數(shù)
    errval = 0  # 初始化誤差
    print('%3s %10s %18s' % ('迭代次數(shù)', '方程近似值', '迭代誤差'))
    for k in range(maxIter):
        x_b = x_n  # 代表x_n
        x_n = 20 / (pow(x_b, 2) + 2 * x_b + 10)  # 寫出第一種迭代函數(shù)格式
        errval = abs(fh.evalf(subs={x: x_n}))  # 第k次迭代誤差的大小
        print('%3d %22.15f %22.15f' % (k + 1, x_n, errval))  # 分別輸出迭代次數(shù)，方程的近似根以及迭代誤差
        if errval <= eps:
            break
    if k + 1 <= maxIter - 1:
        print('方程在滿足精度' + str(eps) + '的條件下，近似解為：'
              + str(x_n) + ',誤差為：' + str(errval))
    else:
        print('不動點迭代法求解方程的根，已經(jīng)達(dá)到最大迭代次數(shù)，也可能不收斂或精度過高...')
    return None


if __name__ == '__main__':
    fh = fun()
    plot(fh)
    x0 = float(input('請輸入迭代初值：'))  # input函數(shù)總是以字符串的形式返回
    eps = float(input('請輸入誤差精度要求：'))  # 方程解的精度要求是近似解與真值之間的誤差
    maxIter = int(input('請輸入最大迭代次數(shù)：'))  # 方程一般會迭代無數(shù)次，必須定義其迭代的次數(shù)，以求收斂
    fixedPoint(x0, eps, maxIter)
    print('方程為：', '%30s' % (str(fun())))
    print('迭代函數(shù)為：', '%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庫中的函數(shù)，pow()需要調(diào)用符號庫中的函數(shù)，即sympy.pow()
    或者用from sympy import * 來導(dǎo)入sympy中的所有函數(shù)，之后直接寫pow()函數(shù)即
    是默認(rèn)的sympy中的函數(shù)。
    """
    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  # 返回函數(shù)的值


def fixedPoint(x0, eps, maxIter):
    """
    不動點迭代法的作用是，求解非線性方程的根，采用逐步逼近的方法進(jìn)行計算
    :param x0:迭代初值
    :param eps:誤差精度要求
    :param maxIter:最大迭代次數(shù)
    :return:返回值為None
    """
    x = symbols('x')
    fh = fun()
    x_n = x0  # 定義初值
    k = 0  # 初始化迭代次數(shù)
    errval = 0  # 初始化誤差
    print('%3s %10s %18s' % ('迭代次數(shù)', '方程近似值', '迭代誤差'))
    for k in range(maxIter):
        x_b = x_n  # 代表x_n
        x_n = 2 - (pow(x_b, 3) + 2 * pow(x_b, 2))/10  # 寫出第二種迭代函數(shù)格式
        errval = abs(fh.evalf(subs={x: x_n}))  # 第k次迭代誤差的大小
        print('%3d %22.15f %22.15f' % (k + 1, x_n, errval))  # 分別輸出迭代次數(shù)，方程的近似根以及迭代誤差
        if errval <= eps:
            break
    if k + 1 <= maxIter - 1:
        print('方程在滿足精度' + str(eps) + '的條件下，近似解為：'
              + str(x_n) + ',誤差為：' + str(errval))
    else:
        print('不動點迭代法求解方程的根，已經(jīng)達(dá)到最大迭代次數(shù)，也可能不收斂或精度過高...')
    return None


if __name__ == '__main__':
    fh = fun()
    plot(fh)
    x0 = float(input('請輸入迭代初值：'))  # input函數(shù)總是以字符串的形式返回
    eps = float(input('請輸入誤差精度要求：'))  # 方程解的精度要求是近似解與真值之間的誤差
    maxIter = int(input('請輸入最大迭代次數(shù)：'))  # 方程一般會迭代無數(shù)次，必須定義其迭代的次數(shù)，以求收斂
    fixedPoint(x0, eps, maxIter)
    print('方程為：', '%30s' % (str(fun())))
    print('迭代函數(shù)為：', '%20s' % ('2 - (x**3 + 2 * x**2)/10'))
image.png
image.png
image.png
image.png

（3） 對迭代格式(1)和(2)分別用Aitken法進(jìn)行加速，并對結(jié)果進(jìn)行分析。
迭代格式
第（1）種迭代格式Aitken加速法代碼詳解：

from sympy import *


def Aitken(x0, eps, iterNum, Phi):
    xk_1 = x0
    errval = 0
    print('%3s %10s %18s' % ('迭代次數(shù)', '方程近似值', '迭代誤差'))
    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函數(shù)總是以字符串形式返回
    eps = float(input('請輸入迭代誤差精度要求：'))  # 方程誤差精度要求為迭代近似值與真值之間的差值
    iterNum = int(input('請輸入最大迭代次數(shù)：'))  # 最大迭代次數(shù)限制了方程的收斂
    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' % ('迭代次數(shù)', '方程近似值', '迭代誤差'))
    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函數(shù)總是以字符串形式返回
    eps = float(input('請輸入迭代誤差精度要求：'))  # 方程誤差精度要求為迭代近似值與真值之間的差值
    iterNum = int(input('請輸入最大迭代次數(shù)：'))  # 最大迭代次數(shù)限制了方程的收斂
    Phi(x)
    Aitken(x0, eps, iterNum, Phi)
image.png

牛頓迭代法代碼詳解：

from sympy import *


def fun():
    """
    求解方程的符號定義
    :return:方程
    error:這里不能用math庫中的函數(shù)，pow()需要調(diào)用符號庫中的函數(shù)，即sympy.pow()
    或者用from sympy import * 來導(dǎo)入sympy中的所有函數(shù)，之后直接寫pow()函數(shù)即
    是默認(rèn)的sympy中的函數(shù)。
    """
    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  # 返回函數(shù)的值


def diffFun():
    """
    求解方程的一階導(dǎo)數(shù)
    :return: 一階導(dǎo)數(shù)
    """
    return fun().diff()


def newtonIterative(x0, eps, maxIter):
    """
    牛頓迭代函數(shù)的作用，求解非線性方程的根，采用逼近法的思想
    :param x0:迭代初值
    :param eps:誤差精度要求
    :param maxIter:最大迭代次數(shù)
    :return:返回值為None
    """
    # x1 = x0 - f(x0)/f'(x0)==>x1 = x2 - f(x1)/f'(x1)...
    x = symbols('x')
    fh = fun()  # 引用方程
    dfh = diffFun()  # 引用方程的一階導(dǎo)數(shù)
    x_n = x0
    k = 0
    errval = 0
    print('%3s %10s %18s' % ('迭代次數(shù)', '方程近似值', '迭代誤差'))
    # 利用牛頓迭代法的思想逐步逼近精確解
    for k in range(maxIter):
        x_b = x_n  # 代表x(n)
        fx = fh.evalf(subs={x: x_b})  # 方程在x_n處的數(shù)值
        dfx = dfh.evalf(subs={x: x_b})  # 一階導(dǎo)數(shù)方程在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('牛頓迭代法求解數(shù)值逼近，已經(jīng)達(dá)到最大迭代次數(shù)，也可能不收斂或精度過高...')
    # print(x0, eps, maxIter)
    return None


if __name__ == '__main__':
    fh = fun()
    dfh = diffFun()
    plot(fh)
    plot(dfh)
    x0 = float(input('請輸入迭代初值：'))  # input函數(shù)總是以字符串的形式返回
    eps = float(input('請輸入誤差精度要求：'))  # 方程解的精度要求是近似解與真值之間的誤差
    maxIter = int(input('請輸入最大迭代次數(shù)：'))  # 方程一般會迭代無數(shù)次，必須定義其迭代的次數(shù)，以求收斂
    newtonIterative(x0, eps, maxIter)
    print('方程為：', '%30s' % (str(fun())))
    print('方程的導(dǎo)數(shù)為：' '%22s' % (str(diffFun())))
image.png
image.png
image.png
image.png

5 結(jié)果分析與討論

圖6 方程的曲線

由方程的曲線可以看出，此方程的根在1附近只有一個根，因此迭代初值取1可以降低迭代次數(shù)，節(jié)省迭代時間。
（1）從不動點迭代法與牛頓迭代法的結(jié)果對比可以看出，在精度1e-15高精度的情況下，方程很難在100次迭代中得到近似解。但是牛頓迭代法對于方程的收斂速度高于不動點迭代法。
（2）在兩種不同的迭代格式下，從運行結(jié)果來看：第（1）種迭代格式收斂，但由于精度過高，達(dá)到最大迭代次數(shù)；第（2）種迭代格式是發(fā)散的，不論精度，都不會收斂。
（3）分別對兩種迭代格式用Aitken法進(jìn)行加速，結(jié)果可以看出兩種迭代格式在1e-15的高精度下，均處于收斂的情況。不同的是第（1）種迭代格式在4此迭代后收斂，而第（2）種迭代格式卻需要5次迭代才能收斂。但通過Aitken法加速之后，出現(xiàn)了以下兩種情況：
1）由于精度過高引起的達(dá)到最大迭代次數(shù)的問題得以解決，函數(shù)收斂；
2）由于迭代格式的發(fā)散問題導(dǎo)致無法求解方程，在加速的情況下，也可以很快收斂，進(jìn)一步說明了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]《數(shù)值分析》歐陽潔等。

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容