¿Puedo aplicar la rutina Gauss-Hermite sólo con un intervalo infinito o puedo transformar el intervalo?

Question

Versión corta de mi pregunta

¿Puedo aplicar la rutina Gauss-Hermite sólo con un intervalo infinito o puedo transformar el intervalo?

Versión larga (razón por la que pregunto)

Estoy interesado en resolver ecuaciones integrales numéricamente y una parte de esto, si lo he entendido bien, es el uso de rutinas de integración numérica. Así que me aventuré en este campo e implementé con éxito el punto medio, el Simpson y la rutina de Gauss-Legendre.

Aquí está mi implementación de la rutina Gauss-Legendre en R (dependiendo del paquete statmod ):

gauss_legendre <- function(f, upper, lower, nr_int, ...) {

  #1. Compute the nodes; those are the roots of the legendre-polynom
  leg_poly <- gauss.quad(nr_int)
  x <- leg_poly$nodes

  #2. Compute the weights as follows: 2/((1-x^2)Ableitung(P(x))^2)
  w <- leg_poly$weights

  #3. Compute the integral
  x_transformed <- ((x+1)*(upper - lower)/2 + lower)
  (upper - lower)/2*sum(w * mapply(f, x_transformed, ...))

}

Así que, básicamente, las dos primeras partes me dan los números tabulados aquí y la tercera parte sólo computa:

$$\int_a^b f(x)dx = \frac{b-a}{2}\sum\limits_{i=1}^n \omega_i f\left( \frac{b-a}{2}x_i + \frac{a + b}{2} \right)$$

Comprobemos que funciona (nótese que integrate es sólo la función incorporada en R):

gauss_legendre(function(x) {x^3 + log(x) + sqrt(x)}, 10, 2, 10)
[1] 2528.836
> integrate(function(x) {x^3 + log(x) + sqrt(x)}, 2, 10)
2528.836 with absolute error < 2.8e-11

Como puede ver, esto funciona con un intervalo que va de 2 a 10 transformando el $x$ valores. Dado que el intervalo está bien definido, no veo ningún problema en utilizarlo para cualquier intervalo.

A continuación, implementé la rutina Gauss-Hermite:

gauss_hermite <- function(f, upper, lower, nr_int, ...) {

  #1. Compute the nodes; those are the roots of the Hermite-polynom
  her_poly <- gauss.quad(nr_int, "hermite")
  x <- her_poly$nodes

  #2. Compute the weights according to the function; attention: this already includes the term exp(x^2)
  w <- her_poly$weights * exp(x^2)

  #3. Compute the integral
  sum(w * mapply(f, x))

}

Esto funciona, por ejemplo, al calcular la integral de la función de densidad de la distribución normal estándar:

gauss_hermite(dnorm, -Inf, Inf, 100)
1

Sin embargo, no obtengo ningún resultado útil cuando el intervalo no va de $[-\infty; \infty]$ . Por otro lado, en ninguno de los materiales que he leído se menciona que el enfoque de Gauss-Hermite esté limitado a ese caso (aunque he leído que es útil en relación con las variables aleatorias normales). Así que mi pregunta (perdón por el largo prólogo): ¿puedo aplicar la rutina de Gauss-Hermite sólo con un intervalo infinito o puedo transformar el intervalo?

Como probablemente habrás adivinado por la pregunta, no soy un matemático, así que cualquier ayuda con esto sería genial.

Answer 1

Sí, hay transformaciones que transforman integrales finitas en integrales doblemente infinitas, como la $\tanh$ -sustitución

$$\int_a^b f(x) \mathrm dx=\frac{b-a}{2}\int_{-\infty}^\infty f\left(\frac{a+b}{2}+\frac{b-a}{2}\tanh\,u\right)\operatorname{sech}^2 u \,\mathrm du$$

y sí, se puede aplicar la cuadratura de Gauss-Hermite a cualquier integral doblemente infinita convenientemente transformada:

$$\int_{-\infty}^\infty f(v)\mathrm dv\approx\sum_{k=1}^n w_k \exp(x_k^2)f(x_k)$$

PERO

Gauss-Hermite no funciona muy bien si las funciones que se integran no son de la forma $\exp(-x^2)f(x)$ , donde $f(x)$ es una función bien aproximada por un polinomio. Recuerde que la idea de la cuadratura gaussiana en general es "factorizar" el comportamiento rebelde en sus integrados, y mantener ese comportamiento a los nodos y pesos de esa regla de cuadratura. Eso supuestamente te deja con un resto que se comporta mucho mejor para la cuadratura que el integrando original. En efecto, tratar de usar Gauss-Hermite para integrar cosas que no tienen el factor $\exp(-x^2)$ es como intentar usar un cepillo de dientes para fregar un inodoro: conseguirás terminar, pero te llevará un tiempo.

Para entenderlo mejor, consideremos las dos integrales

$$\begin{align*} \mathcal I_1&=\int_{-\infty}^\infty \operatorname{sech}\,u;\mathrm du=\pi\\ \mathcal I_2&=\int_{-\infty}^\infty \exp(-u^2)\operatorname{sech}\,u\;\mathrm du\approx 1.4790611714495759 \end{align*}$$

Aquí hay unos cuantos Mathematica rutinas para generar los nodos y pesos para la cuadratura Gauss-Hermite:

(* Golub-Welsch algorithm *)
golubWelsch[d_?VectorQ, e_?VectorQ] := 
 Transpose[
  MapAt[(First[e] Map[First, #]^2) &, 
   Eigensystem[
    SparseArray[{Band[{1, 1}] -> d, Band[{1, 2}] -> Sqrt[Rest[e]], 
      Band[{2, 1}] -> Sqrt[Rest[e]]}, {Length[d], Length[d]}]], {2}]]

(* generate nodes and weights for Gauss-Hermite quadrature *)
ghq[n_Integer, prec_: MachinePrecision] := 
 Transpose[
  Sort[golubWelsch[ConstantArray[0, n], 
    N[Prepend[Range[n - 1]/2, Sqrt[Pi]], prec]]]]

Aquí hay una breve tabla del error relativo en el uso de $n$ -La cuadratura de Gauss-Hermite para la evaluación numérica de estas dos integrales (realizada con Mathematica ):

$$\begin{array}{c|c|c} n&-\log_{10}\left|\frac{\text{estimate}}{\mathcal I_1}-1\right| &-\log_{10}\left|\frac{\text{estimate}}{\mathcal I_2}-1\right|\\\hline 2&0.58153&1.3066\\ 5&1.0804&2.6953\\ 10&1.6518&4.3750\\ 15&2.0926&5.6971\\ 20&2.4643&6.8235\\ 25&2.7917&7.8217\\ \end{array}$$

Claramente, la convergencia de la cuadratura de Gauss-Hermite para $\mathcal I_1$ es bastante pobre en comparación con la convergencia relativamente más rápida de $\mathcal I_2$ . A veces, tendrás suerte y encontrarás una función en la que Gauss-Hermite funcione bien aunque no tenga un $\exp(-x^2)$ factor, pero esas cosas no son tan comunes.