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}];

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}];

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\)]\)"}]

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.

enter image description here

(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.

hq1 squared

  • $\begingroup$ As far as I know the Integral/ 'ha1'is correct. May be integration domain need to be changed. Can you suggest any way to remove the singularity? $\endgroup$
    – jazz1001
    Feb 1 '15 at 14:47
  • $\begingroup$ 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? $\endgroup$
    – bbgodfrey
    Feb 1 '15 at 16:20
  • $\begingroup$ yeah, we can choose the normal branch cut. $\endgroup$
    – jazz1001
    Feb 1 '15 at 20:14
  • $\begingroup$ I used Cauchy Principal value-> true to take care of the singularity. Can you please comment on removing the singularities with any other faster option. $\endgroup$
    – jazz1001
    Feb 18 '15 at 22:50
  • $\begingroup$ I do not believe that PrincipalValue -> True does take care of the singularities in this case. $\endgroup$
    – bbgodfrey
    Feb 19 '15 at 2:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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