Problems with NDSolve and stiffness

I am trying to solve an ODE in chemical kinetics:

\begin{align*} \frac{\mathrm d[x]}{\mathrm dt} &= -k_1 [x][y]\\ \frac{\mathrm d[y]}{\mathrm dt} &= k_1 [x][y] - k_3[y] \end{align*}

My solution seems to drop below zero for some reason. This cannot be the case since x, y are concentrations. Is this due to stiffness perhaps? Can I do better?

Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];
k1 = 1*10^13;
k2 = 1*10^6;
k3 = 2000;
s = NDSolve[{x'[t] == -k1*x[t]*y[t], y'[t] == -k3*y[t] + k1*x[t]*y[t],
x == 1*10^(-8), y == 1*10^(-14)}, {x, y}, {t, 0, 0.001}];
Plot[Evaluate[{x[t], y[t]} /. s], {t, 0, 0.0003},
PlotStyle -> {{AbsoluteThickness,
RGBColor[0, 0, 0]}, {AbsoluteThickness,
RGBColor[.7, 0, 0]}, {AbsoluteThickness, RGBColor[0, .7, 0]}},
Axes -> False, Frame -> True,
PlotLabel ->
StyleForm[A StyleForm[" B*", FontColor -> RGBColor[.7, 0, 0]] ,
FontFamily -> "Helvetica", FontWeight -> "Bold"] ,
PlotRange -> {{0, 0.00030}, {1*10^-16, 1.5*10^-8}}] • I think there is nothing wrong, your initial conditions are very close to $(0,0)$ which is one of the equilibrium points that have the y-axis as a (local) stable manifold. Actually the x-axis looks like an attractor for the solutions. – Spawn1701D Apr 21 '13 at 8:42
• Not the range your choice of initial conditions are the reason why you get this plot. For a moment imagine for some $t$ $y$ becomes $0$, what happens to your system? It becomes $\hat x=0,\ \hat y=0$ and as consequence the solution from that point on is stuck to the $x$-axis. Thats what happens to your case the solutions near $y=0$ are attracted to the $x$-axis and eventually fall to it. – Spawn1701D Apr 21 '13 at 17:49
• I think it goes below zero since the initial condition is so low. This is due to stiffness right? – l3win Apr 21 '13 at 23:37

If you increase the working precision this will work:

Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];
k1 = 1*10^13;
k2 = 1*10^6;
k3 = 2000;
s = NDSolve[{x'[t] == -k1*x[t]*y[t], y'[t] == -k3*y[t] + k1*x[t]*y[t],
x == 1*10^(-8), y == 1*10^(-14)}, {x, y}, {t, 0, 1/1000},
WorkingPrecision -> 20];
Plot[Evaluate[{x[t], y[t]} /. s], {t, 0, 0.0003},
PlotStyle -> {{AbsoluteThickness,
RGBColor[0, 0, 0]}, {AbsoluteThickness,
RGBColor[.7, 0, 0]}, {AbsoluteThickness, RGBColor[0, .7, 0]}},
Axes -> False, Frame -> True,
PlotLabel ->
StyleForm[A StyleForm[" B*", FontColor -> RGBColor[.7, 0, 0]],
FontFamily -> "Helvetica", FontWeight -> "Bold"], PlotRange -> All] • I wonder if an appropriate Method setting might be much more profitable... – J. M. is away Apr 22 '13 at 12:28
• I have seen the same problem as the OP before, and increasing WorkingPrecision is usually the answer! I've never tried playing around with other Method though – Sosi Apr 22 '13 at 12:34
• I meant that in this case, switching to arbitrary precision did the trick, but maybe one could choose to tweak method options instead of increasing the internal precision in computations. Like, maybe using "StiffnessSwitching"... – J. M. is away Apr 22 '13 at 12:39
• @l3win, and how long did you wait? – user21 Apr 22 '13 at 16:26
• Would it help if you blindly switched to the LSODA method? I generally noticed that the LSODA is quite forgiving and pushes stiff equations towards a "solution" better than the BDF at the very least. – dearN Apr 22 '13 at 17:32

I do not think the problem is stiffness; nor is it, arguably, scaling. The values of the solution are on the order of 10^-14 to 10^-8. The default setting for AccuracyGoal is around 8, which means that steps with an error of 10^-8 or less are acceptable. And they shouldn't be accepted with solutions that are so small.

[If you are unsure how precision and accuracy goals work in Mathematica, please see the appendix below first.] A setting AccuracyGoal -> ag needs to be high enough so that

$$\|{\bf X(t)}\| <\mskip-4mu< 10^{-ag}$$

I would suggest setting AccuracyGoal to be about 8 digits above -Log10[size], where size is the smallest value that should be considered significantly different from zero. (One might set it to Infinity, even.)

The above rule of thumb suggests

AccuracyGoal -> 22


based on the initial condition y == 1*10^(-14). In fact setting it to 16 (based on the norm of the initial conditions being approximately $10^{-8}$) or Infinity both work, too.

s = NDSolve[{x'[t] == -k1*x[t]*y[t], y'[t] == -k3*y[t] + k1*x[t]*y[t],
x == 1*10^(-8), y == 1*10^(-14)}, {x, y}, {t, 0, 0.001},
AccuracyGoal -> 22
]; Why do scaling and increasing WorkingPrecision work?

In @Goofy's answer, the rescaled x[t] and y[t] have values on the order of 10 to nearly 100000. So the default AccuracyGoal of around 8 is sufficiently high. One might say that the problem is poorly scaled with respect to the accuracy goal.

Using WorkingPrecision -> 20 seems to be barely sufficient. It raises AccuracyGoal to 10. It raises PrecisionGoal to 10, too, but with a solution whose norm is around 10^-8 or less, such a PrecisionGoal should be negligible compared to the effect of AccuracyGoal. One can check that through most of integration interval, WorkingPrecision -> 20 produces a solution that has 1 to 3 digits of precision, which is in line with an error of 10^-10 in a solution of norm 10^-8. (This was estimated by comparing it with a solution computed with WorkingPrecision -> 40 and assuming this solution was highly accurate.)

FWIW, the AccuracyGoal -> 22 solution above has 6-8 digits of precision throughout most of the integration interval. The worst precision in all methods (except the OP's) is around where the x and y components cross near t == 0.00014.

Appendix: Error norms in Mathematica.

The above discussion assumes the reader is familiar with how error is compared with PrecisionGoal -> pg and AccuracyGoal -> ag. NDSolve uses NDSolveScaledVectorNorm, which is discussed in Norms in NDSolve. The approach is also discussed in the documentation for PrecisionGoal and AccuracyGoal. I give a brief explanation.

Let $\bf E$ be the estimated error in the solution ${\bf X}$. If the error $\bf E$ is less than $10^{-ag} + 10^{-pg}\,\|{\bf X}\|$, then the result is accepted; if not, then in NDSolve the step size is adjusted or a warning or error is generated if the step size cannot be adjusted.

When $10^{-ag} >\mskip-4mu> 10^{-pg}\,\|{\bf X}\|$, $$10^{-ag} + 10^{-pg}\,\|{\bf X}\| \approx 10^{-ag}$$ implies that if the error indicates that ${\bf X}$ is correct to at least $ag$ digits past the decimal point, the step is accepted. When $10^{-pg}\,\|{\bf X}\| >\mskip-4mu> 10^{-ag}$, $$10^{-ag} + 10^{-pg}\,\|{\bf X}\| \approx 10^{-pg}\,\|{\bf X}\|$$ implies that if the error indicates that ${\bf X}$ is correct to at least the leading $pg$ digits, the step is accepted. Thus when $$10^{-ag} >\mskip-4mu> 10^{-pg}\,\|{\bf X}\|$$ which usually happens when $\bf X$ is small, the errors that are nearly equal to $10^{-ag}$ are acceptable. Hence, solution values of this magnitude or less are being treated as roughly equivalent to zero. This is the issue in the OP.

• Interesting. To make your conclusion more convictive: WorkingPrecision -> 20, PrecisionGoal -> 8, AccuracyGoal -> 8 results in the erroneous graph. – xzczd Feb 26 '17 at 5:18
• @xzczd Yes, that's a good example. With PrecisionGoal and AccuracyGoal set to $MachinePrecision/2, you get a nearly indistinguishable solution (relative difference approx. 100 *$MachineEpsilon). It shows the working precision is a relatively negligible factor. – Michael E2 Feb 26 '17 at 14:05

Here's an approach to scaling the problem, first suggested by @BurgerAutomata, which does not require the units. Just replace $x$ and $y$ by $\alpha x$ and $\beta y$ and solve for $\alpha$ and $\beta$ so that the coefficients are not too different in magnitude. In fact one can get the coefficients to be $1$, pretty much by inspection.

k1 = 1*10^13;
k2 = 1*10^6;
k3 = 2000;
α = k3/k1; β = 1/k1;
s = NDSolve[{α*x'[t] == -α*β*k1*x[t]*y[t],
β*y'[t] == -k3*β*y[t] + α*β*k1*x[t]*y[t],
x == 1*10^(-8)/α, y == 1*10^(-14)/β}, {x, y}, {t, 0, 0.001}];
Plot[Evaluate[{α*x[t], β*y[t]} /. s], {t, 0, 0.0003},
PlotStyle -> {{AbsoluteThickness,
RGBColor[0, 0, 0]}, {AbsoluteThickness,
RGBColor[.7, 0, 0]}, {AbsoluteThickness, RGBColor[0, .7, 0]}},
Axes -> False, Frame -> True,
PlotLabel ->
StyleForm[A StyleForm[" B*", FontColor -> RGBColor[.7, 0, 0]],
FontFamily -> "Helvetica", FontWeight -> "Bold"],
PlotRange -> {{0, 0.00030}, {1*10^-16, 1.5*10^-8}}] • Oh finally here comes an example! If I understand it correctly, "scaling" generally doesn't require one to know the unit? – xzczd Feb 24 '17 at 3:49
• @xzczd Scaling in a mathematical-numerical sense has to do with rescaling the variables/coefficients so that calculations can be done in the floating-point format with greater precision. From a mathematical point of view, the numbers matter and not the units. From a scientific point of view, I suppose, we don't know what $x$ and $y$ represent, but the original quantities are represented by the scaled $\alpha x$ and $\beta y$. (If scaling were done internally by a numerical algorithm, why should it care about the units the user has in mind?) – Goofy Feb 24 '17 at 11:46

As prompted by one of the comments, stiffness switching did help some. I usually go with Method->LSODA for stiff equations but tried Method -> {"StiffnessSwitching", Method -> {"ExplicitRungeKutta", Automatic}} straight out of the help manual. A little late, this answer is!

Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];
k1 = 1*10^13;
k2 = 1*10^6;
k3 = 2000;
s = NDSolve[{x'[t] == -k1*x[t]*y[t], y'[t] == -k3*y[t] + k1*x[t]*y[t],
x == 1*10^(-8), y == 1*10^(-14)}, {x, y}, {t, 0, 0.001},
Method -> {"StiffnessSwitching",
Method -> {"ExplicitRungeKutta", Automatic}}];
Plot[Evaluate[{x[t], y[t]} /. s], {t, 0, 0.0003},
PlotStyle -> {{AbsoluteThickness,
RGBColor[0, 0, 0]}, {AbsoluteThickness,
RGBColor[.7, 0, 0]}, {AbsoluteThickness, RGBColor[0, .7, 0]}},
Axes -> False, Frame -> True,
PlotLabel ->
StyleForm[A StyleForm[" B*", FontColor -> RGBColor[.7, 0, 0]],
FontFamily -> "Helvetica", FontWeight -> "Bold"],
PlotRange -> {{0, 0.00030}, {1*10^-16, 1.5*10^-8}}] • I'm sorry, but this result is incorrect, just Show` it with the one obtained by e.g. user21. – xzczd Feb 24 '17 at 3:49

It's usually much better in my experience to scale the variables. You get no underflow this way, and stiffness is partially ameliorated, although if a system is stiff it will probably stay that way. As an example, try to be as a-dimensional as possible. Take some main unit/molecule/entity/whatever magnitude as your reference scale, and then work everything around it so that the variables have no units. Usually you will obtain stuff that is by the majority much closer to 1. In the field Molecular Dynamics there are plenty of references on scaling this way.

• Do you mind fixing OP's code in this way? In fact it's not the first time I see this method suggested, but I never find a specific example in this site. – xzczd Mar 5 '14 at 2:03
• I'd love to try, but I need the units used in the problem. – BurgerAutomata Mar 5 '14 at 10:48