Czy przekształca

15

Słyszałem anegdotycznie, że gdy ktoś próbuje liczbowo wykonać całkę formy

0f(x)J0(x)dx

z gładkim i dobrze się zachowującym (np. nie jest silnie oscylacyjny, niesingularny itp.), wówczas poprawność przepisania go jakofa(x)

1π0π0fa(x)sałata(xgrzechθ)rexreθ

i najpierw wykonaj całkę wewnętrzną numerycznie. Nie widzę powodu, dla którego powinienem oczekiwać, że to zadziała, ale z drugiej strony dokładność metody numerycznej jest rzadko oczywista.

Oczywiście wiem, że najlepszym sposobem na to jest użycie metody zoptymalizowanej dla takich całek oscylacyjnych, ale ze względu na ciekawość załóżmy, że ograniczę się do zastosowania reguły kwadraturowej. Czy ktoś może potwierdzić lub obalić, że dokonanie tej transformacji poprawia dokładność całki? I / lub wskaż mi źródło, które to wyjaśni?

David Z
źródło
1
Zintegrowany ponad ... Jest to jedna z integralnych definicji funkcji Bessela. 0θπ
David Z
4
Zatem twoje pytanie brzmi: biorąc pod uwagę ogólne równania kwadraturowe punktowe Q N [ ] na [ 0 , ) i Q N π [ ] na [ 0 , π ] , to Q N M [ fN.QN.[][0,)QπN.[][0,π] gorszy lub lepszy niż Q M π [ Q N [ f ( x )QN.M.[fajot0] . QπM.[QN.[fa(x)sałata(xgrzechθ)]]
Stefano M
@StefanoM tak, zgadza się.
David Z
FWIW, jedną z najbardziej wydajnych metod oceny funkcji Bessela rzędu zerowego jest reguła trapezoidalna, która jest znana z tego, że daje bardzo dokładne wyniki przy całkowaniu okresowych całek w jednym okresie (nawet lepszym niż zwykły standard, kwadratura Gaussa). Więc: może pomóc, może nie.
JM

Odpowiedzi:

3

Nie sądzę, żeby to miało jakąkolwiek różnicę. Musisz wybrać wystarczająco wysoką kwadraturę dla całki powyżej , aby była ona równa funkcji Bessela J 0 . Wybrałem kolejność 20 w poniższym przykładzie, ale zawsze musisz dokonać zbieżności w odniesieniu do dokładnej funkcji i przedziału, przez który integrujesz. Następnie dokonałem zbieżności z n , rzędem kwadratury Gaussa całki nad x . Wybrałem f ( x ) = e - x x 2 i używam domeny [ 0 , x max ] , możesz zmienić x maxθjot0nxfa(x)=mi-xx2)[0,xmax]xmaxponiżej. Mam:

 n      direct         rewritten
 1  0.770878284949  0.770878284949
 2  0.304480978430  0.304480978430
 3  0.356922151260  0.356922151260
 4  0.362576361509  0.362576361509
 5  0.362316789057  0.362316789057
 6  0.362314010897  0.362314010897
 7  0.362314071949  0.362314071949
 8  0.362314072182  0.362314072182
 9  0.362314072179  0.362314072179
10  0.362314072179  0.362314072179

n=9 obie całki są w pełni zbieżne do 12 cyfr znaczących.

Oto kod:

from scipy.integrate import fixed_quad
from scipy.special import jn
from numpy import exp, pi, sin, cos, array

def gauss(f, a, b, n):
    """Gauss quadrature"""
    return fixed_quad(f, a, b, n=n)[0]

def f(x):
    """Function f(x) to integrate"""
    return exp(-x) * x**2

xmax = 3.

print " n      direct         rewritten"
for n in range(1, 20):
    def inner(theta_array):
        return array([gauss(lambda x: f(x) * cos(x*sin(theta)), 0, xmax, n)
            for theta in theta_array])
    direct = gauss(lambda x: f(x) * jn(0, x), 0, xmax, n)
    rewritten = gauss(inner, 0, pi, 20) / pi
    print "%2d  %.12f  %.12f" % (n, direct, rewritten)

xmax[0,]f(x)rewritten = gauss(inner, 0, pi, 20) / pi

Ondřej Čertík
źródło
Podejrzewam, że masz rację, moje własne testy wykazały podobne wyniki.
David Z