double integral in R

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

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