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.

Why does the following approximation work very well for small theta's but not for large theta's?


$$ I = \int_0^\infty x (1 + \theta x)^{-2 \left(\frac{1}{\theta} + 1 \right)} \, \mathrm{d}x $$

Exact evaluation using integration by part

$$ I = \frac{1}{2(\theta + 2)} $$

Approximation using Gauss-Laguerre quadrature $$ I \approx \sum_{k=1}^N w_k \exp(x_k) \, x_k (1 + \theta x_k)^{-2 \left(\frac{1}{\theta} + 1 \right)} $$ where the $x_k$'s and the $w_k$'s can be found, e.g., here.

Comparison (N = 100)

  • theta = 0.5806:

    true = 0.1937534 and approx = 0.1937522

  • theta = 15:

    true = 0.02941176 and approx = 0.01662866


my R-code

#------ package ------
library(gaussquad)
#---------------------

#------ f(x) ------
f <- function(x, theta)
{
  x * (1 + theta * x)^(-2 * ((1 / theta) + 1))
}
#------------------

#------ xk and wk ------
laguerre <- glaguerre.quadrature.rules(n=100, alpha=0, normalized=FALSE)[[100]]
x <- laguerre[, 1]
w <- laguerre[, 2]
#-----------------------

#------ approx ------
approx <- function(theta, x, w)
{
  sum(w * exp(x) * f(x=x, theta=theta))
}
#--------------------

> approx(theta=0.5806, x=x, w=w)
[1] 0.1937522
> approx(theta=15, x=x, w=w)
[1] 0.01662866
share|improve this question
    
Sorry, your integral is not showing up in the post. –  Ron Gordon Apr 22 '13 at 10:35
    
@RonGordon: Should be fixed now. Thanks! –  Marco Apr 22 '13 at 10:38
    
Sorry again, I still cannot see it. –  Ron Gordon Apr 22 '13 at 10:40
    
@RonGordon: Sorry for the trouble! Hope it is ok now –  Marco Apr 22 '13 at 11:14
    
You should recall that Gauss-Laguerre quadrature works best if the "infinite extent" portion of your integral behaves very much like $\exp(-x)$. This for instance is the reason why Gauss-Laguerre is not very useful for computing integrals like $\int_0^\infty \frac{\mathrm dx}{1+x^2}$. –  J. M. Apr 22 '13 at 12:36
show 7 more comments

Your Answer

 
discard

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

Browse other questions tagged or ask your own question.