时隔一个月后水的一篇博客（

I(lnx)=∫12lnxdxI(ln x)=\int_{1}^{2}lnxdx，取ε=10−4,h=1\varepsilon=10^{-4},h=1，试用Romberg公式计算积分直到∣Rk,k−Rk−1,k−1∣<ε|R_{k,k}-R_{k-1,k-1}| \lt \varepsilon时停止，并做出Romberg函数积分表

代码如下（不保证正确，发现错误还请联系我，请勿无脑抄袭）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
import numpy as np
import math

t = np.empty(100)
s = np.empty(100)


def f(x):
    return math.log(x)


def romberg(a, b, e, func, stdAns):
    R0 = np.empty(100)
    R1 = np.empty(100)
    h = b-a
    R0[0] = (func(a)+func(b))*h/2
    t[1] = abs(R0[0]-stdAns)
    print("{:.15f}".format(R0[0]))
    temp = 1
    k = 1
    while True:
        k += 1
        temp *= 2
        xs = np.linspace(a, b, temp+1)
        iter = (func(x) for x in xs[1::2])
        R1[0] = (R0[0]+h*np.sum(np.fromiter(iter, float)))/2
        t[k] = abs(R1[0]-stdAns)
        print("{:.15f}".format(R1[0]), end=" ")
        h /= 2
        temp2 = 4.0
        for j in range(1, k):
            R1[j] = R1[j-1]+(R1[j-1]-R0[j-1])/(temp2-1)
            print("{:.15f}".format(R1[j]), end=" ")
            temp2 *= 4
        s[k] = abs(R1[1]-stdAns)
        print("")
        if abs(R1[k-1]-R0[k-2]) < e:
            return k
        R0 = R1.copy()


if __name__ == '__main__':
    stdAns = math.log(4)-1
    print("R:")
    k = romberg(1, 2, 1e-4, f, stdAns)
    print("T:")
    for i in range(1, k+1):
        print("N = {}: error = {:.15f}".format(
            int(math.pow(2, i-1)), t[i]), end="")
        if i > 1:
            order = math.log2(t[i-1]/t[i])
            print(" order = {:.15f}".format(order))
        else:
            print("")
    print("S:")
    for i in range(2, k+1):
        print("N = {}: error = {:.15f}".format(
            int(math.pow(2, i-1)), s[i]), end="")
        if i > 2:
            order = math.log2(s[i-1]/s[i])
            print(" order = {:.15f}".format(order))
        else:
            print("")