First I want to solve for $n_{0}$ for a fixed value of $n$, lets say $n=1$ $$ n= n_{0}+ \dfrac{1}{2}\int_{-1/2}^{1/2}dq\left(\dfrac{e_{q}+Un_{0}}{\hbar\omega}-1\right) $$ where $e_{q}=2[1-cos(2\pi q)] $ and $\hbar\omega=(e_{q}^{2}+2Un_{0}e_{q})^{1/2} $

Then I want to plot between $U$ and $n_{0}$. The following is the code written to get it

Clear[eq, hq, ha, vv];
eq[q1_] := (2*(1 - Cos[2*Pi*q1]));
hq[q1_, U_, n0_] := ((eq[q1])^2 + 2*U*n0*(eq[q1]))^(1/2);
ha[q1_, n0_, U_] := (((eq[q1]) + (U*n0))/hq[q1, U, n0]) - 1;
vv[n0_, U_] := n0 + 0.5*Integrate[ha[q1, n0, U], {q1, -0.5, 0.5}];

n[U_] := FindRoot[vv[n0, U] == 1, {n0, 0.1}];

Plot[n[U], {U, 1, 20}, PlotRange -> {{0, 20}, {0, 1}}, 
AxesOrigin -> {0, 0}, Exclusions -> hq == 0, AxesLabel -> {"U", "\!\(\*SubscriptBox[\(n\), \(0\)]\)"}]

the plot should follow similar fashion like that of $y=1/x$ but I am not getting that rather what I get is straight line and a huge bundle of errors.

Is the code written correctly and what should I change to get correct results? Any help is appreciated.

  • $\begingroup$ What kind of errors? The first error I get is "Integrate::idiv: Integral...does not converge on {-0.5,0.5}. >>", just as I did with mathematica.stackexchange.com/questions/73734/… -- you really ought to address the errors produce by your code in your question. $\endgroup$
    – Michael E2
    Feb 17 '15 at 12:29
  • $\begingroup$ Same problem here: mathematica.stackexchange.com/questions/72946/… $\endgroup$
    – Michael E2
    Feb 17 '15 at 12:31
  • 1
    $\begingroup$ your integral is singular. $\endgroup$
    – george2079
    Feb 17 '15 at 14:11
  • $\begingroup$ @MichaelE2: Please remove the hold, since using the Cauchy principal value method, the problem can actually be solved (by circumventing the obvious singularity), as I showed in my answer. $\endgroup$
    – Jinxed
    Feb 19 '15 at 8:23
  • $\begingroup$ @Jinxed, Integrate is not the most reliable witness. As I explained in a comment to a related question, the integrand at q == 0 is asymptotically proportional to 1/Abs[q]; so it does not have a Cauchy principal value. If you disagree, you can edit the question to add your explanation of why the integral converges. Editing the question causes the question to go into the review queue. Then the community will vote on whether to reopen the question. $\endgroup$
    – Michael E2
    Feb 19 '15 at 12:59

If you are prepared to rely on Cauchy principal values, you can "circumvent" the singularity of your integral around the origin like so:

eq = 2 (1 - Cos[2 \[Pi] q]);
hw = Sqrt[eq^2 + 2 U n0 eq];
eqn = (eq + U n0)/hw - 1;
res = n0 + 1/2 Integrate[eqn, {q, -1/2, 1/2}, PrincipalValue -> True];

After some minutes, you will be presented with the following result in res:

 n0 + 1/4 (-2 + (
     8 Sqrt[1/(n0 U)] ArcTan[Sqrt[2]/Sqrt[n0 U]] + 
      Sqrt[2] (-Log[1 + 2/(n0 U)] + Log[4/(n0 U)] + Log[n0 U]))/(
     2 \[Pi] Sqrt[1/(n0 U)])), 
 Im[1/(n0 U)] != 0 || Re[1/(n0 U)] >= -(1/2)]

and can find an instance using e.g. N@FindInstance[res == 1 && U \[Element] Integers, {n0, U}], {{n0 -> 1.19646, U -> 2.}}. This first idea of suitable values is important, since Reduce fails to produce a result, even after hours of computing.

Starting from that instance, you can do an exemplary ContourPlot for res==1:

contour = 
    n0 + 1/4 (-2+(8 Sqrt[1/(n0 U)] ArcTan[Sqrt[2]/Sqrt[n0 U]] + 
      Sqrt[2] (-Log[1 + 2/(n0 U)] + Log[4/(n0 U)] + 
         Log[n0 U]))/(2 \[Pi] Sqrt[1/(n0 U)]))==1.,
    {U,0.1,8.}, {n0,1.1,1.2}]


If you need the contour's coordinates, use:

pts=Cases[Normal@contour, Line[x_] :> x, Infinity][[1]]//Reverse;

and then possibly do a Fit:

fit = Fit[pts, {1, x^-2, x, x^2, Log[x]}, x]
(* 1.20802 + 0.000132045/x^2 - 0.0245595 x + 0.000654567 x^2 + 
   0.0502606 Log[x] *)

which matches the contour quite well, as you may find.

  • $\begingroup$ The code is running past 2 hours, still no output. Do you know how to speed it up? $\endgroup$
    – jazz1001
    Feb 17 '15 at 19:57
  • $\begingroup$ More memory and a faster machine. :) It took me some time also, but not two hours or even more. $\endgroup$
    – Jinxed
    Feb 18 '15 at 6:39
  • $\begingroup$ This is not what I want. I want to plot between U and n0, so using FindInstance doesn't seems a good option. What I want is U in terms of n0. $\endgroup$
    – jazz1001
    Feb 18 '15 at 22:29
  • $\begingroup$ Thanks for your kind response. This makes me feel so inclined to try to help you further. :| $\endgroup$
    – Jinxed
    Feb 18 '15 at 22:36
  • $\begingroup$ hey I didn't meant that way. Well my machine just completed running the code after such a long time :P, that's why I noticed it. Anyways, any help/comment is appreciated. $\endgroup$
    – jazz1001
    Feb 18 '15 at 22:40

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