Message Boards

WOLFRAM COMMUNITY

12081 Views

12 Replies

2 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Problem with Numerical Integration of Aspirin Synthesis Reaction

Robin Hartley

Posted 9 years ago

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

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

POSTED BY: Robin Hartley

12 Replies

Sort By:

Kay Herbert

Kay Herbert, Chief scientist, Sloan Valve Co.

Posted 9 years ago

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 BY: Kay Herbert

Robin Hartley

Posted 9 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

POSTED BY: Robin Hartley

David Reiss

David Reiss, Scientific Arts

Posted 9 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.

POSTED BY: David Reiss

David Reiss

David Reiss, Scientific Arts

Posted 9 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}

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 BY: David Reiss

Robin Hartley

Posted 9 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 BY: Robin Hartley

David Keith

Posted 9 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 BY: David Keith

Robin Hartley

Posted 9 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 BY: Robin Hartley

David Keith

Posted 9 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 BY: David Keith

Robin Hartley

Posted 9 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 BY: Robin Hartley

David Keith

Posted 9 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 BY: David Keith

Robin Hartley

Posted 9 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 BY: Robin Hartley

David Keith

Posted 9 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 BY: David Keith

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback