double integral in R

Question

I'm wondering how to code that takes double integrals in R. I already referred two similar questions.

calculating double integrals in R quickly

double integration in R with additional argument

But I'm still confused how I can get my question from those answers. My question is following.

I would like to code this calculations in R.

From my hand and Wolfram alpha calculation, it becomes 16826.4. I know how to take a integral if both integrals are from exact numbers using adaptIntegrate(). But I'm not sure how to do in my case. Could you guys help me? Thank you so much in advance.

Jeremyjerez

Thank you for the comment. Yes, this is just an example. My true function that I want to take a double integral is more complicated. And I guess it's not easy to calculate it in analytical way. – Jeremyjerez

Answer 1

Let me start with the code and then step through to explain it.

InnerFunc = function(x) { x + 0.805 }
InnerIntegral = function(y) { sapply(y, 
    function(z) { integrate(InnerFunc, 15, z)$value }) }
integrate(InnerIntegral , 15, 50)
16826.4 with absolute error < 1.9e-10

The first line is very easy. We just need the function f(x) = x + 0.805 to be able to compute the inner integral.

The second step is the only thing that is tricky. It seems natural to compute the inner integral with a simpler expression function(z) { integrate(InnerFunc, 15, z)$value } and just integrate it. The problem with that is that integrate expects a vectorized function. You should be able to give it a list of values and it will return a list of values. This simple form of the first integral just works for one value at a time. That is why we need the sapply so that we can pass in a list of values and get back a list of values (the first definite integral).

Once we have this vectorized function for the inner integral, we can just pass that to integrate to get the answer.

Later Simplification
While the above sapply method worked, it is more natural to use the function Vectorize like this.

InnerFunc = function(x) { x + 0.805 }
InnerIntegral = Vectorize(function(y) { integrate(InnerFunc, 15, y)$value})
integrate(InnerIntegral , 15, 50)
16826.4 with absolute error < 1.9e-10

Forgiveness

Thank you so much @Forgiveness with detail explanation. It worked! – Jeremyjerez

@Matt I see your new question. Can you answer my comment there? What triple integral are you trying to compute? – Forgiveness

@Forgiveness Thank you for you answer, it was quite helpful. I also have to calculate some double integral. Just for being clear by taking the example of the OP, do we agree that we can write : (15 <=) x <= y <= 50 ? Thanks – Simmon

@MbrMbr Yes, this integral is over that triangular region. – Forgiveness

@Forgiveness Thanks a lot, I was pretty sure It was the case but I had some doubts for a moment. – Simmon

Answer 2

The domain of integration is a simplex with vertices (15,15), (50,15) and (15,50). Use the SimplicialCubature package:

> library(SimplicialCubature)
> S = cbind(c(15,15),c(50,15),c(15,50))
> adaptIntegrateSimplex(function(v) v[1]+0.805, S)
$integral
[1] 16826.4

$estAbsError
[1] 1.68264e-08

Perplexity

Answer 3

You can actually use the pracma::integral2() function in order to solve this integral, but some work to define it as dydx is needed:

f <- function(x,y) x+0.805
xmin <- 15; xmax<-50; 
ymin<-function(x) x; ymax<-50;
integral2(f,xmin,xmax,ymin,ymax)
$Q
[1] 16826.4

$error
[1] 1.818989e-12

This since:

More info here

En