Transformar

15

Ouvi de maneira anedótica que quando alguém está tentando fazer numericamente uma integral do formulário

0f(x)J0(x)dx

com suave e bem-comportado (por exemplo, não seja altamente oscilatório, não singular, etc.), ajudará a precisão a reescrevê-lo comof(x)

1π0π0f(x)cos(xsinθ)dxdθ

e execute a integral interna numericamente primeiro. Não vejo nenhum motivo para esperar que isso funcione, mas, novamente, a precisão de um método numérico raramente é óbvia.

É claro que sei que a melhor maneira de fazê-lo é usar um método otimizado para integrais oscilatórias como este, mas, por uma questão de curiosidade, suponha que eu me limite a usar alguma regra de quadratura. Alguém pode confirmar ou refutar que fazer essa transformação tende a melhorar a precisão da integral? E / ou me aponte para uma fonte que a explique?

David Z
fonte
1
Integrado acima de ... É uma das definições integrais da função Bessel. 0θπ
David Z
4
Então a pergunta é: Dado genérico -ponto quadratura fórmulas Q N [ ] em [ 0 , ) e Q N π [ ] em [ 0 , π ] , é Q N M [ fNQN[][0,)QπN[][0 0,π] pior ou melhor que Q M π [ Q N [ f ( x )QNM[fJ0 0] . QπM[QN[f(x)porque(xpecadoθ)]]
21813 Stefano M
@StefanoM sim, está certo.
David Z
O FWIW, um dos métodos mais eficientes para avaliar a função de Bessel de ordem zero, é a regra trapezoidal, conhecida por fornecer resultados muito precisos ao integrar integrandos periódicos durante um período (ainda melhor que o padrão usual, quadratura gaussiana). Então: pode ajudar, pode não.
JM 14/05

Respostas:

3

Eu não acho que isso faça alguma diferença. Você deve escolher uma quadratura alta o suficiente para a integral sobre para que seja igual à função de Bessel J 0 . Eu escolhi a ordem 20 no exemplo abaixo, mas você sempre precisa fazer convergência com relação à função e intervalo exatos nos quais você integra. Então fiz convergência com n , a ordem da quadratura gaussiana da integral sobre x . Eu escolhi f ( x ) = e - x x 2 e uso o domínio [ 0 , x max ] , você pode alterar x maxθJ0 0nxf(x)=e-xx2[0 0,xmax]xmaxabaixo. Eu tenho:

 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

Aqui está o código:

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 0,]f(x)rewritten = gauss(inner, 0, pi, 20) / pi

Ondřej Čertík
fonte
Eu suspeito que você esteja certo, meus próprios testes mostraram resultados semelhantes.
David Z