龙贝格积分公式（数值积分）

问题描述

利用龙贝格积分公式计算函数f(x)=(x^2+x+1)cos(x)，在区间[0, pi/2]范围内的定积分近似值。

输入形式

在屏幕上龙贝格积分表行的最大值。

输出形式

龙贝格积分表。

样例输入

6

样例输出

[[0.78539816 0. 0. 0. ]

[1.72681266 2.04061749 0. 0. ]

[1.96053417 2.03844134 2.03829626 0. ]

[2.01879395 2.03821388 2.03819871 2.03819716]

[2.03334734 2.03819847 2.03819745 2.03819743]

[2.03698495 2.03819749 2.03819743 2.03819743]]

样例说明

输入：龙贝格积分表行的最大值n为6。输出：6行4列的龙贝格积分表。

代码

import numpy as np
import math


def Input():
    n = input()
    return 0, math.pi / 2, int(n)


def f(x):
    return (x ** 2 + x + 1) * math.cos(x)


def romber(a, b, n, tol):
    m = 1
    h = b - a
    err = 1
    j = 1
    R = np.zeros((n, n), dtype=np.double)
    R[0, 0] = h * (f(a) + f(b)) / 2
    while err > tol and j < n:
        h /= 2
        s = 0
        for p in range(1, m+1):
            x = a + h * (2 * p - 1)
            s += f(x)
        R[j, 0] = R[j-1, 0] / 2 + h * s
        m *= 2
        for k in range(1, j+1):
            R[j, k] = R[j, k-1] + (R[j, k-1] - R[j-1, k-1]) / (4 ** k - 1)

        err = abs(R[j-1, j-1]-R[j, j])
        j += 1

    return R


def out(x):
    print(x)


def main():
    a, b, n = Input()
    tol = np.finfo(np.float32).eps
    R = romber(a, b, n, tol)
    out(R[:, :4])


if __name__ == '__main__':
    main()