第二类反常积分是值积分区间包含奇异点(singular points)。常规计算方法是将积分积分区间在奇异点内收,然后按照定积分来处理,再将计算结果取极限。如图1所示:
如图2所示:a为奇异点,在a附近划分一个子区间 ,长度为 ,将其一分为二,中点为 ,右半部分参与积分计算,再将左半部分 二等分,中点为 ,右半部分参与积分计算,如此下去,直到误差满足要求。
python代码如下:
import math
### 第二类反常积分数值分析
### y = 1/sqrt(x)
### 积分区间(0, 1]
def Func(x):
return 1/ math.sqrt(x)
def Improp2(Func, a, b, eps = 1e-6):
###
### a为区间的左端点,是奇异点
###子区间积分时,还要调用自适应梯形公式,这里可以任选方法。比如辛普森公式
h0 = 0.1e0 * (b-a) #子区间长度
s = 0e0
h = h0
s1 = 1e0
while(math.fabs(s1) > eps * math.fabs(s)):
h *= 0.5 # 半区间
x = a + h
s1 = AdaptiveTrapzCtrl(Func,x,x+h,eps)
s += s1
a += h0
s += AdaptiveTrapzCtrl(Func, a, b, eps )
return s
### 自适应梯形公式求积分
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 = Improp2(Func, 0, 1, eps = 1e-6)
print(s)
第二类反常积分的数值算法大致思路就是在奇异点附近划分一个子区间,将这个子区间二等分,将其中之一积分,剩下的再二等分,将其中之一积分,如此下去,不断扩展积分区间,若扩展前后的积分的相对误差满足要求,则停止计算。
★★★★往期相关★★★★