solve ODE with divergencies

Question

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$.

Answer 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["DifferentialEquations`InterpolatingFunctionAnatomy`"];
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: