[R] Fast evaluation of functions in 3D domains

[David Winsemius]

On Mar 25, 2015, at 8:45 AM, David Winsemius wrote:

> 
> On Mar 25, 2015, at 4:55 AM, <Lluis.Hurtado at uv.es> <Lluis.Hurtado at uv.es> wrote:
> 
>> Dear all,
>> 
>> Finally I have tried three different options to integrate a function in a 3D volume. Here 
>> I show a test example. The volume is the box [0,100] x [0,100] x [0,100] and the 
>> function is
>> 
>> nfw(d) = 4/((d/5)*(1+(d/5))^2)

snipped quite a bit.

> 
> This is the same as one gets with considering this to be a radial symmetric function and transforming the problem to one dimension using the infinitesimal transformation dV = 4*pi*dr:

I meant to type:
  dV = 4*pi*r^2*dr:

> 
>> f <- function(x) 4*pi*x^2*4.0/((x/5)*(1+(x/5))^2)

But function definition was correct. 

>> integrate(f, 0, 40 , subdivisions=10000)
> 8220.516 with absolute error < 0.0012
> 
> I wasn't able to achieve convergence to that value using adaptIntegrate using a Boolean variation on Duncan's suggestion:
> 
> model <- function(x) 
>   {
>      d <- sqrt((x[1]-50)^2 + (x[2]-50)^2 + (x[3]-40)^2)
>      r <- 4.0/((d/5)*(1+(d/5))^2)
>      dinside <- (d <= 40)*r
>    }
>> adaptIntegrate(model, rep(0,3), rep(100,3), tol=1e-3, maxEval=1000000)
> $integral
> [1] 8199.022
> 
> $error
> [1] 33.33091
> 
> $functionEvaluations
> [1] 1000065
> 
> $returnCode
> [1] 0
> 
> Without the maxEval the method exceeded my patience. If the real problem is radial symmetric (and it might be noting your domain of investigation) this may offer speed and accuracy advantages.
> 
> --
> David.
snipped

David Winsemius
Alameda, CA, USA