1. 程式人生 > >貝爾斯托法求根

貝爾斯托法求根

先上總結,緊跟程式碼,後上講解

總結

之前總結不完整,完整的以後敬上

 

Python程式碼(重修訂)

注:程式碼中的下標正好與上面推論中的相反

import cmath
import math


def f(x):
    return x ** 5 - 3.5 * x ** 4 + 2.75 * x ** 3 + 2.125 * x ** 2 - 3.875 * x + 1.25


rate = 0.01

r = s = -1

a = [1, -3.5, 2.75, 2.125, -3.875, 1.25]
b = []
c = []

def deal() :
    for i in range(len(b)) :
        a[i] = b[i]
    del a[len(a) - 2 : len(a)]

flag= 0

for k in range(9999):

    size = len(a)

    del b[0: len(b)]
    b.append(a[0])
    b.append(a[1] + r * b[0])

    for i in range(2, size):
        b.append(a[i] + r * b[i - 1] + s * b[i - 2])

    del c[0: len(c)]
    c.append(b[0])
    c.append(b[1] + r * c[0])

    for i in range(2, size - 1):
        c.append(b[i] + r * c[i - 1] + s * c[i - 2])


    R = (b[size - 2] * c[size - 3] - b[size - 1] * c[size - 4]) / (c[size - 2] * c[size - 4] - c[size - 3] * c[size - 3])
    S = (b[size - 1] * c[size - 3] - b[size - 2] * c[size - 2]) / (c[size - 2] * c[size - 4] - c[size - 3] * c[size - 3])

    r += R
    s += S

    er = math.fabs(R / r)
    es = math.fabs(S / s)




    if er <= rate and es <= rate:

        x = (r + cmath.sqrt(r * r + 4 * s)) / 2
        print(x)

        x = (r - cmath.sqrt(r * r + 4 * s)) / 2
        print(x)


        if flag == 1 :

            print(-b[1] / b[0])

            break

        deal()
        flag = 1


輸出結果

與答案正好匹配

 

講解