[R] How to double integrate a function in R
Hans W Borchers
hwborchers at googlemail.com
Sat Jul 27 08:32:19 CEST 2013
Tiago V. Pereira <tiago.pereira <at> mbe.bio.br> writes:
> I am trying to double integrate the following expression:
>
> # expression
> (1/(2*pi))*exp(-y2/2)*sqrt((y1/(y2-y1)))
>
> for y2>y1>0.
>
> I am trying the following approach
>
> # first attempt
>
> library(cubature)
> fun <- function(x) { (1/(2*pi))*exp(-x[2]/2)*sqrt((x[1]/(x[2]-x[1])))}
> adaptIntegrate(fun, lower = c(0,0), upper =c(5, 6), tol=1e-8)
>
> However, I don't know how to constrain the integration so that y2>y1>0.
>
> Any ideas?
> Tiago
You could use integral2() in package 'pracma'. It implements the
"TwoD" algorithm and has the following properties:
(1) The boundaries of the second variable y can be functions of the first
variable x;
(2) it can handle singularities on the boundaries (to a certain extent).
> library(pracma)
> fun <- function(y1, y2) (1/(2*pi))*exp(-y2/2)*sqrt((y1/(y2-y1)))
> integral2(fun, 0, 5, function(x) x, 6, singular=TRUE)
$Q
[1] 0.7706771
$error
[1] 7.890093e-11
The relative error is a bit optimistic, the absolute error here is < 0.5e-6.
The computation time is 0.025 seconds.
Hans Werner
More information about the R-help
mailing list