# Find solution of nonlinear ODE in terms of JacobiCN

I am trying to find a specific solution for this differential equation:

$-\frac{1}{2}\frac{d^2}{dx^2}\psi(x)-2k \; \psi(x)^3 + \frac{1}{2}k^2\; \psi(x)=0$

MMA gives me a solution in the form of a JacobiSN function that, after setting the second argument equal to one by choosing the right constant of integration, can be rewritten as a Tanh function. However, I know from a ref that also

$\psi(x)=\sqrt{\frac{k}{2}}\; \text{sech}({kx})$

solves the ode. How to make MMA show this solution? I know that this could be obtained by choosing a JacobiCN instead of a JacobiSN, but how to implement this?

The point is that I want to find the solutions of the 2D version of this equation, and MMA only finds the Tanh solution just like in 1D. How to know if there's a Sech solution (or similar) in 2D?

So I guess the general question is: is there a way to "guide" MMA when it chooses the solutions for ODEs/PDEs, for example in the above case by making it choose a JacobiCN?

EDIT

The Tanh solution I find is the following:

$\frac{\sqrt{k}}{2} \text{tanh}\left(i\frac{k|x|}{\sqrt{2}}\right)$

Here is the code:

The equation:

sol1D=DSolve[-2 k ψ[x]^3 - D[ψ[x], {x, 2}]/2 == -(1/2) k^2 ψ[x], {ψ[x]}, {x}]


The output of the above is a JacobiSN function, it's pretty big and not insightful so it makes no sense to paste it here. Now, JacobiSN[u,m] simplifies to Tanh for m=1. Important to note: JacobiSN never simplifies to Sech, but JacobiCN does. Here is what the solution becomes when you adjust the constant to get the second argument to 1 (and I choose the other constant to be zero for simplicity):

In:

Simplify[sol1D[[1]] /. C[1] -> -(k^3/8) /. C[2] -> 0, x ∈ Reals && k > 0]


Out:

 ψ[x] -> 1/2 I Sqrt[k] Tan[(k Abs[x])/Sqrt[2]]}


Here below I check the Sech solution:

In:

sechSol[x_] := (Sqrt[k] Sech[k x])/Sqrt[2]
Simplify[-1/2 (D[sechSol[x], {x, 2}]) - 2 k (sechSol[x])^3] == -1/2 k^2  sechSol[x]


Out: True

As expected. Now the question is how to make this sechSol show up.

2nd EDIT

For anyone who wants to know: I found analytically that also the 2D case has a Sech solution. In two dimensions the equation reads:

$-\frac{1}{2}\frac{d^2}{dx^2}\psi(x,y)-\frac{1}{2}\frac{d^2}{dy^2}\psi(x,y)- \psi(x,y)^3 + \frac{1}{2}k^2\; \psi(x,y)=0$

The coefficient in front of the nonlinear term changes (dimensionality arguments). This equation is solved by:

$\psi(x,y)=k\;\text{sech}(\frac{k}{\sqrt{2}}(x+y))$.

• (1) You'll get more nibbles if you bait the hook: Copyable code (what you tried, for instance), would make it easy for others to cut, paste & check. (2) Traditionally your equation is called an ODE, not a PDE. You may find this this meta Q&A helpful – Michael E2 Nov 10 '15 at 15:00
• Precisely what solution did you obtain? – bbgodfrey Nov 11 '15 at 5:21
• @bbgodfrey see the edit – user50473 Nov 13 '15 at 9:45
• @MichaelE2 thanks, it should work now – user50473 Nov 13 '15 at 11:04
• @ user50473: there is a difference in your equations for 1D and 2D: the factor in front of psi^3 is 2k for D=1 and 1 for D=2. In fact, your 2D Sech expression solves the 2D equation only with the factor 1, whereas the 1D Sech expression derived from the 2D case by letting y=0 requires the factor 2k to be a soluton. So could you please clarify. – Dr. Wolfgang Hintze Nov 16 '15 at 9:42

Abstract

(Please see section 3.1 for a short answer to the original question of the OP.)

We present here the complete solution of the problem of the OP, which consists of analyzing and solving the ODE numerically and symbolically. It turns out that the essential parameter for the classification of the solutions is the energy w of the equivalent mechanical problem.

The main problem of the OP was to reconcile the basic solution given by DSolve in terms of just the function JacobiSN containing two constants of integration C[1] and C[2] with the explicitly real valued solution which, in its most compact form, is given in terms of the function JacobiCN.

We show that, as expected, the reconciliation is possible by an appropriate choice of the integration constants, but we have found no way to let Mathematica perform the necessary transformations of the Jacobi elliptic functions automatically, therefore these have been done manually.

This answer is organized as follows: we start (section 1) with a qualitative analysis of the solution types in terms of an analogue mechanical problem, then (section 2) we provide the numerical and the exact symbolical solution, and derive (section 3) the exact solution from the basic form, and accompany this with other approaches, such as the integration of the energy equation. In section 4 we provide additional information and discuss the results.

Section 1: reformulation, qualitative picture of solutions

The ODE to be solved is

eq0 = lhs == 0;


where

$$\text{lhs}=\frac{1}{2} k^2 \psi (z)-2 k \psi (z)^3-\frac{\psi ''(z)}{2}$$

First we simplify things by making the ODE dimensionless.

eq = Simplify[
0 == lhs /. {ψ[z] -> x[z k] Sqrt[k/2], ψ''[z] ->
k^2 Sqrt[k/2] x''[z k]} /. z -> t/k, k > 0]

(*
Out[54]= 2 x[t]^3 + (x^′′)[t] == x[t]
*)


It will be convenient in the following to adopt the wording of classical mechanics. The ODE eq. describes the 1D motion of a particle subject to a conservative force with the potential

U = - x^2 + x^4


The energy

w = x'[t]^2-x[t]^2+ x[t]^4


is conserved in time. Actually, the true energy is w/2, but we still keep calling w the energy.

Now let us look for the qualitative picture of the solutions. For a given energy w the motion of the particle is limited to certain intervals in x determined by the obvious requirement x'[t]^2 >= 0 as shown in the graph, in which the characteristic points are given by

solx = x /. Solve[w == -x^2 + x^4, x];
{xA, xB} = {solx[[3]], solx[[4]]};
{xC, xM, xD} = {solx[[3]], solx[[2]], solx[[4]]};
{xE, xF} = {solx[[2]], solx[[4]]};
xG = 1/Sqrt[2];
(* the cell to produce the graph has been closed to save space *)


We have to distinguish three regions of w:

Region I (w > 0)
the point oscillates between the symmetric points A and B symmetrically about x = 0

Region II (-1/4 <= w < 0)
the point oscillates between E and F (or symmetrically in the region x<0, not shown here), i.e. x is never 0. For w = -1/4 the point remains at rest at G

Region III (w = 0)
the point for x > 0 will move between D and M towards M (after having been reflexted in D, if it started out with positive velocity) in the manner of a loopswing moving slowly to its upper position at which it comes to a rest only in infinite time (this is the "soliton" solution Sech)

Furthermore we can state that the oscillations will be increasingly harmonic either for w>>1 or for w close to -1/4.

Section 2: symbolical and numerical solution for the three regions

The symbolic expressions will be derived in section 3.

Region I

xI[t_] = Sqrt[(1 + Sqrt[1 + 4 w])/2]
JacobiCN[t (1 + 4 w)^(1/4), 1/2 (1 + 1/Sqrt[1 + 4 w])]; (* w > 0 *)

With[{ww = 1/2, tt = 20},
x0 = xB /. w -> ww;
x0p = Sqrt[ww + x0^2 - x0^4];
xn = x[t] /. NDSolve[eq && (x[0] == x0) && (x'[0] == x0p), x[t], {t, 0, tt}];
Plot[{xn, 0.1 + xI[t] /. w -> ww}, {t, 0, tt},
PlotRange -> {1.1 xA /. w -> ww, 1.1 xB /. w -> ww},
PlotLabel ->
"Numerical solution of ODE eq: region I\nEnergy w = " <>
ToString[ww, InputForm] <>
"\nInitial position x0 = xB\nblue curve = numerical solution\nyellow \
curve = symbolic solution (slightly shifted)", AxesLabel -> {"t", "x(t)"}]]


Region II

xII[t_] = Sqrt[(1 + Sqrt[1 + 4 w])/2]
JacobiDN[(t (1 + Sqrt[1 + 4 w])^(1/4))/2^(1/4), (1 + 4 w - Sqrt[1 + 4 w])/(
2 w)];(* -1/4 < w < 0 *)

With[{ww = -1/8, tt = 10},
x0 = xF /. w -> ww;
x0p = Sqrt[ww + x0^2 - x0^4];
xn = x[t] /. NDSolve[eq && (x[0] == x0) && (x'[0] == x0p), x[t], {t, 0, tt}];
Plot[{xn, xII[t] /. w -> ww}, {t, 0, tt}, PlotRange -> {-1.5, 1.5},
PlotLabel ->
"Solution of ODE eq: region II\nEnergy w = " <> ToString[ww, InputForm] <>
"\nInitial position x0 = xF\nblue curve = numerical solution\nyellow \
curve = symbolic solution", AxesLabel -> {"t", "x(t)"}]]


Region III

We obtain

xIII[t_] = Sech[t]; (* w = 0 *)


as a limit of region I or of region II.

Region I -> region III

xI[t] /. w -> 0

(* Out[312]= Sech[t] *)

With[{ww = 10^(-5), tt = 20},
x0 = xF /. w -> ww ;
x0p = Sqrt[ww + x0^2 - x0^4];
xn = x[t] /. NDSolve[eq && (x[0] == x0) && (x'[0] == x0p), x[t], {t, 0, tt}];
Plot[{xn, xIII[t] + 0.1}, {t, 0, tt}, PlotRange -> {-1.5, 1.5},
PlotLabel ->
"Numerical solution of ODE eq: region I close to region III\nEnergy w = " <>
ToString[ww, InputForm] <>
"\nInitial position x0 = xB\nblue curve = numerical solution\nyellow \
curve = symbolic solution (slightly shifted upwards)",
AxesLabel -> {"t", "x(t)"}]]


Region II -> region III

JacobiDN @@ Limit[List @@ (xII[t]/Sqrt[((1 + Sqrt[1 + 4 w])/2)]), w -> 0]

(* Out[313]= Sech[t] *)

With[{ww = -10^(-5), tt = 20},
x0 = xB /. w -> ww ;
x0p = Sqrt[ww + x0^2 - x0^4];
xn = x[t] /. NDSolve[eq && (x[0] == x0) && (x'[0] == x0p), x[t], {t, 0, tt}];
Plot[{xn, xIII[t] + 0.1}, {t, 0, tt}, PlotRange -> {-1.5, 1.5},
PlotLabel ->
"Numerical solution of ODE eq: region II close to region III\nEnergy w = " \
<> ToString[ww, InputForm] <>
"\nInitial position x0 = xF\nblue curve = numerical solution (= symbolic \
solution)\nyellow curve = symbolic solution (slightly shifted upwards)",
AxesLabel -> {"t", "x(t)"}]]


Section 3: Exact symbolic solution

Section 3.1: Deriving the exact solutions from the basic solution

Here we show how the solution by DSolve[eq] can be transformed, by chosing apropriate integration constants, to a simple unified form in terms of the function JacobiCN. This was the original task requested by the OP.

We repeat the ODE to be solved

eq

(*
Out[5]= 2 x[t]^3 + (x^\[Prime]\[Prime])[t] == x[t]
*)


DSolve[eq] gives two solutions which differ only in the sign. Hence it is sufficient to consider just one of them. Also, without loss of generality we rename the constants of integration.

Mathematica finds immediately

x[t] /. DSolve[eq, x[t], t][[1]] /. {C[1] -> w, C[2] -> c} // Simplify


During evaluation of In[456]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>

$$- i \frac{ \text{sn}\left(\frac{\sqrt{-(c+t)^2 \left(\sqrt{4 w+1}+1\right)}}{\sqrt{2}}|\frac{1-\sqrt{4 w+1}}{\sqrt{4 w+1}+1}\right)}{\sqrt{2} \sqrt{\frac{1}{\sqrt{4 w+1}-1}}}$$

The designation of the constant C1 is not accidental; in fact it is the energy w of the mechanical problem, in other words, it is the first integral of eq.

After a manual simplification in the factors we define the basic soution as

x0[t_, w_] := -I Sqrt[(-1 + Sqrt[1 + 4 w])/2]
JacobiSN[I (c + t) Sqrt[ ((1 + Sqrt[1 + 4 w])/2)], (1 - Sqrt[1 + 4 w])/(
1 + Sqrt[1 + 4 w])]


Now it turns out that making the following choice for the integration constant c

cw = 1/(1 + 4 w)^(1/4) EllipticK[1/2 (1 + 1/Sqrt[1 + 4 w])];

repc = c -> cw;


x[t_, w_] :=
Sqrt[(1 + Sqrt[1 + 4 w])/2 ]
JacobiCN[t (1 + 4 w)^(1/4), 1/2 (1 + 1/Sqrt[1 + 4 w])]


$$x(t,w)=\sqrt{\frac{1}{2} \left(\sqrt{4 w+1}+1\right)} \text{cn}\left(t \sqrt[4]{4 w+1}|\frac{1}{2} \left(1+\frac{1}{\sqrt{4 w+1}}\right)\right)$$

It came as a surprise to me and is gratifying that this expression for the solution of the ODE is valid for all physically meaningful values of w, i.e. for w>-1/4. Hence the distinction of the expressions for the exact solutions according to regions, as was done in section 2, is not necessary.

The "unifying" idea is to place the origin of time (t=0) so that the motion starts always at the maximum attainable value of x for given w, i.e. with velocity x'[0] zero.

In the limit of w->0 we find the expected Sech

x[t, 0]

(* Out[5]= Sech[t] *)


Now we plot the solution for varying parameter w, and distingush the two regions w>0 and w<0 for better viewing:

Plot3D[x[t, w], {t, -10, 10}, {w, 0, 1}, PlotRange -> {-1.2, 1.2},
PlotLabel -> "Solution of ODE x(t,w) for w>=0",
AxesLabel -> {"t", "w", "x"}]


Plot3D[x[t, w], {t, -10, 10}, {w, -1/4, 0}, PlotRange -> {-1.2, 1.2},
PlotLabel ->
"Solution of ODE x(t,w) for -1/4 \[LessEqual] w \[LessEqual] 0",
AxesLabel -> {"t", "w", "x"}]


Unfortunately, Mathematica seems to be not familiar enough with the transformation properties of the Jacobi elliptic functions to be able to (Full)Simplify the expression x0 after the replacement of c automatically.

Therefore we shall show how to perform the steps manually indicating the relevant formulas of the "Handbook of Mathematical Functions" by Abramovich and Stegun, June 1964, called AS in what follows.

Step 0: we start out with sn(u | m) where m<0

Step 1: removal of explicitly imaginary expressions

AS 16.20

sn (I u | m) = I sc (u | n), n = 1-m > 1


Step 2: transformation to have second parameter in the "natural" region between 0 and 1

AS 16.11 and 16.3

sc = sn/cn

sn (u | n) -> n^(-1/2) sn (u n^(1/2) | 1/n)
cn (u | n) -> dn (u n^(1/2) | 1/n)

sc (u | n) -> n^(-1/2) sn (u n^(1/2) | 1/n) / dn (u n^(1/2) | 1/n)
= n^(-1/2) sd (u n^(1/2) | 1/n)


Step 3: translation of the argument by the complete elliptical integral K to have

x'(0) = 0


AS 16.8

n^(-1/2) sd (u n^(1/2) + K (1/n) n^(1/2) | 1/n)
->
n^(-1/2) (1 - 1/n)^(-1/2) cn (u n^(1/2) | 1/n)


Section 3.2: Integration of the energy equation

An alternative approach to find the exact solution is the integration of the energy equation (the first integral of the original ODE) given by

eqw = x'[t]^2 == w + x[t]^2 - x[t]^4;


Section 3.2.1 The case w = 0

We start with the most simple case of zero energy (w = 0).

As we are looking for the soliton-solution which must have x'[0]==0.

But, unfortunately, DSolve can't deal with neither the initial condition x'[0]==0 or the equivalent x[0]==1.

Hence we DSolve without initial conditions and obtain

solw0 = x[t] /. DSolve[(eqw /. w -> 0), x[t], t]

(*
Out[351]= {-I Coth[t - C[1]] Sqrt[1 - Tanh[t - C[1]]^2],
I Coth[t - C[1]] Sqrt[1 - Tanh[t - C[1]]^2], -I Coth[t + C[1]] Sqrt[
1 - Tanh[t + C[1]]^2], I Coth[t + C[1]] Sqrt[1 - Tanh[t + C[1]]^2]}
*)


It is sufficient to take just one of these and define (letting also C[1]->z)

xw0[t_] = solw0[[1]] /. C[1] -> z

(*
Out[387]= -I Coth[t - z] Sqrt[1 - Tanh[t - z]^2]
*)


This looks like a completely imaginary solution. It took me some time to get this reconciled with the expected Sech. But wait! This need not be imaginary if we allow z to be complex, as we shall see here:

Let's try to formally calculate z from the initial condition x'[0] == 0

xw0p = D[xw0[t], t] /. t -> 0 // Simplify

(*
Out[384]= I Coth[z]^2 Sqrt[Sech[z]^2]
*)

Reduce[0 == xw0p, z ]

(* Out[385]= False *)


Bad luck! OK, obviously the square root prevents Reduce[] from giving results. Hence we try the equivalent condition

Reduce[0 == xw0p^2, z, GeneratedParameters -> K ]

(*
Out[386]= K[1] ∈ Integers && z == (I π)/2 + I π K[1]
*)


Hence we have found that z must be purely imaginary.

Taking z -> I π/2 gives

xw1[t_] = Simplify[xw0[t] /. z -> I π/2, t > 0]

(*
Out[394]= Sech[t]
*)


Section 3.2.2 The energy equation for $w\neq 0$

The energy integral, which gives t + const, can be done explicitly. Here is the antidervative

FullSimplify[Integrate[1/Sqrt[w + x^2 - x^4], x], w > 0]

(*
Out[158]= -((
I Sqrt[2] EllipticF[
I ArcSinh[(Sqrt[2] x)/Sqrt[-1 + Sqrt[1 + 4 w]]], (-1 - 2 w + Sqrt[
1 + 4 w])/(2 w)])/Sqrt[1 + Sqrt[1 + 4 w]])
*)


For -1/4 < w < 0 we have the same result but now I ArcSinh is replaced by ArcSin:

(*
Out[158]= -((
I Sqrt[2] EllipticF[
ArcSin[(Sqrt[2] x)/Sqrt[-1 + Sqrt[1 + 4 w]]], (-1 - 2 w + Sqrt[
1 + 4 w])/(2 w)])/Sqrt[1 + Sqrt[1 + 4 w]])
*)


Section 4: Generalizations and discussion

Section 4.1: General potential with "confinement"

Consider a potential repelling at $x = 0$ in the simplest forU = x^2, and going to +infinity for|x| >> 1.

So the particle is "globally" confined between the "walls" of the potential on both side of the x-axis.

First, locally, the equation of motion close to $x = 0$ is

eqx0 = x''[t] == x[t];


and the solution is

x[t] /. DSolve[eqx0, x[t], t]

(* Out[104]= {E^t C[1] + E^-t C[2]} *)


Letting C[1]->0 to keep the solution finite for t -> [Infinity] we obtain

xx0[t_] = a Exp[-t]

(* Out[111]= a E^-t *)


The (positive) constant "a" is the (negative) velocity of the point at t = 0.

a /. Solve[v0 == xx0'[0], a][[1]]

(* Out[120]= -v0 *)


Hence we have the following scenario:

The point starts at x = x1 = |v0| with a negative velocity v0 towards the repelling potential at x = 0. It will approach the summit of the potential exponentially for t ->[Infinity].

This is the typical behaviour of a loop swing and of situations leading to a sech(t) solution.

Section 4.2: The soliton solution for $x''=x-x^{2 n+1}$

Now let us consider the special confinement of the form x^(2k+1)

The energy integral for w = 0 can be calculated immediately:

fi = Integrate[1/Sqrt[u^2 - u^(2 k + 2)], {u, x, 1},
Assumptions -> {k \[Element] Integers, k > 0, 0 < x < 1}]

(* Out[6]= ArcTanh[Sqrt[1 - x^(2 k)]]/k *)


Inverting gives

sol = x /. FullSimplify[Solve[t == fi, x]] // Quiet

(* Out[7]= {(-Sqrt[Sech[k t]^2])^(1/k), (Sech[k t]^2)^((1/2)/k)} *)


of which the positive solution is adequate.

Hence we have found that for a positive integer n the equation

x''[t] == x[t] - x[t]^(2n+1)


has the solution

x[t] = Sech[n t]^(1/n)

k = 1 -> Sech[t]
k = 2 -> Sech[2t]^(1/2)
...


For t -> [Infinity] we find

Exp[-t] Limit[Exp[t] Sech[k t]^(1/k), t -> \[Infinity], Assumptions -> k > 0]

(* Out[27]= 2^(1/k) E^-t *)


Which is the exponential decay expected close to x = 0.

In the limit of large n we have

Limit[Sech[n t]^(1/n), n -> \[Infinity], Assumptions -> t > 0]

(* Out[] E^-t *)


n->[Infinity] means that the potential has hard walls at x = [PlusMinus] 1. Inside the box there is just the repelling potential from x = 0.

Hence the limit of large n again leads to the expected result, the simple exponential decay.

Section 4.3: The solution of $\nabla ^2u=u-u^3$ in 2D and 3D

For higher geometrical dimensions we confine ourselves to the case of radial symmetry.

Hence for general dimension D we consider

$$\frac{(D-1) u'(r)}{r}+u''(r)=u(r)-2 u(r)^3$$

It turns out that the case D = 2 is sufficient to grasp the general picture also for D > 2.

As before we translate the problem to the mechanical analogue and write the equation of motion as

eq = x''[t] + 1/t x'[t] == x[t] - 2 x[t]^3;


The new feature is that now due to the term with the first derivative we have friction, and this dissipative process forces the "particle" to come to a rest at t = oo.

Applying the usual procedure to derive energy conservation, i.e multiplying by x'[t] we can write the equation as

$$\frac{\partial w(t)}{\partial t}=-\frac{2 x'(t)^2}{t}$$

where the conservative part of the energy is

$$w(t)=x'(t)^2+x(t)^4-x(t)^2$$

We see that the quantity w decreases from its initial value to some finite value with the passage of time. The asymptotic value is given by x'' = x' = 0 which, from the equation of motion leads to xf = 1/Sqrt2and hence wf = -1/4.

We shall now solve the equation numerically. At first sight it might seem that the factor 1/t causes problems. But this is not the case. We shall, in the numeric treatment, start with a very small time [Epsilon] instead of zero.

The solution will be displayed in two forms: the usual dependence of position as a function of time and in a non-standard form which I like to call energy trajectory. The latter has been used in the discussion of the case D=1 above. But there the trajectories were simply oscillations between turning points. Here we see the non-trivial dissipative history of the particle in the energy picture.

The normal history is that after starting with the given energy the particle may make some oscillations in the big energy container but then "decide" to take the left or the right where it comes to rest asymptotically, hereby oscillating in the trough.

However with a suitable choice of the initial conditions the particle comes to rest at x = 0. Here asymptotically the particle moves "uphill" like a loop swing (as discussd earlier). The space time diagram is then a one-sided "soliton".

The following code provides the solution in this form. We give examples of two sets of parameters here. The reader is invited to study more cases.

The oscillating case

s = 10; (* scale factor for plots *)
\[Epsilon] = 10^-2;
tmax = 100;
x0 = 1.2(*1.5959735*);
v0 = -1; (*-1*)
ww[t_] := xx'[t]^2 - xx[t]^2 + xx[t]^4
w0 = v0^2 - x0^2 + x0^4; (* initial energy *)

(* numerical solution *)
xx[t_] = x[t] /.
NDSolve[(x''[t] + x'[t]/t == x[t] - 2 x[t]^3) && x[\[Epsilon]] == x0 &&
x'[\[Epsilon]] == v0, x[t], {t, \[Epsilon], tmax}][[1]];

(* plot x(t) *)
Plot[{-1/Sqrt[2], 1/Sqrt[2], xx[t]}, {t, \[Epsilon], tmax},
AxesLabel -> {"t", "x"},
PlotLabel ->
"Solution of \!$$\*SuperscriptBox[\(x$$, $$\[Prime]\[Prime]$$,\n\
MultilineFunction->None]\)(t)+\!$$\*FractionBox[\(\*SuperscriptBox[\"x\", \"\ \[Prime]\",\nMultilineFunction->None](t)$$, $$t$$]\)\[LongEqual]x(t)-2 \
x(t\!$$\*SuperscriptBox[\()$$, $$3$$]\), x(0) = " <> ToString[x0, InputForm] <>
", v(0) = " <> ToString[v0, InputForm] <>
"\nPosition x as a function of t\n",
PlotStyle -> {{Black, Thin, Dashed}, {Black, Thin, Dashed}, Red},
PlotRange -> {-1.1, 1.1 (* x0 *)}]

(* container plot for energy trajectory *)
s = 10;
xmax = x0; xmin = -xmax;
p0 = Plot[s (-x^2 + x^4), {x, xmin, xmax},
AxesLabel -> {"x", "w (red), U (blue)"},
PlotLabel ->
"Solution of \!$$\*SuperscriptBox[\(x$$, $$\[Prime]\[Prime]$$,\n\
MultilineFunction->None]\)(t)+\!$$\*FractionBox[\(\*SuperscriptBox[\"x\", \"\ \[Prime]\",\nMultilineFunction->None](t)$$, $$t$$]\)\[LongEqual]x(t)-2 \
x(t\!$$\*SuperscriptBox[\()$$, $$3$$]\), x(0) = " <> ToString[x0, InputForm] <>
", v(0) = " <> ToString[v0, InputForm] <> ", w(0) = " <>
ToString[w0, InputForm] <>
"\nEnergy trajectory: w(t) = \!$$\*SuperscriptBox[\(v$$, \
$$2$$]\)-\!$$\*SuperscriptBox[\(x$$, $$2$$]\)+\!$$\*SuperscriptBox[\(x$$, $$4\$$]\) as a function of x(t)\n"];

(* plot energy trajectory *)
p1 = ParametricPlot[{xx[t], s ww[t]}, {t, \[Epsilon], tmax},
PlotStyle -> {Thin, Red}];
Show[{p0, p1}]


The "soliton" case

By varying the initial position x0 we arrive at the picture where the particle moves to x = 0.

s = 10; (* scale factor for plots *)
\[Epsilon] = 10^-2;
tmax = 10 (* = 30 *);
x0 = 1.5959735;
v0 = -1; (*-1*)
ww[t_] := xx'[t]^2 - xx[t]^2 + xx[t]^4
w0 = v0^2 - x0^2 + x0^4; (* initial energy *)

(* numerical solution *)
xx[t_] = x[t] /.
NDSolve[(x''[t] + x'[t]/t == x[t] - 2 x[t]^3) && x[\[Epsilon]] == x0 &&
x'[\[Epsilon]] == v0, x[t], {t, \[Epsilon], tmax}][[1]];

(* plot x(t) *)
Plot[{-1/Sqrt[2], 1/Sqrt[2], xx[t]}, {t, \[Epsilon], tmax},
AxesLabel -> {"t", "x"},
PlotLabel ->
"Solution of \!$$\*SuperscriptBox[\(x$$, $$\[Prime]\[Prime]$$,\n\
MultilineFunction->None]\)(t)+\!$$\*FractionBox[\(\*SuperscriptBox[\"x\", \"\ \[Prime]\",\nMultilineFunction->None](t)$$, $$t$$]\)\[LongEqual]x(t)-2 \
x(t\!$$\*SuperscriptBox[\()$$, $$3$$]\), x(0) = " <> ToString[x0, InputForm] <>
", v(0) = " <> ToString[v0, InputForm] <>
"\nPosition x as a function of t\n",
PlotStyle -> {{Black, Thin, Dashed}, {Black, Thin, Dashed}, Red},
PlotRange -> {-2, 1.1 x0}]

(* container plot for energy trajectory *)
s = 10;
xmax = x0; xmin = -xmax;
p0 = Plot[s (-x^2 + x^4), {x, xmin, xmax},
AxesLabel -> {"x", "w (red), U (blue)"},
PlotLabel ->
"Solution of \!$$\*SuperscriptBox[\(x$$, $$\[Prime]\[Prime]$$,\n\
MultilineFunction->None]\)(t)+\!$$\*FractionBox[\(\*SuperscriptBox[\"x\", \"\ \[Prime]\",\nMultilineFunction->None](t)$$, $$t$$]\)\[LongEqual]x(t)-2 \
x(t\!$$\*SuperscriptBox[\()$$, $$3$$]\), x(0) = " <> ToString[x0, InputForm] <>
", v(0) = " <> ToString[v0, InputForm] <> ", w(0) = " <>
ToString[w0, InputForm] <>
"\nEnergy trajectory: w(t) = \!$$\*SuperscriptBox[\(v$$, \
$$2$$]\)-\!$$\*SuperscriptBox[\(x$$, $$2$$]\)+\!$$\*SuperscriptBox[\(x$$, $$4\$$]\) as a function of x(t)\n"];

(* plot energy trajectory *)
p1 = ParametricPlot[{xx[t], s ww[t]}, {t, \[Epsilon], tmax},
PlotStyle -> {Thin, Red}];
Show[{p0, p1}]
`

But we have cheated a bit: due to numerical instability the particle finally falls into the trough and oscillates about the asymptotic position. You can see this by increasing tmax to say 30.

Section 4.4: Discussion

Some remarks concerning the soliton solution should be made.

We have seen that non oscillating solutions occur if the particle climbes the potential hill at t = 0. Actually this is the tail of the asymptotic behaviour.

We can call the finite non zero part of this structure "soliton", although it must be noted that the position of the particle does not remain finite for t = 0 except for the case D = 1.

Another point is the symmetric structure such as in the Sech[t]. This appears if we include negative times in our consideration. If we would do this for D > 1 we would also have the symmetry.

Concluding we have shown that the mechanical analogue of our ODE can be a helpful device to understand the character of the Solutions. We also have extended this approach to disspative systems by introducing the concept of an energy trajectory depending on the time as a parameter.

• There ought to be a Textbook badge. :) Sorry, can't upvote more than once. Cheers! – Michael E2 Nov 25 '15 at 21:13
• @Michael E2: thank you very much. It was great fun for me to take the opportunity to delve into Jacobi elliptic functions, elliptic functions, inverse functions, non linear differential equations, and I have still some points to make here. Don't worry, not too many ;-) By the way, I have found out (tacitly) that, despite of its name, a non linear Schrödinger equation has nothing to do with quantum mechanics. – Dr. Wolfgang Hintze Nov 25 '15 at 22:57
• @Dr.WolfgangHintze I am so waiting for section 4!!! By the way, the ref in question is this article: J. Zittartz and J. S. Langer Phys. Rev. 148, 741, 1966 where also the different coefficient of the 2D problem is clarified. I approve the Textbook badge. – user50473 Dec 1 '15 at 16:18
• @user50473: oh yes, you are right, and my bad concious is there to remind me. I ask politely for some more patience. – Dr. Wolfgang Hintze Dec 1 '15 at 20:21
• @Dr.WolfgangHintze I was reading again all this, just wondering whether you abandoned the project. Would really love to see how this ends. – user50473 Mar 1 '16 at 17:09