# Imperfect Bose gas nonlinear Schroedinger type ODE numerical solution

Posted 8 years ago
12899 Views
|
37 Replies
|
10 Total Likes
|
 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 = 0where 0 <= x <= +inf and u[x=0]=0 ; u[x=+inf]=1.How can I do? Thank you very much.
37 Replies
Sort By:
Posted 8 years ago
 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[] = 0That 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 8 years ago
 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}]; pd=PadeApproximant[series,{x,0,{pdO,pdO}}]; stepData={\[Alpha],\[Delta]=Tanh[10/Sqrt[2]] - pd/. x->10}; {\[Delta], pd /. x -> \[Xi]} ] Example 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. Monitor[ fr = FindRoot[ difference[\[Alpha], x, {100, 200}][[1]], {\[Alpha], 0.5831, 0.5832}, WorkingPrecision -> 60, PrecisionGoal -> 16, Method -> "Brent", Evaluated -> False];, stepData] fr {\[Alpha] -> 0.583189634418924769412413009673771234659851269345508491782915}The rational approximate solution. uSol = difference[\[Alpha] /. fr, x, {100, 200}][[2]]; Plot[Evaluate[uSol], {x, 0, 10}, PlotRange -> All] Check right boundary condition (uSol /. x -> 10) - Tanh[10/Sqrt[2]] -2.6496232607434627305318392508171048229727621206337120044345649125837*10^-27Check 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] Notebook is attached. Attachments:
Posted 8 years ago
 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 8 years ago
 Ah, ok! When deriving the new domain, I had a typo and used $t = e^{-x}$ instead of $x = e^{-t}$...
Posted 8 years ago
 So Chip what I have to chance in your code? Thanks a lot
Posted 8 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  0. sol /. x -> 10^15  1. Plot[Evaluate[sol], {x, 0, 10}, PlotRange -> All] 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. Plot[Evaluate[{ 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] 
Posted 8 years ago
 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]  Series[t[x], {x, 0, 1}] Series[t[x], {x, Infinity, 1}] Has anyone tried this change of variable to be at ease using BC now at (0,1) using NDSolve ?
Posted 8 years ago
 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 8 years ago
 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 8 years ago
 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 8 years ago
 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 8 years ago
 They used the so called "Chebyshev method", but I don't know how it works.
Posted 8 years ago
 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 8 years ago
 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 8 years ago
 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])[ t]) 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) Attachments:
Posted 8 years ago
 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 8 years ago
 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 8 years ago
 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] 
Posted 8 years ago
 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 8 years ago
 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]", MultilineFunction->None]$$[x])) 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 8 years ago
 ..and here the solution for different values of the 1st derivative ( Range[ 10 ] / 10. )
Posted 8 years ago
 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 8 years ago
 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"] 
Posted 8 years ago
 Luca, Nice, you avoid the stiffness issue by using a shooting method! How did you find the initial value of 1.11?
Posted 8 years ago
 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 8 years ago
 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 8 years ago
 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 8 years ago
 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}] 
Posted 8 years ago
 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)which satisfies your boundary conditions
Posted 8 years ago
 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 8 years ago
 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 Attachments:
Posted 8 years ago
 Thanks a lot Kay!
Posted 8 years ago
 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 8 years ago
 keep in mind, my approximation is very rough and not accurate. I'm sure there is a way to get a better solution
Posted 8 years ago
 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 8 years ago
 DSolve is unlikely to solve equations with nonlinear terms in it
Posted 8 years ago
 Can you tell more about the equation? if you transform with t=Log[x] you getu''[t]+u'[t]-u+E^(2 t) (u-u^3) = 0 with u[-inf]=0 and u[inf]=1clearly 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