La transformation de aide-t-elle à l'intégration numérique?

J'ai entendu de façon anecdotique que lorsque l'on essaie de faire numériquement une intégrale de la forme

\int_{0}^{\infty} F (X) J_{0} (X) ré X

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

avec lisse et bien comporté (par exemple pas lui-même très oscillatoire, non singulier, etc.), alors il aidera la précision à le réécrire comme $f(x)$

\frac{1}{π} \int_{0}^{π} \int_{0}^{\infty} F (X) \cos (X péché θ) ré X ré θ

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

et effectuer d'abord l'intégrale intérieure numériquement. Je ne vois aucune raison pour que je m'attende à ce que cela fonctionne, mais là encore, la précision d'une méthode numérique est rarement évidente.

Bien sûr, je sais que la meilleure façon de le faire est d'utiliser une méthode optimisée pour les intégrales oscillantes comme celle-ci, mais par curiosité, supposons que je me limite à utiliser une règle de quadrature. Quelqu'un peut-il confirmer ou réfuter que la réalisation de cette transformation tend à améliorer la précision de l'intégrale? Et / ou me pointer vers une source qui l'explique?

quadrature special-functions David Z
la source

Intégré sur ... C'est l'une des définitions intégrales de la fonction Bessel.

0 \leq θ \leq π

$0\le\theta \le\pi$

David Z

Votre question est donc la suivante: étant donné les formules génériques de quadrature à points sur et sur , est pire ou meilleur que .

N

$N$

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

$Q_\infty^N[\cdot]$

[0, \infty)

$[0,\infty)$

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

$Q_\pi^N[\cdot]$

[0, π]

$[0,\pi]$

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 oui, c'est vrai.

David Z

FWIW, l'une des méthodes les plus efficaces pour évaluer la fonction de Bessel d'ordre zéro est la règle trapézoïdale, qui est bien connue pour donner des résultats très précis lors de l'intégration d'intégrandes périodiques sur une période (encore mieux que la norme habituelle, la quadrature gaussienne). Donc: ça pourrait aider, ça pourrait ne pas l'être.

Réponses:

Je ne pense pas que cela fasse une différence. Vous devez choisir une quadrature suffisamment élevée pour l'intégrale sur afin qu'elle soit égale à la fonction de Bessel . J'ai choisi l'ordre 20 dans l'exemple ci-dessous, mais vous devez toujours faire une convergence en ce qui concerne la fonction exacte et l'intervalle sur lequel vous intégrez. Ensuite, j'ai fait la convergence avec , l'ordre de la quadrature gaussienne de l'intégrale sur . J'ai choisi et j'utilise le domaine , vous pouvez changer $\theta$ $J_0$ $n$ $x$ $f(x) = e^{-x} x^2$ $[0, x_\text{max}]$ $x_\text{max}$ au dessous de. J'ai eu:

 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

Comme vous pouvez le voir, pour deux intégrales sont entièrement convergées vers 12 chiffres significatifs. $n=9$

Voici le code:

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)

Vous pouvez jouer avec cela vous-même, simplement changer xmax, peut-être pourriez-vous avoir besoin de diviser l'intervalle en éléments et d'intégrer élément par élément. Vous pouvez également modifier la fonction . Assurez-vous que vous faites toujours converger l'intégrale , c'est-à-dire commencez par un ordre faible et continuez à l'augmenter jusqu'à ce que les résultats imprimés cessent de changer. $[0, \infty]$ f(x)rewritten = gauss(inner, 0, pi, 20) / pi

Ondřej Čertík
la source

Je suppose que vous avez raison, mes propres tests ont montré des résultats similaires.

David Z