# Solving a stiff nonlinear ODE system

The system I am trying to solve is simple, but looks pretty stiff and I have unsuccessfully tried to solve it with StiffnessSwitching. It is the following one:

x0=.001;
y0=10^64;
const=3*10^9;
h[a_] := Sqrt[y[a] + x[a]];
Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];
system = {
a*y'[a] + 3*y[a] == -const*y[a]/h[a],
a*x'[a] + 4*x[a] == const*y[a]/h[a],
x == x0,
y == y0
};

system2 = ReplacePart[system /. {x -> (F[#]^2 &), y -> (G[#]^2 &)}, {{3} -> F == Sqrt[x0], {4} -> G == Sqrt[y0]}];
{nds1, nds2} = NDSolveValue[system2, {F, G}, {a, 1, 10^7}, Method -> "StiffnessSwitching", MaxSteps -> 10^6, WorkingPrecision -> MachinePrecision];


where system2 has been defined to ensure the positiveness of the functions. But here I receive the stiffness message:

NDSolveValue::ndsz: At a == 1., step size is effectively zero; singularity or stiff system suspected.

I have also tried with ExplicitRungeKutta, but nothing happened. It is clear that the initial condition y0 is extremely high, but I need it. Trying with a lower y0=10^16 it works, but only if I do not include const.

Any hint?

• Evaluating system2[[1;;2]] at a == 1 yields {F' -> 4.74342*10^42, G' -> -1.5*10^32}. Mathematica probably concludes that these enormous values indicate stiffness. Try rescaling your variables to obtain smaller coefficients. By the way, you do not need the code containing Needs. Sep 2, 2019 at 19:14
• Which coefficients are you referring to?
– Lele
Sep 4, 2019 at 14:43
• y0 in particular. Sep 4, 2019 at 14:48
• Yes, with lower y0 it works, as I said. But, physically, I need that high value!
– Lele
Sep 4, 2019 at 15:27
• Try rescaling your variables so that these numbers are not so large. Sep 4, 2019 at 15:47

The computation can be performed as follows. First, solve for {x, y} instead of {F, G}, because the ODEs are simpler, and then obtain {F, G} by taking the square roots of {x, y}. Even then,

NDSolveValue[system, {x, y}, {a, 1, 10^7}]


immediately fails, claiming that the "step size is effectively zero; singularity or stiff system suspected". Perhaps, NDSolve concludes this, because the second equation in system is approximately a*x'[a] + 4*x[a] == 3 10^41. Rescaling dependent variables to eliminate such large constants. i.e., y by y0 and x by const Sqrt[y0], probably would allow the computation to proceed. In addition, changing the independent variable from a to Log[a] reduces the enormous domain of integration to a more manageable size. As it turns out, this last transformation alone is sufficient to obtain the solution. With b == Log[a], system becomes

system1 = {y'[b] + 3*y[b] == -const*y[b]/h[b], x'[b] + 4*x[b] == const*y[b]/h[b],
x == x0, y == y0};


which is, incidentally, autonomous. Then,

NDSolveValue[system1, {x[b], y[b]}, {b, 0, Log[10^7]}];
LogPlot[%, {b, 0, Log[10^7]}, ImageSize -> Large,
AxesLabel -> {b, "{x,y}"}, LabelStyle -> {Bold, Black, 15}]
Plot[First@%%, {b, 0, 2}, PlotRange -> All, ImageSize -> Large,
AxesLabel -> {b, x}, LabelStyle -> {Bold, Black, 15}]  The second plot, a blowup of x near the origin, shows that it grows rapidly but linearly there. Why transforming only the independent variable is sufficient to obtain a solution is unclear. Perhaps, NDSolve uses a different computational method for autonomous ODEs.

• Great! Thank you! It was much easier than I had thought... last question: how do I use the solutions from NDSolveValue? I tried with {xsol,ysol}=NDSolveValue[system1, {[b], y[b]}, {b, 0, Log[10^7]}];, but when I evaluate one of them, for example, as xsol it returns only InterpolatingFunction[...][b]. I would like to use those solutions to define other functions and so on... thank you!
– Lele
Sep 5, 2019 at 10:27
• Try xsol /. b -> 10. Alternatively, start with {xsol, ysol} = NDSolveValue[system1, {x, y}, {b, 0, Log[10^7]}]; (Look at the output from NDSolveValue using {x[b], y[b]} vs {x, y}` to see the difference.) Sep 5, 2019 at 11:50