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 (:=).
:=
(SetDelayed
) instead of=
. $\endgroup$x = 0.379441
. I suspect everything is correct. $\endgroup$