# Negative numerical result when integrating a strictly non-negative integrand

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]


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]


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.

• 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 4 '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 5 '16 at 0:50
• 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 5 '16 at 2:16
• 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 5 '16 at 5: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 5 '16 at 13:55

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


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]


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


## 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:

• I'd prolly run //PiecewiseExpand //Simplify on res3 before feeding it to the integration, if I were the OP. – J. M. is away Jun 5 '16 at 13:57
• It is a good idea. I did not post the expanded version of res3 because that would be harder to evaluate for correctness of its definition. – Anton Antonov Jun 5 '16 at 14:01