变步长梯形公式数值积分的Python程序

基本代码：

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 = 1 #步长,就是将(a,b)区间分为多少个块
a = 8
b = 30
#快速迭代
I = T(2*n,a,b) + (T(2*n,a,b) - T(n,a,b))/3
print(I)

结果：

11061.3355350810
11065.7163277322

迭代至特定误差

让我们看看变步长梯形公式比普通梯形公式快多少，相对误差小于10^-6，迭代停止。
普通梯形公式的效果见我之前的文章
梯形公式的数值积分的Python程序

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
    #快速迭代
    I = T(2*n,a,b) + (T(2*n,a,b) - T(n,a,b))/3
    error = abs(truth-I)/truth
    if error <= 10**(-6):
        print('I为：',I)
        print('n为：',n)
        break

结果：

11061.3355350810
I为： 11061.3434684075
n为： 5