16

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.

enter image description here

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.

G5W
  • 36,531
  • 10
  • 47
  • 80
moggirio
  • 161
  • 1
  • 1
  • 4
  • 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. – moggirio Apr 03 '17 at 20:35

2 Answers2

20

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
G5W
  • 36,531
  • 10
  • 47
  • 80
  • Thank you so much @G5W with detail explanation. It worked! – moggirio Apr 04 '17 at 07:43
  • @Matt I see your [new question](https://stackoverflow.com/q/44551816/4752675). Can you answer my comment there? What triple integral are you trying to compute? – G5W Jun 14 '17 at 18:55
  • @G5W 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 – Mbr Mbr Nov 19 '21 at 14:41
  • 1
    @MbrMbr Yes, this integral is over that triangular region. – G5W Nov 19 '21 at 14:47
  • @G5W Thanks a lot, I was pretty sure It was the case but I had some doubts for a moment. – Mbr Mbr Nov 19 '21 at 14:53
4

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
Stéphane Laurent
  • 75,186
  • 15
  • 119
  • 225