2013-10-09 12 views
10

jestem rozwiązywania integralną numerycznie przy użyciu Pythona:Python: Znajdź główną wartość stanowi integralną numerycznie

enter image description here

gdzie a (x) może przyjmować dowolną wartość; dodatni, ujemny, wewnątrz lub na zewnątrz [-1; 1] i eta jest nieskończenie małą ilością dodatnią. Istnieje druga zewnętrzna całka który zmienia wartość (x)

próbuję rozwiązać ten problem za pomocą Sokhotski–Plemelj theorem: enter image description here

Jednak wiąże się to z ustalenia wartości do zasady, których nie mogę znajdź dowolną metodę w python. Wiem, że jest on zaimplementowany w Matlab, ale czy ktoś wie o bibliotece lub innym sposobie określania wartości głównej w pythonie (jeśli istnieje podstawowa wartość)?

+0

Jak zaimplementować go w MATLAB? – kyle

+0

W języku MATLAB symboliczna integracja "int" może obsłużyć wartości główne: http://se.mathworks.com/help/symbolic/int.html W przeciwnym razie całka numeryczna "integralna" może również obsługiwać osobliwości w punktach końcowych. Więc możesz podzielić całkę na dwie dokładnie dodać osobliwość, a następnie dodać dwa wyniki: http://se.mathworks.com/help/matlab/ref/integral.html?searchHighlight=integral –

Odpowiedz

5

Możesz użyć sympy do bezpośredniej oceny integralności. Jej część rzeczywistą z ETA> 0 jest wartością główny:

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)) 

symboliczny przybornik Matlaba int daje ten sam wynik, oczywiście (nie jestem świadomy innych odpowiednich narzędzi w Matlab dla tego proszę --- określić, czy znasz konkretny).

Pytasz o obliczenia numeryczne wartości głównej. Odpowiedź jest taka, że ​​jeśli masz tylko jedną funkcję, której analityczną formę lub zachowanie nie znasz, generalnie niemożliwe jest obliczenie ich numerycznie. Musisz wiedzieć, na przykład, gdzie znajdują się bieguny integralności i jaka jest ich kolejność.

Jeśli z drugiej strony wiem, swoją integralną ma postać f(y)/(y - y_0), scipy.integrate.quad może obliczyć wartość podstawową dla ciebie, na przykład:

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 

Zobacz here szczegóły.