Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
1  
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

2 Answers 2

up vote 1 down vote accepted

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

share|improve this answer
    
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.

share|improve this answer
    
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 at 2:03
    
I'd love to try, but I need the units used in the problem. –  BurgerAutomata Mar 5 at 10:48

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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