数值积分|自适应辛普森积分公式

fem178

发布于 2020-06-09 03:43:55

4K00

代码可运行

文章被收录于专栏：数值分析与有限元编程数值分析与有限元编程

运行总次数：0

代码可运行

在

数值积分| 辛普森公式

提到，辛普森积分最简单的形式是

\int_a^b f(x) dx \approx \frac{h}{3}(y_0+4y_1+y_2)

也就是说至少要三个积分点，两个积分子区间。所以，自适应辛普森积分公式要从S1起步,即

S_1=\frac{h_1}{3}[f(a)+4f(a+h_1)+f(b)] \quad (1)

将

(1)

式与自适应梯形公式

T_0=h[\frac{f(a)}{2}+\frac{f(b)}{2}]

T_1=h_1[\frac{f(a)}{2}+f(a+h_1)+\frac{f(b)}{2}]

比较，可得

S_1=\frac{4T_1-T_0}{3}\quad (2)

S_2=\frac{h_2}{3}[f(a)+4f(a+h_2)+2f(a+2h_2)+4f(a+3h_2)+f(b)]=\frac{4T_2-T_1}{3} \quad (3)

由此可以得到递推式

S_k=\frac{4T_k-T_{k-1}}{3}\qquad (4)

若以

\epsilon

表示前后两次计算结果的相对误差，即

|S_k-s_{k-1}| < \epsilon|S_k| \qquad (5)

若满足要求，则停止计算。python代码


import math
###自适应辛普森公式求积分
### y = 1/( 1+x^2 )


def Func(x):
    return 1/( 1+pow(x,2) )
    

def AdaptiveSimpsonCtrl(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)
        
        s = (4*t - t0)/3
        
        if (k > 1):
            if (math.fabs(s-s0) <= eps * math.fabs(s)): break
            if (math.fabs(s) <= eps and math.fabs(s) <= math.fabs(s-s0)): break
            
        h *= 0.5
        n *= 2
        s0 = s
        t0 = t
        
    if (k >= kmax): print("已经达到最大迭代步数!")
    
    return s
    

S = AdaptiveSimpsonCtrl(Func, 0, 1, eps = 1e-6)
print(S)