Численное интегрирование с особенностями в python (главное значение) - PullRequest
0 голосов
/ 08 октября 2018

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

from scipy.integrate import quad
import numpy as np
def fun(x):
    return 1./(1-x**2)
quad(fun, -2, 2, points=[-1, 1])

Это приводит к IntegrationWarning и возвращает значение около 0.4.

Полюса для функции [-1,1].Ответ должен быть примерно 1,09 (рассчитано с помощью ручки и бумаги).

Ответы [ 2 ]

0 голосов
/ 08 октября 2018

Опцию weight='cauchy' можно использовать для эффективного вычисления основного значения расходящихся интегралов, как этот.Это означает, что функция, предоставленная для quad, будет неявно умножена на 1/(x-wvar), поэтому настройте эту функцию соответствующим образом (умножьте ее на x-wvar, где wvar - точка сингулярности).

i1 = quad(lambda x: -1./(x+1), 0, 2, weight='cauchy', wvar=1)[0]
i2 = quad(lambda x: -1./(x-1), -2, 0, weight='cauchy', wvar=-1)[0]
result = i1 + i2

результат 1.0986122886681091.

С помощью такой простой функции вы также можете выполнять символическую интеграцию с SymPy:

from sympy import symbols, integrate
x = symbols('x')
f = integrate(1/(1-x**2), x)
result = (f.subs(x, 2) - f.subs(x, -2)).evalf()

Результат: 1.09861228866811.Без evalf() это было бы log(3).

0 голосов
/ 08 октября 2018

Я также не смог заставить его работать с оригинальной функцией.Я придумал это, чтобы оценить главное значение в scipy:

def principal_value(func, a, b, poles, eps=10**(-6)):
    #edges
    res = quad(func,a,poles[0]-eps)[0]+quad(func,poles[-1]+eps,b)[0]
    #inner part
    for i in range(len(poles)-1):
        res += quad(func, poles[i]+eps, poles[i+1]-eps)[0]
    return res

Где func - ваш дескриптор функции, a и b - ограничения, poles - список полюсов.и eps это то, как близко вы хотите приблизиться к полюсам.Вы можете сделать eps меньше и меньше, чтобы получить лучший результат, но, возможно, sympy будет лучше для такой проблемы.

С этой функцией и стандартом eps я получаю 1.0986112886023367 в результате, чтопочти так же, как дает Вольфрамальфа.

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...