I am using various numerical methods to solve a system of differential equations. These equations model chemical reactions. To be specific, the ODEs are sets of equations for monitoring the change in the concentration of chemical that participate in reactions -- chained irreversible reactions and chained consecutive equilibrium reactions. Due to rounding and inherit errors in the methods, I discovered that the solutions violate the "conservation of matter law" because the sum of the concentrations (or weights to be more exact) of the chemicals change with time!

My first fix to this problem is to readjust the concentration of one of the chemicals (say the main chemical that starts the reaction) to preserve the conservation of matter law. The conservation of matter law represents a constraint on the ODE solutions, which admittedly is unusual.

Your opinion? and suggestions? Is any fix corrupting the "math" and rendering the numerical solution invalid?

Namir

It's quite likely that one runs into conservation law violations when trying to (numerically) simulate complex chained reactions with simple ODE integrators (like the Newton model). Especially for long computations or when it comes to oscillating reactions one will notice inconsistencies in total mass (but also charge, energy, etc.).

One possible (quick-and-dirty) solution would be some kind of intermediate normalization. But better use a predictor-corrector integrator (e.g. Adams-Bashforth-Moulton) instead.

Here are two further links regarding chemical kinetics simulations:

http://chemistry.illinoisstate.edu/stand...12-S06.pdf
http://www.springer.com/cda/content/docu...306-c1.pdf
Would it be possible at all for you to post the equations here? Or at least the reaction mechanism you derived the equations from?

This might be the case in which you might need to add an algebraic equation to your system of differential equations to check if things are right.

Marcio

For the consecutive reversible reactions:

A <--> B <--> C

The equations are:

d[A]/dt = -k1[A] + k'1[B]

d[B]/dt = k1[A] - (k'1 + k2) [B] + k'2[C]

d{c]/dt = k2[B] - k`2[C]

Where [A], [B], and [C] are the concentrations for chemicals A, B, and C, respectively. The forward reaction rates constants are k1 and k2. The reverse reaction rate constants are k'1 and k'1.

As you can see the ODEs are not complicated. The analytical solutions are complicated. I am trying to find a general method where the analytical solution for the differential equations is way too complicated.

Hm. The equations are very simple. So the problem is that the methods you've been using to solve them numerically are returning unsatisfactory results?

If you are not using a canned method, you might try deriving a method using 3 steps so that A is updated on steps 1(mod 3), B on steps 2(mod 3), and C on steps 0(mod 3). This type of overlapping update has helped in dynamic equations. Of course, (as I didn't derive anything), the error terms would have to be worked out.

(01-22-2016 07:41 PM)ttw Wrote: [ -> ]If you are not using a canned method, you might try deriving a method using 3 steps so that A is updated on steps 1(mod 3), B on steps 2(mod 3), and C on steps 0(mod 3). This type of overlapping update has helped in dynamic equations. Of course, (as I didn't derive anything), the error terms would have to be worked out.

I am using Runge-Kutta 4th order, Runge-Kutta Fehlberg, Cash-Karp, and Dormandâ€“Prince ODE methods. I divide each interval into 10 sub-intervals. I then apply the ODE methods on the sub-intervals. At the end of each iteration of a sub-interval I re-calculate the concentration of one of the chemical, say compound A, as equal to the

initial sum of concentrations (for all chemicals) minus the sum of the

current concentrations of the

other chemicals. This approach maintains the conservation of matter law.

I will do some more testing using reactions with known analytical solutions vs using the ODE methods. This kind of comparison will be very telling.

Namir

Namir,

In your example equations, if you add the three together, you get d(A+B+C)/dt=0, which means A+B+C=M, M being a constant throughout the solution. I suppose this is your conservation of mass, actually embedded in the equations. If I am not missing something, you could use this equation to eliminate one of the unknowns, say A, from the remaining two equations, which you could then solve for B and C with the numerical method of your choice, at each step calculating A from the above conservation relation. Would not this work?

Gordon

(01-27-2016 02:45 AM)gestrickland Wrote: [ -> ]In your example equations, if you add the three together, you get d(A+B+C)/dt=0, which means A+B+C=M, M being a constant throughout the solution. I suppose this is your conservation of mass, actually embedded in the equations. If I am not missing something, you could use this equation to eliminate one of the unknowns, say A,

Isn"t more or less the same as :

(01-22-2016 08:08 PM)Namir Wrote: [ -> ]At the end of each iteration of a sub-interval I re-calculate the concentration of one of the chemical, say compound A, as equal to the initial sum of concentrations (for all chemicals) minus the sum of the current concentrations of the other chemicals. This approach maintains the conservation of matter law.

but maybe less efficient?

(01-27-2016 02:45 AM)gestrickland Wrote: [ -> ]Namir,

In your example equations, if you add the three together, you get d(A+B+C)/dt=0, which means A+B+C=M, M being a constant throughout the solution. I suppose this is your conservation of mass, actually embedded in the equations. If I am not missing something, you could use this equation to eliminate one of the unknowns, say A, from the remaining two equations, which you could then solve for B and C with the numerical method of your choice, at each step calculating A from the above conservation relation. Would not this work?

Gordon

Good point! I can replace [A] with = M - [B] - [C]. I think it will work and make the solution a bit simpler. I can still do my optimization for a nonlinear curve fit by minimizing the square root of:

Sum of ([A(i)]obs - [A(i)]calc)^2 + ([B(i)]obs - [B(i)]calc)^2 + ([C(i)]obs - [C(i)]calc)^2 for i=1 to Number of Data points

Will give it a try!

Thanks,

Namir

Namir

May I suggest the following reference material?

1. Differential Equation Analysis in Biomedical Science & Engineering, ODE applications w/ R

2. Differential Equation Analysis in Biomedical Science & Engineering, PDE applications w/ R

BEST!

SlideRule

(01-28-2016 01:31 PM)SlideRule Wrote: [ -> ]Namir

May I suggest the following reference material?

1. Differential Equation Analysis in Biomedical Science & Engineering, ODE applications w/ R

2. Differential Equation Analysis in Biomedical Science & Engineering, PDE applications w/ R

BEST!

SlideRule

Thanks for the referral. I ordered the ODE version of the book!

Namir