Message Boards Message Boards


Problem with Numerical Integration of Aspirin Synthesis Reaction

Posted 5 years ago
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},
    {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, 
12 Replies

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:


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

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!


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!

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

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.

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.

Hi, Taking a look at your equations, your problem is the fractional power terms (...)^d and especially (...)^c . If the arguments in the parenthesis become negative you can get complex answers. If you don't expect complex answers, then something is wrong with the arguments, which could be the signs of the terms or the initial conditions. Cheers, Kay

Posted 5 years ago

Dear Kay,

You're absolutely right, thank you very much! With this information I have now managed to integrate the system of equations and am very happy with the solution which has resulted.

Thank you kindly for having a look at this, I'm very grateful!

Best regards, Robin

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract