Hello,
I am attempting to numerically integrate a set of kinetic equations to describe the reaction for aspirin synthesis. Unfortunately, when I run the below code, the integration appears to take place however the graph is empty with only the axes displayed. As I am new to Mathematica, I was hoping someone with a bit of experience may be kind enough to point out where I'm going wrong. Just for info, I am running the whole piece of code at once.
Many thanks,
Robin
ClearAll[CSA, CAA, CASA, CHA, CASAA, CH2O, V, t]
Reaction = NDSolve[{
CAA'[t] == -(k1 * CSA[t] * CAA[t]) - (k2 * CASA[t] * CAA[t]) - 0.5 *(k4 * CAA[t] * CH2O[t]) - CAA[t] ((V'[t])/(V[t])),
CHA'[t] == (k1 * CSA[t] * CAA[t]) + (k2 * CASA[t] * CAA[t]) + (k3 * CASAA[t] * CH2O[t]) + (k4 * CAA[t] * CH2O[t]) - CHA[t] (V'[t]/V[t]),
CSA'[t] == (kd *(CSatSA - CSA[t])^d) - (k1 * CSA[t] * CAA[t]) - CSA[t] ((V'[t])/(V[t])),
CASAA'[t] == (k2 * CASA[t] * CAA[t]) - (k3 * CASAA[t] * CH2O[t]) - CASAA[t] ((V'[t])/(V[t])),
CH2O'[t] == - (k3 * CASAA[t] * CH2O[t]) - 0.5 (k4 * CAA[t] * CH2O[t]) - CH2O[t] ((V'[t])/(V[t])) + CInH2O (f/V[t]),
CASA'[t] == (k1 * CSA[t] * CAA[t]) - (k2 * CASA[t] * CAA[t]) + (k3 * CASAA[t] * CH2O[t])
- (kc *(CASA[t] - CSatASA)^c) - CASA[t]((V'[t])/(V[t])),
V'[t] == V[t]*(vASA ((k1 * CSA[t] * CAA[t]) - (k2 * CASA[t] * CAA[t]) + (k3 * CASAA[t] * CH2O[t]) - (kc *(CASA[t] - CSatASA)^c)) +
vH2O (-(k3 * CASAA[t] * CH2O[t]) - 0.5 (k4 * CAA[t] * CH2O[t]) + CInH2O (f/V[t]))
+ vASAA ((k2 * CASA[t] * CAA[t]) - (k3 * CASAA[t] * CH2O[t])) + vSA ((kd *(CSatSA - CSA[t])^d) - (k1 * CSA[t] * CAA[t])) +
vHA ((k1 * CSA[t] * CAA[t]) + (k2 * CASA[t] * CAA[t]) + (k3 * CASAA[t] * CH2O[t]) + (k4 * CAA[t] * CH2O[t]))
+ vAA (-(k1 * CSA[t] * CAA[t]) - (k2 * CASA[t] * CAA[t]) - 0.5 *(k4 * CAA[t] * CH2O[t]))),
CAA[0] == 10.57,
CHA[0] == 0.0,
CSA [0] == 0.19,
CASAA [0] == 0.0,
CH2O [0] == 0.0,
CASA [0] == 0.0,
V [0] == 120}
/. {k1 -> 0.0337, k2 -> 0.43, k3 -> 1525, k4 -> 85, kd -> 4.405, kc -> 1.09, c -> 1.34, d -> 1.90, CSatASA -> 1.52,
CSatSA -> 2.19, vASA -> 0.014, vSA -> 0.014, vASAA -> 0.014, vHA -> 0.014, vAA -> 0.014, vH2O -> 0.014, CInH20 -> 55.5, f -> 0},
{CSA, CASA, CASAA, CHA, CH2O, CAA, V},
{t, 0, 200}]
myplot = Plot[
{Evaluate[CSA[t] /. Reaction], Evaluate[CASA[t] /. Reaction],
Evaluate[CASAA[t] /. Reaction], Evaluate[CHA[t] /. Reaction],
Evaluate[CH2O[t] /. Reaction], Evaluate[CAA[t] /. Reaction],
Evaluate[V[t] /. Reaction]},
{t, 0, 200},
PlotPoints -> 500,
PlotRange -> {{0, 200}, {0, 0.2}},
AxesLabel -> {time, conc},
LabelStyle -> Directive[Black, Medium, "Palatino"],
PlotStyle -> {{Red, Thickness[0.008]}, {Green,
Thickness[0.008]}, {Blue, Thickness[0.008]}, {Black,
Thickness[0.008]}, {Yellow, Thickness[0.008]}, {Orange,
Thickness[0.008]}}]