# Integrating a function over a surface integral

From a first principles bandstructure calculation I get an energy scalar field in three dimensions $E(x,y,z)$. It's now easy to plot a constant energy (contour)-surface for dedicated values $\epsilon=E(x,y,z)$ by using ListContourPlot3D. Lets assume the following example scalar field (e.g. the case of a graviatation potential field):

data = Table[1/Sqrt[x^2 + y^2 + z^2], {x, -1, 1, 0.1}, {y, -1, 1, 0.1}, {z, -1, 1,0.1}]
/. ComplexInfinity -> 10^16 // Quiet;


I exchanged the infinite value at $x=y=z=0$ with a finite one to prevent any error messages coming from the visualisation routines. Now the iso-surfaces can be easily plotted for e.g. $\epsilon=2.0$

ListContourPlot3D[data, Contours -> {2.0}, DataRange -> {{-1, 1}, {-1, 1}, {-1, 1}}]


Coming back to the bandstructure calculation, the energy surfaces are much more complicated and cannot be described analytically for the general case as can be seen from the following plot of an energy surface for Silicon:

An important electronic property, the DOS (Density of States) is calculated through the integral $$\oint \limits_{E(x,y,z)=\epsilon}\frac{dS}{\lvert\nabla E(x,y,z)\rvert}.$$

and additional integrals of interest have the general form

$$\oint \limits_{E(x,y,z)=\epsilon}\frac{dS}{\lvert\nabla(\Delta E(x,y,z))\rvert}A(x,y,z)$$

with $A(x,y,z)$ as some additional scalar function (e.g. for effective mass tensor, $A(x,y,z)$ is the tensor product of the gradients of the energy).

Now for the question. How to most efficiently implement such integrals on numerically specified data (see above) over iso-surfaces in Mathematica?

EDIT: Trying out the proposal by Mastrok I run into convergence problems for e.g. $\epsilon$ = 2.0 of the integration (most probably due to the singularity in the coordinate center and the additional derivatives on the interpolation function when using Gauss' Theorem):

e = ListInterpolation[data, {{-1, 1}, {-1, 1}, {-1, 1}}];
dive = Simplify[Div[#/Norm[#]^2& @ Grad[e[x, y, z], {x, y, z}], {x, y, z}] /.
Derivative[1][Abs][x_] :> x/Abs[x] /. {x->#1, y->#2, z->#3} &];
NIntegrate[Piecewise[{{dive[x, y, z], e[x, y, z] < 2.0}}, 0] ,
{x, -1, 1}, {y, -1, 1}, {z, -1, 1}
]

(*Output*)
NIntegrate::slwcon: Numerical integration converging too slowly;
suspect one of the following: singularity, value of the integration is 0,
highly oscillatory integrand, or WorkingPrecision too small.

• If I'm not mistaken, the integral is equal to $\displaystyle\frac{\mathrm d}{\mathrm d\epsilon}\iiint_{E(x,y,z)\le\epsilon}\mathrm dV$, is that right? It might be easier to compute that instead, but I'm not sure. – Rahul Apr 23 '15 at 6:07
• Could you add to the Question a plot of a typical complicated energy surface over which you wish to integrate? – bbgodfrey Apr 23 '15 at 6:12
• The literal implementation Integrate[1/Norm[Grad[e[x, y, z], {x, y, z}]], {x, y, z} ∈ ImplicitRegion[e[x, y, z] == ϵ, {x, y, z}]] seems to work when $E$ is given analytically, but not when it is an InterpolatingFunction constructed from the data. – Rahul Apr 23 '15 at 6:32
• @bbgodfrey: Example added – Rainer Apr 23 '15 at 6:33
• @Rahul, I did not check if the Gauss theorem yields an easier representation through a volume integral. To my view one would just exchange the problem of integrating over a complicated surface by the problem of integrating over a volume integral with bounds of similiar complexity. – Rainer Apr 23 '15 at 6:34

First note that,

$$\oint_{E(\vec{x})=\epsilon}\frac{dS}{|\nabla E(\vec{x})|} =\oint_{E(\vec{x})=\epsilon}\frac{\nabla E(\vec{x})\cdot d\vec{S}}{|\nabla E(\vec{x})|^2} = \int_{E(\vec{x})\leq\epsilon} \nabla\cdot\left( \frac{\nabla E(\vec{x})}{|\nabla E(\vec{x})|^2}\right) dV$$

By now, the volume integral can be evaluated by

Integrate[f Boole[E[x,y,z]<epsilon],{x,xmin,xmax},{y,ymin,ymax},{z,zmin,zmax}]


where $f$ is the integrant and

 Boole[E[x,y,z]<epsilon]


forces the integrant to be 0 outside the region. Therefore, you need to have large enough {x,xmin,xmax},{y,ymin,ymax},{z,zmin,zmax} to cover your whole domain.

• Note also that if you really want to do surface integral numerically, you can use $\delta$-function to impose the surface to be integrated. For the domain problem, since you are dealing with lattice (I guess), then it is much easier since you just need to set the domain to be one lattice size. – mastrok Apr 23 '15 at 6:38
• thanks for your comprehensive answer. However when adopting your method to my example above, the numerical integration does not converge. I assume that the singularity in the coordinate center is causing this problem (see update in my question on this topic). – Rainer Apr 26 '15 at 8:36

By using the output of ListContourPlot3D and the new Mathematica 10.0 feature DiscretizeGraphics, one can nicely generate a meshed contour region which is suitable for NIntegrate. We can show this for the above example for energy contours from 1.0 to 2.0 with a step width of 0.1:

Monitor[Table[
e = ListInterpolation[data, {{-1, 1}, {-1, 1}, {-1, 1}}];
f = Simplify[1/Norm[Grad[e[x, y, z], {x, y, z}]]];
plot = ListContourPlot3D[data, Contours -> {ϵ},
DataRange -> {{-1, 1}, {-1, 1}, {-1, 1}}];
R = DiscretizeGraphics[Normal[plot /. (Lighting -> _) :> Lighting -> Automatic]];
NIntegrate[f, {x, y, z} ∈ R], ϵ, 1.0, 2.0, 0.1}], ϵ]
, {ϵ, 1.0, 2.0, 0.1}], ϵ]
(*Output*)
{12.5221, 8.55496, 6.03428, 4.38027, 3.25017, 2.46631, 1.90372, 1.49193, 1.18722, 0.952876, 0.777326}


The result is very close to the analytical solution. The differences are coming from the low sampling rate of the original data generation:

anasol = Table[
Integrate[
Simplify[1/Norm[Grad[1/Sqrt[x^2 + y^2 + z^2], {x, y, z}]]],
{x, y, z} ∈
ImplicitRegion[x^2 + y^2 + z^2 == (1/ϵ)^2, {x, y, z}]],
{ϵ, 1.0, 2.0, 0.1}]
(*Output*)
{12.5664, 8.583, 6.06017, 4.39984, 3.27113, 2.48225, 1.91748, 1.50458, 1.19707, 0.964263, 0.785398}