Tiago,

Se for possível calcular o valor da função em qualquer ponto, então talvez você possa utilizar integração por Monte Carlo. Por exemplo:

# Cálculo de volume (1/8 de esfera de raio r)

r<-1

ux<-runif(2e4,0,r)

uy<-runif(2e4,0,r)

uz<-runif(2e4,0,r)

id<-ux^2+uy^2+uz^2<=r^2

mean(id)*r^3        # Aproximado

4*pi*r^3/(3*8)      # Exato

# Vizualização

np=200

x<-seq(0,r,l=np)

z<-outer(x,x,function(x,y) sqrt(r^2-x^2-y^2))

library(rgl)

persp3d(x,x,z)

points3d(ux,uy,uz,pch=20,col=ifelse(id==TRUE,4,2))

Att.,

Rubem