# Problem with Numerical Integration of Aspirin Synthesis Reaction

Posted 5 years ago
7553 Views
|
12 Replies
|
2 Total Likes
|
 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]}}] 
12 Replies
Sort By:
Posted 5 years ago
 I don't know the details of your differential equation, but it appears that, for some reason, the results are not strictly real, for example:I n[87]:= ({CSA[t], CASA[t], CASAA[t], CHA[t], CH2O[t], CAA[t], V[t]} /. Reaction[[1]]) /. t -> 10 Out[87]= {2.14012 + 0.0397231 I, 1.69264 + 0.23278 I, 7.36331 + 0.656905 I, 8.99929 - 0.0635311 I, 0. + 0. I, -0.0623169 - 0.0743696 I, 141.893 + 2.18947 I} 
Posted 5 years ago
 Hi David, thank you for your response.That is very unusual, the reactions are supposed to be of the form dCSA/dt = r1 -r2 + r3 where dCSA is the rate of change of the concentration of SA and r1 is a rate expression of the form k1Conc(species x)Conc(species y). I'm not entirely sure why that isn't creating a real solution.Thanks again for your help
Posted 5 years ago
 Hi Robin and David,In the attached notebook the calculation indicates that in the solution for the derivatives at t=0, given the initial conditions, the derivatives of both V and CASA have complex values. Is there perhaps an error in the equations?Kind regards, David Attachments:
Posted 5 years ago
 Dear David and David,I must thank you for such a comprehensive answer, I'm absolutely delighted! I was completely unsure of what the problem was however now you've pointed it out I can go away and double check the equations. I did derive them so I'll have to have a good run through them again.Another great benefit is that you've given me a really simple way of solving the equations, the syntax I was using before was very clunky (as I'm sure you've noticed!) so I can use this much simpler version from now on.Kindest regards and thanks for your help, especially on Christmas Eve!Robin
Posted 5 years ago
 You're very welcome, Robin. Here is a second notebook. Rather that calculating the first step, I just obtained the derivatives by solution for them, and graphed them by differentiating the solution for the concentrations. Several derivatives start with imaginary components.Merry Christmas to you and David! Attachments:
Posted 5 years ago
 Dear David, Thank you very much for the extra file and I hope you have enjoyed your Christmas! Having looked into it, I have found many forum posts from people with a similar result from integration. Apparently this is a property of NDSolve that takes into account the possibility of a complex differential equation as the general solution and therefore, as a general solution, returns complex numbers. As you have been so helpful in sharing your knowledge so far, I was wondering if you would know how to confine the integration to the real domain in the attached file? Many thanks, Robin Attachments:
Posted 5 years ago
 Hi Robin,I don't have an answer for this. We do have existence and uniqueness theorems for linear ODEs. But this in clearly nonlinear, where things are far more complicated. It's been a long time since my classes in analysis, and even then nonlinear systems were a special topic. Maybe one of the mathematicians in the group will speak up, but I don't know of a way to test this system for real solutions. And I don't see domain options to NDSolve which could guide it to a real solution, if one exists.Best regards, David
Posted 5 years ago
 Hi David,That's no trouble, I was just wondering if there was just a quick 'assumptions -> ' method which I had missed out on however that appears not to be the case! As this has somewhat deviated from the original question I may open up a new discussion in the Maths section. It appears to be a common problem so hopefully that may prove useful to several people. I cannot thank you enough for volunteering your time and knowledge, keep up the good work David!Best regards, Robin
Posted 5 years ago
 Hi Robin. I saw your new post. I'm glad you did that -- I am also very interested in this. I hope one of the WRI mathematicians will comment.
Posted 5 years ago
 Yes, a very quick look at the equations does not show me that I'd expect complex results, but I'd have to dig in deeper.