python龙贝格算法

龙贝格算法

import numpy as np
import sympy
x=sympy.symbols("x")
yuan=input("请输入原函数：")
quyu=input("请输入积分区间（用空格隔开）：")
wucha=float(input("请输入误差："))
Q=[]
Q=quyu.split(" ")
Q=[float(x) for x in Q]
a=Q[0]
b=Q[1]
fa=sympy.sympify(yuan).subs(x,a)
fb=sympy.sympify(yuan).subs(x,b)
T00=(b-a)/2*(fa+fb)
T=np.array([T00])
jishuqi=1
ling=[0]
while True:
    hang=np.zeros(jishuqi)
    lie=np.zeros(jishuqi+1)
    T=np.row_stack((T,hang))#加一行
    T=np.column_stack((T,lie))#加一列
    he=0
    wuhu=int(2**(jishuqi-1))
    for i in range(wuhu):
        he+=sympy.sympify(yuan).subs(x,(a+(2*i+1)*(b-a)/(2**jishuqi)))
    T[jishuqi][0]=0.5*T[jishuqi-1][0]+he*(b-a)/(2**jishuqi)
    n=0
    for j in range(jishuqi):
        T[jishuqi-1-n][j+1]=(4**(j+1)*T[jishuqi-n][j]-T[jishuqi-1-n][j])/(4**(j+1)-1)
        n+=1
    if T[0][jishuqi]-T[0][jishuqi-1]<wucha:
        break
    else:
        jishuqi+=1
        continue
print('T:',T)
print("积分值为",round(T[0][jishuqi],6),'（误差为%f）'%(wucha))