无穷区间的积分又称第一类反常积分。常规计算方法是将积分上限
视为常数,然后按照定积分来处理,再将计算结果取极限。如图1所示:
2 算法实现
第一类反常积分的数值算法大致思路就是不断扩展积分区间,若扩展前后的积分的相对误差满足要求,则停止计算。
如图2所示,计算反常积分
时,先计算
,再计算
,然后计算
,
若
的相对误差满足要求,则停止计算。
python代码如下:
import math
### 第一类反常积分(无穷区间)数值分析
### y = 1/( x^2 )
### 积分区间[1,+inf)
def Func(x):
return 1/ pow(x,2)
def Improp1(Func, a, inf, eps = 1e-6):
### [a,+inf) (inf > 0时) 或者 (-inf,a](inf < 0),
### inf为积分收敛时区间的右(左)端点
### 子区间积分时,还要调用自适应梯形公式,这里可以任选方法。比如辛普森公式
h = 1e0 #子区间长度
x = a
if (inf < 0): h = -h
kstep = 0
s1 = 1e0
s = 0
while(math.fabs(s1) > eps*math.fabs(s) or math.fabs(Func(x)) > eps):
s1 = AdaptiveTrapzCtrl(Func,x,x+math.fabs(h),eps) #计算区间[x,x+h]积分
s += s1 #扩展区间
x += h
kstep += kstep
if(kstep > 1e7): break #若迭代次数过多,则停止
inf = x
return (s, inf)
### 自适应梯形公式求积分
def AdaptiveTrapzCtrl(Func, a, b, eps = 1e-6):
kmax = 9000 #最大迭代步数
h = b-a # 积分区间
n = 1
t0 = 0.5*h*(Func(a) + Func(b))
for k in range(1,kmax+1):
sumf = 0
for i in range(1,n+1):
sumf += Func(a+(i-0.5)*h)
t = 0.5*(t0 + h*sumf)
if (k > 1):
if (math.fabs(t-t0) <= eps*math.fabs(t)): break
if (math.fabs(t) <= eps and math.fabs(t) <= math.fabs(t-t0)): break
h *= 0.5
n *= 2
t0 = t
if (k >= kmax): print("已经达到最大迭代步数!")
return t
s = Improp1(Func, 1, 10, eps = 1e-6)
print(s)
输出结果
PS:无穷区间的长度并不是有几十万甚至几十亿那么长。本例的无穷区间在收敛时为1001.