How to solve a stiff nonlinear second-order ode?

Question

I want to solve a nonlinear second order ODE as an initial value problem (IVP). Unfortunately, for certain values of first derivative at starting point I get a stiff-problem. I tried to overtake this issue looking at Mathematica help, but those advices don't work in my case. My equation is the following:

$$ \gamma(\lambda) \lambda_{zz} + \beta(\lambda) \lambda_z^2 + \tau(\lambda) = T_M $$

where

$$ \varphi(\lambda) = a_0 + a_1 \lambda^2 + a_2 \lambda^2 + a_3 \lambda^3 + a_4 \lambda^4 \\ \tau(\lambda) = \varphi'(\lambda) \\ \gamma(\lambda) = - \frac{h_0^2}{12} \frac{\tau(\lambda)}{\lambda^5} \\ \beta(\lambda)=\frac{1}{2} \gamma'(\lambda) $$

H0 = 1;
par = {
       a0 -> 2.0377638272727268,
       a1 ->  -7.105521894545453,
       a2 -> 9.234000147272726,
       a3 -> -5.302489919999999,
       a4 -> 1.136247839999999
      };

φ[λ_] := a0 + a1 λ + a2 λ^2 + a3 λ^3 + a4 λ^4 /. par;
τ[λ_] := φ'[λ];
γ[λ_] := -(H0^2/12) τ[λ]/λ^5;
β[λ_] := 1/2 γ'[λ];

TM = 0.0059225557575715015;

Eq = γ[λ[x]] D[λ[x], {x,2}] + β[λ[x]] D[λ[x],x]^2 + τ[λ[x]];

xL = 0;
xR = 1;

ym = 1.1666666666666892;
yp = 1;

SYSf = Flatten[NDSolve[{
                Eq == TM,
                λ[0] == ym,
                λ'[0] == yp
               }, λ[x], {x, xL, xR}]
              ]

$ym$ is the value of function in $x = 0$, and $yp$ the first derivative. If I set $yp = 2$ I get the following error:

NDSolve::ndsz: At x == -0.311892, step size is effectively zero; singularity or stiff system suspected. >> NDSolve::ndsz: At x == 0.37944200653091864`, step size is effectively zero; singularity or stiff system suspected. >>

So I tried some advices written in the help:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];

SYSf2 = Flatten[NDSolve[
       {Eq == TM,
        λ[0] == ym,
        λ'[0] == yp},
        λ[x], {x, xL, xR}, Method -> "StiffnessSwitching"]
       ]

or

SYSf3 = Flatten[NDSolve[
     {Eq == TM,
      λ[0] == ym,
      λ'[0] == yp},
      λ[x], {x, xL, xR}, 
      Method -> {"StiffnessSwitching", Method -> {"ExplicitRungeKutta", Automatic}}]
     ]

but with no success. What can I do in order to set $yp \ge 2$?

Petrus

--- EDIT ---

I changed function definitions using SetDelay operator (:=).

Answer 1

I think if you look at the surface defined by the differential equation in λ, λ', λ'' space, you can see that the solution when yp = 2 goes off to infinity. Somewhere in the interval 1.95 < yp < 2, there is a point where the integral curve transitions from being closed to going off to infinity. Three integral curves are plotted on the surface below for yp = 1., 1.95, 2. (blue, green, red, resp.). The other lines indicate the flow of the DE along the surface.

Clear[soln];
soln[yp_] := soln[yp] = λ /. 
   First @ NDSolve[{Eq == TM, λ[0] == ym, λ'[0] == yp}, λ, {x, xL, xR}]

db = First[c /. Solve[Eq == TM /. {λ[x] -> a, λ'[x] -> b, λ''[x] -> c}, c]];
Clear[vf];
vf[a1_] := vf[a1] = First[ (* flow lines at offset a = a1 *)
     StreamPlot[
      {b, db},
      {a, 1, 2}, {b, -5, 10},
      StreamScale -> Full, StreamStyle -> "Line", 
      StreamPoints -> Tuples[{Range[a1, 2, 0.2], Range[-5, 10, 0.75]}]]
     ] /. {{a0_Real, b0_Real} :> {a0, b0, Quiet@Check[db /. {a -> a0, b -> b0}, 10.^6]}};

phasespace = ContourPlot3D[ (* the surface of the DE *)
   Eq == TM /. {λ[x] -> a, λ'[x] -> b, λ''[x] -> c} // Evaluate,
   {a, 1, 2}, {b, -5, 10}, {c, -270, 270},
   Mesh -> None, ContourStyle -> Opacity[0.7], 
   AxesLabel -> {λ[x], λ'[x], λ''[x]}
   ];

λplot[yp_, color_] := (* plot a solution *)
 ParametricPlot3D[{λ[x], λ'[x], λ''[x]} /. λ -> soln[yp] // Evaluate, 
  Evaluate[Prepend[First@soln[yp]["Domain"], x]], PlotStyle -> {Thick, color}];

Show[
 phasespace,
 Graphics3D[{vf[1.03], vf[1.1], vf[1.17]}],
 λplot[1., Blue], λplot[2., Red], λplot[1.95, Darker@Green]
 ]

[The idea for plotting the flow lines on the surface was used in this answer. For some reason StreamPlot would not plot many stream lines, so I combined multiple plots with different stream lines.]

In terms of plotting the solutions λ vs. x, we see seemingly periodic solutions for yp equal to 1 and 1.95. But for yp -> 2., the solution seems to approach an asymptote at approximately x -> 0.379442. Since the integration can't be carried past an asymptote, NDSolve complains. That does not mean there was an error. In this case, it is just a notification that the interval {x, 0, 1} cannot be achieved.

GraphicsGrid[MapThread[
  Plot[Evaluate[soln[#][x]], 
    Evaluate[Prepend[First@soln[#]["Domain"], x]], PlotRange -> #2, 
    PlotLabel -> HoldForm[ym -> #], PlotStyle -> #3] &,
  {{{1., 1.95}, {2., 2.}},
   {{{1, 1.3}, {1, 1.3}}, {Automatic, {0, 100}}},
   {{Blue, Darker@Green}, {Red, Red}}},
  2]
 ]

The change from periodic solutions to one that approaches infinity corresponds to the unstable equilibrium point one sees above in the phase portrait where the red and green curves diverge. The equilibrium point may be found thus:

equilibrium = FindRoot[{
  Eq == TM /. {λ[x] -> a, λ'[x] -> b, λ''[x] -> c},
  b == 0, db == 0},
 {{a, 1.4}, {b, 0}, {c, 0}}]
(* {a -> 1.30794, b -> 0., c -> 1.64351*10^-19} *)

The separatrix may be found approximately by finding an initial point close to the equilibrium point such that the vector field defined by the differential equation points either toward or away from the equilibrium point; that is, the vector should be parallel to the displacement. To get this to work, I needed to use arbitrary precision (tutorial/ArbitraryPrecisionNumbers).

sep = Block[{a0, b0, c0},
  {a0, b0, c0} = {a, b, c} /. equilibrium;
  FindRoot[
   SetPrecision[
    {Eq == TM /. {λ[x] -> a, λ'[x] -> b, λ''[x] -> c},
     b (b - b0) == db (a - a0), (b - b0)^2 + (a - a0)^2 == 10^-6},
    20],
   {{a, 1.3}, {b, -10^-6}, {c, 0.}},
   WorkingPrecision -> 20]
  ]
(* {a -> 1.3079170871961404005,
    b -> 0.00099964512736502608624, 
    c -> -0.037512757244442424374}  *)

Redefining soln to use arbitrary precision numbers, we can use the values above as initial values to solve for the separatrix:

Clear[soln];
soln[yp_] := soln[yp] =
   λ /. First@NDSolve[
      SetPrecision[{Eq == TM, λ[0] == ym, λ'[0] == yp}, 20],
      λ, {x, xL, xR},
      WorkingPrecision -> 20];
sepIF = Block[{ym = a /. sep, xL = -3, xR = 3},
  soln[b /. sep]]

It may be added to the phase space plot above in similar manner, this time with an Orange style. One can see that is slightly off, in that it does not approach the equilibrium but gets very close and keeps on going.