3
$\begingroup$

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["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
system = {
           a*y'[a] + 3*y[a] == -const*y[a]/h[a],
           a*x'[a] + 4*x[a] == const*y[a]/h[a],
           x[1] == x0,
           y[1] == y0
         };

system2 = ReplacePart[system /. {x -> (F[#]^2 &), y -> (G[#]^2 &)}, {{3} -> F[1] == Sqrt[x0], {4} -> G[1] == 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?

Improve this question
$\endgroup$
5
  • $\begingroup$ Evaluating system2[[1;;2]] at a == 1 yields {F'[1] -> 4.74342*10^42, G'[1] -> -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. $\endgroup$
    – bbgodfrey
    Sep 2, 2019 at 19:14
  • $\begingroup$ Which coefficients are you referring to? $\endgroup$
    – Lele
    Sep 4, 2019 at 14:43
  • $\begingroup$ y0 in particular. $\endgroup$
    – bbgodfrey
    Sep 4, 2019 at 14:48
  • $\begingroup$ Yes, with lower y0 it works, as I said. But, physically, I need that high value! $\endgroup$
    – Lele
    Sep 4, 2019 at 15:27
  • $\begingroup$ Try rescaling your variables so that these numbers are not so large. $\endgroup$
    – bbgodfrey
    Sep 4, 2019 at 15:47

1 Answer 1

Reset to default
5
$\begingroup$

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[0] == x0, y[0] == 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}]

enter image description here

enter image description here

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.

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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[10] it returns only InterpolatingFunction[...][b][10]. I would like to use those solutions to define other functions and so on... thank you! $\endgroup$
    – Lele
    Sep 5, 2019 at 10:27
  • 1
    $\begingroup$ 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.) $\endgroup$
    – bbgodfrey
    Sep 5, 2019 at 11:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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