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?
system2[[1;;2]]
ata == 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 containingNeeds
. $\endgroup$y0
in particular. $\endgroup$