Unexpected results from NDSolve

Question

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

(* 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[[1]],
   Derivative[1, 0][c2][radius, t] == influx[[2]],
   Derivative[1, 0][c3][radius, t] == influx[[3]],
   Derivative[1, 0][c4][radius, t] == influx[[4]],
   Derivative[1, 0][c5][radius, t] == influx[[5]]
   };

(* 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[[1]],
   c2[r, 0] == cinit[[2]],
   c3[r, 0] == cinit[[3]],
   c4[r, 0] == cinit[[4]],
   c5[r, 0] == cinit[[5]]
   };

(* RDEs to be solved *) 
diffeqn = {
   dc[[1]]*(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]]*(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[[3]]*(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[[4]]*(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[[5]]*(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]"}]

Answer 1

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