0
$\begingroup$

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.

Improve this question
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – N0va
    Oct 21, 2020 at 13:37
  • $\begingroup$ 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). >> $\endgroup$
    – TsuKi
    Oct 21, 2020 at 17:29

3 Answers 3

Reset to default
1
$\begingroup$ 
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}]]]

enter image description here

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

enter image description here

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

(* True *)

Plotting,

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

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – TsuKi
    Oct 22, 2020 at 18:05
1
$\begingroup$

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.

Improve this answer
$\endgroup$
0
$\begingroup$

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

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]]]
volume2[J_?NumberQ]:=NIntegrate[dAreadx[J,L1,L2,x],{x,0,1},{L1,-1,1},{L2,-1,1}]
volume2Analytical=Integrate[dAreadx[J,L1,L2,x],{x,0,1},{L1,-1,1},{L2,-1,1}]

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

which results (after a few seconds of computation time) in

volume2

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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Dear @N0va: I have tried your code, it works fine and I get the given plot, thank you. If you don't mind a further question, is it possible to extend this code to a two-variable integral? Something like this: $\endgroup$
    – TsuKi
    Oct 22, 2020 at 16:17
  • $\begingroup$ 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! $\endgroup$
    – TsuKi
    Oct 22, 2020 at 16:19
  • $\begingroup$ 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 NIntegrates involved. The solution presented in the answer circumvents those issues by means of analytic and/or unnested numerical integration. $\endgroup$
    – N0va
    Oct 22, 2020 at 19:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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