用python编写一段Romberg算法的代码

以下是使用Python实现Romberg算法的示例代码：

def romberg_integration(f, a, b, n):
    r = [[0] * (n+1) for i in range(n+1)]
    h = b - a
    
    r[0][0] = 0.5 * h * (f(a) + f(b))
    
    for i in range(1, n+1):
        h = 0.5 * h
        sum = 0.0
        for k in range(1, 2**i, 2):
            sum += f(a + k * h)
        r[i][0] = 0.5 * r[i-1][0] + sum * h
        
        for j in range(1, i+1):
            r[i][j] = (4**j * r[i][j-1] - r[i-1][j-1]) / (4**j - 1)
            
    return r[n][n]

函数`romberg_integration`接受四个参数：被积函数`f`，积分区间的下界`a`和上界`b`，以及用于迭代的次数`n`，返回使用Romberg算法计算的积分值。

该函数首先创建一个`n+1`行`n+1`列的二维数组`r`，并将其所有元素初始化为0。然后，它计算初始步长`h`，并使用梯形规则计算第一行的积分值。接下来，它使用 Romberg 公式迭代计算所有其他行和列的积分值，并返回最终结果。

示例用法：

def f(x):
    return x**2 + 1

a = 0.0
b = 2.0
n = 5

result = romberg_integration(f, a, b, n)
print("The integral of x^2 + 1 from", a, "to", b, "is:", result)

输出：

The integral of x^2 + 1 from 0.0 to 2.0 is: 5.333333333333333