Python: Найти главное значение интеграла численно

Question

Я решении интегрального численно с помощью Python:

где а (х) может принимать любое значение; положительный, отрицательный, внутри или вне [-1; 1], а эта - бесконечно малая положительная величина. Существует второй внешний интеграл который изменяет значение (х)

Я пытаюсь решить эту проблему с помощью Sokhotski–Plemelj theorem:

Однако это предполагает определение принципа значения, которое я не могу найти любой метод в python. Я знаю, что это реализовано в Matlab, но знает ли кто-либо из библиотеки или каким-либо другим способом определения основного значения в python (если существует принципиальное значение)?

Answer 1

Вы можете использовать sympy для непосредственного вычисления интеграла. Его действительная часть с эта-> 0 главное значение:

from sympy import * 
x, y, eta = symbols('x y eta', real=True) 
re(integrate(1/(x - y + I*eta), (x, -1, 1))).simplify().subs({eta: 0}) 
# -> log(Abs(-y + 1)/Abs(y + 1)) 

символический набор инструментов от Matlab int дает тот же результат, конечно (я не знаю других соответствующих инструментов в Matlab для этого --- пожалуйста укажите, знаете ли вы конкретный).

Вы задали вопрос о численном вычислении основного значения. Ответ заключается в том, что если у вас есть функция f(y), аналитическая форма или поведение которой вы не знаете, их вообще невозможно вычислить численно. Вам нужно знать такие вещи, как, например, полюсы подынтегрального выражения и какой порядок они есть.

Если с другой стороны, знаете, что ваш интеграл имеет вид f(y)/(y - y_0), scipy.integrate.quad может вычислить главное значение для вас, например:

import numpy as np 
from scipy import integrate, special 

# P \int_{-1}^1 dx 1/(x - wvar) * (1 + sin(x)) 
print(integrate.quad(lambda x: 1 + np.sin(x), -1, 1, weight='cauchy', wvar=0)) 
# -> (1.8921661407343657, 2.426947531830592e-13) 

# Check against known result 
print(2*special.sici(1)[0]) 
# -> 1.89216614073 

См here подробности.