使用for循环的辛普森规则(数值积分)

问题描述 投票:1回答:1

[我正在尝试使用for循环在python中编写Simpson规则,但我不断收到断言错误,无法找到原因。

def integrate_numeric(xmin, xmax, N):
    ''' 
    Numerical integral of f from xmin to xmax using Simpson's rule with 
        N panels.
    '''
    xsum = 0
    msum = 0
    h = (xmax-xmin)//N

    for i in range(0, N):
        xsum += f(xmin + i*h)
        print (xsum)

    for i in range(0,N-1):
        msum += f(xmin + (h/2) + i*h)    
        print (msum)

    I = (h/6) * (f(xmin) + 4*(msum) + 2*(xsum) + f(xmax))
    return I

f:

def f(x):
    '''Function equivalent to x^2 sin(x).'''
    return (x**2) * numpy.sin(x)

样本:

assert numpy.isclose(integrate_numeric(xmin=0, xmax=4, N=50), 1.096591)
python python-3.x loops numerical-methods simpsons-rule
1个回答
1
投票

这是您的代码的固定版本:

import numpy

def integrate_numeric(xmin, xmax, N):
    '''                                                                                                                               
    Numerical integral of f from xmin to xmax using Simpson's rule with                                                               
        N panels.                                                                                                                     
    '''
    xsum = 0
    msum = 0
    h = (xmax-xmin)/N

    for i in range(1, N):
        xsum += f(xmin + i*h)
        print(xsum)

    for i in range(0, N):
        msum += f(xmin + (h/2) + i*h)
        print(msum)

    I = (h/6) * (f(xmin) + 4*msum + 2*xsum + f(xmax))
    return I


def f(x):
    '''Function equivalent to x^2 sin(x).'''
    return (x**2) * numpy.sin(x)


assert numpy.isclose(integrate_numeric(xmin=0, xmax=4, N=50), 1.096591)

注意:

  • 两个for循环中的范围已更改:我们希望第一个for循环以xmin + h的步长从xmin + (N-1)*hh(因此N-1总值),并且从xmin + h/2xmin + (N-1)*h + h/2的秒循环,以h的步长(N总值)。
  • [在最终计算中,无需将f应用于msumxsum:这些值已经是f值的和。我们仍然需要评估f的唯一位置是xminxmax。 (注意:此问题已在问题编辑中得到解决。)
  • h = (xmax-xmin)//N行必须为h = (xmax-xmin)/N。您只想在这里进行常规划分,而不是楼层划分。这可能是您最初得到零的原因:h本来是0
© www.soinside.com 2019 - 2024. All rights reserved.