龙贝格求积算法python实现
import numpy as np
def trapezoid(a, b, n, func):
"""
复化梯形公式求函数func在区间[a,b]上的积分值
n是等分的区间数目
"""
x = np.linspace(a, b, num=n + 1)
y = func(x)
h = (b - a) / (2 * n)
return h * (y[0] + 2 * np.sum(y[1:-1]) + y[-1])
def romberg(a, b, func, eps, max_iter=100):
"""
龙贝格求积算法
求函数func在区间[a,b]上的积分值
eps:误差
max_iter:最大迭代上限，应该大于等于2
"""
previous = [trapezoid(a, b, 1, func)]
for i in range(1, max_iter):
current = [trapezoid(a, b, 2 ** i, func)]
for k in range(1, i + 1):
tmp = ((4 ** k) * current[k - 1] - previous[k - 1]) / (4 ** k - 1)
current.append(tmp)
if np.abs(current[-1] - previous[-1]) < eps:
return current[-1]
previous = current
return previous[-1]
if __name__ == '__main__':
# x^2+x^3+exp(x)在[0,2]上的积分值：13.055722765711728
print(
romberg(0, 2, lambda x: x ** 2 + x ** 3 + np.exp(x), 1e-5)
)
