Optimising the code for better performance, when using NIntegrate to solve a double/triple integral

Clear[eq11, hq1, ha1, vv32];
eq11[q1_?NumericQ,q2_?NumericQ] := (2*(1 - Cos[2*Pi*q1] + Cos[2*Pi*q2]));
hq1[q1_?NumericQ, q2_?NumericQ, U_?NumericQ,n0_?NumericQ] :=
((eq11[q1, q2])^2 + 2*U*n0(eq11[q1, q2]))^(1/2);
ha1[q1_?NumericQ, q2_?NumericQ, n0_?NumericQ,U_?NumericQ] :=
(((eq11[q1, q2]) + (U*n0))/hq1[q1, q2, U, n0]) - 1;
vv32[n0_?NumericQ, U_?NumericQ] :=
n0 + 0.5*NIntegrate[ha1[q1, q2, n0, U], {q1, -0.5, 0.5}, {q2, -0.5, 0.5}];

Clear[n];
n[U_?NumericQ] := FindRoot[vv32[n0, U] == 0.5, {n0, .1}][[1, 2]];

Quiet@Plot[n[U], {U, 1, 20}, PlotRange -> {{0, 20}, {0, 1}},AxesOrigin -> {0, 0},
AxesLabel -> {"U", "\!$$\*SubscriptBox[\(n$$, $$0$$]\)"}]


When I am solving for 1-d it runs in about 5 secs, but for 2-d (the one above) it takes a lot of time and still running. Can anyone help me out with the problem?

For 1-d case, this the code:

Clear[eq11, hq1, ha1, vv32];
eq11[q1_?NumericQ] := (2*(1 - Cos[2*Pi*q1]));
hq1[q1_?NumericQ, U_?NumericQ, n0_?NumericQ] := ((eq11[q1])^2 + 2*U*n0 (eq11[q1]))^(1/2);
ha1[q1_?NumericQ, n0_?NumericQ, U_?NumericQ] := (((eq11[q1]) + (U*n0))/hq1[q1, U, n0]) - 1;
vv32[n0_?NumericQ, U_?NumericQ] := n0 + 0.5*NIntegrate[ha1[q1, n0, U], {q1, -0.5, 0.5}];

Clear[n];
n[U_?NumericQ] := FindRoot[vv32[n0, U] == 1, {n0, .1}][[1, 2]];

Quiet@Plot[n[U], {U, 1, 20}, PlotRange -> {{0, 20}, {0, 0.5}},
AxesOrigin -> {0, 0},
AxesLabel -> {"U", "\!$$\*SubscriptBox[\(n$$, $$0$$]\)"}]

• What code do you use for the 1-d case? – bbgodfrey Feb 1 '15 at 13:34
• @bbgodfrey I edited the question with the 1-d case's code – jazz1001 Feb 1 '15 at 14:11
• – Michael E2 Feb 1 '15 at 14:56

2_D Integration is slow (and incorrect), because the integrand, ha1 is singular in parts of the integration domain for n0 = 0.1 (and elsewhere), as can be seen in the plot below for U = 5. Perhaps, your equations or domain of integration need to be corrected. In general, using Quiet is unwise, when debugging a code.

(The plot is ragged at the singularity, because too few PlotPoints were used in creating it.)

ha1 is singular where hq1 = 0. The ContourPlot below of hq1^2 indicates that this is true on four semicircles. Note that hq1^2 < 0 in the darkest region. Thus, hq1 and ha1 are complex there.
• Changing the domain to exclude the singularity is, of course, easy to do, but the choice for doing so should be based on the actual mathematical problem you are trying to solve. I notice that hq1 involves Sqrt, which has a branch cut. Is it possible that a different branch than the normal one should be chosen? – bbgodfrey Feb 1 '15 at 16:20
• I do not believe that PrincipalValue -> True does take care of the singularities in this case. – bbgodfrey Feb 19 '15 at 2:02