89267

Solving an integral in R gives error “The integral is probably divergent”

Question:

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

<img alt="enter image description here" class="b-lazy" data-src="https://i.stack.imgur.com/Bka7l.png" data-original="https://i.stack.imgur.com/Bka7l.png" src="https://etrip.eimg.top/images/2019/05/07/timg.gif" />

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.

Answer1:

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) } <hr /><h2>Edit (facepalm moment)</h2>

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

Recommend

  • ValueError when using if commands in function
  • why sympy count_ops() fail on this result from integration?
  • Google Street View container inside Shiny Application
  • Apply kurtosis to a distribution in python
  • Elmah not logging 404 (missing files / images)
  • Login failed for user ''
  • Verbose log abbriviations meaning in SVC, scikit-learn
  • double precision error when converting to scientific notation
  • Can I disable IE compatibility mode only for content within a ?
  • web shop (shopping cart) on google app engine
  • Can't run Appium tests on iOS 10 on real device
  • What's an elegant way of accessing parent controller's member from child controller?
  • OAuth and the YouTube API
  • How can I create a plugin mechanism that calls functions only when the plugin is available?
  • RabbitMQ java client stops consuming messages
  • Using SWIG with a build system [closed]
  • Adding directive inside the directive programatically
  • Vigenere cipher not working
  • PHP multiple file uploads
  • Using same constraints in multiple classes
  • Upload file that is in the cpan database
  • How to pass nginx proxy url for socket
  • Error processing multiple files
  • Android cannot disable cut copy paste
  • Servlet stops working on Tomcat server after some hits or time
  • hide missing dates from x-axis ggplot2
  • blade.php method outputting it's result to the form
  • Thread 1: EXC_BAD_ACCESS (code =1 address = 0x0)
  • Query to find the duplicates between the name and number in table
  • JSON response opens as a file, but I can't access it with JavaScript
  • Avoid links criss cross / overlap in d3.js using force layout
  • Cannot connect to cassandra from Spark
  • Regex thinks I'm nesting, but I'm not
  • Bug in WPF DataGrid
  • TFS: Get latest causes slow project reloading
  • Possible to stop flickering java tooltip in heavyweight mode?
  • Javascript Callbacks with Object constructor
  • In LanguageTool, how do you create a dictionary and use it for spell checking?
  • How to make Safari send if-modified-since header?
  • JTable with a ScrollPane misbehaving