On Mar 25, 2015, at 8:45 AM, David Winsemius wrote:
>
> On Mar 25, 2015, at 4:55 AM, <lluis.hurt...@uv.es> <lluis.hurt...@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
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
