Numerical Triple Integration in R

Question

Is it possible to do triple integration in R without using the cubature package?

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

For example, to code this triple integral:

I tried

InnerMostFunc = function(v) { v + y^2 }
InnerMostIntegral = function(w) { sapply(w, 
   function(x) { integrate(InnerMostFunc, 1, 2)$value }) }
InnerIntegral = function(y) { sapply(y, 
   function(z){integrate(InnerMostIntegral, 1, 2)$value }) }
integrate(InnerIntegral, 0, 1)

I am having trouble reading what triple integral you are trying to compute. Can you provide the formula? Perhaps using an image of the formula like in the post that you cited. — G5W, Jun 14 '17 at 18:53
A triple integral requires three variables (and three separate ranges). Your integrations only have one variable, although I'm guessing you are intending to integrate over a square cube. — IRTFM, Jun 14 '17 at 19:08
So you cannot expect R to deliver that symbolic result, but you could get a result if you specified a value for `k`. — IRTFM, Jun 14 '17 at 19:10
@42 - ah yes I was trying to integrate over a square cube. I keep finding examples with repeated ranges though (like the most recent one I added). You really can't do a triple integral with repeated ranges? — maddie, Jun 14 '17 at 19:16
You would need to specify the ranges in some 3d structure and use them to specify the inner function calls to integrate. The new version of the question is going to require a function capable of dealing with partial derivatives and I don't think that `integrate is there yet. Could be wrong.... discarding my draft answer since new version too complex for my strategy. — IRTFM, Jun 14 '17 at 19:21
@42 so there is essentially no way for me to do triple integration using the the model I found in the other post? — maddie, Jun 14 '17 at 19:30
@G52 I've tried your solution from the other post and it works for cube functions? — maddie, Jun 14 '17 at 19:57
integrate(function(z) { sapply(z, function(z){ integrate(function(y) { sapply(y, function(y) { integrate(function(x) myfun(x,y,z), llim, ulim)$value }) }, llim, ulim)$value }) }, llim, ulim) — maddie, Jun 14 '17 at 19:57
@GW5 it's basically integrating where y and z is kept as a constant via in the inner function, and then y is applied as a variable in the first outer function keeping only z as constant and it hits the outer application loop which applies the last layer — maddie, Jun 14 '17 at 19:59
I think that what you said is right. It will work for triple integral where the inner function has only one variable, but the limits depend on the others. Would it be helpful if I wrote up an example as an answer? — G5W, Jun 14 '17 at 20:10

G5W · Accepted Answer · 2017-06-14T22:01:17.690

Further down is the extension of the cited previous post that was asked for in the comments. But this top part will show how to compute the integral given in the updated post.

This is harder to write than the type of integral below, because the functions must be completely nested. You can't separate them out as in the second example, so the single expression is rather long and hard to read. Here is the code to compute the requested integral.

integrate(Vectorize(function(x) { 
    integrate(Vectorize(function(y) { 
        integrate(function(z) { x^2 + y*z }, 1, 2)$value }), 1,2)$value }), 0,1)
2.583333 with absolute error < 2.9e-14

Notice that it computes the correct answer 31/12. The cited source of the integral incorrectly gives the answer as 31/4.

Extension of referenced previous post

Here is an example of extending the process from the previous post to a triple integral. The example that I will compute is simple to do analytically, so we can check that we are getting the correct answer. I will use:

In order to make it easier to follow, I am breaking the code up into several steps.

InnerFunc = function(x) { x + 1 }
InnerIntegral = Vectorize(function(y) { integrate(InnerFunc, 0, y)$value})
Integral2     = Vectorize(function(z) { integrate(InnerIntegral , 0, z)$value})
(Tripleintegral = integrate(Integral2 , 0, 4))
21.33333 with absolute error < 2.4e-13

This extends the previous example, but does not cover the type of integral in the new statement of the problem.

Bravo! And with no changes of variables. I didn't think that was possible in such a compact form. But it did need explicit ranges. — IRTFM, Jun 15 '17 at 00:27
Color me "puzzled". If the upper limit of the innermost integral had been 5, you wouldn't needed to have changed anything? — IRTFM, Jun 15 '17 at 00:42

score 3 · Answer 2 · answered Jan 03 '18 at 09:17

3

For such an integral, the cubature package is the way to go.

> library(cubature)
> hcubature(function(v) v[1]^2+v[2]*v[3], c(0,1,1), c(1,2,2))
$integral
[1] 2.583333

answered Jan 03 '18 at 09:17

Stéphane Laurent

75,186
15
119
225

Numerical Triple Integration in R

2 Answers2

Linked