Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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[\[Eta]_ /; \[Eta] > 0] = 1.;
f1[\[Eta]_ /; \[Eta] == 0] = 2.;
f2[\[Eta]_ /; \[Eta] > 0] = 1./\[Eta];
f2[\[Eta]_ /; \[Eta] == 0] = 0.;
f3[\[Eta]_ /; \[Eta] > 0] = 4./\[Eta]^2 f[\[Eta]];
f3[\[Eta]_ /; \[Eta] == 0] = 4.;

NDSolve[{f1[\[Eta]] \!\(
\*SubscriptBox[\(\[DifferentialD]\), \(\[Eta], \[Eta]\)]\(f[\[Eta]]\)\
\) + f2[\[Eta]] \!\(
\*SubscriptBox[\(\[DifferentialD]\), \(\[Eta]\)]\(f[\[Eta]]\)\) + 
f[\[Eta]] - f3[\[Eta]] - f[\[Eta]]^3 == 0, f[0] == 0, f[10] == 1},
f[\[Eta]], {\[Eta], 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
add comment

1 Answer 1

up vote 10 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/\[Eta] D[\[Eta] D[f[\[Eta]],\[Eta]],\[Eta]]+(1-s^2/\[Eta]^2)f[\[Eta]]-f[\[Eta]]^3==0//Expand//Collect[#,f[\[Eta]]]&

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

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Module[{min\[Eta],max\[Eta]},
  Manipulate[
    sol = Quiet@NDSolve[Join[{ode/.s->s0},{f[1]==f1,f'[1]==fd1}],f,{\[Eta],0,10}];
    {min\[Eta], max\[Eta]} = First@InterpolatingFunctionDomain[f/.sol[[1]]];
    Plot[f[\[Eta]] /. sol, {\[Eta], min\[Eta], max\[Eta]}, 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)]&) /. \[Eta] -> 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:

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[(\[Eta]^2-1)/(\[Eta]^2+1)]/.sol2,f[\[Eta]]/.sol},{\[Eta],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 -> (\[Eta]^2-1)/(\[Eta]^2+1) //Simplify

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

Plot[
  {h[(\[Eta]^2 - 1)/(\[Eta]^2 + 1)] /. sol2,
   \[Eta]^2/(\[Eta]^2 + 1)}, {\[Eta], 0, 20},
  PlotRange -> {0, 1}, AxesLabel -> {\[Eta]}, 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 singularities...like 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
    
This was meant more as a starting point for your own experimenting. When you vary the initial values with the sliders you will see how sensitively the $\eta \approx 10$ part reacts to the initial conditions. Since the ODE for this problem is stiff, numerical shooting methods always tend to become unstable further from the start point. There are other approaches to solve this problem, e.g. taking advantage of the cylindrical nature of the problem by using a Hankel transform or trying a series solutions with BesselJ[s,\[Eta]] functions. –  Thies Heidecke Jul 9 '13 at 13:39
    
great! @ThiesHeidecke, could you please detail on the Hankel/Bessel approaches? –  Andrei Jul 9 '13 at 19:02
show 4 more comments

Your Answer

 
discard

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.