数值积分|中点法则(Midpoint Rule)

fem178

发布于 2020-07-09 11:14:08

10.4K0

发布于 2020-07-09 11:14:08

黎曼(Riemann)对定积分的定义是：积分区间划分为无数子区间，子区间内任意一点的函数值乘以子区间的长度得到一个矩形面积，然后将这些矩形面积累加起来可以得到积分值。中点法则(Midpoint Rule)是取子区间的中点的函数值作为矩形的高，如图所示

对于定积分

I = \int_1^2 {1\over x}dx

，将积分区间

[1,2]

划分成4个等长子区间

[1,{5\over 4}],[{5\over 4}, {3\over 2}],[{3\over 2}, {7\over 4}],[{7\over 4},2]

,每个子区间中点的函数值为

{8\over 9},{8\over 11},{8\over 13},{8\over 15}

. 可以得到定积分的近似值为：

I_4 = {1\over 4}[{8\over 9}+{8\over 11}+{8\over 13}+{8\over 15}]=0.69129189...

如果划分为8个子区间，可得到

I_8 = {1\over 8}[{16\over 17}+{16\over 19}+{16\over 21}+{16\over 23}+{16\over 25}+{16\over 27}+{16\over 29}+{16\over 31}]=0.69126605...

精确值是0.69314781... 子区间划分越多，误差就越小。

为了配合编程实现，可以通过将子区间的数量增加三倍来设计算法，而不是前面的辛普森，龙贝格公式增加两倍子区间。如图所示，空心圆圈表示必须计算的新值，而实心圆圈表示前一个迭代步计算的值。

M_0= h_0f(a+{h_0\over 2} ), h_0=b-a

M_1= h_1[f(a+{h_1\over 2})+f(a+{3h_1\over 2})+f(a+{5h_1\over 2} )], h_1=h_0/3

\begin{split} M_2&=h_2[f(a+{h_2\over 2})+f(a+{3h_2\over 2})+f(a+{5h_2\over 2} )+f(a+{7h_2\over 2} )+ \\ &f(a+{9h_2\over 2} )+f(a+{11h_2\over 2} )+f(a+{13h_2\over 2} )+f(a+{15h_2\over 2} )+f(a+{17h_2\over 2} )]\\ &\qquad h_2=h_1/3 \\ \end{split}

可得到递推式：

M_0= h_0f(a+{h_0\over 2} ), h_0=b-a， n_0=1

\begin{split} M_k & = {1\over 3} M_{k-1} + {1\over 3} {h_{k-1}\sum_{i=1}^{n_k-1}[f(a+(i-5/6)h_{k-1} )+f(a+(i-1/6)h_{k-1} )] }\\ &\quad h_k=h_{k-1}/3,n_k=3n_{k-1} ,k=1,2,... \\\end{split}

[算例]用中点法则求

I = \int_0^1 {4\over {1+x^2}}dx

import math


### Midpoint ruler求积分
### y = 4/( 1+x^2 )


def Func(x):
    return 4/( 1+pow(x,2) )
    
    
def MidPointIntegrate(Func, a, b, eps = 1e-6):
    kmax = 5000   # 最大迭代步数                                                       
    f1p6 = 1./6.
    f5p6 = 5./6.
    
    h = b-a
    n = 1
    
    m0 = h * Func(a+0.5*h)
    
    for k in range(1,kmax+1):
        sumf = 0
        for i in range(1,n+1): 
            sumf += Func(a+(i-f5p6)*h) + Func(a+(i-f1p6)*h)
            
        m = (m0 + h*sumf)/3
        if (math.fabs(m - m0) <= eps * math.fabs(m)): break
        
        h /= 3
        n *= 3
        
        m0 = m
        
    if (k >= kmax): print("已经达到最大迭代步数!")                                           
        
    return m
    
    
S = MidPointIntegrate(Func, 0, 1, eps = 1e-6)
print(S)

结果为

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2020-07-04，如有侵权请联系 cloudcommunity@tencent.com 删除

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

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

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

登录后参与评论

0 条评论

热度

数值积分|中点法则(Midpoint Rule)

数值积分|中点法则(Midpoint Rule)

社区

活动

资源

关于

腾讯云开发者

热门产品

热门推荐

更多推荐