import numpy as np

# 被积因变量

def fun(x):

f = pow(x, 1.5)

return f

# 求梯形值

def T2n(a, b, n, Tn):

h = (b - a) / n # 步长

sum = 0.

for k in range(n):

sum += fun(a + (k + 0.5) * h)

T2n = Tn / 2. + sum * h / 2.

return T2n

# 求加快值

def romberg(max, a, b, epsilon): # max为计划最大度数,a、b为积分下、下限,epsilon为控制缺点

tm = np.zeros(max, dtype=float) # 第m行元素

tm1 = np.zeros(max, dtype=float) # 第m+1行元素

tm[0] = (b - a) * (fun(a) + fun(b)) / 2. # 初始值

print(tm)

k = 0

err = 1

while (err > epsilon and k < max - 1): # 当缺点小于预订缺点,或计划度数大于最大度数时遏止

n = 2 ** k

m = 1

tm1[0] = T2n(a, b, n, tm[0])

while (err > epsilon and m <= (k + 1)): # 当缺点小于预订缺点,或本行十足计划结束后遏止

tm1[m] = tm1[m - 1] + (tm1[m - 1] - tm[m - 1]) / (4. ** m - 1)

result = tm1[m]

err1 = abs(tm1[m] - tm[m - 1]) # 对角元素差

err2 = abs(tm1[m] - tm1[m - 1]) # 同业前后两元素差

err = min(err1, err2)

m += 1

tm = np.copy(tm1) # 下移一条龙

k += 1

print(tm)

return result

if __name__ == '__main__':

f1 = romberg(6, 0, 1, 1.e-10)

print(f1)