1

Although there are similar posts here Negative integral of a positive function , Positive integrand giving negative answer ,

my case is even more surprising. I have a very simple function

Clear[y]
y[q_, w_] := Module[{res},
 res = If[w > q (2 - q), 0, w];
 If[w > Abs[q (2 - q)] && w < q (2 + q), res = 1 - 1/4 (q - w/q)^2];
 If[w < 0, res = 0];
 res
 ]

It looks like this

g1 = Plot[y[0.25, w], {w, 0, 4}, PlotRange -> All]

enter image description here

Numerical and analytical integrals are equal

Integrate[y[0.25, w], {w, 0, 4}]
NIntegrate[y[0.25, w], {w, 0, 4}]
(*0.0957031*)

Now let us shift the function

g2 = Plot[y[0.25, w-2], {w, 0, 4}, PlotRange -> All]

enter image description here

and integrate again in the same limits

Integrate[y[0.25, w - 2], {w, 0, 4}]
NIntegrate[y[0.25, w - 2], {w, 0, 4}]
(*-1.9043*)

The result is now negative and not equal to the previous value! What's up?

In this particular case i am not so much interested in the origin of such peculiar result, but rather in the practical prescription on how to deal with it.

Community
  • 1
yarchik
  • 18,202
  • 2
  • 28
  • 66
  • I think the result 0.0957031 is also wrong. Redefine your function as Clear[y];y[q_?NumericQ, w_?NumericQ] := ... and run both NIntegrate[y[0.25, w], {w, 0, 4}] and NIntegrate[y[0.25, w - 2], {w, 0, 4}]. – march Jun 04 '16 at 20:04
  • @march I see that I get the same result (0.124349) for the w and w-2 cases by forcing the arguments to be numeric, the question is why does that work? – Jack LaVigne Jun 05 '16 at 00:50
  • 1
    I think your function doesn't generally do what you intend. Try y[a,b]. It's strange to be using a module for this. Less strange, but more troublesome, is using If instead of Piecewise. – John Doty Jun 05 '16 at 02:16
  • 2
    Yes, what @JohnDoty said. Certainly Integrate is not going to work correctly, because y doesn't evaluate correctly symbolically, and Integrate manipulates symbolic expressions. NIntegrate is also failing probably because it does symbolic pre-processing before it does the numerical integration, and so the function doesn't evaluate correctly as an integrand. – march Jun 05 '16 at 05:35
  • @march Thank you for your suggestion. I have also explored another route by expressing everything in terms of the "UnitStep" function (this is in line suggestions of @ Jack Lavigne ). So I have solutions in both, numerical and analytical, ways. Should I delete this question as too simple? Can one, actually, learn anything from my question? – yarchik Jun 05 '16 at 13:55

1 Answers1

4

There is some confusion using Plot. Compare the plot of the non-evaluated y in the question with the plot of the evaluated y:

Grid[{{Plot[y[0.25, w - 2], {w, 0, 4}, PlotRange -> All, 
    ImageSize -> Medium],
   Plot[Evaluate@y[0.25, w - 2], {w, 0, 4}, PlotRange -> All, 
    ImageSize -> Medium]}}]

enter image description here

We can see that going through the algorithmic steps of defining the value y with symbolic arguments produces an expression that would give different values over the domain of y if using numerical arguments.

In this particular case i am not so much interested in the origin of such peculiar result, but rather in the practical prescription on how to deal with it.

Since in the question both Integrate and NIntegrate are used let us redefine y with Piecewise.

res1 = Piecewise[{{0, w > q (2 - q)}}, w]
res2 = Piecewise[{{1 - 1/4 (q - w/q)^2, 
    w > Abs[q (2 - q)] && w < q (2 + q)}}, res1]
res3 = Piecewise[{{0, w < 0}}, res2]

enter image description here

With that definition we get positive results for the integrals in the question.

Block[{q = 0.25},
 NIntegrate[res3, {w, 0, 4}]
]

(* 0.124349 *)

Integrate[res3, {w, 0, 4}]

enter image description here

Using IntegrationMonitor

Using the option IntegrationMonitor can help in these kind of situations. In this case we clearly see that the integrant is not what it is expected:

enter image description here

Community
  • 1
Anton Antonov
  • 37,787
  • 3
  • 100
  • 178