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;
leads to the expression
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.
