Message Boards Message Boards

37 Replies
10 Total Likes
View groups...
Share this post:

Imperfect Bose gas nonlinear Schroedinger type ODE numerical solution

I have to find and plot a numerical solution for this second order differential equation:

u''[x] + (u'[x]/x) - (u[x]/(x^2)) + u[x] - u[x]^3 = 0

where 0 <= x <= +inf and u[x=0]=0 ; u[x=+inf]=1.

How can I do? Thank you very much.

POSTED BY: Mattia Udina
37 Replies
Posted 9 years ago

I attacked this problem by making the substitution $x = e^{-t}$. When subbing for the independent variable, we need to make careful use of the chain rule to express $u(x), u'(x), u''(x)$ in terms of $u(t), u'(t), u''(t)$.

ode = u''[x]+(u'[x]/x)-(u[x]/(x^2))+u[x]-u[x]^3 == 0;
sub = Exp[-t];

newode = ode /. {
    u[x] -> u[t],
    u'[x] -> u'[t]/D[sub, t],
    u''[x] -> (u''[t] - D[sub, {t, 2}] u'[t]/D[sub, t])/D[sub, t]^2,
    x -> sub

u[t] - E^(2 t) u[t] - u[t]^3 - E^(2 t) u'[t] + E^(2 t) (u'[t] + u''[t]) == 0

sol = NDSolveValue[{newode, u[1] == 0, u[0] == 1}, u[t], t] /. t -> Exp[-x];
sol /. x -> 0


sol /. x -> 10^15


Plot[Evaluate[sol], {x, 0, 10}, PlotRange -> All]

enter image description here

Here's my issue though... If I plug this solution into the original ODE, I don't get zero, so I think I did something wrong but I can't figure out where. In fact I am having the same issue with Hans Dolhaine's solution he posted above (his is sol2 in the plot below).

Maybe someone here can figure out where I went wrong? I feel it's a simple oversight.

   D[sol, {x, 2}] + (D[sol, x]/x) - (sol/x^2) + sol - sol^3,
   D[sol2, {x, 2}] + (D[sol2, x]/x) - (sol2/x^2) + sol2 - sol2^3
 }], {x, 4, 20}, PlotRange -> All]

enter image description here

POSTED BY: Greg Hurst

I like this approach: It is sufficient to try various initial conditions for u'[eps] to get the correct boundary condition u[Infinity]=1. For instance, below there is the code to get a reasonable plot with s=0.2. By the way, it is the plot of the quantum vortex (in a constant background) of the modified Gross-Pitaevskii equation.

s = 0.2 ; 
eps = 10^-4;

eq = u''[x] + (u'[x]/x) - (u[x]/(x^2)) + u[x] - u[x]^3 - s*D[x*D[u[x]^2, x], x]*(u[x]/x) == 0 ; 

sol = u /. NDSolve[{eq, u[eps] == 0, u'[eps] == 1.11}, u, {x, eps, 5}][[1]]

Plot[sol[x], {x, eps, 5}, PlotTheme -> "Detailed"]

enter image description here

POSTED BY: Luca Salasnich

It is very easy to show (and everyone knows) that Nonlinear Schrodinger in the simplest case has Tanh kink solution:

f[x_] := Tanh[x/Sqrt[2]];
u''[x] + u[x] - u[x]^3 /. u -> f // FullSimplify

Out[] = 0

That is your equation at asymptotic 'x -> + Infinity' because

Limit[u'[x]/x - u[x]/x^2,x->0] == 0

And it is clear that solution around zero can be linearized. Now that's being said Micahel Trott had offered a proper way of finding this solution that you can see below in his reply to tis comment.

POSTED BY: Vitaliy Kaurov

High-order series solution of the ODE around x=0 followed by a Pade approximation. difference calculates the difference of the approximation to the right boundary value. It also returns the actual Pade approximation.

difference[\[Alpha]In_Real,\[Xi]_,{prec_, pdO_}] := 
Module[{\[Alpha]=SetPrecision[\[Alpha]In, prec],diff,sol,series,pd,u,x,\[Delta]},  
series=\[Alpha] x -\[Alpha]/8  x^3 ;
Do[diff=u''[x]+u'[x]/x-u[x]/x^2+u[x]-u[x]^3 /.
           u->(Function@@{x,series+C x^M+O[x]^(M+1)});
 series=series  + Together[C x^M /. First[Solve[diff[[3,1]]==0,C]]],{M,5,2pdO+5,2}];
stepData={\[Alpha],\[Delta]=Tanh[10/Sqrt[2]] - pd/. x->10};
{\[Delta], pd /. x -> \[Xi]}


difference[0.58318 , x, {20, 20}]

{0.2435578857088556, (0.58318000000000003 x + 0.27407625855333046 x^3 + 0.051256728918896911 x^5 + 0.0048322482238297771 x^7 + 0.00023632082879286407 x^9 + 5.1729410403285544*10^-6 x^11 + 8.2185843137751949*10^-9 x^13 - 1.06931812684189022*10^-9 x^15 - 2.2863733027673254*10^-12 x^17 + 7.6295697669084428*10^-14 x^19)/(1 + 0.59496854925294154 x^2 + 0.142883728435766338 x^4 + 0.0176772849277092897 x^6 + 0.00118228825938525016 x^8 + 0.000039832791262778028 x^10 + 4.5345802881727914*10^-7 x^12 - 5.5226863238549098*10^-9 x^14 - 9.7254616067217217*10^-11 x^16 + 6.0544145050658359*10^-13 x^18 + 3.9887426264821583*10^-15 x^20)}

Find slope at x=0 to any desired precision (by increasing prec and pdO). Say 15 digits.

 fr = FindRoot[
    difference[\[Alpha], x, {100, 200}][[1]], {\[Alpha], 0.5831, 0.5832}, WorkingPrecision -> 60, PrecisionGoal -> 16, 
    Method -> "Brent", Evaluated -> False];,


{\[Alpha] -> 0.583189634418924769412413009673771234659851269345508491782915}

The rational approximate solution.

uSol = difference[\[Alpha] /. fr, x, {100, 200}][[2]];
Plot[Evaluate[uSol], {x, 0, 10}, PlotRange -> All]

enter image description here

Check right boundary condition

(uSol /. x -> 10) - Tanh[10/Sqrt[2]]


Check ODE. It is fulfilled to 20 digits everywhere.

odeDiff = u''[x] + u'[x]/x - u[x]/x^2 + u[x] - u[x]^3 /. u -> Function[x, Evaluate[uSol]];

Quiet[ Plot[Log10[Abs[odeDiff]], {x, 0, 10}, WorkingPrecision -> 200, 
  GridLines -> {None, {-20}}, PlotRange -> {Full, {All, 0}}], Plot::precw]

enter image description here

Notebook is attached.

POSTED BY: Michael Trott

Redoing my math actually in Mathematica, I found I had one term wrong. The final approximation should be:

x/(1/2 (1 + Sqrt[5]) + x)

Attached my notebook

POSTED BY: Kay Herbert

Can you tell more about the equation? if you transform with t=Log[x] you get

u''[t]+u'[t]-u+E^(2 t) (u-u^3) = 0 with u[-inf]=0 and u[inf]=1

clearly as t goes to infinity the two right hand terms dominate u->1 as t goes to -infinity the left hand terms dominate
u->A E^((-(1/2) + Sqrt[5]/2) t) where A is a constant.
However, I don't think there is an easy numerical solution because of the infinite domain

POSTED BY: Kay Herbert

If you assume that the solution is u = x^(Sqrt[5]/2 - 1/2)/(B + x^(Sqrt[5]/2 - 1/2))

and the equation is exact at x=1 then B -> (-3 + Sqrt[5] - Sqrt[38 - 14 Sqrt[5]])/(2 (-3 + Sqrt[5])) or about 2.19353 or approximation u=x^0.618/(2.1935+x^0.618)

enter image description here

which satisfies your boundary conditions

POSTED BY: Kay Herbert

Thanks Kay. The equation gives the vortex solution to the GPE, and it represents a wave function for dilute BECs. The domain has not to be necessarily infinite, but big enough. I also have to generalize the method to another equation with is similar, but has a parameter. How have you find the solution u = x^(Sqrt[5]/2 - 1/2)/(B + x^(Sqrt[5]/2 - 1/2)) and the plot (which is exactly what I expected)? What code should I use? (I'm not an expert with Mathematica). Thanks a lot.

POSTED BY: Mattia Udina

Hello Mattia,

as it is a diff.eq. of second order we need to have information about the first derivative as well. And the division by x leads to problems at x == 0. But we could try to start in the vicinity of 0:

eq = u''[x] + (u'[x]/x) - (u[x]/(x^2)) + u[x] - u[x]^3 == 0;
eps = 10^-4;
sol = u /. 
  NDSolve[{eq, u[eps] == 0, u'[eps] == .2}, u, {x, eps, 5}][[1]]
Plot[sol[x], {x, eps, 3}]

enter image description here

POSTED BY: Hans Dolhaine

..and here the solution for different values of the 1st derivative ( Range[ 10 ] / 10. )

enter image description here

POSTED BY: Hans Dolhaine

Dear Hans, unfortunately the behavior of the first derivative is unknown. This is a typical Schroedinger problem, in which only BC's are known. Your solution seems to be correct near 0, but the function is monotone and most reach asymptotically 1 at +inf.

POSTED BY: Mattia Udina

Thanks a lot Kay!

POSTED BY: Mattia Udina

Dear Mattia,

a really interesting problem, and I see the point, but what are BC's? And I am somewhat astonished that a quantum-mechanical function does not vanish for x -> infinity.

But besides that I tried a "brute force" solution:

Your original equation was

eq1 = u''[x] + (u'[x]/x) - (u[x]/(x^2)) + u[x] - u[x]^3;

Now I looked for a function (hopefully) converging towards 1 for x -> Infinity. This should workk if v does not grow faster than exp.

uu[x_] := 1 - v[x] Exp[-x]

(Later I tried Exp[ - x^2 ] as well )

This gave a modified differential equation for v[ x ]

 In[3]:= eq2 = x^2 FullSimplify[eq1 /. u -> uu]

  Out[3]= E^(-3 x) (E^(2 x) (1 + x + x^2) v[x] - 3 E^x x^2 v[x]^2 + 
     x^2 v[x]^3 - E^(2 x) (E^x + (1 - 2 x) x 
  \!\(\*SuperscriptBox["v", "\[Prime]",
  MultilineFunction->None]\)[x] + x^2 
  \!\(\*SuperscriptBox["v", "\[Prime]\[Prime]",

Which I tried to solve numerically

eps = 10^-3;
sol = v /. 
  NDSolve[{eq2 == 0, v[eps] == 1, v'[eps] == 0}, v, {x, eps, 5}, 
    Method -> "ImplicitRungeKutta"][[1]]

But even after about 50 minutes of evaluation no result was given back.

Regards Hans

POSTED BY: Hans Dolhaine

The other equation is the same with another term and a parameter: $$u''[x] + (u'[x]/x) - (u[x]/(x^2)) + u[x] - u[x]^3 - b*(1/x)*D[x*D[u[x]^2]*u[x]=0$$. I hope to be able to get a numerical solution varying b using your idea.

POSTED BY: Mattia Udina

Well, not monotone for small values of x, but......

eq1 = u''[x] + (u'[x]/x) - (u[x]/(x^2)) + u[x] - u[x]^3;
uu[x_] := 1 - v[x] Exp[-x^1.3];
eq2 = x^2 FullSimplify[eq1 /. u -> uu];
eq2a = FullSimplify[Solve[eq2 == 0, v''[x]][[1, 1]]] /. Rule -> Equal
eps = 10^-3;
sol = v /. 
  NDSolve[{eq2a, v[eps] == 1, v'[eps] == -.1}, v, {x, eps, 50}, 
    Method -> "ExplicitRungeKutta"][[1]]
Plot[1 - Exp[-x^2] sol[x], {x, eps, 10}, PlotRange -> All]

enter image description here

POSTED BY: Hans Dolhaine

Dear Hans, it's a very peculiar situation. But it is possible to show that for a large system with a fix number of particles, the state of lowest energy is a uniform distribution. In this calculation, the effect of the boundaries may be neglected. You can consider the system as "unbounded". This is a n-body wave function, for n identical particles (an imperfect Bose gas). Ginzburg and Pitaevskii solved numerically that equation in 1958, I don't know how, and showed that the solution is a monotone curve...

POSTED BY: Mattia Udina

Hans, Interesting approach. I tried to map the infinite domain to a finite one i.e. something like x1=x/(1+x), but the resulting ode doesn't seem to solve with NDSolve.

POSTED BY: Kay Herbert

keep in mind, my approximation is very rough and not accurate. I'm sure there is a way to get a better solution

POSTED BY: Kay Herbert

Unfortunately there is a mistake in the last plot-command, working with Exp[ - x^2 ] instead of Exp[ - x^1.3 ] which was used to solve the differential equation. With the last setting the approach doesn't work.

But another question: the differentials in the original equation are obviously the laplacian in cylindical coordinates. The function part is

u[x] - u[x]/x^2 - u[x]^3

and this may be written as

-u[x] (u[x] - Sqrt[-1 + x^2]/x) (u[x] + Sqrt[-1 + x^2]/x)

What does that mean for the allowed values of x? If u is real then x > 1 ?

POSTED BY: Hans Dolhaine

The new equation is:

ode1[u_] := 
  D[u, {x, 2}] + D[u, x]/x - u/x^2 + u - u^3 - 
    D[x*D[u^2, x], x]*u/x  == 0;

After the transformation and considering only dominating terms, I tried to use:

DSolve[u1''[t] - u1[t] - u1[t]*(u1'[t])^2 - 2*u1''[t]*(u1[t])^2 == 0, 
 u1[t], t]

but I had no output...

POSTED BY: Mattia Udina

DSolve is unlikely to solve equations with nonlinear terms in it

POSTED BY: Kay Herbert

Folks I have not check this to the end but I have a question to all of you. One of the tight spots here is the BC at Infinity. I've seen someone mention mapping (0, Infinity) -> (0, 1). Have you looked into it? There is a beautiful algebraic mapping for the positive half-plane to (0,1):

t[x_] := 1 /(1 + x^(-1/2));
Plot[t[x], {x, 0, 10}, PlotRange -> All]

enter image description here

Series[t[x], {x, 0, 1}]
Series[t[x], {x, Infinity, 1}]

enter image description here

Has anyone tried this change of variable to be at ease using BC now at (0,1) using NDSolve ?

POSTED BY: Vitaliy Kaurov

Dear Vitaliy, this equation gives the vortex solution of the GPE for ultra cold imperfect Bose gases. I have no idea about the behavior of the first derivate but I know that the wave function u[x] is u[x=0]=0 and u[x=+inf]=1, and it's monotone. I don't think that domain has necessary to be infinitive but big enough. The problem is that I don't know how to use NDSolve knowing only BC's. The second problem is that I have to solve a second equation for the modified GPE, which is: $$u?[x]+(u?[x]/x)?(u[x]/(x2))+u[x]?u[x]3?b?(1/x)?D[x?D[u[x]2]?u[x]=0$$ depending on the parameter b, and to compare the two plots.

POSTED BY: Mattia Udina

You do not have to repeat the problem to me, I've read it. The idea that I am suggesting is that the statement

"I don't think that domain has necessary to be infinitive but big enough"

leads to an approximation of BC at Infinity - nothing wrong with that. I am just saying potentially there is a mapping of domain that would allow numerical routine to run on the full interval.

NDSolve does not take Infinity as a spec for BC. So one has to trick out, use asymptotic, approximations. Nothing wrong with that - I am just suggesting a potential way to avoid that and rely completely on numerical routine via easy BC obtained from domain mapping.

This mapping could be wrong too because of potential issues that could appear during solution. Maybe other mapping can be tried then. Just an idea.

POSTED BY: Vitaliy Kaurov

But can I use NDSolve to solve the equations with u[x=0] = 0 and u[x=1000] = 1 (for example)? How? I'm not an expert in Mathematica. Thanks a lot.

POSTED BY: Mattia Udina

Chip, That's the same transformation that I used (i used Exp[t] instead of Exp[-t]) see above. Problem is that now your domain in t goes from -infinity to infinity (for t) instead from 0 to infinity (for x). That's why your answer isn't right as you solve for t in the interval from 0 to 1.

POSTED BY: Kay Herbert
Posted 9 years ago

Ah, ok! When deriving the new domain, I had a typo and used $t = e^{-x}$ instead of $x = e^{-t}$...

POSTED BY: Greg Hurst

So Chip what I have to chance in your code? Thanks a lot

POSTED BY: Mattia Udina

Vitaly, Your transform is similar to the one I tried t=x/(1+x) yours: t=Sqrt[x]/(1+Sqrt[x]) . Any reason Sqrt[x}? In any case, NDSolve was unable to solve the transformed equation. I think the reason is that the equation becomes very stiff near the boundaries. Although I'm sure there is a way, I didn't really spend a lot of time on it.
Mattia mentioned that the equation was solved before, maybe it is worth looking at some literature to see how the problem has been handled in the past.

POSTED BY: Kay Herbert

They used the so called "Chebyshev method", but I don't know how it works.

POSTED BY: Mattia Udina

Actually, Kay, I like your transform much better. Could you post your code or attach a notebook or the transformed ODE and NDSolve errors? We could try to figure out a custom method in NDSolve that would handle stiffness.

POSTED BY: Vitaliy Kaurov

You can attach PDF papers to your posts and also hyperlink to the papers online. Why don't you edit your top post and point people to all essential information?

POSTED BY: Vitaliy Kaurov

Vitaly, Ok here it is:

In[201]:= ode1[u_] := D[u, {x, 2}] + D[u, x]/x - u/x^2 + u - u^3;

In[202]:= pdc[ode1[u[x]]]

Out[202]//TraditionalForm= -(u(x)/x^2)+\[PartialD]^2u(x)/\[PartialD]x^2-u(x)^3+u(x)+\[PartialD]u(x)/\[PartialD]x/x

transform equation

In[203]:= teq1 = Expand[ode1[u[x/(1 + x)]]] /. x/(1 + x) -> t /. x -> t/(1 - t)

Out[203]= u[t] - ((1 - t)^2 u[t])/t^2 - u[t]^3 + (
 2 t Derivative[1][u][t])/((1 - t) (1 + t/(1 - t))^3) - (
 3 Derivative[1][u][t])/(1 + t/(1 - t))^2 + ((1 - t) Derivative[1][u][t])/(
 t (1 + t/(1 - t))) + (
 t^2 (u^\[Prime]\[Prime])[t])/((1 - t)^2 (1 + t/(1 - t))^4) - (
 2 t (u^\[Prime]\[Prime])[t])/((1 - t) (1 + t/(1 - t))^3) + (
  u^\[Prime]\[Prime])[t]/(1 + t/(1 - t))^2

In[204]:= teq2 = Simplify[%] t^2

Out[204]= (-1 + 2 t) u[t] - 
 t^2 u[t]^3 + (-1 + 
    t)^3 t ((-1 + 2 t) Derivative[1][u][t] + (-1 + t) t (u^\[Prime]\[Prime])[

In[205]:= Factor[teq2]

Out[205]= -u[t] + 2 t u[t] - t^2 u[t]^3 + t Derivative[1][u][t] - 
 5 t^2 Derivative[1][u][t] + 9 t^3 Derivative[1][u][t] - 
 7 t^4 Derivative[1][u][t] + 2 t^5 Derivative[1][u][t] + 
 t^2 (u^\[Prime]\[Prime])[t] - 4 t^3 (u^\[Prime]\[Prime])[t] + 
 6 t^4 (u^\[Prime]\[Prime])[t] - 4 t^5 (u^\[Prime]\[Prime])[t] + 
 t^6 (u^\[Prime]\[Prime])[t]

In[206]:= eps = .1;

In[207]:= NDSolve[{teq2 == 0, u[0] == 0, u[1.] == 1.}, u[t], {t, eps, 1.}]

During evaluation of In[207]:= Power::infy: Infinite expression 1/0.^2 encountered. >>

During evaluation of In[207]:= Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

During evaluation of In[207]:= Power::infy: Infinite expression 1/0.^2 encountered. >>

During evaluation of In[207]:= Power::infy: Infinite expression 1/0.^2 encountered. >>

During evaluation of In[207]:= General::stop: Further output of Power::infy will be suppressed during this calculation. >>

During evaluation of In[207]:= Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

During evaluation of In[207]:= Infinity::indet: Indeterminate expression 0. ComplexInfinity encountered. >>

During evaluation of In[207]:= General::stop: Further output of Infinity::indet will be suppressed during this calculation. >>

During evaluation of In[207]:= NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.`. >>

Out[207]= NDSolve[{(-1 + 2 t) u[t] - 
    t^2 u[t]^3 + (-1 + 
       t)^3 t ((-1 + 2 t) Derivative[1][u][t] + (-1 + t) t (
         u^\[Prime]\[Prime])[t]) == 0, u[0] == 0, u[1.] == 1.}, 
 u[t], {t, 0.1, 1.}]

This blows up because of division by zero at the low end. I can fix that by using the low end asymptotic term as t->0 t u'==u:

In[217]:= NDSolve[{teq2 == 0, t u'[eps] == u[eps], u[1.] == 1.}, 
 u[t], {t, eps, 1.}]

During evaluation of In[217]:= NDSolve::ndsz: At t == 0.9999112864801112`, step size is effectively zero; singularity or stiff system suspected. >>

Out[217]= NDSolve[{(-1 + 2 t) u[t] - 
    t^2 u[t]^3 + (-1 + 
       t)^3 t ((-1 + 2 t) Derivative[1][u][t] + (-1 + t) t (
         u^\[Prime]\[Prime])[t]) == 0, 
  t Derivative[1][u][0.1] == u[0.1], u[1.] == 1.}, u[t], {t, 0.1, 1.}]

Now it blows up near t=1 because the differential equation turns into an algebraic equation (u=u^3) via an apparent stiffness at 1 (coefficients of the derivatives start to vanish)

POSTED BY: Kay Herbert

Luca, Nice, you avoid the stiffness issue by using a shooting method! How did you find the initial value of 1.11?

POSTED BY: Kay Herbert

I tried several values for the derivative. Obviously one can improve my choice u'[eps] = 1.11 but I do not think that the improved plot will be very different.

Clearly there is ONLY ONE exact value of u'[eps] to obtain that u[x] to 1 for x \to Infinity, but it is not necessary to find the EXACT value to produce a reasonable plot (for a BSc thesis).

POSTED BY: Luca Salasnich

Yes, that sounds good. But what if you integtate to, say, x = 10? At least with my system the solution "explodes".

And it shows a peculiar behaviour around 1.11. Just have a try with 1.105 and 1.108 and so on....

POSTED BY: Hans Dolhaine

Mattia is working on it. Notice that with s > 1/2 the equartion becomes nasty (with s = 0 one has the familiar Gross-Pitaevskii equation for a quantized vortex).

POSTED BY: Luca Salasnich
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract