6
$\begingroup$

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["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
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}}]

ode

$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$
    – Spawn1701D
    Apr 21, 2013 at 8:42
  • $\begingroup$ 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. $\endgroup$
    – Spawn1701D
    Apr 21, 2013 at 17:49
  • $\begingroup$ I think it goes below zero since the initial condition is so low. This is due to stiffness right? $\endgroup$
    – l3win
    Apr 21, 2013 at 23:37

5 Answers 5

4
$\begingroup$

If you increase the working precision this will work:

Needs["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
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]

enter image description here

$\endgroup$
9
  • $\begingroup$ I wonder if an appropriate Method setting might be much more profitable... $\endgroup$ Apr 22, 2013 at 12:28
  • 1
    $\begingroup$ 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 $\endgroup$
    – Sos
    Apr 22, 2013 at 12:34
  • 1
    $\begingroup$ 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"... $\endgroup$ Apr 22, 2013 at 12:39
  • 1
    $\begingroup$ @l3win, and how long did you wait? $\endgroup$
    – user21
    Apr 22, 2013 at 16:26
  • 1
    $\begingroup$ 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. $\endgroup$
    – dearN
    Apr 22, 2013 at 17:32
8
$\begingroup$

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[0] == 1*10^(-14). In fact setting it to 16 (based on the norm of the initial conditions being approximately $\scriptstyle 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[0] == 1*10^(-8), y[0] == 1*10^(-14)}, {x, y}, {t, 0, 0.001},
   AccuracyGoal -> 22
   ];

Mathematica graphics

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 NDSolve`ScaledVectorNorm, 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.

$\endgroup$
2
  • 2
    $\begingroup$ Interesting. To make your conclusion more convictive: WorkingPrecision -> 20, PrecisionGoal -> 8, AccuracyGoal -> 8 results in the erroneous graph. $\endgroup$
    – xzczd
    Feb 26, 2017 at 5:18
  • $\begingroup$ @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. $\endgroup$
    – Michael E2
    Feb 26, 2017 at 14:05
3
$\begingroup$

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[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}}]

enter image description here

$\endgroup$
2
  • $\begingroup$ Oh finally here comes an example! If I understand it correctly, "scaling" generally doesn't require one to know the unit? $\endgroup$
    – xzczd
    Feb 24, 2017 at 3:49
  • 2
    $\begingroup$ @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?) $\endgroup$
    – Goofy
    Feb 24, 2017 at 11:46
1
$\begingroup$

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["DifferentialEquations`NDSolveProblems`"];
Needs["DifferentialEquations`NDSolveUtilities`"];
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}, 
   Method -> {"StiffnessSwitching", 
     Method -> {"ExplicitRungeKutta", Automatic}}];
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}}]

enter image description here

$\endgroup$
1
  • $\begingroup$ I'm sorry, but this result is incorrect, just Show it with the one obtained by e.g. user21. $\endgroup$
    – xzczd
    Feb 24, 2017 at 3:49
0
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$
    – xzczd
    Mar 5, 2014 at 2:03
  • $\begingroup$ I'd love to try, but I need the units used in the problem. $\endgroup$ Mar 5, 2014 at 10:48

Your Answer

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

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