前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >数值积分|第二类反常积分

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

作者头像
fem178
发布2020-07-20 16:15:59
2.3K0
发布2020-07-20 16:15:59
举报
文章被收录于专栏:数值分析与有限元编程

1 概述

第二类反常积分是值积分区间包含奇异点(singular points)。常规计算方法是将积分积分区间在奇异点内收,然后按照定积分来处理,再将计算结果取极限。如图1所示:

2 算法实现

python代码如下:

代码语言:javascript
复制

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)

输出结果:

第二类反常积分的数值算法大致思路就是在奇异点附近划分一个子区间,将这个子区间二等分,将其中之一积分,剩下的再二等分,将其中之一积分,如此下去,不断扩展积分区间,若扩展前后的积分的相对误差满足要求,则停止计算。

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2020-07-16,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 数值分析与有限元编程 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1 概述
  • 2 算法实现
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档