4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Have you tried replacing HeavisideTheta[] with UnitStep[]? $\endgroup$ – J. M. will be back soon Oct 9 '16 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$ – Rahul Oct 10 '16 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$ – Rahul Oct 10 '16 at 5:30
  • $\begingroup$ Really? Things like this are not doable numerically? No one knows how to pull this off? $\endgroup$ – user16320 Oct 11 '16 at 21:59
4
$\begingroup$

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[] - %%

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

$\endgroup$
0
$\begingroup$

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]

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}];
Show[ListPlot[DoS], 
  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

reg

enter image description here

has an outlier and plateaus on the right.

$\endgroup$
  • 1
    $\begingroup$ Use NIntegrate with Method -> "LocalAdaptive". This evaluates without any warnings. $\endgroup$ – Julien Kluge Oct 10 '16 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 '16 at 9:20
-1
$\begingroup$

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]]"}]
$\endgroup$
  • $\begingroup$ I don't think the output (obtained in 543 sec) is very pleasing: imgur.com/a/jLwfQ $\endgroup$ – corey979 Oct 10 '16 at 9:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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