Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

Short version of my question

Can I only apply the Gauss-Hermite routine with an infinite interval or can I transform the interval?

Long version (reason I'm asking)

I am interested in solving integral equations numerically and one part of this, if I understood it correctly, is the usage of numerical integration routines. So I ventured into this field and implemented the midpoint, the Simpson, and the Gauss-Legendre routine successfully.

Here is my implementation of the Gauss-Legendre routine in R (depending on the package 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, ...))

}

So basically, the first two parts give me the numbers tabulated here and the third part just computes:

$$\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)$$

Let's check that it works (note that integrate is just the built-in function in 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

As you can see, this works with an interval ranging from 2 to 10 by transforming the $x$ values. Given that the interval is well defined, I don't see any problem using that for any interval.

Next, I implemented the Gauss-Hermite routine:

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))

}

This works for instance when computing the integral of the density function of the standard normal distribution:

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

However, I don't get any useful results when the interval doesn't go from $[-\infty; \infty]$. On the other side, in none of the material I read was it mentioned that the Gauss-Hermite approach is limited to that case (I read that it's useful in connection with Normal random variables though). So my question (sorry for the long prolog): can I only apply the Gauss-Hermite routine with an infinite interval or can I transform the interval?

As you probably guessed from the question, I'm not a mathematician, so any help with this would be great.

share|improve this question
    
My phone is a poor platform for this, but I would wager that it might be possible to use an approach similar to generalized polynomial chaos to manage this. The connection between the Gauss-Hermite quadrature and normal random variables is exactly the Gauss-Hermite polynomial chaos. –  Arkamis Aug 8 '12 at 23:36

1 Answer 1

up vote 1 down vote accepted

Yes, there are transformations that transform finite integrals to doubly infinite integrals, like the $\tanh$-substitution

$$\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$$

and yes, you can apply Gauss-Hermite quadrature to any suitably transformed doubly infinite integral:

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

BUT

Gauss-Hermite doesn't work very well if the functions you are integrating are not of the form $\exp(-x^2)f(x)$, where $f(x)$ is a function that is well approximated by a polynomial. Remember that the idea of Gaussian quadrature in general is to "factor out" unruly behavior in your integrands, and keep that behavior to the nodes and weights of that quadrature rule. That supposedly leaves you with a remainder that is much more well-behaved for quadrature than the original integrand. In effect, trying to use Gauss-Hermite to integrate things that don't have the factor $\exp(-x^2)$ is like trying to use a toothbrush to scrub a toilet: you'll manage to finish, but it'll take you a while.

To drive the point home, consider the two integrals

$$\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*}$$

Here are a few short Mathematica routines for generating the nodes and weights for Gauss-Hermite quadrature:

(* 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]]]]

Here's a short table of the relative error in using $n$-point Gauss-Hermite quadrature for numerically evaluating these two integrals (done with 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}$$

Clearly, the convergence of Gauss-Hermite quadrature for $\mathcal I_1$ is rather shabby compared to the relatively quicker convergence for $\mathcal I_2$. Sometimes, you'll get lucky and find a function where Gauss-Hermite performs well even if it does not have an explicit $\exp(-x^2)$ factor, but those things aren't that common.

share|improve this answer
    
To quote you: Remember that the idea of Gaussian quadrature in general is to "factor out" unruly behavior in your integrands, and keep that behavior to the nodes and weights of that quadrature rule. To be honest, I didn't remember that because wherever I looked, they only showed how the different Gauss quadrature worked, but didn't really try to explain why there are different methods and when to apply a specific one. So thanks a lot for the great answer in which you explained when Gauss-Hermite works well and when not. That helps me a lot. –  Christoph_J Aug 9 '12 at 6:33
1  
Yes, I suppose a lot of the usual textbooks have been remiss in that regard. Anyway, that's the point of being interested in Gaussian quadrature in general: "if the integrand has an $\exp(-x)$ factor and you are integrating over $[0,\infty)$, use Gauss-Laguerre; if you are integrating over $(-1,1)$ and your integrand has a $\dfrac1{\sqrt{1-x^2}}$ factor, use Gauss-Chebyshev..." and so on. –  J. M. Aug 9 '12 at 6:40
    
If the "misbehaving factor" of your integrand does not correspond to a well-known system of orthogonal polynomials, then you have more work to do. Gautschi's book deals with that more complicated situation. –  J. M. Aug 9 '12 at 6:42

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.