Numerical Triple Integration in R

Question

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

based on the answer in this post

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)

Answer 1

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.