Python: Finden Hauptwert eines integralen numerisch

Question

ich das Integral numerisch mit Python zu lösen:

wo a (x) kann auf jeden beliebigen Wert annehmen; positiv, negativ, innerhalb oder außerhalb des [-1; 1] und eta ist eine infinitesimale positive Größe. Es gibt ein zweites äußerte Integral ich ändert den Wert eines (x)

Ich versuche dies die Sokhotski–Plemelj theorem mit zu lösen:

aber dies beinhaltet den Grundsatz Wert zu bestimmen, was ich nicht kann Finde irgendeine Methode in Python. Ich weiß, dass es in Matlab implementiert ist, aber kennt jemand eine Bibliothek oder eine andere Möglichkeit, den Hauptwert in Python zu bestimmen (wenn ein Hauptwert existiert)?

Answer 1

Mit sympy können Sie das Integral direkt auswerten. Ihr eigentlicher Teil mit eta-> 0 ist der Hauptwert:

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 symbolische Toolbox int gibt Ihnen das gleiche Ergebnis, natürlich (ich bin mir nicht bewusst andere relevante Tools in Matlab für diesen --- bitte geben Sie an, wenn Sie einen bestimmten kennen).

Sie haben nach der numerischen Berechnung eines Hauptwerts gefragt. Die Antwort dort ist, dass wenn Sie nur eine Funktion f(y) haben, deren analytische Form oder Verhalten Sie nicht kennen, ist es im Allgemeinen unmöglich, sie numerisch zu berechnen. Sie müssen wissen, wie die Pole des Integranden sind und in welcher Reihenfolge sie sind.

Wenn Sie auf der anderen Seite kennen Ihre Integral der Form f(y)/(y - y_0) ist, kann scipy.integrate.quad den Hauptwert berechnen für Sie, zum Beispiel:

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 für Details.