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)
上一篇:贪心算法
下一篇:Selenium常用方法