代码
例子
下面的代码的意图是改变积分区域的梯形块数,得到不同精度的积分值。
以对函数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
tra_result = 0
for i in range(n):
tra_result += 1/2*h*(f(a+i*h)+f(a+(i+1)*h)) #梯形积分算法
print(tra_result) #梯形积分值
结果:
11061.3355350810
11069.5835420617
下面来探讨,积分区域内,需要划分多少块,梯形积分算法计算出来的积分值才能与真值之间的截断误差小于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
tra_result = 0
for i in range(n):
tra_result += 1/2*h*(f(a+i*h)+f(a+(i+1)*h)) #梯形积分算法
R = abs(truth - tra_result)
if R <= 10**(-1):
print(f'需要将积分区域划分{n}块,截断误差<0.1')
break
结果:
需要将积分区域划分91块,截断误差<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
tra_result = 0
for i in range(n):
tra_result += 1/2*h*(f(a+i*h)+f(a+(i+1)*h)) #梯形积分算法
error = abs((truth - tra_result)/truth)
if error <= 10**(-6):
print(f'需要将积分区域划分{n}块')
break
结果:
需要将积分区域划分274块
可以对比一下Simpson积分算法:抛物线公式即辛卜生(Simpson)公式的数值积分的Python程序