Использование правила Симпсона для циклов (численное интегрирование) - PullRequest
Новая
       73

Использование правила Симпсона для циклов (численное интегрирование)

1 голос
/ 27 октября 2019

Я пытаюсь закодировать правило Симпсона в python, используя циклы for, и получаю ошибку подтверждения и не могу понять, почему. 

def integrate_numeric(xmin, xmax, N):
    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):
    return (x**2) * numpy.sin(x)

sample:

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

1 Ответ

2 голосов
/ 27 октября 2019

Вот исправленная версия вашего кода: 

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)*h с шагом h (т. е. N-1 суммарные значения), а второй цикл для перехода с xmin + h/2 до xmin + (N-1)*h + h/2 с шагом h (N итоговые значения).
  • В конечном вычислении нет необходимости применять f к msum и xsum: эти значения уже являются суммами значений f. Единственные места, которые нам еще нужно оценить f, это xmin и xmax. (Примечание: это уже было исправлено при редактировании вопроса.)
  • Строка h = (xmax-xmin)//N должна быть h = (xmax-xmin)/N. Вы просто хотите регулярное разделение здесь, а не разделение по полу. Вероятно, это причина появления нулей, которые вы получали изначально: h было бы 0.
...