Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
    
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 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 at 19:16

2 Answers 2

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) *)
share|improve this answer

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.

share|improve this answer
    
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 at 23:14

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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