Diverging solution to coupled second order ODEs from NDSolve

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?

• Could you rework your code a bit such that it can be copy&pasted? – user21 Apr 23 '15 at 16:12
• @user21 I have updated the code, it should be copy-pastable now. – Funzies Apr 24 '15 at 8:00
• 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... – Albert Retey Jun 5 '15 at 9:37
• 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... – Albert Retey Jun 5 '15 at 9:52
• @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... – Albert Retey Jun 5 '15 at 15:46

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

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]


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.

Method -> {"Shooting", "StartingInitialConditions" -> {r0  u'[r0] == (m + 1/2) u[r0],

increases the speed of BdGplot1[2, 20, 2] by a factor of about three.
• 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. – bbgodfrey Jun 8 '15 at 13:55
• @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. – bbgodfrey Jun 9 '15 at 0:11
• 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. – Funzies Jun 9 '15 at 7:31