# NIntegrate returns zero for non zero integrand, 4d Integration

I am trying to calculate the following integral.

sigma1 = 10.0; sigma2 = 5.0; delta = 0.5;
t[x1_, y1_, x_, y_] := 100*HeavisideLambda[sigma1^-1*(x - x1), sigma2^-1*(y - y1)];
B2[x1_, y1_, x_, y_] := HeavisideTheta[(delta/2)^2 - (x - x1)^2, (delta/2)^2 - (y - y1)^2];
trans[x1_, y1_, x2_, y2_] :=
NIntegrate[B2[x1, y1, xz, yz]*t[xz, yz, xp, yp]*
(B2[x2, y2, xz, yz] - B2[x2, y2, xp, yp]),
{xp, x2 - 2.0*sigma1, x2 + 2.0*sigma1},
{yp, y2 - 2.0*sigma2, y2 + 2.0*sigma2},
{xz, x1 - 0.5*delta, x1 + 0.5*delta},
{yz, y1 - 0.5*delta, y1 + 0.5*delta},
WorkingPrecision -> 12, AccuracyGoal -> 8, MinRecursion -> 8, MaxRecursion -> 100];


I am interested in the value of the integral for the following inputs:

trans[0, 0, delta, 0]


Here is my problem: for values of delta greater than 0.5 (I tried 1.0, 0.9, 0.8, 07, 0.6, 0.515), the result is a negative number, and it takes some time for Mathematica to come up with the result. For any value of delta smaller than 0.5, Mathematica immediately returns 0, and doesn't give any hints about what is wrong.

This is a part of a refinement study and I need to choose smaller and smaller values for delta. Do you know how I can make this work?

-
works for me: trans[0,0,1/4,0] -> 624.99 . You had some strange character instead of a semicolon after delta=0.5 by the way. – george2079 Jun 17 '14 at 17:40
yes that works for me too, but what i need is trans[0,0,delta,0]. Does that work for you? – amirdel Jun 17 '14 at 19:16

The Heaviside functions are essentially piecewise functions, and NIntegrate knows how to handle Piecewise functions but not Heaviside functions. In particular, it will analyze the domain of Piecewise functions and adjust its sampling accordingly. Here are two rules for conversion, ignoring boundary points which won't affect the integral anyway:

heaviside2piecewise = {
HeavisideTheta -> (Piecewise[{{1, #1 > 0 && #2 > 0}}, 0] &),
HeavisideLambda -> (Piecewise[{{#1 + 1, -1 < #1 < 0}, {1 - #1, 0 <= #1 < 1}}, 0] *
Piecewise[{{#2 + 1, -1 < #2 < 0}, {1 - #2, 0 <= #2 < 1}}, 0] &)};


Then we can apply them to the integrand:

trans[x1_, y1_, x2_, y2_] := NIntegrate[
B2[x1, y1, xz, yz] * t[xz, yz, xp, yp] * (B2[x2, y2, xz, yz] - B2[x2, y2, xp, yp]) /.
heaviside2piecewise,
{xp, x2 - 2.0*sigma1, x2 + 2.0*sigma1},
{yp, y2 - 2.0*sigma2, y2 + 2.0*sigma2},
{xz, x1 - 0.5*delta, x1 + 0.5*delta},
{yz, y1 - 0.5*delta, y1 + 0.5*delta}]


It evaluates rather quickly, too:

trans[0, 0, delta, 0]
(* -5.73958 *)

Needs["GeneralUtilities"]; (* V10 *)
trans[0, 0, delta, 0] // AccurateTiming
(* 0.0757289 *)


Exact solution

In fact, this integral may be solved exactly, with exact values for the parameters.

transExact[x1_, y1_, x2_, y2_] := Integrate[
B2[x1, y1, xz, yz] * t[xz, yz, xp, yp] * (B2[x2, y2, xz, yz] - B2[x2, y2, xp, yp]) /.
heaviside2piecewise,
{xp, x2 - 2*sigma1, x2 + 2*sigma1},
{yp, y2 - 2*sigma2, y2 + 2*sigma2},
{xz, x1 - 1/2*delta, x1 + 1/2*delta},
{yz, y1 - 1/2*delta, y1 + 1/2*delta}];

transExact[0, 0, delta, 0]
(* -(551/96) *)

-

ah, you are right for that condition: B2[x2, y2, xz, yz] and B2[x2, y2, xp, yp] are returning 0 for every sample point in the domain. This is not so much an answer but I thought i may be useful to show how to use EvaluationMonitorand Reap/Sow to get at this:

 trans[x1_, y1_, x2_, y2_] :=
NIntegrate[B2[x1, y1, xz, yz]*t[xz, yz, xp, yp]*
((b21 = B2[x2, y2, xz, yz]) - (b22 = B2[x2, y2, xp, yp])),
{xp, x2 - 2.0*sigma1, x2 + 2.0*sigma1},
{yp, y2 - 2.0*sigma2, y2 + 2.0*sigma2},
{xz, x1 - 0.5*delta, x1 + 0.5*delta},
{yz, y1 - 0.5*delta, y1 + 0.5*delta},
WorkingPrecision -> 12, AccuracyGoal -> 8, MinRecursion -> 8,
MaxRecursion -> 100, EvaluationMonitor :> Sow[{b21, b22}]];

(First@Last@Reap[trans[0, 0, delta, 0]])


result 14,000 {0,0}'s

I didn't get so far as studying why this happens, but use of the same the symbol x1 for different things makes it a little confusing to read.

## edit

looking at this a little deeper, i suspect you may have a situation where the integrand is zero everywhere except for a small region. NIntegrate samples over a pretty coarse grid, doesn't happen to hit any non-zero points and stops. Try some of the other integration Methods.( MinRecursion` doesnt seem to do anything BTW ) You might also see if you can id the nonzero region and tighten the integration limits.

-
I think the integral limits are already as tight as possible, Based on the shape of "t" and "B2". I also think that NIntegrate doesn't sample point where it should, but I don't know how to make it do that. Do you have any specific method in mind? – amirdel Jun 17 '14 at 23:14