0

I am trying to carry out a 1D numerical integration of the probability2DBox_free_t function here below:

prob1Dbox<-function(invL, t, invtau, x0, x, n_lim) {
  
  c =  pi * (pi/4) * (t * invtau)
  res = 0
  for(n in 1:n_lim){
    res = res + (exp(-1 * (n * n) * c) * cos((n * pi * x) * invL) * cos((n * pi * x0) * invL))
  }
  return(invL + ((2 * invL) * res))
}


probability2DBox_free_t<-function( t,  x0,  xt_pos,  y0,  yt_pos,  invL, invtau, n_lim,  De,  mut_rate) {
  
  temp_pbx = prob1Dbox(invL, t, invtau, x0, xt_pos, n_lim)
  temp_pby = prob1Dbox(invL, t, invtau, y0, yt_pos, n_lim)
  return((1 / (2 * De)) * (temp_pbx * temp_pby) * exp(-2 * t * mut_rate))
}

Fixing De=25, invL = 0.09090909, x0=2, xt_pos=5, y0=7, yt_pos=5, invtau = 0.0008264463, n_lim = 1000, when I try to integrate the function using mut_rate = 7e-3, I do get:

0.001997469 with absolute error < 1.8e-05

whereas when I try to integrate the function using mut_rate = 7e-8, I do get:

1180.58 with absolute error < 0.0072

The integration result for mut_rate = 7e-8 is not as expected (values between 0 and 1), as also point out by the absolute error. Any suggestion on how I can get convergence and decrease error?

I have been suggested to use Sparse grid methods, but the only R package I found that do that (SparseGrid R package) is not very clear on how to specify the domain of integration.

CafféSospeso
  • 1,101
  • 3
  • 11
  • 28
  • Your function `prob1Dbox` produces negative values and so your function `probability2DBox_free_t` also produces negative values. when I integrate `probability2DBox_free_t` between 0 and 1, I get a (small) negative value (not 1180.58). – G5W Jun 27 '22 at 23:59
  • Sorry, I was not clear in the question. I do integrate between 0 and Infinity, not between 0 and 1. The results I expect from the integrand should be between 0 and 1. – CafféSospeso Jun 28 '22 at 10:18

0 Answers0