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.
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
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.
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}]
dgl.bmp (350 KB)
..and here the solution for different values of the 1st derivative ( Range[ 10 ] / 10. )
5868dgl.bmp (339 KB)
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.
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
Thanks a lot Kay!
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
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.
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]
7406dgl.bmp (337 KB)
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…
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.
keep in mind, my approximation is very rough and not accurate. I’m sure there is a way to get a better solution
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 ?
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…
DSolve is unlikely to solve equations with nonlinear terms in it
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 ?
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.