# 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[0] == 1*10^(-8), y[0] == 1*10^(-14)}, {x, y}, {t, 0, 0.001}];
Plot[Evaluate[{x[t], y[t]} /. s], {t, 0, 0.0003},
PlotStyle -> {{AbsoluteThickness[2],
RGBColor[0, 0, 0]}, {AbsoluteThickness[2],
RGBColor[.7, 0, 0]}, {AbsoluteThickness[2], 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[0] == 1*10^(-8), y[0] == 1*10^(-14)}, {x, y}, {t, 0, 1/1000},
WorkingPrecision -> 20];
Plot[Evaluate[{x[t], y[t]} /. s], {t, 0, 0.0003},
PlotStyle -> {{AbsoluteThickness[2],
RGBColor[0, 0, 0]}, {AbsoluteThickness[2],
RGBColor[.7, 0, 0]}, {AbsoluteThickness[2], 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. 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
@J.M., sorry I had the wrong code in the past buffer. Now the option is there. I am not quite sure I can follow you with the Method option. WorkingPrecision is something different than a method option. Could explain a bit what you mean. – user21 Apr 22 '13 at 12:38
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. Apr 22 '13 at 12:39
@J.M., if look at the results from other methods like "Adams", "BDF", they go completely south. This looks like a badly conditioned stiff systems. I think arbitrary precision is the way to go. But I'd most interested in other routes. – user21 Apr 22 '13 at 12: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