Unexpected results from NDSolve

I am trying to solve a stiff reaction diffusion system with NDSolve. However, it does not produce the expected results.

My problem is a spherical cell with 5 different species, of which only one can flow across the cell membrane. Initially the cell contains nothing of species 4 and 5 (only #5 can cross the membrane), and the idea is that the external concentration of #5 should increase instantaneously from 0 to outConc.

However, when having this step-like increase the problem is not solved correctly by NDSolve, so a fix using a delay factor of (1-Exp[t^2]) for the increase of the outer concentration. But this does not help all the way. Though the system has been solved now, and it looks rather good, there appear some negative values of the concentrations at small t:s.

In an attempt to locate these, the "EventLocator" method was used, but this in turn gave an error about complex values. Which should not be there. Any ideas how to get a stable solution from the stiff problem?

NDSolve::evre: The value of the event function at t = 7.450580596923828`*^-6 was not a real number. The event will be ignored in steps where it does not evaluate to real numbers at both ends.

My code so far:

permcoeff = 9.5*10^(-3); (* Membrane permeability of ButH [dm/s] *)
Btot = 2.84106*10^(-2); (* Total concentration of cell buffer [mol/L*)

Kbeq =  44.0945*10^(-9); (* Equilibrium constant for cell buffer [mol/L]*)
Kbuteq = 10^(-4.82); (* Eq constant for ButH [mol/L] *)
pH0 = 7.2; (* Initial pH-value *)
outConc = 20*10^(-3) ; (* Conc of ButH outside cell [mol/L] *)
radius = 5*10^(-3); (* [dm] *)
tmax = 50; (* [s] *)
delta = 10^(-10); (* Min value for radius [dm]*)
epsilon = 5*10^(-2); (* Distance from membrane of electrode [dm] *)

(* Coefficients of diffusion [dm^2/s] *)
dc = {
10^(-5),
6.5*10^(-8),
6.5*10^(-8),
2*9.2*10^(-8),
2*9.2*10^(-8)
};

(* Reaction rates
k3f - [H][B] - > [HB]
k3b - [HB] - > [H]+[B]
k4f - [H][But] - > [ButH]
k4b - [ButH] - > [H]+[But]
*)
k3f = 10^(10);
k3b = k3f*Kbeq;
k4f = k3f;
k4b = k4f * Kbuteq;

(* Concentrations
c1 - [H+]
c2 - [B-]
c3 - [HB]
c4 - [But-]
c5 - [ButH]
*)

(* Flux at outer boundary. Only ButH enters the cell *)
influx = {
0,
0,
0,
0,
-(1 - Exp[-t]) (c5[radius, t] - outConc)*permcoeff/dc[]
};

(* Initial concentrations [mol/L] *)
cinit = {
10^(-pH0),
Kbeq*Btot/(10^(-pH0) + Kbeq),
10^(-pH0)*Btot/(10^(-pH0) + Kbeq),
0,
0
};

(* BC at membrane *)
bcMembrane = {
Derivative[1, 0][c1][radius, t] == influx[],
Derivative[1, 0][c2][radius, t] == influx[],
Derivative[1, 0][c3][radius, t] == influx[],
Derivative[1, 0][c4][radius, t] == influx[],
Derivative[1, 0][c5][radius, t] == influx[]
};

(* BC at center *)
bcCenter = {
Derivative[1, 0][c1][delta, t] == 0,
Derivative[1, 0][c2][delta, t] == 0,
Derivative[1, 0][c3][delta, t] == 0,
Derivative[1, 0][c4][delta, t] == 0,
Derivative[1, 0][c5][delta, t] == 0
};

(* Initial conditions *)
initcond = {
c1[r, 0] == cinit[],
c2[r, 0] == cinit[],
c3[r, 0] == cinit[],
c4[r, 0] == cinit[],
c5[r, 0] == cinit[]
};

(* RDEs to be solved *)
diffeqn = {
dc[]*(2*D[c1[r, t], r]/r + D[c1[r, t], r, r]) -
k3f*c1[r, t]*c2[r, t] + k3b*c3[r, t] - k4f*c1[r, t]*c4[r, t] +
k4b*c5[r, t] == D[c1[r, t], t],
dc[]*(2*D[c2[r, t], r]/r + D[c2[r, t], r, r]) -
k3f*c1[r, t]*c2[r, t] + k3b*c3[r, t] == D[c2[r, t], t],
dc[]*(2*D[c3[r, t], r]/r + D[c3[r, t], r, r]) +
k3f*c1[r, t]*c2[r, t] - k3b*c3[r, t] == D[c3[r, t], t],
dc[]*(2*D[c4[r, t], r]/r + D[c4[r, t], r, r]) -
k4f*c1[r, t]*c4[r, t] + k4b*c5[r, t] == D[c4[r, t], t],
dc[]*(2*D[c5[r, t], r]/r + D[c5[r, t], r, r]) +
k4f*c1[r, t]*c4[r, t] - k4b*c5[r, t] == D[c5[r, t], t]
};

sol = NDSolve[
Join[diffeqn, bcMembrane, bcCenter, initcond], {c1, c2, c3, c4,
c5}, {r, delta, radius}, {t, 0, tmax},
Method -> {"EventLocator", "Event" -> {c2[r, t], c2[r, t]},
"Direction" -> {1, -1}},  PrecisionGoal -> 100,
MaxStepSize -> {0.001, 0.01}, MaxSteps -> 10^7];

Plot3D[Evaluate[c1[r, t] /. sol], {t, 0, tmax}, {r, delta, radius},
PlotPoints -> 50, PlotLabel -> "[H+]",
AxesLabel -> {"Time [s]", "Distance from center [dm]", "[mol/L]"}]
• a) in your "Event" should it be c2[r,t] instead of c2[t,r]?b) When I add StepMonitor :> Print[c2[r,t]] it appears that r is not getting set, leading to a non-real number (rather than complex) Neither of these comments help with your problem, though. What was the original equation? Maybe that can be debugged and all this EventLocator stuff can be avoided. :-) – Eric Brown Oct 8 '12 at 13:33
• Thank you, I had totally messed up the order of r,t. However as you say the problem still remains after changing them to the right positions. – Daniel Oct 15 '12 at 7:52

My understanding of what the problem is that there is an abberation at a certain {r,t}, for example in c2, that has negative values for concentration. When I used a finer grid, then that went away:

sol = NDSolve[
Join[diffeqn, bcMembrane, bcCenter, initcond],
{c1, c2, c3, c4, c5}, {r, delta, radius}, {t, 0, tmax},
Method -> {"MethodOfLines", "Method" -> "BDF",
"SpatialDiscretization" -> {"TensorProductGrid", "MinPoints" -> 500}}];  • I found about 10-15% savings in run time with LSODA instead of BDF. Interesting, innit? :P – dearN Oct 8 '12 at 18:47
• LSODA rocks. It has some nice stiff/nonstiff switching capabilities that seem to "just work" for a lot of problems. It goes on land and sea. – Eric Brown Oct 8 '12 at 19:02
• This person has a nifty table comparing several (non)stiff solvers. – dearN Oct 8 '12 at 21:50
• Thank you Eric, I had previously tried to experiment with the maxStepSize for both the radius and the time, however, decreasing these enough to achieve a good solution requires too much computation time. With your solution everything seems to work just fine! – Daniel Oct 15 '12 at 8:06