# How to solve a stiff nonlinear second-order ode?

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["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];

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 (:=).

-
You should define functions with := (SetDelayed) instead of =. – shrx Dec 10 '13 at 23:43
@shrx : Thanks for your suggestion. I tried to do what you suggested (i.e. φ[λ_]:= a0 + a1 λ + a2 λ^2 + a3 λ^3 + a4 λ^4 /. par), but the result does not change – Petrus Dec 11 '13 at 11:13
I think your solution blows up (approaches infinity) near x = 0.379441. I suspect everything is correct. – Michael E2 Dec 11 '13 at 14:55

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.

-
Thanks very much for your answer. Anyway, I'm new in numerical solution of ode. What does it mean practically for me? I apologize in advance if my doubt is trivial – Petrus Dec 12 '13 at 11:31
@Petrus Your doubt is fine. Practically it means that the solution is accurate, the stiffness message is a warning not an error, the solution approaches an asymptote near 0.397, and thus the domain of the solution cannot be extended beyond that point (for mathematical reasons). And, practically, you will have to interpret what it means for the problem you're working on. I added a more extensive analysis for fun. – Michael E2 Dec 12 '13 at 14:43
@ Michael thanks again for your answer, it is extremely accurate! – Petrus Dec 13 '13 at 23:53
@ Michael:if you have time, should you have a look at the "original" problem? mathematica.stackexchange.com/questions/38798/… – Petrus Dec 14 '13 at 12:46