Czy przekształca

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

\int_{0}^{\infty} f (x) J_{0} (x) d x

$\int_0^\infty f(x) J_0(x)\,\mathrm{d}x$

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

\frac{1}{π} \int_{0}^{π} \int_{0}^{\infty} fa (x) sałata (x grzech θ) re x re θ

$\frac{1}{\pi}\int_0^\pi \int_0^\infty f(x) \cos(x\sin\theta) \,\mathrm{d}x\,\mathrm{d}\theta$

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?

quadrature special-functions David Z
źródło

Zintegrowany ponad

... Jest to jedna z integralnych definicji funkcji Bessela.

0 \leq θ \leq π

$0\le\theta \le\pi$

David Z

Zatem twoje pytanie brzmi: biorąc pod uwagę ogólne równania kwadraturowe

punktowe

, to

N

$N$

Q_{\infty}^{N} [\cdot]

$Q_\infty^N[\cdot]$

[0, \infty)

$[0,\infty)$

Q_{π}^{N} [\cdot]

$Q_\pi^N[\cdot]$

[0, π]

$[0,\pi]$

gorszy lub lepszy niż

Q_{\infty}^{N \cdot M} [f J_{0}]

$Q^{N\cdot M }_\infty[f\,J_0]$

Q_{π}^{M} [Q_{\infty}^{N} [f (x) \cos (x \sin θ)]]

$Q^M_\pi[Q_\infty^N[f(x)\,\cos(x\sin\theta)]]$

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.

Odpowiedzi:

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 . 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 , rzędem kwadratury Gaussa całki nad . Wybrałem i używam domeny , możesz zmienić $\theta$ $J_0$ $n$ $x$ $f(x) = e^{-x} x^2$ $[0, x_\text{max}]$ $x_\text{max}$ poniż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, \infty]$ 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