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

Я пытаюсь закодировать правило Симпсона на 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)

образец:

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

python-3.x python loops numerical-methods simpsons-rule
person Oliver Moore    schedule 27.10.2019    source источник
comment
вы можете включить f и образец ввода, который дает сбой?   -  person Christian Sloper    schedule 27.10.2019
comment
@ChristianSloper да, конечно   -  person Oliver Moore    schedule 27.10.2019
comment
Почему вы выполняете только 4 оценки функций?   -  person Mark Dickinson    schedule 27.10.2019
comment
@MarkDickinson, что ты имеешь в виду?   -  person Oliver Moore    schedule 27.10.2019
comment
У вас N = 1, диапазон (0, N) - всего 1 шаг   -  person Christian Sloper    schedule 27.10.2019
comment
@OliverMoore: Вы делаете только четыре вызова функций f(something) в предпоследней строке. Для правила Симпсона с N + 1 очками (N делений) вы должны коллировать f N + 1 раз. Вы хотите суммировать f значений в ваших циклах. Вместо этого вы суммируете аргументы с этими вызовами (значения x).   -  person Mark Dickinson    schedule 27.10.2019
comment
@MarkDickinson вы имеете в виду, что для этого должно быть: xsum + = f (xmin + i * h)   -  person Oliver Moore    schedule 27.10.2019
comment
@OliverMoore: Да, именно так.   -  person Mark Dickinson    schedule 27.10.2019
comment
@MarkDickinson, когда я меняю код на этот, я все равно получаю 0   -  person Oliver Moore    schedule 27.10.2019
comment
@OliverMoore: Да; это не единственная проблема с кодом. Вы хотите (xmax-xmin)/N вместо (xmax-xmin)//N. Могут быть и другие проблемы.   -  person Mark Dickinson    schedule 27.10.2019
comment
@MarkDickinson ok Спасибо, я теперь не получаю только 0 результатов, но я все еще получаю ту же ошибку утверждения.   -  person Oliver Moore    schedule 27.10.2019

Ответы (1)


arrow_upward
2
arrow_downward

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

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 общих значений), а второй цикл for переходил с 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.
person Mark Dickinson    schedule 27.10.2019