Decreasing computation time for numerical plot

My code is as follows:

M = 1.876;
m = 0.9389;
q = Sqrt[Q2 + ν^2];
E3[p3_] := Sqrt[p3^2 + m^2]
E4[p3_] := Sqrt[p3^2 + q^2 - 2*p3*q*Cos[θ] + m^2]
p3zero1 := p3 /.FullSimplify[Solve[M + ν - E3[p3] - E4[p3] == 0, p3]][[1]]
p3zero2 := p3 /.FullSimplify[Solve[M + ν - E3[p3] - E4[p3] == 0, p3]][[2]]
f[p3_] := M + ν - E3[p3] - E4[p3]
h[p3_] :=
Abs[-(p3/Sqrt[m^2 + p3^2]) -
(2 p3 - 2 Sqrt[Q2 + ν^2] Cos[θ]) /
(2 Sqrt[m^2 + p3^2 + Q2 + ν^2 -2 p3 Sqrt[Q2 + ν^2] Cos[θ]])]
AngularIntegrand[p1_, p2_] :=
Piecewise[{
{p1^2/(M^2 (M - 2*E3[p1])^2) 1/(E3[p1] E4[p1])
HeavisideTheta[p3zero1]/h[p3zero1] +
p2^2/(M^2 (M - 2*E3[p2])^2) 1/(E3[p2] E4[p2])
HeavisideTheta[p3zero2]/h[p3zero2],
Im[p3zero1] == 0 && Im[p3zero2] == 0},
{0, Im[p3zero1] != 0 || Im[p3zero2] != 0 }}]
Plot3D[
NIntegrate[AngularIntegrand[p3zero1, p3zero2]*Sin[θ], {θ, 0, π}],
{ν, 0, 3}, {Q2, 0, 6}]


The computation time on my laptop is approximately 4 hours, but I was hoping to decrease this to less than an hour.

• Your code is missing several definitions without which the numerical computation can't proceed. – Jens Sep 28 '15 at 21:17
• My apologies. I have added the necessary definitions. – T-Ray Sep 28 '15 at 21:46

M = 1.876 // Rationalize
m = 0.9389 // Rationalize;
q = Sqrt[Q2 + ν^2];
E3[p3_] = Sqrt[p3^2 + m^2];
E4[p3_] = Sqrt[p3^2 + q^2 - 2*p3*q*Cos[θ] + m^2];


For p3zero1 and p3zero2 use Set rather than SetDelayed so that the Solve and FullSimplify are done just once and not for each call.

p3zero1 = FullSimplify[
p3 /. Solve[M + ν - E3[p3] - E4[p3] == 0, p3][[1]]];

p3zero2 = FullSimplify[
p3 /. Solve[M + ν - E3[p3] - E4[p3] == 0, p3][[2]]];

f[p3_] = M + ν - E3[p3] - E4[p3];

h[p3_] = Abs[-(p3/Sqrt[m^2 + p3^2]) - (2 p3 -
2 Sqrt[Q2 + ν^2] Cos[θ])/(2 Sqrt[
m^2 + p3^2 + Q2 + ν^2 - 2 p3 Sqrt[Q2 + ν^2] Cos[θ]])];


Since AngularIntegrand is called by NIntegrate its definition should be restricted to numeric arguments.

AngularIntegrand[p1_?NumericQ, p2_?NumericQ] :=
Piecewise[{{p1^2/(M^2 (M - 2*E3[p1])^2) 1/(E3[p1] E4[p1]) HeavisideTheta[
p3zero1]/h[p3zero1] +
p2^2/(M^2 (M - 2*E3[p2])^2) 1/(E3[p2] E4[p2]) HeavisideTheta[p3zero2]/
h[p3zero2], Im[p3zero1] == 0 && Im[p3zero2] == 0}, {0,
Im[p3zero1] != 0 || Im[p3zero2] != 0}}]

Plot3D[
NIntegrate[
AngularIntegrand[p3zero1, p3zero2]*Sin[θ],
{θ, 0, π},
MinRecursion -> 4],
{ν, 0, 3}, {Q2, 0, 6},
ClippingStyle -> None,
ImageSize -> 360] //
AbsoluteTiming // Column


NIntegrate::izero: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option. If value of integral may be 0, specify a finite value for the AccuracyGoal option. >>

The Plot3D took less than three minutes on my MacBook Pro.

• It looks great. I had a few questions as I'm still learning however. 1.) You rationalized the values for M and m, was this to further decrease the computation time? 2.) Did my use of setdelay cause the variable to be called and reset at every use? If you don't have time, I understand. Thanks either way! – T-Ray Sep 29 '15 at 2:16
• @JohnTerry - 1.) I tend to rationalize constants so that I don't have to be concerned about whether I specify higher precision in subsequent numerical algorithms (e.g., NIntegrate) or graphics (e.g., Plot3D) and it can sometimes help with Simplify and efficiency. 2.) The right hand side (RHS) of SetDelayed is evaluated each time the left hand side is used. Consequently, if the RHS is time-consuming, you only want to evaluate if necessary. Your Solve and FullSimplify need only be done once so Set is appropriate. Whereas SetDelayed is necessary with AngularIntegrand. – Bob Hanlon Sep 29 '15 at 2:48
• Passing the option Method->{Automatic,“SymbolicProcessing“->0} to NIntegrate might also give additional improvement in terms of time. – Lukas Sep 29 '15 at 7:07

I am able to gain ~25% speedup compared to the method suggested by @BobHanlon by slightly modifying two things

• make AngularIntegrand remember its values (see point 5 of this tutorial) so that they are not unnecessarily recomputed for same input values
• Pass the option Method->{Automatic,"SymbolicProcessing->0} to the NIntegrate in Plot3D to prevent it from preprocess the integrand symbolically

That is, the last part of the code now looks like this:

AngularIntegrand[p1_?NumericQ, p2_?NumericQ] :=
AngularIntegrand[p1, p2] =
Piecewise[{{p1^2/(M^2 (M - 2*E3[p1])^2) 1/(E3[p1] E4[p1])HeavisideTheta[p3zero1]/h[p3zero1] + p2^2/(M^2 (M - 2*E3[p2])^2) 1/(E3[p2] E4[p2]) HeavisideTheta[
p3zero2]/h[p3zero2], Im[p3zero1] == 0 && Im[p3zero2] == 0}, {0, Im[p3zero1] != 0 || Im[p3zero2] != 0}}];

Plot3D[NIntegrate[AngularIntegrand[p3zero1, p3zero2]*Sin[\[Theta]], {\[Theta],0, \[Pi]}, MinRecursion -> 4, Method -> {Automatic, "SymbolicProcessing" -> 0}], {\[Nu], 0, 3}, {Q2, 0, 6}, ClippingStyle -> None, ImageSize -> 360] //AbsoluteTiming //Column


which also results in the warning given by BobHanlon. However compared to his method (~166s on my machine), my modifications let it finish after ~125s. So an improvement by ~25%.