Numerical Triple Integration in R

Numerical Triple Integration in R

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.

Triple integration using R when each inner integral depends on the others

I think, that basing on the suggestion from the post you have attached the key is to add the Vectorize function and we come up with something like this:

# Example function 
phi <- function(t){t}

integrate(Vectorize(function(x) { 
  integrate(Vectorize(function(y) { 
    integrate(function(z) { 
      phi(x)*phi(y)*phi(z)   ## function you want to calculate
    }, -10, 6-x-y)$value     ## 2nd level interval
   }), -5, 3-x)$value        ## 1st level interval
}), -5,4)                    ## top interval

# 16480.2 with absolute error < 1.8e-10

the above approach is quite general, so you can test on other functions. (The answer is the same as provided by WolframAlpha)

How to solve non-numeric argument.. error in numerical integration?

The first integral must be parameterized by y and z and the second by z. Then we can perform the final integration.

int1 <- Vectorize(function(y, z) integrate(fxyz, 0, 5, y = y, z = z)$value)
int2 <- Vectorize(function(z) integrate(int1, 0, 3, z = z)$value)
integrate(function(z) log(z) * int2(z), 2, 6)$value
## [1] 2071.71

---

Related Topics