数值积分|第一类反常积分

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)