抛物线公式即辛卜生(Simpson)公式的数值积分的Python程序

代码

例子

下面的代码的意图是改变积分区域的抛物面块数，得到不同精度的积分值。
以对函数f(t) = 2000×ln(140000/(140000-2100×t))-9.8×t 进行积分为例，积分区域为[8,30]

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)  #真值

n = 10 #步长,就是将(a,b)区间分为多少个块
a = 8
b = 30
h = (b-a)/n
par_result = 0
for i in range(n):
    par_result += 1/6*h*(f(a+i*h)+4*f((a+i*h+a+(i+1)*h)/2)+f(a+(i+1)*h))  #抛物线积分算法(Simpson)
print(par_result)  #Simpson算法积分值

结果：

11061.3355350810
11061.3360331828

下面来探讨，积分区域内，需要划分多少块，抛物线积分算法计算出来的积分值才能与真值之间的截断误差小于0.1。

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()  #真值

a = 8
b = 30
n = 0 #步长,就是将(a,b)区间分为多少个块
while True:
    n += 1
    h = (b-a)/n
    par_result = 0
    for i in range(n):
        par_result += 1/6*h*(f(a+i*h)+4*f((a+i*h+a+(i+1)*h)/2)+f(a+(i+1)*h))  #抛物线积分算法(Simpson)
    R = abs(truth - par_result)
    if R <= 10**(-1):
        print(f'需要将积分区域划分{n}块,截断误差<0.1')
        break

结果：

需要将积分区域划分3块,截断误差<0.1

最后我们来看看真值和抛物线积分算法的积分值误差小于10^(-6)，需要将积分区域划分多少块？

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()  #真值

a = 8
b = 30
n = 0 #步长,就是将(a,b)区间分为多少个块
while True:
    n += 1
    h = (b-a)/n
    par_result = 0
    for i in range(n):
        par_result += 1/6*h*(f(a+i*h)+4*f((a+i*h+a+(i+1)*h)/2)+f(a+(i+1)*h))  #抛物线积分算法(Simpson)
    error = abs((truth - par_result) / truth)
    if error <= 10 ** (-6):
        print(f'需要将积分区域划分{n}块')
        break

结果：

需要将积分区域划分5块

对比一下梯形公式的数值积分梯形公式的数值积分的Python程序，Simpson积分算法简直完胜。