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

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.

  • $\begingroup$ Have you tried replacing HeavisideTheta[] with UnitStep[]? $\endgroup$ Oct 9, 2016 at 23:56
  • $\begingroup$ 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$. $\endgroup$
    – user484
    Oct 10, 2016 at 5:26
  • 1
    $\begingroup$ 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. $\endgroup$
    – user484
    Oct 10, 2016 at 5:30
  • $\begingroup$ Really? Things like this are not doable numerically? No one knows how to pull this off? $\endgroup$
    – user16320
    Oct 11, 2016 at 21:59

3 Answers 3


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

 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.


Plot[n[omega, 1/100], {omega, 0.1, 1.7}, PlotPoints -> 15, 
 MaxRecursion -> 1, PlotRange -> All, GridLines -> {{1, Sqrt[2]}, None}]
SessionTime[] - %%

Mathematica graphics

(*  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:

Mathematica graphics


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, 
      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]

enter image description here

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}];
  Plot[g'[x], {x, 0, 1.7}, PlotStyle -> Red], Frame -> True]

enter image description here

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

enter image description here

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]

enter image description here

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]

enter image description here

You can smoothen the derivative afterwards:

ListPlot[MovingAverage[DoS, 5], Frame -> True]

enter image description here

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]

enter image description here

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


enter image description here

has an outlier and plateaus on the right.

  • 1
    $\begingroup$ Use NIntegrate with Method -> "LocalAdaptive". This evaluates without any warnings. $\endgroup$ Oct 10, 2016 at 1:06
  • $\begingroup$ This method causes evaluation to stall at some points, like omega = 1...but yeah, it goes on and on without warning ^^ $\endgroup$
    – user16320
    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]:

ComputenumericalDerivative[f_, x_] := 
 Block[{xn, xn1, vars, guess, dx = 0.0001, dfdx, derivative},
  xn = x;
  dfdx = (f[xn + dx] - f[xn])/dx;
  {xn, dfdx}
f[x_] := x^2 + Sin[x]
ComputenumericalDerivative[f, 1]
D[f[x], x] /. x -> 1.
f[x_, y_, z_] := Sqrt[Sin[x/2.]^2. + Sin[y/2.]^2. + Sin[z/2.]^2.];
funcomega[omega_] := 
  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]]"}]
  • $\begingroup$ I don't think the output (obtained in 543 sec) is very pleasing: imgur.com/a/jLwfQ $\endgroup$
    – corey979
    Oct 10, 2016 at 9:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.