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

Leave a Comment