关于使用python完成龙贝格算法

Question

关于使用python完成龙贝格算法

a=0             # 积分下限
b=1             # 积分上限
eps=10**-5      # 精度
T=[]            # 复化梯形序列
S=[]            # Simpson序列
C=[]            # Cotes序列
R=[]            # Romberg序列
def f(x):    # 被积函数
    y=4/(1+x**2)
    return y
def Romberg(a,b,eps,f):
    h = b - a
    T.append(h * (f(a) + f(b)) / 2)     # 定义T[0]
    ep=eps+1
    m=0
    while(ep>=eps):
        m=m+1
        t=0                           # 每一次循环都会将t初始化
        for i in range(2**(m-1)-1):
            t=t+f(a+(2*(i+1)-1)*h/2**m)*h/2**m
        t=t+T[-1]/2                   #  T[-1]指的是数组T中倒数第一个元素
        T.append(t)
        if m>=1:
            S.append((4/3)*T[-1]-(1/3)*T[-2])
        if m>=2:
            C.append((16/15)*S[-1]-(1/15)*S[-2])
        if m>=3:
            R.append((64/63)*C[-1]-(1/63)*C[-2])
        if m>4:
            ep=abs(10*(R[-1]-R[-2]))
Romberg(a,b,eps,f)
print(T)
print(S)
print(C)
print(R)

print("积分结果为：{:.5f}".format(R[-1]))

运行结果
[3.0, 1.5, 1.6911764705882353, 2.1358026537831245, 2.506292558631274, 2.759621739265887, 2.9189057910568406, 3.0145118477143518, 3.07021171660127, 3.1019889340333107, 3.1198359194773158, 3.1297372868164106, 3.1351765796633675, 3.1381404486809648, 3.1397444739929528, 3.140607526927665, 3.1410695722537456, 3.141315854025983, 3.141446624386673, 3.1415158242843315, 3.1415523315867686, 3.1415715389135075, 3.141581619414385, 3.14158689808331, 3.141589656627041, 3.1415910955038555]
[1.0, 1.7549019607843137, 2.2840113815147545, 2.629789193580657, 2.844064799477424, 2.9720004749871585, 3.0463805332668548, 3.0887783395635755, 3.1125813398439908, 3.125784914625317, 3.1330377425961085, 3.1369896772790193, 3.139128405020164, 3.140279149096948, 3.140895211239236, 3.141223587362439, 3.1413979479500616, 3.1414902145069035, 3.141538890916884, 3.141564500687581, 3.141577941355754, 3.141584979581344, 3.141588657639619, 3.141590576141618, 3.14159157512946]
[1.8052287581699347, 2.3192853428967837, 2.652841047718384, 2.858349839870542, 2.9805295200211406, 3.0513392038188347, 3.0916048599833568, 3.1141682065293517, 3.1266651529440725, 3.133521264460828, 3.137253139591213, 3.139270986869574, 3.1403558653687336, 3.1409362820487217, 3.1412454791039854, 3.1414095719892363, 3.141496365610693, 3.141542136010883, 3.1415662080056275, 3.141578837400299, 3.1415854487963832, 3.141588902843504, 3.141590704041751, 3.1415916417286494]
[2.3274449712257814, 2.658135582715552, 2.8616118841904172, 2.982468880023531, 3.0524631670537183, 3.092243997382793, 3.114526354887224, 3.1268635171728776, 3.13363009162776, 3.1373123757043935, 3.1393030161914526, 3.1403730856623713, 3.140945495011896, 3.141250386993751, 3.1414121766382084, 3.1414977432872235, 3.141542862525171, 3.1415665901007817, 3.141579037866881, 3.141585553739178, 3.1415889576696485, 3.1415907326321992, 3.1415916566125683]
积分结果为：3.14159

疑惑：
运行结果应该是

为什么数组T[]，S[]，C[]，R[]，中会出现那么多其他结果

Answer 1

题主的问题恐怕不是数组元素多而是计算错误，比如，T数组的第2个元素应该是3.1，题主计算出来却是1.5。下面是我写的一个算法，请测试。

def romberg(f, a, b, eps):
    T = [(f(a)+f(b))/2]
    for i in range(4):
        k, m = pow(2, len(T)), 0
        for j in range(1, k, 2):
            m += f(j/k)
        T.append(T[-1]/2 + m/k)
    S = [4*T[i]/3 - T[i-1]/3 for i in range(1, len(T))]
    C = [16*S[i]/15 - S[i-1]/15 for i in range(1, len(S))]
    R = [64*C[i]/63 - C[i-1]/63 for i in range(1, len(C))]
    while abs(R[-1] - R[-2]) > eps:
        k, m = pow(2, len(T)), 0
        for j in range(1, k, 2):
            m += f(j/k)
        T.append(T[-1]/2 + m/k)
        S.append(4*T[-1]/3 - T[-2]/3)
        C.append(16*S[-1]/15 - S[-2]/15)
        R.append(64*C[-1]/63 - C[-2]/63)
    print('T = ', T)
    print('S = ', S)
    print('C = ', C)
    print('R = ', R)
    return R[-1]

romberg(lambda x:4/(1+x*x), 0, 1, 1e-5)
T =  [3.0, 3.1, 3.131176470588235, 3.138988494491089, 3.140941612041389]
S =  [3.1333333333333337, 3.1415686274509804, 3.141592502458707, 3.1415926512248227]
C =  [3.1421176470588232, 3.141594094125889, 3.1415926611425635]
R =  [3.141585783761874, 3.1415926383967965]
3.1415926383967965
romberg(lambda x:4/(1+x*x), 0, 1, 1e-10)
T =  [3.0, 3.1, 3.131176470588235, 3.138988494491089, 3.140941612041389, 3.1414298931749745, 3.141551963485656]
S =  [3.1333333333333337, 3.1415686274509804, 3.141592502458707, 3.1415926512248227, 3.141592653552836, 3.1415926535892167]
C =  [3.1421176470588232, 3.141594094125889, 3.1415926611425635, 3.1415926537080368, 3.1415926535916423]
R =  [3.141585783761874, 3.1415926383967965, 3.1415926535900285, 3.141592653589795]
3.141592653589795