chemistry – NDSolve: InterpolatingFunction of Chemical Kinetics System Differs From Expected (AKA – Why Are My Plots Squiggly?)

I’m trying to use Mathematica (12.0.0) to model radiation chemical kinetics using a known reaction set, in this case it’s water radiolysis.

It solves the set of equations I throw at it without any errors, however the InterpolatingFunction it’s throwing out for the result is different from expected. The magnitudes in some cases is a bit off and most notably for a number of the products I’ve specified they are very irregular, being a bit “squiggly” where they should be smooth curves and straight lines.

I’ve put the same reactions/rate constants into other software (FACSIMILE) to show what outputs I should expect. One example below:

H (Mathematica)

Plot of the concentration of the hydrogen radical over time as generated by NDSolve in Mathematica 12.0.0!

H (Facsimile)

Plot of the concentration of the hydrogen radical over time as generated by FACSIMILE and plotted in Excel!

As you can see there are 5 peaks over the plot in the Mathematica plot whereas the FACSIMILE plot is smooth. My question is what is causing this behaivour and how can I resolve it?

Code below:

ClearAll("Global'*")

doserate = 10;
time = 60;
geaqa = 2.6*1.0364*10^-7;
gh = 0.6*1.0364*10^-7;
goh = 2.7*1.0364*10^-7;
gh2 = 0.45*1.0364*10^-7;
gh2o2 = 0.7*1.0364*10^-7;

kgeaqa = geaqa*doserate;
kgh = gh*doserate;
kgoh = goh*doserate;
kgh2 = gh2*doserate;
kgh2o2 = gh2o2*doserate;

kw2 = 7.26*10^9;
kw3 = 5*10^9;
kw4 = 4.8*10^9;
kw5 = 3.4*10^10;
kw6 = 3.5*10^10;
kw7 = 1.1*10^10;
kw8 = 1.4*10^10;
kw9 = 2.3*10^10;
kw10 = 1.3*10^10;
kw11 = 1.3*10^10;
kw12 = 3.6*10^7;
kw13 = 1.3*10^10;
kw14 = 1.13*10^10;
kw14a = 1.14*10^10;
kw15 = 1.13*10^10;
kw16 = 2.9*10^7;
kw17 = 1.1*10^10;
kw18 = 8.8*10^9;
kw19 = 8.4*10^5;
kw20 = 1*10^8;
kw21 = 3*10^-1;
kw22 = 5*10^5;
kw22a = 2.3*10^-7;
kw23 = 1.17*10^-3;
kw23b = 1.18*10^11;
kw24 = 8.9*10^-2;
kw24b = 4.78*10^10;
kw25 = 1.27*10^10;
kw25b = 1.4*10^6;
kw26 = 8.9*10^-2;
kw26b = 4.78*10^10;
kw27 = 1.27*10^10;
kw27b = 1.4*10^6;
kw28 = 4.78*10^10;
kw28b = 7.35*10^5;
kw29 = 1.27*10^10;
kw29b = 1.63*10^-1;
kw30 = 5.83;
kw30b = 2.1*10^10;
kw31 = 2.44*10^7;
kw31b = 1.74*10^1;
kw32 = 4.58*10^-5;
kw32b = 3.95*10^7;


rw2 = kw2*eaqa(t)^2*h2o(t)^2;
rw3 = kw3*h(t)^2;
rw4 = kw4*oh(t)^2;
rw5 = kw5*eaqa(t)*h(t)*h2o(t);
rw6 = kw6*eaqa(t)*oh(t);
rw7 = kw7*h(t)*oh(t);
rw8 = kw8*eaqa(t)*h2o2(t);
rw9 = kw9*eaqa(t)*o2(t);
rw10 = kw10*eaqa(t)*o2a(t)*h2o(t);
rw11 = kw11*eaqa(t)*ho2(t);
rw12 = kw12*h(t)*h2o2(t);
rw13 = kw13*h(t)*o2(t);
rw14 = kw14*h(t)*ho2(t);
rw14a = kw14a*h(t)*ho2(t);
rw15 = kw15*h(t)*o2a(t);
rw16 = kw16*oh(t)*h2o2(t);
rw17 = kw17*oh(t)*o2a(t);
rw18 = kw18*oh(t)*ho2(t);
rw19 = kw19*ho2(t)^2;
rw20 = kw20*o2a(t)*ho2(t)*h2o(t);
rw21 = kw21*o2a(t)^2*h2o(t)^2;
rw22 = kw22*h2o2(t);
rw22a = kw22a*h2o2(t);
rw23 = kw23*h2o(t);
rw23b = kw23b*hc(t)*oha(t);
rw24 = kw24*h2o2(t);
rw24b = kw24b*ho2a(t)*hc(t);
rw25 = kw25*h2o2(t)*oha(t);
rw25b = kw25b*ho2a(t)*h2o(t);
rw26 = kw26*oh(t);
rw26b = kw26b*hc(t)*oa(t);
rw27 = kw27*oh(t)*oha(t);
rw27b = kw27b*oa(t)*h2o(t);
rw28 = kw28*ho2(t);
rw28b = kw28b*hc(t)*o2a(t);
rw29 = kw29*ho2(t)*oha(t);
rw29b = kw29b*o2a(t)*h2o(t);
rw30 = kw30*h(t);
rw30b = kw30b*eaqa(t)*hc(t);
rw31 = kw31*h(t)*oha(t);
rw31b = kw31b*eaqa(t)*h2o(t);
rw32 = kw32*h(t)*h2o(t);
rw32b = kw32b*h2(t)*oh(t);

watersolver = NDSolve({eaqa'(t) == kgeaqa + rw30 + rw31 - rw2 - rw5 - rw6 - rw8 - rw9 - rw10 - rw11 - rw30b - rw31b,
h2o'(t) == rw7 + rw12 + rw16 + rw18 + rw22 + rw23b + rw25 + rw27 + rw29 + rw31 + rw32b - rw2 - rw5 - rw10 - rw20 - rw21 - rw23 - rw25b - rw27b - rw29b - rw31b - rw32,
h2'(t) == kgh2 + rw2 + rw3 + rw5 + rw32 - rw32b,
oha'(t) == rw2 + rw5 + rw6 + rw8 + rw10 + rw17 + rw20 + rw21 + rw23 + rw25b + rw27b + rw29b + rw31b - rw23b - rw25 - rw27 - rw29 - rw31 ,
oh'(t) == kgoh + rw8 + rw12 + rw14a + rw22a + rw26b + rw27b + rw32 - rw4 - rw6 - rw7 - rw16 - rw17 - rw18 - rw26 - rw27 - rw32b,
h2o2'(t) == kgh2o2 + rw4 + rw10 + rw14 + rw19 + rw20 + rw21 + rw24b + rw25b - rw8 - rw12 - rw16 - rw22 - rw22a - rw24 - rw25,
h'(t) == kgh + rw30b + rw31b + rw32b - rw3 - rw5 - rw7 - rw12 - rw13 - rw14 - rw14a - rw15 - rw30 - rw31 - rw32,
o2a'(t) == rw9 + rw28 + rw29 - rw10 - rw15 - rw17 - rw20 - rw21 - rw28b - rw29b,
ho2a'(t) == rw11 + rw15 + rw24 + rw25 - rw24b - rw25b,
hc'(t) == rw23 + rw24 + rw26 + rw28 + rw30 - rw23b - rw24b - rw26b - rw28b - rw30b,
oa'(t) == rw26 + rw27 - rw26b - rw27b,
o2'(t) == rw17 + rw18 + rw19 + rw20 + rw21 + rw22 - rw9 - rw13,
ho2'(t) == rw13 + rw16 + rw28b + rw29b - rw11 - rw14 - rw14a - rw18 - rw19 - rw20 - rw28 - rw29,
eaqa(0) == 0,
h2o(0) == 55.5,
h2(0) == 0,
oha(0) == 0,
oh(0) == 0,
h2o2(0) == 0,
h(0) == 0,
o2a(0) == 0,
ho2a(0) == 0,
hc(0) == 0,
oa(0) == 0,
o2(0) == 0,
ho2(0) == 0}, {eaqa, h2o, h2, oha, oh, h2o2, h, o2a, ho2a, hc, oa, o2, ho2}, {t, 0, 1})

I’m fairly new to Mathematica and I’m still trying to learn all its intricacies/limits. Have I missed anything glaring? Is there a better InterpolatingFunction/NDSolve setting for what I’m trying to do? Am I running into some sort of limitation in Mathematica already? Any help on this would be much appreciated.

Many thanks.