4
$\begingroup$

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}]

enter image description here

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.

$\endgroup$
4
$\begingroup$

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}]

enter image description here

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}]

enter image description here

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.