Density of states from numerical integration and then differentiation

Question

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.

Answer 1

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[[1]], 2 ArcSin@xr[[-1]]},
     {y, 2 ArcSin@yr[[1]], 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[2] 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[2]}, 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:

Answer 2

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, 
    NDSolve`FiniteDifferenceDerivative[Derivative[1], 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"][[1]]

Plot[gi'[x], {x, dom[[1]], dom[[2]]}, 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.

Answer 3

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]]"}]