# solve ODE with divergencies

The solutions of the second order differential equation

$$\frac{1}{\eta}\frac{d}{d\eta}\left(\eta \frac{df}{d\eta}\right)+\left(1-\frac{s^2}{\eta^2}\right)f-f^3=0$$

is shown in Fig. 5.2 below, for two different values of $s$.

The result is taken from the book "Bose-Einstein Condensation", by Pitaevskii&Stringari (2003).

As the authors comment, $\lim_{\eta \rightarrow 0} f(\eta) = \eta^{|s|}$, whereas for $\eta \rightarrow \infty$, the limit is 1. This can also be seen from the plot above.

The problems arise due to the divergencies at $\eta=0$. The one in the first term is relatively easy to take care of, but I do not know how to get rid of the divergency in the term $\frac{s^2}{\eta^2}$.

For the $s=2$ case, where $f$ goes to zero as $\eta^2$, this term has the limit 1:

$$\lim_{\eta \rightarrow 0} \frac{4}{\eta^2}f = 4,$$ so I tried to make use of that.

Here is what I came up with so far:

ClearAll[f1, f2, f3];
f1[η_ /; η > 0] = 1.;
f1[η_ /; η == 0] = 2.;
f2[η_ /; η > 0] = 1./η;
f2[η_ /; η == 0] = 0.;
f3[η_ /; η > 0] = (4./η^2) f[η];
f3[η_ /; η == 0] = 4.;
NDSolve[{f1[η] Dt[f[η], η, η] + f2[η] Dt[f[η], η] + f[η] - f3[η] - f[η]^3 == 0,
f[0] == 0, f[10] == 1}, f[η], {η, 0, 10}]


Although this gets rid of the divergencies for the particular case of $s=2$, it results in another series of errors and no result is obtained.

Edit: I am particularly interested in the case $s=1$.

Shooting method with Manipulate

One pragmatic approach of getting a solution for your boundary value problem is just guessing efficiently (which is what most numerical BVP codes do anyway...). A nice way of doing this in Mathematica after setting up our ordinary differential equation

ode=1/η D[η D[f[η], η], η]+(1 - s^2/η^2) f[η] - f[η]^3 == 0//Expand//Collect[#, f[η]]&


is to wrap the NDSolve code in a Manipulate which lets us try different combinations for our initial values neatly with sliders instead of typing them in every time

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];
Module[{minη, maxη},
Manipulate[
sol = Quiet@NDSolve[Join[{ode/.s->s0}, {f[1] ==f1, f'[1] ==fd1}], f, {η, 0, 10}];
{minη, maxη} = First@InterpolatingFunctionDomain[f /. sol[[1]]];
Plot[f[η] /. sol, {η, minη, maxη}, PlotRange->{{0, 10}, {0, 1}}]
,{{f1, 0.520052}, 0.52, 0.5205}, {{fd1,0.410635}, 0.409, 0.412}, {{s0,1}, 0, 3, 1}
]
]


A few things to note: Quiet suppresses the singularity warnings of NDSolve (which we are already aware of) which can be extremely annoying, especially inside a Manipulate. We loaded an extra package to get the domain of our numerical solution, so we can actually get an approximation of the singularities and use them to choose an appropriate PlotRange. This could also be a starting point of doing more fancy automatic optimization of the initial values to maximize the support of our ODE solution.

Changing variables via stereographic projection

To better see the behaviour of $f$ when $\eta$ approaches infinity and also to be able to use other tools that work (only or best) in a finite interval (e.g. fourier series, finite difference relaxation methods) we can transform the ODE to the cylindric domain $z\in (-1,1)$ via change of variables given by a stereographic projection $$z=\frac{\eta^2-1}{\eta^2+1} \Leftrightarrow \eta=\sqrt{\frac{1+z}{1-z}}.$$

We can let Mathematica do all the tedious stuff of working out the correctly transformed partial derivatives and just get rid of the unnecessary denominator:

ode /. f -> (h[(#^2 - 1)/(#^2 + 1)]&) /. η -> Sqrt[(1 + z)/(1 - z)] // Simplify
odez = (1 + z) # & /@ %


Now we can solve the ODE in the new domain and compare it to our original solution in the old domain:

Manipulate[
sol2=Quiet@NDSolve[Join[{odez/. s-> 1}, {h[4/5]== h0, h'[4/5]== hd0}], h, {z, -1, 1}];
GraphicsColumn[
{Plot[{h[z] /. sol2}, {z, -1, 1}, PlotRange -> {0, 1}],
Plot[{h[(η^2 - 1)/(η^2 + 1)]/. sol2, f[η]/. sol}, {η, 0, 20}, PlotRange->{0, 1}]}
]
, {{h0, 0.917477}, 0.90, 0.94}, {{hd0, 0.546211827}, 0.5, 0.6}
]


It's a bit easier to keep an overview of how the solution evolves, although now the equations are a bit more sensitive as we squeezed alot of the $\eta$ to infinity range into the tiny end near $z\approx 1$. But the good news is, now we have all the tools at hand which work in the finite interval (-1,1). I'll just present a quick and cheap one: At first sight the solution in the new domain looks almost linear, so why not use that to get a first approximation to our function?

Plot[{h[z] /. sol2, (1 + z)/2}, {z, -1, 1}, PlotRange -> {0, 1}]


The function $(1+z)/2$ is the easiest approximation one can think of which still satisfies both of our boundary conditions and is not looking too bad actually. Now we can transform this approximate solution back to our original domain:

(1 + z)/2 /. z -> (η^2 - 1)/(η^2 + 1) // Simplify


$$\frac{\eta^2}{\eta^2+1}$$

Plot[
{h[(η^2 - 1)/(η^2 + 1)] /. sol2,
η^2/(η^2 + 1)}, {η, 0, 20},
PlotRange -> {0, 1}, AxesLabel -> {η}, AspectRatio -> 0.4
]


and see that it doesn't compare too badly given its simplicity:

• Practical approach. Like it +1. Does work if the function doesn't have much singularities...like in this case. Jul 8, 2013 at 21:09
• $s$ is not a parameter, it is in fact the quantum number corresponding to the z component of the angular momentum, therefore it can only have values 1,2,3,.. and I am particularly interested in the case $s=1$. Jul 8, 2013 at 21:45
• I also don't understand why your solution has a negative slope at the end of the interval, when it should flatten out towards $f$=1 at $\eta=10$. Jul 9, 2013 at 12:09
• @Andrej i looked into it and it seems the Hankel transform approach is only applicable in case where the function is sufficiently bounded (see domain of definition). Another thing i tried is changing variables via stereographic projection to the interval (-1,1), which makes it easier to see the behaviour of $\eta$ near $\infty$. Jul 10, 2013 at 21:22
• @Andrei Two books that i really liked for their clear and approachable style are 'Numerical Methods for Scientists and Engineers' by R.W. Hamming (his book 'Digital Filters' is also great) and 'Partial Differential Equations for Scientists and Engineers' by Stanley J. Farlow, both published by Dover. I liked that they are not too theoretical, just the amount needed to get going and getting motivated to try out stuff on your own. Jul 11, 2013 at 21:41