Computing Separatrix of Second-Order Nonlinear Autonomous ODE

Question

Numerically solving differential equations of the form

{x''[t] == f[x[t], x'[t]], x[0] == x0, x'[tm] == 0}

where tm is large, typically is quite difficult, if the linearization of f[x[tm], 0] about the solution of f[x[tm], 0] == 0 is a x[tm] with a > 0; i. e., if the solution coincides with the separatrix of the ODE at large tm. This is because even infinitesimal numerical errors grow exponentially there.

Section 2.2 of the Handbook of Exact Solutions for Ordinary Differential Equations uses the transformation x'[t] == w[x] to reduce the order of the ODE by one,

w'[x] w[x] == f[x, w[x]]

Then, true to its name, the book solves several particular cases of f analytically. I wish to know whether solving this equation numerically for autonomous second-order ODEs without analytical solutions and then inverting the transform would be easier than numerically solving the original ODE.

As a test case, consider the equation

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

which is compact but nonetheless nonlinear in both x[t] and x'[t]. The separatrix lies at x == 0. It can, in fact, be solved directly by

tm = 10; 1/2 (1 - Sqrt[5]);
NDSolve[{eq, x[0] == 1/2, x'[tm] == c x[tm]}, x, {t, 0, tm}, 
    Method -> {"Shooting", "StartingInitialConditions" -> {x[0] == 1/2, 
    x'[0] == -30724/100000}}, WorkingPrecision -> 30];
Plot[x[t] /. %, {t, 0, tm}, PlotRange -> All, AxesLabel -> {x, w}]

but only with an excellent initial guess. Even then, increasing tm much beyond ten causes NDSolve to fail.

My question, then, is, can this ODE be solved numerically by means of the transformation above with larger tm and without an excellent initial guess.

Answer 1

The ODE in the question can be solved as follows.

eqw = eq /. {x''[t] -> w'[x] w[x], x'[t] -> w[x], x[t] -> x}
(* w[x] w'[x] == x (1 - x^2) + w[x] + w[x]^2 *)

The separatrix occurs at x == 0, and w[0] == 0 as well. To understand the behavior near x == 0, linearize the equation and assume that w[x] is approximately equal to a x, with a a constant. The linearized equation then is a^2 x == x + a x, and a is determined by

a /. Solve[a^2 x == x + a x, a]
(* {1/2 (1 - Sqrt[5]), 1/2 (1 + Sqrt[5])} *)

If 1 > x0 > 0, the desired expression for a must be negative; conversely, it must be positive for 0 > x0 > -1. Let us choose x0 == 1/2, in which case set

c = First@%
(* 1/2 (1 - Sqrt[5]) *)

Then, the equation for w is solved in about a second by

offset = 10^-10;
sw = Flatten@NDSolve[{eqw, w[offset] == c offset}, w[x], {x, offset, 1/2}, 
    WorkingPrecision -> 30];
Plot[(w[x] /. sw), {x, offset, 1/2}, AxesLabel -> {x, w}]

Finally, the transformation is inverted in a few seconds by

xp = x'[t] == Piecewise[{{(w[x] /. sw /. x -> x[t]), x[t] > offset}, 
    {c, offset > x[t] > 0}}];
Flatten@Quiet@NDSolve[{xp, x[0] == 1/2}, x[t], {t, 0, 40}, WorkingPrecision -> 30];
Plot[x[t] /. %, {t, 0, 40}, PlotRange -> All, AxesLabel -> {t, x}]

In effect, the approach just described decomposes one second-order boundary-value problem into two first-order initial-value problems. For completeness, note that x[t] == offset occurs at

t /. FindRoot[(x[t] /. %%) == offset, {t, 30}]
(* 36.0493 *)

The last NDSolve integration could have been stopped there, and an approximate analytical solution, offset Exp[c (t - 36.0493)] used for larger t.

Incidentally, although the approach described here could be applied when f depends on x only, the method of the answer to 140328 seems simpler.