For a physics application I am considering the radial Bogoliubov-de Gennes equations in two dimensions. In dimensionless form they read as follows $$ \begin{cases} \bigg[\frac{1/4-m^2}{r^2}+2(E-1)\bigg] u_m(r)-2v_m(r)+u_m''(r)=0 \newline \bigg[\frac{1/4-m^2}{r^2}-2(E+1)\bigg] v_m(r)-2u_m(r)+v_m''(r)=0, \end{cases} $$ where $m = 0,1,2$ and $E$ is the energy, a continuous parameter. In fact, one can show that $$ \begin{cases} u_m(r) = -\frac{2\sqrt{r}J_m(kr)}{2+k^2-2E} \newline v_m(r) = \sqrt{r}J_m(kr) \end{cases} $$ is a solution for $k = \sqrt{2(\sqrt{1+E^2}-1)}$ by explicitly substituting it in the equations.

At some point I want to multiply the off-diagonal terms in the system of equations by a function of the radius $r$. This will not be exactly solvable. Hence, I would like to understand how to numerically solve these equations in Mathematica before I make the equations more difficult.

The code I have up to now, is the following:

BdGplot[E1_, Bc1_, m1_] :=
  Module[{E = E1, Bc = Bc1, m = m1}, {btc = 50,

  k = Sqrt[2]*Sqrt[-1 + Sqrt[1 + E^2]],
  uw = -(2 /(2 + k^2 - 2E)),
  ub[r_] := (2/Pi)^(1/2)*uw*Cos[k*r - m*Pi/2 - Pi/4],
  vb[r_] := (2/Pi)^(1/2)*Cos[k*r - m*Pi/2 - Pi/4],

  s = NDSolve[{((1/4 - m^2)/r^2 + 2 (-1 + E)) u[r] - 
    2 v[r] + u''[r] == 0, -2 u[r] + ((1/4 - m^2)/r^2 - 2(1 + E)) v[r] + v''[r] == 0, u[Bc] == ub[Bc], u'[Bc] == ub'[Bc], v[Bc] == vb[Bc], v'[Bc] == vb'[Bc]}, {u,v}, {r, 1/1000, btc}]};

 Plot[Evaluate[{u[r]/(uw*Sqrt[r]), v[r]/Sqrt[r], 
 BesselJ[m, k*r]} /. s], {r, 0.001, btc}, 
 PlotRange -> {{0, btc}, {-2, 2}}, PlotStyle -> {Red, Blue, Orange}]

The idea is to give boundary conditions at large distance $r$ which are the asymptotic values of the m-th Bessel functions (in the more difficult case, the boundary conditions will be similar). Subsequently I let Mathematica numerically integrate the equations to a small radius $r=0.001$, such that it does not run into trouble at $r=0$.

However, plotting with a command like plotAlpha1 = Manipulate[BdGplot[1.5, Bc, 0], {Bc, 1, 50, 1/10}], one sees that the numerical solutions only match to the exact solution in a certain region. At some point the solutions diverge heavily. Does anybody have any idea how I can resolve this issue?

  • $\begingroup$ Could you rework your code a bit such that it can be copy&pasted? $\endgroup$
    – user21
    Apr 23, 2015 at 16:12
  • $\begingroup$ @user21 I have updated the code, it should be copy-pastable now. $\endgroup$
    – Funzies
    Apr 24, 2015 at 8:00
  • $\begingroup$ I think it might be just numeric precision loss. Try something like WorkingPrecision -> 50 which will improve the situation a lot. You will have to take care that all constants will be given as exact numbers, but that would only mean to change 1.5 to 3/2 for the value in the Manipulate... $\endgroup$ Jun 5, 2015 at 9:37
  • $\begingroup$ I think WorkingPrecision->30 should be good enough, and I forgot to mention that this is of course to be used as an option to NDSolve... $\endgroup$ Jun 5, 2015 at 9:52
  • 1
    $\begingroup$ @xzczd : I'm not sure why you did see other behaviour before, but that is the usual rule: precision will never automatically increase, so if you start out with machine precision numbers only in very limited cases increasing WorkingPrecision higher than MachinePrecision will have any effect, and usually not the desired one (AFAIK). To be on the safe side I do use completely "exact" input if I play around with the precision stuff beyond machine precision and that so far has been free from any surprises... $\endgroup$ Jun 5, 2015 at 15:46

1 Answer 1


As indicated in the comments, a larger WorkingPrecision (and, when necessary, also a larger MaxSteps) greatly improves the results of BdGplot as defined in the Question. For instance, with WorkingPrecision -> 50, MaxSteps -> 20000

BdGplot[2, 20, 2]


enter image description here

However, the code can be modified to still further improve the accuracy of the result:

BdGplot1[e_,  m_] := Module[
  {btc = 50, r0 = 1/1000, k = Sqrt[2]*Sqrt[-1 + Sqrt[1 + e^2]]},
  uw = -2 /(2 + k^2 - 2 e);
  s = Quiet@NDSolve[{((1/4 - m^2)/r^2 + 2 (e - 1)) u[r] - 2 v[r] + u''[r] == 0,
                     ((1/4 - m^2)/r^2 - 2 (e + 1)) v[r] - 2 u[r] + v''[r] == 0,
                     r0 u'[r0] == (m + 1/2) u[r0], r0 v'[r0] == (m + 1/2) v[r0], 
                     u[r0] == uw Sqrt[r0] BesselJ[m, k*r0], u[btc] == uw v[btc]}, 
                     {u, v}, {r, r0, btc}, WorkingPrecision -> 40, MaxSteps -> 20000];
  Plot[Evaluate[{u[r]/(uw*Sqrt[r]), v[r]/Sqrt[r],  BesselJ[m, k*r]} /. s], {r, r0, btc}, 
      PlotRange -> {{0, btc}, {-1, 1}}, PlotStyle -> {Red, Blue, Orange}]]

BdGplot1[2, 2]

enter image description here

although the modified code is a bit slower, and NDSolve complains about the boundary conditions before producing the correct answer. (The three curves, which appear as one in the plot, agree to six significant figures.)

BdGplot1 readily generalizes to the modified calculation mentioned in the question. The boundary conditions r0 u'[r0] == (m + 1/2) u[r0], r0 v'[r0] == (m + 1/2) v[r0] should be unchanged, the boundary condition u[r0] == uw Sqrt[r0] BesselJ[m, k*r0] is arbitrary and was chosen merely for normalization, and the boundary condition u[btc] == uw v[btc] changes insofar as uw may need to be modified, if the asymptotic ratio of u and v changes.


The added option,

Method -> {"Shooting", "StartingInitialConditions" -> {r0  u'[r0] == (m + 1/2) u[r0], 
    r0 v'[r0] == (m + 1/2) v[r0], u[r0] == uw v[r0], v[r0] == Sqrt[r0] BesselJ[m, k*r0]}, 
    Method -> {"StiffnessSwitching"}}

increases the speed of BdGplot1[2, 20, 2] by a factor of about three.

  • $\begingroup$ Great, detailed answer, thanks! One question: I tried changing the boundary condition from 50 to 100. In this case, the solutions start to diverge around r = 60, even if I increase the maximum step size and accuracy. Do you have any clue what I should change in this case? $\endgroup$
    – Funzies
    Jun 8, 2015 at 8:33
  • $\begingroup$ Additionally, for higher energies (from about E = 5) there are some problems as well.. Do you know of any solutions to this? Or perhaps I am just asking too much of Mathematica. $\endgroup$
    – Funzies
    Jun 8, 2015 at 8:41
  • $\begingroup$ Your equations support two solutions, one with the k in the Question, and the other with k = Sqrt[2]*Sqrt[-1 + I*Sqrt[1 + e^2]], The second solution has the character of a Modified Bessel Function for large e, which grows exponentially. For larger e larger values of WorkingPrecision are needed to prevent the second solution from polluting the first. So, for e == 8, WorkingPrecision -> 80 is needed. Also, be sure to include the Method option that I gave in the Addendum. At some point, just increasing WorkingPrecision is not practical, and I have some other ideas. $\endgroup$
    – bbgodfrey
    Jun 8, 2015 at 13:55
  • $\begingroup$ @Funzies You state in your question that you want to multiply the off-diagonal terms by a function of the radius r. Is your goal to obtain numerical expressions for u and v, or obtain some averaged quantity, like the phase shift at large r? Also, do your off-diagonal multipliers differ significantly from 1 over a large range of r, or are they highly localized? I have obtained a solution for the multiplier (1 + 5 Exp[-(r - btc/2)^2]), which is fairly localized. The calculation requires a high WorkingPrecision and, therefore, is slow but produces accurate results. $\endgroup$
    – bbgodfrey
    Jun 9, 2015 at 0:11
  • $\begingroup$ I need to calculate the numerical expressions for u and v for different values of m and E. I'm afraid that the off-diagonal function is horrendous: it is $\exp(i \pm \theta(r))$ for the upper (lower) diagonal, with $\theta(r) = c_0 \log(r)$, where $c_0$ is a negative constant. (It has this weird form because in these type of physics problems the derivative of theta is proportional to the superfluid velocity, which in this case is giving by 1/r). So the function is highly non-localized. However, intuitively I still expect the ordinary BdG results to hold at large r. $\endgroup$
    – Funzies
    Jun 9, 2015 at 7:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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