Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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

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

share|improve this question
up vote 15 down vote accepted

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[η]]&

the ordinal differential equation for our problem

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

Module[{minη, maxη},
    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}

using Manipulate to choose initial values for our BVP

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) # & /@ %

ODE transformed to cylindrical coordinates ODE after getting rid of the unnecessary denominator

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

  sol2=Quiet@NDSolve[Join[{odez/. s-> 1}, {h[4/5]== h0, h'[4/5]== hd0}], h, {z, -1, 1}];
    {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}

comparing the solutions

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}]

our solution in the new domain with a linear approximation

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


  {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:

numerical solution plus simple analytical approximation

share|improve this answer
Practical approach. Like it +1. Does work if the function doesn't have much in this case. – Stefan Jul 8 '13 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$. – Andrei Jul 8 '13 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$. – Andrei Jul 9 '13 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$. – Thies Heidecke Jul 10 '13 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. – Thies Heidecke Jul 11 '13 at 21:41

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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