Cauchy principal value numerical integration

Question

I'm using R's integrate function, but the integrand is giving me some trouble (it has a singularity). Consider this,

distrib <- function(x){
  dnorm(x, mean=700,sd=50)
}

phi <- function(nu, epsilon=20, maxi=Inf){

  fun <- function(nup)  {
    lambdap <- 1e7/nup
    distrib(lambdap) / (nup*(nup - nu))
  }
  # first part: from epsilon to nu - epsilon
  I1 <- integrate(fun, epsilon, nu-epsilon, rel.tol = 1e-8,
                  subdivisions = 200L)
  # second part: from nu + epsilon to Infty
  I2 <- integrate(fun, nu+epsilon, maxi, rel.tol = 1e-8, 
                  subdivisions = 200L)

  I1$value + I2$value 

}


x <- seq(1e7/500, 1e7/1000, length=200)

phi1 <- sapply(x, phi)    
phi1 <- sapply(x, phi, epsilon=5)
#Error in integrate(fun, epsilon, nu - epsilon, rel.tol = 1e-08, subdivisions = #200L) : 
#  the integral is probably divergent

Is there a better way of calculating such principal values in R, than playing with a cutoff parameter and adjusting it manually until it gives a result (not necessarily accurate)?

Answer 1

If you consider the integral on an interval centered on the singularity, you can use a change of variable to symmetrize the integrand: the singularity disappears, as detailed in the paper Numerical evaluation of a generalized Cauchy principal value (A. Nyiri and L. Bayanyi, 1999).

# Principal value of \( \int_{x_0-a}^{x_0+a} \dfrac{ f(x) }{ h(x) - h(x0) } dx \)
pv_integrate <- function( f, h, x0, a, epsilon = 1e-6 ) {
  # Estimate the derivatives of f and h
  f1 <- ( f(x0+epsilon) - f(x0-epsilon) ) / ( 2 * epsilon )
  h1 <- ( h(x0+epsilon) - h(x0-epsilon) ) / ( 2 * epsilon )
  h2 <- ( h(x0+epsilon) + h(x0-epsilon) - 2 * h(x0) ) / epsilon^2 
  # Function to integrate. 
  # We just "symmetrize" the integrand, to get rid of the singularity,
  # but, to be able to integrate it numerically, 
  # we have to provide a value at the singularity.
  # Details: http://mat76.mat.uni-miskolc.hu/~mnotes/downloader.php?article=mmn_16.pdf
  g <- function(u) 
    ifelse( u != 0,
      f( x0 - u ) / ( h(x0 - u) - h(x0) ) + f( x0 + u ) / ( h(x0+u) - h(x0) ),
      2 * f1 / h1 - f(x0) * h2 / h1^2
    )
  integrate( g, 0, a )
}

# Compare with the numeric results in the article
pv_integrate( function(x) 1,      function(x) x,    0, 1  ) # 0
pv_integrate( function(x) 1,      function(x) x^3,  1, .5 ) # -0.3425633
pv_integrate( function(x) exp(x), function(x) x,    0, 1  ) # 2.114502
pv_integrate( function(x) x^2,    function(x) x^4,  1, .5 ) # 0.1318667