8
$\begingroup$

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?

Improve this question
$\endgroup$
6
  • $\begingroup$ Could you rework your code a bit such that it can be copy&pasted? $\endgroup$ Commented Apr 23, 2015 at 16:12
  • $\begingroup$ @user21 I have updated the code, it should be copy-pastable now. $\endgroup$ Commented 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$ Commented 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$ Commented 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$ Commented Jun 5, 2015 at 15:46

1 Answer 1

Reset to default
13
+50
$\begingroup$

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]

yields

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.

Addendum

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.

Improve this answer
$\endgroup$
8
  • $\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$ Commented 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$ Commented 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$ Commented 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$ Commented 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$ Commented Jun 9, 2015 at 7:31

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.