# NIntegrate with a function defined with conditions

Good afternoon,

I am trying to integrate a function that is defined with a conditional. The simplified code is the following:

Area[J_, L_] := Block[ {int},
If[L < 0,
int = L*NIntegrate[Max[J x, x^2], {x, 0, 1}] ,
int = L*NIntegrate[Min[J x, x^2], {x, 0, 1}]
];
Return[Abs[int]]
]

Vol[J_] :=
NIntegrate[Area[J, L], {L, -1, 1}, Method -> "QuasiMonteCarlo",
AccuracyGoal -> 20, PrecisionGoal -> 20]


The Area[J,L] function is working properly, in fact i.e. Area[2,L] can be plotted for a range of L and it is a normal function. However, its integral does not work: Vol[2] does not give any result. I am assumig that I am not calling properly the NIntegrate command.

How can I define the NIntegrate so it works with conditional functions?

Any help is welcome. Thank you.

• You have to prevent symbolic evaluation of Area. I had to use a different symbol name in testing since Area is protected in my version of Mathematica but something like area[J_?NumberQ, L_?NumberQ] :=... does the trick. Using Vol[J_?NumberQ] :=... is also a good a idea. To get a numerical result fast I removed the options for NIntegrate in vol. – N0va Oct 21 '20 at 13:37
• Dear @N0va, thanks for your answer. I haver tried adding "?NumberQ" in the definition of the functions, however the vol function does not work. Mathematica shows the following error: NIntegrate::inumr: The integrand |Int| has evaluated to non-numerical values for all sampling points in the region with boundaries (0 1). >> – TsuKi Oct 21 '20 at 17:29

Clear["Global*"]


Rather than using If, I recommend that you use Piecewise

area[J_, L_] =
Piecewise[{{Abs[L*Integrate[Max[J x, x^2], {x, 0, 1}]], L < 0}},
Abs[L*Integrate[Min[J x, x^2], {x, 0, 1}]]]


Vol[J_] = Integrate[area[J, L], {L, -1, 1}] // Simplify


Vol[-J] == Vol[J] // Simplify

(* True *)


Plotting,

Plot[Vol[J], {J, 0, 5}, AxesOrigin -> {0, 0}]


• Dear @BobHanlon, thanks for your answer. I have tried this code and it works fine for a simple analytical integration, but it does not work with numerical integration. I have presented a simplified code, and the actual one has more difficult functions and more variables, so I think I really need to rely on numerical integration. I will try to see if there is a possibility to implement this solution. – TsuKi Oct 22 '20 at 18:05

Fix area[] (_?NumericQ, single NIntegrate):

ClearAll[area, vol, vol2];

area[J_?NumericQ, L_?NumericQ] :=
NIntegrate[
Abs[L]*If[L < 0,
Max[J x, x^2],
Min[J x, x^2]
],
{x, 0, 1}];


Iterated integral version, with singularity L == 0 (not as efficient as a multiple integral):

vol[J_?NumericQ] := NIntegrate[area[J, L], {L, -1, 0, 1}];

vol[2] // RepeatedTiming

(*  {0.18, 0.666667}  *)


Multiple integral version (requires some tricky code generation, but much faster):

vol2[J_?NumericQ] := Block[{L},
Block[{NIntegrate = Hold},
Block[{NumericQ = True &},
area[J, L] /. DownValues[area]]
] /. Hold[f_, rest___] :> NIntegrate[f, {L, -1, 0, 1}, rest]
];

vol2[2] // RepeatedTiming

(*  {0.0028, 0.666667}  *)


Note: The single NIntegrate in area[] facilitates the code generation in vol2[]; otherwise, it is not strictly necessary.

Regarding my comment under OPs question here a solution using nested numerical integration

area[J_?NumberQ,L_?NumberQ]:=Block[{int},If[L<0,int=L*NIntegrate[Max[J x,x^2],{x,0,1}],int=L*NIntegrate[Min[J x,x^2],{x,0,1}]];Return[Abs[int]]]
volume[J_?NumberQ]:=NIntegrate[area[J,L],{L,-1,1}]
areaBobHanlon[J_,L_]=Piecewise[{{Abs[L*Integrate[Max[J x,x^2],{x,0,1}]],L<0}},Abs[L*Integrate[Min[J x,x^2],{x,0,1}]]];
volumeBobHanlon[J_]=Integrate[areaBobHanlon[J,L],{L,-1,1}]//Simplify;

Table[{J,volume[J]},{J,-5,5,1/2}];
Show[{Plot[volumeBobHanlon[J],{J,-5,5},Frame->True,GridLines->Automatic,FrameLabel->{"J","volume[J]"}],ListPlot[%]}]


for comparison I included the analytical solution presented by Bob Hanlon. The code above should produce this plot

EDIT: Regarding the extension to higher dimensions mentioned in the comments under this answer: yes it is possible to extend this code but the way the code is written makes it highly ineffective. I would strongly suggest using piecewise functions and either integrating analytically (as pointed out by Bob Hanlon) or numerically without nesting:

dArea2dx[J_,L1_,L2_,x_]:=Abs[Piecewise[{{L1*Max[J x,x^2],L1<0}},L1*Min[J x,x^2]]-Piecewise[{{L2*Max[J x^2,x],L2<0}},L2*Min[J x^2,x]]]

Mathematica is quite good in integrating over Piecewise functions with simple values and seems to handle the Min and Max in this problem quite well. If high dimensional integration becomes necessary one should avoid unnecessary nesting of integrals or if this is for some reason unavoidable or desired (this problem does not call for it) one should decrease the PrecisionGoal and AccuracyGoal to speed up computation.
• area[J_?NumberQ, L1_?NumberQ, L2_?NumberQ] := Block[{int1, int2}, If[L1 < 0, int1 = L1*NIntegrate[Max[J x, x^2], {x, 0, 1}], int1 = L1*NIntegrate[Min[J x, x^2], {x, 0, 1}]]; If[L2 < 0, int2 = L2*NIntegrate[Max[J x^2, x], {x, 0, 1}], int2 = L2*NIntegrate[Min[J x^2, x], {x, 0, 1}]]; Return[Abs[int1 - int2]]] volume[J_?NumberQ] := NIntegrate[area[J, L1, L2], {L1, -1, 1}, {L2, -1, 1}] If I plot the area function it shows the correct plot ContourPlot[area[1, L1, L2], {L1, -1, 1}, {L2, -1, 1} However if I evaluate volume[1] the programm gets stuck. Thanks! – TsuKi Oct 22 '20 at 16:19
• I addressed the extension to higher dimensions in an edit of the answer.The code in the comment above is valid but to run it in a feasible time one would need to decrease the PrecisionGoal and AccuracyGoal of the NIntegrate`s involved. The solution presented in the answer circumvents those issues by means of analytic and/or unnested numerical integration. – N0va Oct 22 '20 at 19:44