HFUT计算方法实验：龙贝格算法

import math

def Composite_Trapezoidal(f, a, b, n):
    h = (b - a) / n
    s = f(a) + f(b)
    for i in range(1, n):
        s += 2 * f(a + i * h)
    return s * h / 2

def Composite_Simpson(f, a, b, n):
    h = (b - a) / n
    s = f(a) + f(b)
    for i in range(1, n+1):
        s += 4 * f(a + (2 * i - 1) * h / 2)
    for i in range(1, n):
        s += 2 * f(a + i * h)
    return s * h / 6

def T2n(f, a, b, n):
    h = (b - a) / n
    s=Composite_Trapezoidal(f, a, b, n)*0.500000
    for i in range(1, n+1):
        s+=f(a + (2 * i - 1) * h / 2)*h/2
    return s

def T4n(f, a, b, n):
    h = (b - a) / n
    s =T2n(f, a, b, n)*0.500000
    for i in range(1, 2*n+1):
        s += f(a + (2 * i - 1) * h / 4)*h/4
    return s

def T8n(f, a, b, n):
    h = (b - a) / n
    s =T4n(f, a, b, n)*0.500000
    for i in range(1, 4*n+1):
        s += f(a + (2 * i - 1) * h / 8)*h/8
    return s

def S2n(f, a, b, n):
    return T4n(f, a, b, n)*4.000000/3.000000-T2n(f, a, b, n)*1.000000/3.000000

def S4n(f, a, b, n):
    return T8n(f, a, b, n)*4.000000/3.000000-T4n(f, a, b, n)*1.000000/3.000000

def Cn(f, a, b, n):
    return S2n(f, a, b, n)*16.000000/15.000000-Composite_Simpson(f, a, b, n)*1.000000/15.000000

def C2n(f, a, b, n):
    return S4n(f, a, b, n)*16.000000/15.000000-S2n(f, a, b, n)*1.000000/15.000000

def Rn(f, a, b, n):
    return C2n(f, a, b, n)*64.000000/63.000000-Cn(f, a, b, n)*1.000000/63.000000

def fun(x):
    return float(math.sin(x)/x)

if __name__ == '__main__':
    print("f=sin(x)/x  a=0  b=1  n=2")
    print("Tn\t",Composite_Trapezoidal(fun, 0.000000000001, 1, 2))
    print("T2n\t",T2n(fun, 0.000000000001, 1, 2))
    print("T4n\t",T4n(fun, 0.000000000001, 1, 2))
    print("T8n\t",T8n(fun, 0.000000000001, 1, 2))
    print("Sn\t",Composite_Simpson(fun, 0.000000000001, 1, 2))
    print("S2n\t",S2n(fun, 0.000000000001, 1, 2))
    print("S4n\t",S4n(fun, 0.000000000001, 1, 2))
    print("Cn\t",Cn(fun, 0.000000000001, 1, 2))
    print("C2n\t",C2n(fun, 0.000000000001, 1, 2))
    print("Rn\t",Rn(fun, 0.000000000001, 1, 2))
    m=Rn(fun, 0.000000000001, 1, 2)
    print("积分结果为：",'%.9f'%m)

---

代码用于交流学习，如有不足欢迎斧正！