1

I am trying to solve an integral in R. However, I am getting an error when I am trying to solve for that integral.

The equation that I am trying to solve is as follows:

$$ C_m = \frac{{abs{x}}e^{2x}}{\pi^{1/2}}\int_0^t t^{-3/2}e^{-x^2/t-t}dt $$

enter image description here

The code that I am using is as follows:

a <- seq(from=-10, by=0.5,length=100)

## Create a function to compute integration
Cfun <- function(XX, upper){
  integrand <- function(x)x^(-1.5)*exp((-XX^2/x)-x)
  integrated <- integrate(integrand, lower=0, upper=upper)$value
  (final <- abs(XX)*pi^(-0.5)*exp(2*XX)*integrated) }


b<- sapply(a, Cfun, upper=1)

The error that I am getting is as follows:

Error in integrate(integrand, lower = 0, upper = upper) : 
  the integral is probably divergent

Does this mean I cannot solve the integral ?

Any possible ways to fix this problem will be highly appreciated.

Thanks.

Jd Baba
  • 5,948
  • 18
  • 62
  • 96
  • Mnel, I am sorry. I forgot to edit my integrand. It should be `x^(-1.5)`. I will edit that. – Jd Baba Apr 09 '13 at 03:21
  • It is **always** a good idea to plot your integrand in cases like this. Had you done so, you'd have seen something funny going on near `x=0` , as mnel has pointed out. – Carl Witthoft Apr 09 '13 at 11:41

1 Answers1

3

You could wrap the call to Cfun in a try statement 

# note using `lapply` so errors don't coerce the result to character
b <- lapply(a, function(x,...) try(Cfun(x, ...), silent = TRUE), upper = 1)

If you wanted to replace the errors with NA values and print a warning that the integration threw an error

Cfun <- function(XX, upper){
  integrand <- function(x)x^(-1.5)*exp((-XX^2/x)-x)
  int <- try(integrate(integrand, lower=0, upper=upper), silent = TRUE)
  if(inherits(int ,'try-error')){
    warning(as.vector(int))
    integrated <- NA_real_
  } else {
    integrated <- int$value
  }
  (final <- abs(XX)*pi^(-0.5)*exp(2*XX)*integrated) }

Edit (facepalm moment)

Your error arises because your function is not defined when t=0 (you divide by t within the integrand).

Cfun <- function(XX, upper){
  integrand <- function(x)x^(-1.5)*exp((-XX^2/x)-x)
  # deal with xx=0
  if(isTRUE(all.equal(XX, 0)){
      warning('The integrand is not defined at XX = 0')
      return(NA_real_)
  }
  # deal with other integration errors
  int <- try(integrate(integrand, lower=0, upper=upper), silent = TRUE)
  if(inherits(int ,'try-error')){
    warning(as.vector(int))
    integrated <- NA_real_
  } else {
    integrated <- int$value
  }
  (final <- abs(XX)*pi^(-0.5)*exp(2*XX)*integrated) }
mnel
  • 113,303
  • 27
  • 265
  • 254
  • Thanks for your solution. But I still get the same warning with the code you provided . I don't know what is going on here. – Jd Baba Apr 09 '13 at 03:55
  • The 21st point is always NA no matter what upper limit I use. In your computer, when you plot(a,b) do you get a smooth curve ? – Jd Baba Apr 09 '13 at 04:02
  • I am using R version 2.15.2 . Which version did you try on ? – Jd Baba Apr 09 '13 at 04:05
  • @Jdbaba -- I did get an error, because the integrand divides by t (so t=0 will not be finite), see my edit. – mnel Apr 09 '13 at 04:10
  • mnel, you're correct, but the integrand probably does not diverge: the numerator goes to zero as well, so we need to look at the limit of the integrand as x approaches zero. Of course, since `integrate` is doing finite-element work, we have to make sure it doesn't try to evaluate the integral at `x=0`. L'Hospital's rule, anyone? :-0 – Carl Witthoft Apr 09 '13 at 11:51
  • 1
    In summary, try integrating from `x=1e-8` to `1` and you'll be close enough for government work. – Carl Witthoft Apr 09 '13 at 11:56
  • @CarlWitthoft Thank you for your comment. Yes, changing the x limits did work and didn't throw any errors. I hadn't tried to plot the integrand alone. Now, I see your point on plotting and looking at the results. The L'Hospital rule you suggested also makes sense. Thanks. – Jd Baba Apr 09 '13 at 13:15