龙贝格算法求数值积分的Python程序

Ronberg Integration

分成n等份时辛普森公式的值

from sympy import *

def f(t):
    f = 2000*log(140000/(140000-2100*t))-9.8*t
    return f

x = symbols('x')
truth = integrate(f(x),(x,8,30)).evalf()
print(truth) #真值

def T(n,a,b):
    h = (b-a)/n
    tra_result = 0
    for i in range(n):
        tra_result += 1 / 2 * h * (f(a + i * h) + f(a + (i + 1) * h))  # 梯形积分算法
    return tra_result

a = 8
b = 30
list = []
for k in range(4):
    m = 2**k
    list.append(T(m,a,b))
print(list)

while True:
    new_list = []
    for i in range(len(list)-1):
        R = list[i+1] + (list[i+1]-list[i])/3
        new_list.append(R)
    print(new_list)
    if len(new_list) == 1:
        print(new_list[0])
        break
    else:
        list = new_list

结果：

11061.3355350810
[11868.3481898411, 11266.3742932594, 11112.8206763693, 11074.2212976601]
[11065.7163277322, 11061.6361374059, 11061.3548380903]
[11060.2760739638, 11061.2610716518]
[11061.5894042144]
11061.5894042144

辛普森法的积分值Sn

from sympy import *

def f(t):
    f = 2000*log(140000/(140000-2100*t))-9.8*t
    return f

x = symbols('x')
truth = integrate(f(x),(x,8,30)).evalf()
print(truth) #真值

def T(n,a,b):
    h = (b-a)/n
    tra_result = 0
    for i in range(n):
        tra_result += 1 / 2 * h * (f(a + i * h) + f(a + (i + 1) * h))  # 梯形积分算法
    return tra_result

n = 0 #步长,就是将(a,b)区间分为多少个块
a = 8
b = 30
while True:
    n += 1
    S = T(2*n, a, b) + (T(2*n, a, b) - T(n, a, b)) / 3
    error = abs(truth - S) / truth
    if error <= 10 ** (-6):
        print('S为：', S)
        print('n为：', n)
        break

结果：

11061.3355350810
S为： 11061.3434684075
n为： 5

柯特斯法的积分值Cn

from sympy import *

def f(t):
    f = 2000*log(140000/(140000-2100*t))-9.8*t
    return f

x = symbols('x')
truth = integrate(f(x),(x,8,30)).evalf()
print(truth) #真值

def T(n,a,b):
    h = (b-a)/n
    tra_result = 0
    for i in range(n):
        tra_result += 1 / 2 * h * (f(a + i * h) + f(a + (i + 1) * h))  # 梯形积分算法
    return tra_result

def S(n,a,b):
    I = T(2*n, a, b) + (T(2*n, a, b) - T(n, a, b)) / (4**1-1)
    return I

n = 1 #步长,就是将(a,b)区间分为多少个块
a = 8
b = 30
C = S(2*n, a, b) + (S(2*n, a, b) - S(n, a, b)) / (4**2-1)
print(C)

结果：

11061.3355350810
11061.3641247175

from sympy import *

def f(t):
    f = 2000*log(140000/(140000-2100*t))-9.8*t
    return f

x = symbols('x')
truth = integrate(f(x),(x,8,30)).evalf()
print(truth) #真值

def T(n,a,b):
    h = (b-a)/n
    tra_result = 0
    for i in range(n):
        tra_result += 1 / 2 * h * (f(a + i * h) + f(a + (i + 1) * h))  # 梯形积分算法
    return tra_result

def S(n,a,b):
    I = T(2*n, a, b) + (T(2*n, a, b) - T(n, a, b)) / (4**1-1)
    return I

n = 0 #步长,就是将(a,b)区间分为多少个块
a = 8
b = 30
while True:
    n += 1
    C = S(2*n, a, b) + (S(2*n, a, b) - S(n, a, b)) / (4**2-1)
    error = abs(truth - C) / truth
    if error <= 10 ** (-6):
        print('C为：', C)
        print('n为：', n)
        break

结果：

11061.3355350810
C为： 11061.3360848026
n为： 2

龙贝格（Romberg）公式Rn

from sympy import *

def f(t):
    f = 2000*log(140000/(140000-2100*t))-9.8*t
    return f

x = symbols('x')
truth = integrate(f(x),(x,8,30)).evalf()
print(truth) #真值

def T(n,a,b):
    h = (b-a)/n
    tra_result = 0
    for i in range(n):
        tra_result += 1 / 2 * h * (f(a + i * h) + f(a + (i + 1) * h))  # 梯形积分算法
    return tra_result

def S(n,a,b):
    I = T(2*n, a, b) + (T(2*n, a, b) - T(n, a, b)) / (4**1-1)
    return I

def C(n,a,b):
    C = S(2*n, a, b) + (S(2*n, a, b) - S(n, a, b)) / (4**2-1)
    return C

n = 1 #步长,就是将(a,b)区间分为多少个块
a = 8
b = 30
R = C(2*n, a, b) + (C(2*n, a, b) - C(n, a, b)) / (4**3-1)
print(R)

结果：

11061.3355350810
11061.3356397246

from sympy import *

def f(t):
    f = 2000*log(140000/(140000-2100*t))-9.8*t
    return f

x = symbols('x')
truth = integrate(f(x),(x,8,30)).evalf()
print(truth) #真值

def T(n,a,b):
    h = (b-a)/n
    tra_result = 0
    for i in range(n):
        tra_result += 1 / 2 * h * (f(a + i * h) + f(a + (i + 1) * h))  # 梯形积分算法
    return tra_result

def S(n,a,b):
    I = T(2*n, a, b) + (T(2*n, a, b) - T(n, a, b)) / (4**1-1)
    return I

def C(n,a,b):
    C = S(2*n, a, b) + (S(2*n, a, b) - S(n, a, b)) / (4**2-1)
    return C

n = 0 #步长,就是将(a,b)区间分为多少个块
a = 8
b = 30
while True:
    n += 1
    R = C(2*n, a, b) + (C(2*n, a, b) - C(n, a, b)) / (4**3-1)
    error = abs(truth - R) / truth
    if error <= 10 ** (-6):
        print('R为：', R)
        print('n为：', n)
        break

结果：

11061.3355350810
R为： 11061.3356397246
n为： 1