数值分析-劈因子法（贝尔斯托法）python3实现

"""
贝尔斯托法：与牛顿法有松散关系的一种迭代方法，为求多项式的复根
公式：f5(x)=(x+1)(x-6)(x+7)(x-2)(x+3)
目的：求的该方程式的复根
具体方法：
①多项式迭代除去二次因式，x^2-r*x-s,可以得到复根
  为什么除以二次因式可以得到复根：

②求偏导数
③得到r和s以此改进
  计算误差 r和s
  误差均低于给定的终止条件时可以得到根
④结果可能：
    1.三阶或者以上以r,s为初始值应用
    2.商是二次可利用求根公式得到值
    3.商为一次多项式 由单根公式得到
补充：牛顿法-高阶方程求根存在复根情况
具体方法：
构造复变函数，x是复数
二元实值函数复合定义域以及函数值
对牛顿迭代公式的f(x)和F‘（x）除以因式
得最终表达式

"""
import sys
sys.setrecursionlimit(10000) #例如这里设置为一百万   递归深度不够
#RecursionError: maximum recursion depth exceeded in comparison
#1073741571错误是栈溢出，系统默认的栈空间大小是1MB
print('贝尔斯托法求解多项式f(x)=x^5-3.5*x^4+2.75*x^3+2.125*x^2-3.875*x+1.25的根，使用初始估计r=s=-1并迭代满足条件t=1%')
i = 0 # 循环标记量
#数组存储函数式得各个系数
a = [1.25,-3.875,2.125,2.75,-3.5,1]#依次存放a0,a1,a2...
b = ['b0', 'b1', 'b2', 'b3', 'b4', 'b5']
c = ['c0','c1','c2','c3','c4','c5']
rs=[-1,-1,'r','s']#上一次的rs和这一次rs
rs_0 = ['%r','%s']#误差分析
n = 5 #代表多项式的项数
x = ['x0','x1','x2','x3','x4']
n_0=[3,1]
#求B
def B(r,s):
    global i
    i = 3
    b[5]=round(a[5],5)#bn = an
    b[4]=round(a[4]+r*b[5],5)
    while i!=-1:
        b[i]=round(a[i]+r*b[i+1]+s*b[i+2],5)
        i=i-1
#求C
def C(r,s):
    global i
    i = 3
    c[5]=round(b[5],4)
    c[4]=round(b[4]+r*c[5],4)
    while i!=-1:
        c[i] = round(b[i]+r*c[i+1]+s*c[i+2],5)
        i = i-1

#计算当前r,s存入rs[]
def Analays():
    r = (b[0]*c[3]-b[1]*c[2])/(c[2]*c[2]-c[1]*c[3])+rs[0]
    s = (b[0]*c[2]-b[1]*c[1])/(c[3]*c[1]-c[2]*c[2])+rs[1]
    rs[2] = round(r, 4)
    rs[3] = round(s, 4)

#误差分析,存入rs_0[]
def Error():
    rs_0[0]=round((rs[2]-rs[0])/rs[2],5)
    rs_0[1]=round((rs[3]-rs[1])/rs[3],5)
    rs[0]=rs[2]
    rs[1]=rs[3]


w = 0#题设题解的下标
j = 0#存放答案的下标
count = 0#标记迭代次数
def Result(r,s):#贝尔斯托法计算结果
    global n,w,j,count
    count = count+1
    B(r, s)
    C(r, s)
    Analays()
    Error()
    if (abs(rs_0[0])<0.01)|(abs(rs_0[1])<0.01): #误差满足条件
        print('经过',count,'次迭代得')
        x[j] = (rs[2]+(rs[2]**2+4*rs[3])**(1/2))/2  # 4
        print('x',j,'=',x[j])
        j = j + 1
        x[j] = (rs[2]-(rs[2]**2+4*rs[3])**(1/2))/2  # 3
        print('x',j,'=',x[j])
        j =j+1
        n = n_0[w] # n代表当前计算多项式项数
        w=w+1
        if n == 1:
            x[4] = round((-rs[2]) / rs[3],0)  # 单项式求解
            print('x', 4, '=', x[4])
            print('多项式的根:',x)
        elif (n == 2):
            print('二次项求实根：此题目不需计算')
        elif (n > 2):#传入修正的r,s值  商作为因式迭代
            count = 0
            a[0] =b[2]
            a[1] = b[3]
            a[2] = b[4]
            a[3] = b[5]
            a[4] = 0
            a[5] = 0
            Result(rs[2], rs[3])
    else:#不满足近似条件
        Result(rs[2],rs[3])#用修正的值继续迭代

Result(-1,-1)#调用

结果：
贝尔斯托法求解多项式f(x)=x^5-3.5*x^4+2.75*x^3+2.125*x^2-3.875*x+1.25的根，使用初始估计r=s=-1并迭代满足条件t=1%
经过 4 次迭代得
x 0 = 0.5
x 1 = -1.0
经过 5 次迭代得
x 2 = (1+0.49939963956735095j)
x 3 = (1-0.49939963956735095j)
x 4 = 2.0
多项式的根: [0.5, -1.0, (1+0.49939963956735095j), (1-0.49939963956735095j), 2.0]

Process finished with exit code 0