# Density of states from numerical integration and then differentiation

for some reason I want to calculate following expression:

$$N(\omega) = \int \limits_{-\pi}^\pi \mathrm{d} x \int \limits_{-\pi}^\pi \mathrm{d} y \int \limits_{-\pi}^\pi \mathrm{d} z \; \delta \left( \omega - f(x,y,z) \right)$$ where $f(x,y,z) = \sqrt{\sin^2 \frac{x}{2}+\sin^2 \frac{y}{2}+\sin^2 \frac{z}{2}}$. Of course, function $f$ might be similar to this, just starting from $f = 0$ at $x=y=z=0$ and $f>0$ in the cube $-\pi \leq x,y,z \leq \pi$.

However this is not possible not analytically nor numerically due to integrand being the Dirac delta function. So instead I took another function:

$$G(\omega) = \int \limits_{-\pi}^\pi \mathrm{d} x \int \limits_{-\pi}^\pi \mathrm{d} y \int \limits_{-\pi}^\pi \mathrm{d} z \; \theta \left( \omega - f(x,y,z) \right)$$ and the following holds true: $$N(\omega) = \frac{\mathrm{d}}{\mathrm{d} \omega} G(\omega)$$

So I've setup a Mathematica code to calculate function $G$ (analytically possible only in 1D when $f(x) = \left| \sin \frac{x}{2} \right|$ and I want to see the result in 2D and 3D):

f[x_,y_,z_]:= Sqrt[Sin[x/2]^2+Sin[y/2]^2+Sin[z/2]^2];
G = Table[{omega, NIntegrate[HeavisideTheta[omega - f[x,y,z]], {x, -Pi, Pi}, {y, -Pi, Pi}, {z, -Pi, Pi}]}, {omega, 0, 1.7, 0.01}];


Now the result is pretty function $G$ which can be plotted with a ListPlot, but what I am seeking is its derivative. So I finally calculate $N$:

DoS = Table[{G[[i, 1]], (G[[i + 1, 2]] - G[[i - 1, 2]])/(G[[i + 1, 1]] - G[[i - 1, 1]])}, {i, 2, Length[G] - 1}];
ListPlot[DoS]


This isn't very nice code, but it should work, right? No. Wrong. The result is very jerky function (and making step finer is not helping). What am I doing wrong? This kind of result is unexpected. I expect it to be more or less piecewise continuous function, like the analytical result I can get in 1D. Why it is not smooth in 2D and 3D? How can I get a correct, smooth numerical solution to $N (\omega)$? Thanks.

• Have you tried replacing HeavisideTheta[] with UnitStep[]? Oct 9, 2016 at 23:56
• Another symptom of the problem: NIntegrate[...] thinks $G$ is zero for $\omega<0.5$, even though for small $\omega$ it should be close to the volume of the sphere of radius $\omega/2$.
– user484
Oct 10, 2016 at 5:26
• On the other hand, g[ω_] := RegionMeasure@ImplicitRegion[f[x, y, z] <= ω, {{x, -Pi, Pi}, {y, -Pi, Pi}, {z, -Pi, Pi}}] seems to gives better answers.
– user484
Oct 10, 2016 at 5:30
• Really? Things like this are not doable numerically? No one knows how to pull this off? Oct 11, 2016 at 21:59

Another way to "discretize" the problem is to approximate the delta function in some standard way, such as with a normal distribution or box function. I'll show an approach with a box function. One advantage is that the inner integral is easy to compute, reducing the dimension of the numerical integration by one.

Notes on the set-up: 1. I used the even symmetry of the integrand to reduce the integral to one over the cube $[0,\pi]^3$. 2. The parameter s in n[omega, s] below determines the width of the unit box approximation to the Dirac delta function; the value 1/100 seemed to be a good compromise. 3. I used Reduce to figure out the cases where the singular surface

f[x, y, z] == omega


intersects the cube and set up the limits of integration; the region of integration is restricted to the support of the inner integral. To simplify the work of Reduce, I used a == Sin[x/2], and b == Sin[y/2] to convert the trigonometric problem to a polynomial one. The function toPWCases constructs, from the output of Reduce, the integrals for the different cases arising from the intersection of the singular surface with the sides of the cube for varying omega. 4. Inactivating NIntegrate with Block[] lets us construct the code for n[omega, s], calling Reduce and toPWCases just once to break down the integration into its pieces.

ClearAll[f, n, i0, toPWCases];

f[x_, y_, z_] := Sqrt[Sin[x/2]^2 + Sin[y/2]^2 + Sin[z/2]^2];

Block[{NIntegrate, x, y, z, omega, s},
i0[x_, y_, omega_, s_] = First@Differences[
z /. Last@Solve[omega - f[x, y, z] == #, z] & /@ ({1, -1} s/2)
]/s /. ArcSin[a_] :> ArcSin@Sqrt[Clip[a^2, {0, 1}]] //
PiecewiseExpand;

toPWCases[And[cond_, xr_, yr_]] := {NIntegrate[i0[x, y, omega, s],
{x, 2 ArcSin@xr[], 2 ArcSin@xr[[-1]]},
{y, 2 ArcSin@yr[], 2 ArcSin@yr[[-1]]}], cond} /.
{a -> Sin[x/2], b -> Sin[y/2]};
toPWCases[And[cond_, ranges_Or]] :=
{Total[First@toPWCases[Hold && #] & /@ List @@ ranges], cond};

n[omega_?NumericQ, s_?NumericQ] = 8 Piecewise[
toPWCases /@ List @@ Reduce[
omega > 0 && 0 < a < 1 && 0 < b < 1 && omega^2 - (a^2 + b^2) > 0,
{a, b}]
];
]


The integration is still a bit slow, even with the reduction in dimension. It can be sped up by around 100 sec. by adding the options PrecisionGoal -> 4, AccuracyGoal -> 6 without a degradation of the plot. A look at the output of Reduce, one can see that omega == 1 and omega == Sqrt divides the cases; we see, correspondingly, a change in the behavior of n near those values. Only once does NIntegrate complain about slow convergence, but it does not complain about the error.

Dynamic@omega

SessionTime[];
Plot[n[omega, 1/100], {omega, 0.1, 1.7}, PlotPoints -> 15,
MaxRecursion -> 1, PlotRange -> All, GridLines -> {{1, Sqrt}, None}]
SessionTime[] - %% (*  255.846962  <-- Timing *)


Note MaxRecursion is set to 1 for the sake of speed; a higher setting resolves the locations of the corners in the graph, which have a slight rounding, presumably due to the finite value s = 1/100: G isn't just smooth enough.

f[x_, y_, z_] := Sqrt[Sin[x/2]^2 + Sin[y/2]^2 + Sin[z/2]^2];

G = ParallelTable[{omega,
NIntegrate[
HeavisideTheta[omega - f[x, y, z]], {x, -Pi, Pi}, {y, -Pi,
Pi}, {z, -Pi, Pi}, Method -> "LocalAdaptive"]}, {omega, 0, 1.7,
0.01}]; // AbsoluteTiming


{441.301, Null}

For a first shot consider fitting a polynomial:

FitPolynomial[data_, n_] := Fit[data, Table[x^i, {i, 0, n}], x];
g[x_] := Evaluate@FitPolynomial[G, 10]

Show[ListPlot[G], Plot[g[x], {x, 0, 1.7}, PlotStyle -> Red], Frame -> True] Note, however, that for $\omega<0.41$, $G(\omega)$ is a numerical zero.

DoS = Table[{G[[i,
1]], (G[[i + 1, 2]] - G[[i - 1, 2]])/(G[[i + 1, 1]] -
G[[i - 1, 1]])}, {i, 2, Length[G] - 1}];
Show[ListPlot[DoS],
Plot[g'[x], {x, 0, 1.7}, PlotStyle -> Red], Frame -> True] The finite difference method is implemented with an undocumented function (under Spatial Derivative Approximations):

{g1, g2} = Transpose @ G;
plot17 = ListPlot[#, Frame -> True] & @ Drop[#, -1] & @
Transpose @ {g1,
NDSolveFiniteDifferenceDerivative[Derivative, g1, g2,
DifferenceOrder -> 2]} A polynomial isn't the best choice here, but one try other smoothing methods; try with MovingAverage:

G1 = MovingAverage[G, 5];
ListPlot[G1, Frame -> True] DoS = Table[{G1[[i,
1]], (G1[[i + 1, 2]] - G1[[i - 1, 2]])/(G1[[i + 1, 1]] -
G1[[i - 1, 1]])}, {i, 2, Length[G1] - 1}];
ListPlot[DoS, Frame -> True] You can smoothen the derivative afterwards:

ListPlot[MovingAverage[DoS, 5], Frame -> True] which is a quite pleasing result.

Or smoothing and Interpolation:

gi = Interpolation @ G1

dom = gi["Domain"][]

Plot[gi'[x], {x, dom[], dom[]}, Frame -> True] To display the method suggested by Rahul in the comments:

reg = ListPlot[#, Frame -> True] &@
Table[{x, m[x]}, {x, 0.01, 1.7, 0.01}]; // AbsoluteTiming


{262.235, Null}

Faster, but

reg has an outlier and plateaus on the right.

• Use NIntegrate with Method -> "LocalAdaptive". This evaluates without any warnings. Oct 10, 2016 at 1:06
• This method causes evaluation to stall at some points, like omega = 1...but yeah, it goes on and on without warning ^^ Oct 10, 2016 at 9:20

The function G[omega] is not smooth, but you can use the following numerical differentiation to find N[omega]:

Clear[f]
ComputenumericalDerivative[f_, x_] :=
Block[{xn, xn1, vars, guess, dx = 0.0001, dfdx, derivative},
xn = x;
dfdx = (f[xn + dx] - f[xn])/dx;
{xn, dfdx}
]
(*testing*)
f[x_] := x^2 + Sin[x]
ComputenumericalDerivative[f, 1]
D[f[x], x] /. x -> 1.
Clear[f]
f[x_, y_, z_] := Sqrt[Sin[x/2.]^2. + Sin[y/2.]^2. + Sin[z/2.]^2.];
funcomega[omega_] :=
NIntegrate[
HeavisideTheta[omega - f[x, y, z]], {x, -Pi, Pi}, {y, -Pi,
Pi}, {z, -Pi, Pi}]
func = Table[{omega, funcomega[omega]}, {omega, 0, 1.7, 0.1}];
dfunc = Table[
ComputenumericalDerivative[funcomega, omega], {omega, 0, 1.7, 0.1}]
ListLinePlot[{func, dfunc}, Frame -> True, PlotMarkers -> Automatic,
PlotLegends -> {"G[\[Omega]]", "N[\[Omega]]"}]
`
• I don't think the output (obtained in 543 sec) is very pleasing: imgur.com/a/jLwfQ Oct 10, 2016 at 9:54