Интеграция функции с особенностями с использованием четырехъядерной процедуры scipy

Я использую функцию quad из scipy.integrate v0.19.1 для интеграции функций с сингулярностью, подобной квадратному корню, на каждом конце интервала интегрирования, например,

In [1]: quad(lambda x: 1/sqrt(1-x**2), -1, 1)

(Я использую функцию sqrt из numpy v1.12.0), которая немедленно дает правильный результат pi:

Out[1]: (3.141592653589591, 6.200897573194197e-10)

В соответствии с документацией функции quad ключевое слово points следует использовать для указания местоположения сингулярностей или разрывов подынтегральной функции, но если я укажу точки [1, -1], где указанная выше подынтегральная функция сингулярна, я получу предупреждение и nan в результате:

In [2]: quad(lambda x: 1/sqrt(1-x**2), -1, 1, points=[-1, 1])

RuntimeWarning: divide by zero encountered in double_scalars
IntegrationWarning: Extremely bad integrand behavior occurs at some
points of the integration interval.

Out[2]: (nan, nan)

Может кто-нибудь пояснить, почему quad выдает эти проблемы, если указаны особенности подынтегральной функции, и просто работает, если эти точки не указаны?

EDIT: я думаю, что нашел правильное решение этой проблемы. На случай, если кто-то еще столкнется с подобными проблемами, я хочу быстро поделиться своими выводами:

Если вы хотите интегрировать функцию вида f(x)*g(x) с гладкой функцией f(x) и g(x) = (x-a)**alpha * (b-x)**beta, где a и b — пределы интегрирования, а g(x) имеет особенности в этих пределах, если alpha, beta < 0, то вы должны просто интегрировать f(x), используя g(x) как < em>функция взвешивания. Для подпрограммы quad это возможно с использованием аргументов weight и wvar. С помощью этих аргументов вы также можете обрабатывать сингулярности различных видов и проблематичное колебательное поведение. Весовую функцию g(x), определенную выше, можно использовать, установив weight='alg' и используя wvar=(alpha, beta) для указания показателей степени в g(x).

Поскольку 1/sqrt(1-x**2) = (x+1)**(-1/2) * (1-x)**(-1/2) теперь я могу обрабатывать интеграл следующим образом:

In [1]: quad(lambda x: 1, -1, 1, weight='alg', wvar=(-1/2, -1/2))
Out[1]: (3.1415926535897927, 9.860180640534107e-14)

который дает правильный ответ pi с очень высокой точностью, независимо от того, использую ли я аргумент points=(-1, 1) или нет (который, насколько я теперь понимаю, следует использовать только в том случае, если особенности/разрывы не могут быть обработаны путем выбора подходящего весовая функция).


python scipy quad
person TMueller83    schedule 30.06.2017    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Параметр points предназначен для особенностей/разрывов, возникающих внутри интервала интегрирования. Он не предназначен для конечных точек интервала. Итак, в вашем примере версия без points является правильным подходом. (Я не могу точно определить, что пошло не так, когда конечные точки включены в points, не углубляясь в код FORTRAN, который обертывает SciPy.)

Сравните со следующим примером, где особенности возникают на интервале интегрирования:

>>> quad(lambda x: 1/sqrt(abs(1-x**2)), -2, 2)
(inf, inf)
>>> quad(lambda x: 1/sqrt(abs(1-x**2)), -2, 2, points = [-1, 1])
(5.775508447436837, 7.264979728915932e-10)

Здесь включение points уместно и дает правильный результат, в то время как без points вывод бесполезен.

person Community    schedule 03.07.2017
comment
Привет, Уолт, спасибо за добрый ответ! Я думаю, что это правильный путь. - person TMueller83; 03.07.2017