# 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? Commented Feb 1, 2015 at 13:34
• @bbgodfrey I edited the question with the 1-d case's code Commented Feb 1, 2015 at 14:11
• Commented Feb 1, 2015 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.)