Your opinion on a class of ODE problems

[Namir]

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

[ttw]

I haven't looked at concentration equations but there are some formulations that try to conserve energy throughout long computations. You might look up some papers on these methods. Often they try to use midpoint differences or backward differences rather than forward differences. Perhaps some of the methods used to conserve energy could help in conserving mass.

I don't have access to journals (except for free ones on the net) so I can't check things out easily.

https://en.wikipedia.org/wiki/Leapfrog_integration

https://web.njit.edu/~jiang/math614/dahlquist.pdf

https://en.wikipedia.org/wiki/Numerical_..._equations

http://www.unige.ch/~hairer/preprints/santander.pdf

https://people.maths.ox.ac.uk/suli/nsodes.pdf

http://user.it.uu.se/~perl/ODE_07/NumODEsbook.pdf

These may give you some indication of what method to use.

[Thomas Ritschel]

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

[Namir]

(01-22-2016 05:45 AM)ttw Wrote:  I haven't looked at concentration equations but there are some formulations that try to conserve energy throughout long computations. You might look up some papers on these methods. Often they try to use midpoint differences or backward differences rather than forward differences. Perhaps some of the methods used to conserve energy could help in conserving mass.

I don't have access to journals (except for free ones on the net) so I can't check things out easily.

https://en.wikipedia.org/wiki/Leapfrog_integration

https://web.njit.edu/~jiang/math614/dahlquist.pdf

https://en.wikipedia.org/wiki/Numerical_..._equations

http://www.unige.ch/~hairer/preprints/santander.pdf

https://people.maths.ox.ac.uk/suli/nsodes.pdf

http://user.it.uu.se/~perl/ODE_07/NumODEsbook.pdf

These may give you some indication of what method to use.

Thanks for the equations. You mentioned conservation of energy. My problem is similar to solving a system of ODEs that calculate temperatures of various connected objects, where the calculation "errors" end up violating the laws of conservation of energy.

I am using the system of ODE to calculate concentrations of chemicals used in optimization problems. So far, I have used chemical reactions whose ODE have analytical solutions. So I was fine. I am trying to stretch in dealing with chemical reactions whose ODE have no easy or published analytical solutions.

[Marcio]

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

[Namir]

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.

[Marcio]

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?

[ttw]

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.

[Namir]

(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

[gestrickland]

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

[pascal_meheut]

(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?

[Namir]

(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

[SlideRule]

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

[Namir]

(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