Numerical Integration of the Chemical Rate Equations via a

ARTICLE pubs.acs.org/IECR

Numerical Integration of the Chemical Rate Equations via a Discretized Adomian Decomposition Jarod M. Younker*,† Department of Chemistry, The Pennsylvania State University, University Park, Pennsylvania 16802, United States

bS Supporting Information ABSTRACT: Chemists are frequently interested in rate equations, which are first-order differential equations. Numerical integration of these equations allows the researcher to accurately predict the concentrations of chemical species at any time given the initial conditions. Explicit Runge-Kutta (RK) integration is widely used for solving the rate equations. In this article, Adomian decomposition methods (ADM) are used to obtain the solutions of chemical rate equations. The Adomian method outlined here outperforms high-order RK routines in the arenas of accuracy and truncation error. Additionally, four modifications are introduced that place the Adomian integration on par with RK in terms of speed (a primary reason for which Adomian decomposition methods are currently underemployed). The inclusion of up to the fifth term in the Adomian expansion gives a truncation error of order O(h10). The method as presented yields solutions which are step-size independent in the nonstiff regime. The problem of rapid polynomial divergence is addressed through discretizing the time axis. Performance of the ADM method against an implicit algorithm is also given.

’ INTRODUCTION The construction of a chemical mechanism is one of the first steps a chemist takes following the discovery of a novel reaction. Many tools are employed in order to determine the sequence of intermediates and the rates at which they change. The results of these studies can be succinctly expressed as coupled differential equations, the solutions of which yield the concentrations of the experimentally observed reactants at all times of interest. To uncover the desired rate constants, chemists and biochemists frequently analyze reactions under limiting conditions. These pseudo-first-order, steady-state,1 and quasi-steady-state kinetics2 simplify the differential equations through elimination of one or more degrees of freedom. When these conditions cannot be met, kinetic inversion or global analysis is employed to determine the rates.3-14 These nonlinear regression techniques require solutions of the coupled chemical rate equations. Symbolic integration of the equations can be difficult and is, in most cases, impossible.5-7 Thus, numerical integration is necessary. One popular numerical integration method was originally proposed in 1895 by Runge and later expanded by Kutta.15,16 Fourth-order (RK4) and fourth-/fifth-order adaptive step size (RKF45) RungeKutta routines are frequently used to approximate ordinary differential equations as initial-value problems (IVP),17,18 of which the rate equations are prime examples.3,8,11 An alternative approach for solving differential equations was proposed by Adomian in the 1980s. His solution involves decomposing the nonlinear terms in the differential equation(s) into a series of polynomials.19,20 The integration of this infinite series of polynomials, under known initial or boundary conditions, is called the Adomian Decomposition Method (ADM).21 Truncation of the series yields an analytic approximation to the solution. The strength of this approach is its ability to provide analytical equations that describe the integrated differential equations.22-25 r 2011 American Chemical Society

Current ADM implementations often involve symbolic integration of a large series of polynomials.25-27 This symbolic integration is time-consuming, and as a result, ADM falls behind RK methods in performance. In what follows, we show that symbolic integration is unnecessary for linear first-order differential equations. This and four additional modifications to the ADM method place it on par with RK in terms of the computational time required to solve the chemical rate equations. Last, we avoid polynomial divergence by propagating the solution via sequential IVPs (i.e., the final values for one region are used as initial values for the next). Our initial interest in this method stemmed from kinetic inversion/global analysis problems wherein the RKF45 routine was found to be extremely slow and the RK4 was insufficient. These techniques extract rate constants and/or reaction orders from concentration profiles via nonlinear regression. The regression requires multiple solutions of the rate equations. In this article, we provide an ADM solution to coupled ordinary first-order differential equations. Furthermore, we develop a computational approach to the solutions which allows fast and easy numerical integration. Our decision to bypass analytic solutions eliminates the need to use symbolic solvers such as Matlab or Mathematica. We show that these ADM solutions outperform both RK4 and RKF45 in the arenas of accuracy and truncation error. The method, as presented, yields solutions that are step-size independent in the nonstiff regime. Within the stiff regime, like all explicit methods, the ADM solution requires extremely small step sizes to avoid polynomial divergence.

Received: April 12, 2010 Accepted: February 1, 2011 Revised: January 31, 2011 Published: February 21, 2011 3100

dx.doi.org/10.1021/ie1008647 | Ind. Eng. Chem. Res. 2011, 50, 3100–3109

Industrial & Engineering Chemistry Research

ARTICLE

’ EXPERIMENTAL SECTION

Comparing eqs 3 and 8 shows

Solution of the Chemical Rate Equations. The theory and proof of convergence of ADM have been documented.28 The ADM solution to the coupled differential equations of interest proceeds as follows: given p coupled differential equations where L(t) is the linear differential operator d/(dt),

y ðtÞÞ LðtÞ yi ðtÞ ¼ fi ðy1 ðtÞ, y2 ðtÞ, :::, yp ðtÞÞ ¼ fi ð B l LðtÞ yp ðtÞ ¼ fp ðy1 ðtÞ, y2 ðtÞ, :::, yp ðtÞÞ ¼ fp ð B y ðtÞÞ

¥

Given the need to recursively evaluate the individual terms in the series solution of yi(t), we require that the following relation be true: yði, n þ 1Þ ðtÞ ¼ L-1 ðtÞðAði, nÞ ðyBðtÞÞ þ

yði, 1Þ ðtÞ ¼ L ðtÞ ðAði, 0Þ ðyBðtÞÞ þ

ð2Þ

p

δðn,

The nonlinear term gi(yB(t)) is expanded in an infinite series of polynomials A(i,n)(yB(t)).

The Adomian polynomials A(i,n)(yB(t)) can be obtained via a judicious rearrangement of the more familiar Taylor series expansion of gi(yB(t)) (see Supporting Information). This arrangement allows an efficient method of recursively evaluating the terms in the series solution of yi(t). The remaining term hi(yB(t)) can be written as p linear terms yj(t) weighted by k(i,j): p

hi ðyBðtÞÞ ¼

∑ kði, jÞ yjðtÞ j¼1

ð6Þ

hi ð B y ðtÞÞ ¼

¥

∑ ð ∑ kði, jÞyðj, nÞðtÞÞ n¼0 j¼1

ð7Þ

Substitution of eqs 5 and 7 into eq 2 yields the classical ADM result for linear first-order differential equations (eq 8). It reveals the relationship between eqs 2 and 3 needed to obtain the integrated solutions. As written, it is exact, but in practice, only a finite number of terms can be included. yi ðtÞ ¼ yi ð0Þ þ L-1 ðtÞ

¥

∑ ðlðj, 1Þ þ 2lðj, 2Þ þ :::

j¼1

ð8Þ

lðp, nÞ yðp, nÞ

lð1, 1Þ ! 3 3 3 lðp, nÞ !

1 A ð12Þ

p

∑ ððA^ði, nÞ þ j∑¼ 1 kði, jÞyðj, nÞÞ bnðtÞÞ n¼0

ð13Þ

where bn ðtÞ ¼

tn þ 1 ðn þ 1Þ!

ð14Þ

This removes the need for symbolic integration. In the Supporting Information, we detail how this is also possible with functions that possess a direct polynomial dependence on t. Finally, in vector formulation, the ADM solution is

p

∑ ðAði, nÞð By ðtÞÞ þ j∑¼ 1 kði, jÞyðj, nÞðtÞÞ n¼0

lð1, 1Þ yð1, 1Þ

∂lð1, 1Þ þ ::: þ lðp, nÞ gi ð B y ð0ÞÞ lð1, 1Þ þ ::: þ lð1, nÞ lðp, 1Þ þ ::: þ lðp, nÞ :::∂yp ∂yi

yi ðtÞ ¼ yi ð0Þ þ

p

, , :::, lðp, nÞ ¼ 0

ð2 1Þ

where eq 12 has been adapted from that given by Abbaoui et al.29,30 As is apparent from eq 12, a large number of terms needs to be evaluated in order to generate the polynomials. A fourth-order RK method requires four evaluations of the differential equations at each time point, whereas ADM requires many times that number. The integrated solution of three coupled equations (p = 3, n = 5) may contain up to 108 terms. It can be shown that for the solution of time-independent firstorder equations, one may separate the time dependence from eq 8 by introducing the series bn(t) and the time-independent polynomial Â(i,n) (see Supporting Information for derivation), such that

Substitution of eq 3 for yj(t) yields the following expansion: ¥

lð1, 1Þ , lð1, 2Þ

∑ , :::, l

þ ðn - 1Þlðj, n - 1Þ þ nlðj, nÞ ÞÞ

ð4Þ

ð5Þ

ð11Þ

n

where yi(0) = y(i,0)(t). Upon the basis of his approach, we need to uncover a relationship between L-1(t) fi(yB(t)) and the series given in eq 3. Our first step is to realize that the function fi(yB(t)) can be expressed as a sum of linear, hi(yB(t)), and nonlinear terms, gi(yB(t)):

¥

∑ kði, jÞ yðj, 0Þ ðtÞÞ

j¼1

where A(i,0)(yB(t)) is the first term in the Taylor expansion of gi(yB(t))|t=0 (see Supporting Information). The subsequent A(i,n)(yB(t)) terms are determined recursively: Aði, n6¼ 0Þ ð B y ðtÞÞ ¼

∑ ðAði, nÞð By ðtÞÞÞ n¼0

ð10Þ

p

-1

ð3Þ

fi ðyBðtÞÞ ¼ gi ðyBðtÞÞ þ hi ðyBðtÞÞ

∑ kði, jÞyðj, nÞðtÞÞ j¼1

Solving for y(i,1)(t):

¥

∑ yði, nÞðtÞ ¼ n∑¼ 0 yði, nÞðtÞ n¼1

gi ðyBðtÞÞ ¼

p

ð9Þ

ð1Þ

Solving eq 2 for coupled chemical rate equations is the focus of this paper. Adomian proposed that the solution could be written as an infinite series of polynomials in t:22 yi ðtÞ ¼ yi ð0Þ þ

¥

p

Operating L-1(t) on both sides of eq 1 yields the solutions yi(t) minus a constant yi(0), which is the initial value. yi ðtÞ ¼ yi ð0Þ þ L-1 ðtÞfi ðyBðtÞÞ

¥

∑ yði, nÞ ðtÞ ¼ L-1 ðtÞ n∑¼ 0ðAði, nÞð By ðtÞÞ þ j∑¼ 1 kði, jÞyðj, nÞðtÞÞ n¼1

y ðtÞ ¼ B y 0 þ Y^ bBðtÞ B

ð15Þ

where the matrix Ŷ(i,n)(yB(t)) encompasses the nonlinear and linear terms. This approach eliminates the need to categorize the different 3101

dx.doi.org/10.1021/ie1008647 |Ind. Eng. Chem. Res. 2011, 50, 3100–3109


ARTICLE

terms in the differential equations, making the algorithm much easier to implement computationally. The two approaches can be shown to be equivalent (Supporting Information). In conclusion, we have a time-independent matrix Â(i,n) þ ∑pj = 1 k(i,j) y(j,n) or Ŷ(i,n), which can be used to solve the IVP for any time. This matrix need only be calculated once. However, as the Adomian series is truncated, this matrix is only valid within the y 0, the series region near B y 0. For times outside the neighborhood of B quickly diverges.31 The neighborhood can be expanded by including more terms, but a better approach has been put forth by Ghosh et al.24 They proposed solving multiple IVPs through discretizing the time axis. They successfully applied this approach to single equation, one-dimensional systems but did not provide a solution for coupled differential equations. We solved the case of the multidimensional rate equations by separating the polynomial term coefficients Ŷ(i,n) from the time dimension bB(t). This allows Ŷ(i,n) to be propagated from yi(0) to time tq, where a new IVP is formulated with initial conditions yi(tq); tq is subsequently treated as the new time zero. Earlier, we spoke of four modifications that increase performance for numerical ADM. The first, the separation of variables and the use of a propagator, has been detailed. We introduce two parameters that describe the propagator. The first is the step size h. This parameter describes the interval between the values that are sampled from the continuous ADM solutions yi(t) as well as the step between discrete time values at which the RK methods evaluate yi(t). The discrete time axis for the ADM solutions is governed by tq, which is the length of the discrete time axis or the time after which a new Ŷ(i,n) is evaluated with initial values taken from the final values of the previous solution. For example, if h = 10 ms and tq = 1 s, ADM solutions would be sampled and recorded every 10 ms up until 1 s. Then, the values of the ADM solution at time 1 s would be used as initial values for a new IVP, and values would be sampled up until 2 s every 10 ms and so on. The second modification is a capping differential order. For equation fi(yB(t)) with a second-order dependence on yi(t)’s: y ðtÞÞ ¼ kði, jÞ yi yj fi ð B

ð16Þ

Any third-order derivative of fi(yB(t)) (e.g., (∂ would be zero. This fact significantly decreases the number of terms in the polynomial. The third modification is an adaptive discrete axis. Instead of altering the step-size h (like RKF45), the distance over which the local propagator is valid is changed dynamically. A comparison of the solutions evaluated at tq þ h and h for the prior and current propagators, respectively (these correspond to the same value of t), determines whether the prior propagator can be extended to greater times by increasing tq. A tolerance tol is in place wherein, if the difference between the two evaluations is smaller than some fraction of tol, tq is extended. If it is greater than a fraction of tol, tq is decreased, and the solution of yi(t) is returned to the prior region. The last modification is to accept that derivatives of fi(yB(t)), which are identically zero at t = 0, are zero for every local approximation. Thus, those terms are not evaluated. Consider, for example, the following rate equations where y3(t) has no dependence on y2: 3

dy1 ¼ - k1 ½y1 ½y3 þ k2 ½y2 dt dy2 ¼ k1 ½y1 ½y3 - k2 ½y2 dt dy3 ¼ - k1 ½y1 ½y3 dt

fi(t))/(∂y2i ∂yj))

ð17Þ

Figure 1. Numerical solutions for scheme 201 where y1 is shown in red, y2 in blue, and y3 in black. Analytic solution (solid lines), RK4 h = 0.5 s (dashed lines), ADM5 h = 0.5 s, tq = 10.0 s, all modifications enabled (solid lines). Inset shows differences.

Therefore, (∂y3)/(∂y2) = (∂2y3)/(∂y1∂y2) = ... = 0. This modification may not always be valid. For example, it would fail for dyi ¼ ki ½yi ½yj ¼ fi ðtÞ dt

ð18Þ

where [yj]0 = 0.0. Why would it fail? ∂fi ð0Þ ¼ ki ½yj 0 ∂yi

ð19Þ

The above derivative is identically zero at time zero, but if [yj] is nonzero at t > 0, then (∂fi(0))/(∂yi) would exist. It is necessary to make some remarks about stiff differential equations. A problem is classified as stiff if one or more components undergo rapid changes (decay or growth), while other components simultaneously experience very little or no net change.32 An example familiar to chemists would be a secondorder reaction under pseudo-first-order conditions. In this regime, extremely small step sizes are needed, both for RK and ADM approximations. Stiffness is defined in this paper as the minimum eigenvalue of the first derivative matrix.32 The ADM routines are dependent upon the first derivative of the differential equations, making this value trivial to compute. The kinetic schemes and conclusions presented in the Results and Discussion are nonstiffspecific. We will consider routines that include n = 3-10 Adomian polynomials (i.e., ADM(n)). Runge-Kutta routines RK4 and RKF45 were derived from Linz and Wang’s text.32 The integration tolerance for the Runge-Kutta routines was 1 10-8 . The “correct” integrated equations, with the exception of scheme 201, to which an analytical solution exists, were taken as the RKF45 solutions with an initial step-size of 1 10-3 s (they matched with h = 1 10-4 s). If needed, cubic splines were fit to these “correct” data, and interpolated points were used to quantify differences. Sums of absolute differences (l1) were compared. For kinetic inversion problems, the concentration profile with the smallest step size was interpolated to smaller data sets via a cubic spline. Fits employed an l2 merit function (i.e., square root of the sums of differences squared). Central difference finite derivatives of second order were used throughout.33 All derivatives were calculated via a recursive algorithm capable of determining partial derivatives of any order. Derivatives were calculated on the fly and stored for quick retrieval if needed again. Calculations were run on a 2.4 GHz Intel 3102



ARTICLE

where [y1]0 = 1 10-4 M, [y2]0 = 2 10-4 M, [y3]0 = 0 100 M and k1 = 300.0 s-1 M-1. Stiffness = -0.09. Kinetic scheme 201 has an exact analytical solution.1 Kinetic scheme 212 is defined as

Duo-Core Mac. All classes were written in Cþþ on top of the Template Numerical Toolkit (TNT; available at http:// math.nist.gov/tnt/). Classes can be obtained by writing to the author. Kinetic scheme 201 is defined as dy1 dy2 ¼ ¼ - k1 ½y1 ½y2 dt dt dy3 ¼ k1 ½y1 ½y2 dt

dy1 dt dy2 dt dy3 dt dy4 dt dy5 dt dy6 dt

ð20Þ

Table 1. Runge-Kutta Results for Scheme 201a RKF45

RK4

h

m

τ

l1/m

m

τ

0.1

1001

6.08

3.11

1001

3.22

3.14

0.2

501

2.99

6.18

501

1.61

6.20

0.5

201

1.37

15.04

201

0.66

15.06

1.0

123

0.82

15.17

101

0.36

28.36

2.0

115

0.91

19.37

51

0.18

48.15

l1/m

¼ - k1 ½y1 ½y4 þ k3 ½y3 ½y5 þ k5 ½y2 þ k6 ½y3 ¼ k1 ½y1 ½y4 - k2 ½y2 ½y5 - k5 ½y2 ¼ k2 ½y2 ½y5 - k3 ½y3 ½y5 - k6 ½y3

ð21Þ

¼ - k1 ½y1 ½y4 ¼ - k2 ½y2 ½y5 - k3 ½y3 ½y5 þ k4 ½y6 2 =2 ¼ k2 ½y2 ½y5 þ k3 ½y3 ½y5 - k4 ½y6 2

where [y1]0 = 12 10-6 M, [y2]0 = [y3]0 = [y6]0 = 0 100 M, [y4]0 = 12 10-6 M, [y5]0 = 48 10-6, k1 = 3.8 106 s-1 M-1, k2 = 3.5 103 s-1 M-1, k3 = 3.3 102 s-1 M-1, k4 = 1 107 s-1 M-1, k5 = 0.25 s-1, and k6 = 0.01 s-1. Stiffness = -91.2.

h is the initial step size in s, m is the number of integration points, τ is ̂ the computational duration in ms, and l1/m = (1/m) ∑m i |yi - yi| is scaled by 1E6. a

Table 2. ADM Results for Scheme 201a ADM9 modification

h

m

τ

dddd

0.1

1001

eddd

0.1

1001

eedd

0.1

eeed

0.1

eeee

ADM7 l1/m

τ

81405.26

0.00

489.57

0.00

1001

6.36

1001

4.41

0.1

1001

dddd

0.2

eddd

ADM5

ADM3

l1/m

τ

32821.97

0.00

4706.32

0.00

540.74

0.00

394.39

0.00

251.87

0.00

158.29

0.00

0.15

5.18

0.83

3.57

4.79

2.41

29.72

0.19

3.60

1.13

2.57

6.88

1.79

42.15

2.77

0.19

2.45

1.13

1.81

6.88

1.38

42.15

501

40431.59

0.00

16384.53

0.00

2343.91

0.00

269.56

0.00

0.2

501

244.52

0.00

197.15

0.00

125.87

0.00

78.81

0.00

eedd eeed

0.2 0.2

501 501

5.76 3.79

0.15 0.19

4.56 2.99

0.83 1.12

3.04 2.04

4.78 6.81

1.98 1.49

29.70 41.88

eeee

0.2

501

2.05

0.19

1.71

1.12

1.21

6.81

0.86

41.88

dddd

0.5

201

16181.13

0.00

6495.49

0.00

935.94

0.00

107.25

0.00

eddd

0.5

201

97.11

0.00

78.46

0.00

50.08

0.00

31.43

0.00

eedd

0.5

201

5.22

0.15

4.21

0.83

2.86

4.77

1.75

29.62

eeed

0.5

201

3.26

0.19

2.63

1.12

1.73

6.83

1.12

41.88

eeee

0.5

201

1.68

0.19

1.35

1.12

0.92

6.83

0.66

41.88

dddd eddd

1.0 1.0

101 101

8022.95 48.08

0.00 0.00

3194.76 38.78

0.00 0.00

464.49 24.75

0.00 0.00

52.82 15.41

0.02 0.02

eedd

1.0

101

5.07

0.15

4.09

0.83

2.63

4.76

1.66

29.52

eeed

1.0

101

3.12

0.18

2.51

1.05

1.74

6.44

1.03

40.24

eeee

1.0

101

1.54

0.18

1.32

1.05

0.85

6.44

0.54

40.24

dddd

2.0

51

3932.54

0.00

1568.10

0.00

227.54

0.00

25.84

0.15

dddd

2.0

51

3932.54

0.00

1568.10

0.00

227.54

0.00

25.84

0.15

eddd

2.0

51

23.76

0.00

19.05

0.00

12.17

0.00

7.53

0.15

eedd eeed

2.0 2.0

51 51

5.10 3.04

0.15 0.17

4.10 2.45

0.83 1.02

2.57 1.57

4.75 6.22

1.61 0.99

29.40 39.05

eeee

2.0

51

1.48

0.17

1.32

1.02

0.86

6.22

0.50

39.05

l1/m

τ

l1/m

The first column designates which modifications are enabled (e) and disabled (d): first letter f capping derivative, second f propagator (tq = 10 s), third f adaptive discrete axis (tolerance = 1 10-2), fourth f initial zero derivatives excluded. h is the step size in s, m is the number of ̂ integration points, τ is the computational duration in ms, and l1/m = (1/m) ∑m i |yi - yi| is scaled by 1E7. Results for ADM4, ADM6, ADM8, and ADM10 can be found in the Supporting Information. a

3103



ARTICLE

Figure 2. Left: Scheme 201 total integration time versus initial step size h. Right: l1/m difference versus initial step size h. RK4 (red circles), RKF45 (blue crosses), ADM5 tq = 10.0 s, all modifications enabled (black squares).

Figure 3. Scheme 201 ADM5 modification statistics. Left: Total integration time versus step size h. Right: l1/m difference versus step size h. Modifications which are enabled (e) and disabled (d): first letter f capping derivative, second f propagator (tq = 10 s), third f adaptive discrete axis (tolerance = 1 10-2), fourth f initial zero derivatives excluded. dddd (red circles), eddd (blue crosses), eedd (green up-triangles), eeed (purple downtriangles), and eeee (black squares).

Table 3. Runge-Kutta Results for Scheme 212a RKF45

RK4

h

m

τ

l1/m

m

τ

l1/m

0.001

10001

96.77

0.00

10001

56.43

0.00

0.002

5001

49.19

0.65

5001

28.14

0.86

0.005

2001

20.76

2.64

2001

11.22

2.86

0.01

1001

11.24

6.30

1001

5.6

6.52

h is the initial step size in s, m is the number of integration points, τ is ̂ the computational duration in ms, and l1/m = (1/m) ∑m i |yi - yi| is scaled by 1E8. a

Figure 4. Mechanism for scheme 212 where y1 (ferric enzyme) is shown in red, y2 (compound-I) in blue, and y3 (compound-II) in black. In purple are the reagents y4 (peracetic acid), y5 (ascorbic acid), and y6 (dehydroxyascorbate). The individual rate constants are also defined.

’ RESULTS AND DISCUSSION Scheme 201. Bimolecular Reaction. Two different kinetic schemes are analyzed. Scheme 201 is a second-order reaction to which an analytical solution exists.1 This scheme will be used to determine which approach (RK or ADM) more closely approximates the correct solution. Representative concentration profiles for scheme 201 are shown in Figure 1. It is clear from Tables 1 and 2 and Figure 2 that the ADM approximation is far superior to RK. The ADM solutions are orders of magnitude closer to the actual analytical results, even when in comparison with RK solutions at low step size h. This is a significant increase in accuracy and highlights that Adomian’s theories are not only valid but closer to reality for first-order IVPs than the commonly accepted Runge-Kutta approximations. Select computational times and the average of the absolute errors (l1/m = (1/m)

̂ ∑m i |yi-yi|) are found in the tables where m is the number of integration points, yi is a concentration point from the “correct” reference profile, and ̂yi is the corresponding ADM- or RKevaluated point. The order of the truncation error of the Runge-Kutta routines is O(h4).32 In other words, the truncation error increases as the step size increases, not linearly, but to the fourth power (where h is assumed to be less than one). This is easily seen in the RK4 results where the error grows with increasing h. For the scheme 201 RKF45 solutions, the error likewise increases with step size but reaches a point where h (g 1.0 s) is adapted to allow no further increase in truncation error (see Figure 2). The ADM order of error has been shown to be O(hp(kþ1)), where p is the number of terms taken from the Adomian series and k is the order of the differential operator L (k = 1).24 For ADM5, we obtain O(h10). This orders of magnitude improvement in the error is observed in the insensitivity of the ADM solutions to step size (Figure 2); the error increases 6 orders of magnitude slower than fourth-order RK. As shown in Figure 3, the l1/m error when no modifications are enabled obeys the O(h10) trend 3104



ARTICLE

Table 4. ADM Results for Scheme 212a ADM9

ADM7 l1/m

τ

23432.77

0.64

2621.21

0.64

270.89

0.63

ADM5 l1/m

τ

17739.16

0.64

1984.99

0.63

206.38

0.62

ADM3

modification

h

m

τ

l1/m

τ

l1/m

eddd

0.001

10001

eedd

0.001

10001

10144.36

0.64

5619.39

0.64

1135.76

0.59

632.68

eeed

0.001

10001

0.40

121.60

0.53

69.99

0.48

eeee

0.001

10001

91.98

0.64

65.69

0.63

40.09

eddd

0.002

5001

11721.70

0.64

8870.24

0.64

5075.15

0.54

25.84

0.59

0.64

2814.19

eedd

0.002

5001

2352.05

0.63

1781.36

0.62

1020.13

0.55

566.01

0.63 0.41

eeed

0.002

5001

264.77

0.63

201.18

0.61

116.86

0.51

65.90

0.51

eeee eddd

0.002 0.005

5001 2001

84.19 4686.45

0.64 0.64

58.93 3548.56

0.62 0.64

34.20 2030.44

0.52 0.63

21.38 1124.21

0.63 0.59

eedd

0.005

2001

2352.82

0.63

1779.51

0.62

1019.93

0.55

564.33

0.41

eeed

0.005

2001

267.42

0.63

202.32

0.61

116.20

0.51

64.85

0.49

eeee

0.005

2001

81.16

0.64

56.30

0.62

31.62

0.53

19.17

0.60

*eddd

0.01

1001

2338.53

0.63

1772.11

0.61

1013.83

0.54

561.03

0.43

eeed

0.01

1001

256.51

0.63

193.92

0.61

111.13

0.51

61.78

0.50

eeee

0.01

1001

76.67

0.64

53.07

0.62

29.21

0.53

17.42

0.62

The first column designates which modifications are enabled (e) and disabled (d): first letter f capping derivative, second f propagator (tq = 0.01 s), third f adaptive discrete axis (tolerance = 1 10-4), fourth f initial zero derivatives excluded. h is the step size in s, m is the number of integration ̂ points, τ is the computational duration in ms, and l1/m = (1/m) ∑m i |yi - yi| is scaled by 1E8. *eddd and eedd are equivalent for h = 0.01 s. Results for ADM4, ADM6, ADM8, and ADM10 can be found in the Supporting Information. a

(red circles). When the capping derivative is enabled, there is no increase in the error. This is expected as the modification is exact. However, the use of the propagator does significantly impact the overall magnitude of the error (green up-triangles). Likewise, enabling the adaptive discrete axis and excluding all terms with initial zero derivatives increases the error and places it an order of magnitude lower than that of the low h RK methods. Of interest is the insensitivity of the method to h when the modifications are enabled. The truncation error is constant over the step sizes considered. We feel that this independence, in connection with the computational speedup we will discuss next, is a fair trade for the loss in accuracy; the final ADM accuracy is better than RK. In modern computing, time is an important consideration. How much CPU time does the solution require? The user’s desired accuracy has an inverse relation to time. Four (five) evaluations of the rate equations are needed for each RK4 (RKF45) step.32 Therefore, the computational duration increases as four (five) times the number of steps or 4m (5m). For ADM5, the integrand requires the evaluation of 108 terms at each time step (108m, p = 3). How can we decrease computational cost? First, we limit the number of times the rate equations are evaluated. With the exception of n = 0, all other evaluations of fi(t) are needed to generate derivatives. Knowing that fi(yB(t)) has order dependence d on yi(t)’s, all derivatives greater than d will be zero. By excluding these terms (“capping”), we obtain a substantial gain in speed without losing accuracy (see Figure 3). Nevertheless, ADM is slower than RK. Our second modification eliminates the need to evaluate at every discrete point by taking advantage of the continuous nature of ADM. Separation of the time dimension allows us to use a local approximation at time zero and then propagate through t to other points of interest. This modification needs to be tempered as the series quickly diverges when the propagator leaves the neighborhood of the local approximation.31 The neighborhood can be expanded by including more terms, but a

Figure 5. Numerical solutions for scheme 212 where y1 is shown in red, y2 in blue, and y3 in black. Reference (solid lines), RKF45 h = 0.005 s (dashed lines), ADM5 h = 0.005 s, tq = 0.01 s, all modifications enabled (solid lines). Inset shows differences.

better approach has been put forth by Ghosh et al.24 They proposed solving multiple IVPs on adjacent discrete axes. Thus, the propagator is only evaluated at time zero of each discretized axis. This second modification puts ADM on par with RKF45 in terms of speed. A further speedup is observed if we allow the algorithm to adapt the discrete axis dynamically. A final modification decreases the number of derivatives that need to be calculated by accepting that if a derivative is identically zero for the first local propagator, it is zero throughout. This final modification may not be valid for all cases but works well here. RK methods scale linearly with m. In other words, if one doubles the number of time points, the time required to perform the integration likewise doubles. An analysis of ADM5 for this set of coupled differential equations reveals a scaling of 0.2m. For a small number of points, ADM5 trails behind RK4/RKF45. However, with larger data sets, our modified approach is faster (see Figure 2). Scheme 212. Chloroperoxidase. Now that we have shown that ADM solutions have a smaller error and can be run on par 3105



ARTICLE

Figure 6. Left: Scheme 212 total integration time versus initial step size h. Right: l1/m difference versus initial step size h. RK4 (red circles), RKF45 (blue crosses), ADM5 tq = 0.01 s, all modifications enabled (black squares).

Figure 7. Scheme 212 ADM5 modification statistics. Left: Total integration time versus step size h. Right: l1/m difference versus step size h. Modifications which are enabled (e) and disabled (d): first letter f capping derivative, second f propagator (tq = 0.01 s), third f adaptive discrete axis (tolerance = 1 10-4), fourth f initial zero derivatives excluded. eddd (blue crosses), eedd (green up-triangles), eeed (purple down-triangles), and eeee (black squares).

Table 5. Runge-Kutta Optimization Results for Scheme 212 0.001a

0.2

3989.9 466.3 0.211 0.048

6

40.2

241.2

0.002a 0.005a

0.2 0.1

3185.6 381.9 0.203 0.026 3418.5 444.9 0.209 0.034

5 5

20.0 8.0

100.2 40.2

0.01a

0.0

3747.6 497.8 0.190 0.048

5

4.1

20.4

0.001b

0.2

3989.0 466.2 0.211 0.048

6

6.0

36.0

0.002b

0.2

3185.5 381.9 0.202 0.026

5

3.0

15.0

0.005b

0.0

3418.2 444.9 0.209 0.034

5

1.2

6.0

0.01b

0.0

3747.6 497.9 0.191 0.048

5

0.6

3.0

standard deviation

k3

k5

k6

τ

l2

average

k2

steps τ/step

h

3585.2 447.7 0.203 0.039 328.4

45.3 0.009 0.010

a

RKF45 and. b RK4 results. h is the initial step size in s, l2 = ̂ 2 1/2 is scaled by 1E8, ki are the optimized rate constants, (∑m i (yi - yi) ) steps is the number of iterations in the optimization, τ/step is the average computational duration in s per iteration, and τ is the total optimization time in s. To include cubic spline interpolation of data m = 1001, add 56.6 s to τ for RKF45.

with RKF45, we turn our attention to the problem of interest: complex enzymatic reactions. Scheme 212 contains both firstand second-order rates and is seen in the reaction of peroxidases.34 The model was derived from work done on chloroperoxidase.35 The oxidation and reduction of chloroperoxidase is important as a model system for the study of P450 chemistry.36 Lambeir et al. unraveled the mechanism and determined the rate constants for several of these chloroperoxidase reactions.35 Using their reported numbers, we have constructed Scheme 212. A resting state ferric enzyme y1 is activated by peracetic acid y4 (a peroxide substitute) to a radical cationic species called compound-I (y2). Compound-I

is reduced (here, by ascorbic acid, y5) to yield what is known as compound-II (y3), which decays back to a ferric enzyme via first- and second-order reactions. Compound-I is also competitively decaying to ferric enzyme in a first-order manner. Dehydroxyascorbate y6 undergoes a fast disproportionation step, wherein one y6 is reduced by a second, yielding one equivalent of y5. We will assign y1, y2, and y3 as spectroscopically observable enzymatic species and assume y4, y5, and y6 are unobservable reactants. Select computational times and the average of the absolute errors (l1/m) are found in the tables where m is the number of integration points. The mechanism for scheme 212 is illustrated in Figure 4. The reference profiles used to determine differences are defined as the RKF45 h = 0.001 s solution. As such, the l1/m value is biased in favor of RKF45. The l1/m values in Tables 3 and 4 are relative, and the corresponding l1/m results for an ADM reference would be þ0.64 10-8 for RK and -0.64 10-8 for ADM. Figure 5 shows a resulting integrand for a moderate step size. Although appearing to be good for both methods, the early time RK points possess a 4-7% error (see inset). Truncation error increases as O(h4) for RK methods (see Figure 6). There is no gain in accuracy by employing the adaptive step size of RKF45. In contrast, the error of the ADM solutions increases several orders of magnitude slower, being independent of the step sizes considered. Both methods require step sizes h e 22 ms in order for the integrands to converge to a viable solution; the problem is quite stiff (stiffness = -91.2). The ADM step size h, like RK methods, is constrained according to the minimum eigenvalue λ of the Jacobian:37 h < 2=jλjmin 3106

ð22Þ



ARTICLE

Table 6. ADM5 Optimization Results for Scheme 212a h

l2

k2

k3

k5

k6

steps

τ/step

τ

eddd

0.001

101.4

3502.2

253.5

0.250

0.014

15

1013.2

15198.0

eedd

0.001

70.4

3502.0

253.1

0.250

0.014

15

106.5

1597.0

eeed

0.001

59.2

3501.8

250.3

0.250

0.014

15

12.4

186.0

eeee

0.001

59.1

3500.0

253.0

0.250

0.014

16

4.2

68.0

eddd

0.002

71.9

3503.3

253.7

0.250

0.014

13

504.9

6564.0

modification

eedd

0.002

49.1

3502.0

253.4

0.250

0.014

14

121.4

1700.0

eeed

0.002

42.4

3501.3

254.4

0.250

0.013

14

11.7

164.0

eeee eddd

0.002 0.005

43.6 45.2

3503.2 3504.0

254.1 254.7

0.250 0.250

0.014 0.014

15 13

3.8 204.5

53.0 2658.0

eedd

0.005

32.4

3505.6

289.0

0.250

0.012

14

115.4

1615.0

eeed

0.005

27.5

3501.9

253.4

0.250

0.014

15

11.3

170.0

eeee

0.005

29.6

3504.0

289.4

0.250

0.012

14

3.2

45.0

*eddd

0.01

19.8

3501.7

254.3

0.250

0.013

12

102.0

1224.0

eeed

0.01

18.6

3500.3

251.8

0.250

0.014

12

11.5

138.0

eeee

0.01

18.4

3499.7

252.9

0.250

0.014

16

4.1

49.0

3502.2 12.7

258.1 3.8

0.250 0.0003

0.014 0.0007

average standard deviation

The first column designates which modifications are enabled (e) and disabled (d): first letter f capping derivative, second f propagator (tq = 0.01 s), third f adaptive discrete axis (tolerance = 1 10-4), fourth f initial zero derivatives excluded. h is the initial step size in s, l2 = (∑m i (yi ̂yi)2)1/2 is scaled by 1E8, ki are the optimized rate constants, steps is the number of iterations in the optimization, τ/step is the average computational duration in s per iteration, and τ is the total optimization time in s. *eddd and eedd are equivalent for h = 0.01 s. a

An analysis of the integration times for small h finds ADM outperforming RK when our modifications are enabled; the tradeoff in accuracy is small (Figure 7). For this problem, ADM5 scales as 0.1m. Whereas the RKF45 cannot adapt its discrete step size, the continuous nature of the local ADM approximation yields an impressive speedup for large m. When the h = 0.001 s results are compared, ADM5 with all modifications enabled is 1.4 (2.4) times faster than RK4 (RKF45). Both RK and ADM are explicit or forward-integraton methods. In order to determine y(t þ 2), y(t þ 1) must be evaluated first. A more recent development is implicit methods wherein a nonlinear algebraic system is solved at each step.38,39 Implicit methods are better equipped to handle stiff equations, the tradeoff being the need to calculate the Jacobian matrix, as well solve the resulting Newton equations iteratively. The chloroperoxidase problem was solved with the IDA suite from the SUNDIALS package.39 IDA integration is a variable-order, variable-coefficient backward differentiation formula method. In summary, for m = 101, 201, 501, 1001, 2001, 5001, and 10001, the number of Jacobian evaluations and nonlinear iterations was 19 and 315, respectively (residual tolerance 1 10-10). For m = 1001 (10001), the computational duration was 17 (55) ms. The l1/m error is 0.34 10-8 for all m. Scheme 212. Kinetic Inversion. Kinetic inversion problems involve determining rate constants and reactions orders from concentration profiles.4,6-8,11,12,14 This minimization problem is solved using numerical optimization methods (i.e., nonlinear regression). Initial concentration profile guesses were produced from scheme 212 with k1 = 3.8 106 s-1 M-1, k2 = 1.0 103 s-1 M-1, k3 = 1.0 102 M-1 s-1, k4 = 1 107 M-1 s-1, k5 = 0.0 s-1, and k6 = 0.00 s-1. Reaction orders are known. We minimized the square root of the sum-squared differences (l2) via an interior-point algorithm that employs both Newton-Raphson and LevenbergMarquardt steps.40-42 Rate constants were constrained to be positive. The rate of formation of compound-I, k1, was constrained

to be exact (k1 can be determined independent of the other rates via exclusion of ascorbic acid), as was that of the disproportionation of dehydroxyascorbate, k4, which is orders of magnitude faster than the enzymatic rates. Tables 5 and 6 give the results employing RK and ADM5 integration. As noted earlier, the RK solutions for scheme 212 are only good at low h. If we compare the ADM5 solution with the lowest h RKF45 solution, we observe a speedup of ∼10 per iterative step in the optimization. This speedup does not include point interpolation, which further slows down the RKF45 routine. Interpolation is necessary because m will probably not equal m0 , where m0 is the size of the data to be fitted. Also worth noting is the impact of the O(h4) truncation error on the RK solutions. The data points are fit well, but the RK approximations yield inaccurate rate constants. The fitted rate constants have large standard deviations, up to 26% of the average rate and 30% off the actual value. However, the ADM5 O(h10) error allows rate constants to be determined accurately independent of step size. Furthermore, as ADM yields continuous solutions, there is no need for point interpolation. The largest error (in essence the only error) is found for k3 (error = 22%), which is the reduction of compound-II to ferric enzyme by ascorbic acid. All of the other constants are reproduced extremely well. Fitting against an ADM reference yields rates within the range of those shown in Table 6. The standard deviations are small, even for k2. As such, the results are independent of h (unlike RK). They also validate the speedup modifications needed to get ADM on par with RK. Caveats. The ability to propagate Ŷ(i,n) through time, thereby limiting the number of times the differential equations must be calculated, yields a significant speedup in processing time. An important caveat is that Ŷ(i,n) can only be propagated within the local neighborhood of the approximation. Time propagation outside of this region will yield an incorrect result. Luckily, 3107



ARTICLE

accurate and have a very low truncation error. Employing our modifications places numerical ADM at or beyond RK routines in terms of speed with little tradeoff in accuracy. However, ADM is an explicit method and as such is ill-suited to solve stiff equations. Further, systems involving a large number of species are problematic given the large number of possible terms in the Adomian polynomial.

’ ASSOCIATED CONTENT

bS Figure 8. Integration of scheme 201 with h = 1.0 s (tq = 20.0 s), where y1 is shown in red, y2 in blue, and y3 in black. Reference (dotted lines), ADM5 (solid lines). Cusps result at tq = 20 s as Ŷ(i,n) leaves the local neighborhood of the approximation and the integrand diverges.

the local neighborhood is easily identified following the calculation. Sharp points arise at the junctions of the discretized axis (see Figure 8). Although continuous, the derivatives at these points do not exist. These discontinuities inform us that the integration failed and a smaller propagation step is needed. A similar fail-safe does not exist with RK methods. This is the basis for adaptive tq. A comparison of the solutions evaluated at tq þ h and h for the prior and current propagators, respectively (these correspond to the same value of t), determines whether the prior propagator can be extended to greater times by increasing tq or whether the solutions are invalid and a smaller value of tq is needed. Within the current literature, it is common to find coupled rate equations involving a dozen or more species.12,43-47 This begs the question as to whether ADM is applicable to such problems. The increase in the number of differential equations is primarily reflected in A(i,n)(yB(t)). For example, when n = 5, we find that the number of possible terms is 7, 36, 108, 252, 506, 918, 1547, and 2464 for p = 1-8. As such, iterating over the nested sums in eq 12 is time-consuming and becomes prohibitive for coupled equations of large p. The two examples utilized in this article are demonstrative of simple isothermal reactions involving common bimolecular reactions. Equations involving noninteger reaction orders are fully amenable to our ADM approach. However, the use of the capping differential order speedup will become null. Equations involving time dependence within nonpolynomial terms (e.g., exponentials, trigonometric functions, logarithms) are not amenable to our numeric treatment of the ADM method, as it depends on the ability to separate the time dimension from the rest of the equation. For such problems, it will be necessary to utilize a symbolic integrator. Lastly, we must address the problem of stiff differential equations. The ADM approximation will not eliminate the need to use extremely small steps sizes at t near zero. However, the ability to extend the discrete axis dynamically provides us with a way to increase the efficiency of the algorithm at larger t. For stiff problems, an implicit method is the recommended option.

’ CONCLUSIONS In conclusion, a computationally simple solution for solving the coupled chemical rate equations is given as a truncated series of Adomian polynomials. The ADM routines detailed here outperform the commonly employed RK routines for the solution of such equations. The ADM solutions are extremely

Supporting Information. Derivations connecting the Adomian polynomials to the Taylor seriers. Separation of the time dimension from the Adomian polynomials. Derivations equating the classical and modified Adomian polynomials. ADM4, ADM6, ADM8, and ADM10 results for Schemes 201 and 212. This information is available free of charge via the Internet at http://pubs.acs.org/.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Present Addresses †

Oak Ridge Associated Universities, 1 Bethel Valley Road, Oak Ridge TN, 37831-6164. email: [email protected]

’ ACKNOWLEDGMENT We would like to thank Dr. Michael T. Green of the PSU Department of Chemistry and Serge Baliff of the PSU Department of Mathematics for their reading of the manuscript and input. J.M.Y. is a PSU Academic Computing Fellow. ’ REFERENCES (1) Atkins, P.; de Paula, J. Physical Chemistry, 7th ed.; Oxford University Press: New York, 2002; pp 862-898. (2) Li, B.; Shen, Y.; Li, B. Quasi-steady-state laws in enzyme kinetics. J. Phys. Chem. A 2008, 112, 2311–2321. (3) Gay, I. D. On the numerical integration of rate equations. J. Am. Chem. Soc. 1964, 86, 2747–2748. (4) Zimmerle, C. T.; Frieden, C. Analysis of progress curves by simulations generated by numerical integration. Biochem. J. 1989, 258, 381–387. (5) Topham, C. Computer simulations of the kinetics of irreversible enzyme inhibition by an unstable inhibitor. Biochem. J. 1986, 240, 817–820. (6) Duggleby, R. G. Analysis of progress curves for enzyme-catalyzed reactions: application to unstable enzymes, coupled reactions and transient-state kinetics. Biochim. Biophys. Acta 1994, 1205, 268–74. (7) Duggleby, R. G. Analysis of enzyme progress curves by nonlinear regression. Methods Enzymol. 1995, 249, 61–90. (8) Stojan, J. Analysis of progress curves in an acetylcholinesterase reaction: A numerical integration treatment. J. Chem. Inf. Comput. Sci. 1997, 37, 1025–1027. (9) Goudar, C. T.; Sonnad, J. R.; Duggleby, R. G. Parameter estimation using a direct solution of the integrated Michaelis-Menten equation. Biochim. Biophys. Acta 1999, 1429, 377–383. (10) Golicnik, M. On a nonelementary progress curve equation and its application in enzyme kinetics. J. Chem. Inf. Comput. Sci. 2002, 42, 157–161. (11) Liu, G.; Mi, Z.; Wang, L.; Zhang, X. Kinetics of dicyclopentadiene hydrogenation over Pd/Al2O3 Catalyst. Ind. Eng. Chem. Res. 2005, 44, 3846–3851. 3108


Industrial & Engineering Chemistry Research (12) Tang, W.; Zhang, L.; Linninger, A. A.; Tranter, R. S.; Brezinsky, K. Solving kinetic inversion problems via a physically bounded GaussNewton (PGN) method. Ind. Eng. Chem. Res. 2005, 44, 3626–3637. (13) Varon, R.; Garcia-Moreno, M.; Masia-Perez, J.; Garca-Molina, F.; Garca-Canovas, F.; Arias, E.; Arribas, E.; Garca-Sevilla, F. An alternative analysis of enzyme systems based on the whole reaction time: evaluation of the kinetic parameters and initial enzyme concentration. J. Math. Chem. 2007, 42, 789–813. (14) Maeder, M.; Zuberb€uhler, A. D. Nonlinear least-squares fitting of multivariate absorption data. Anal. Chem. 1990, 62, 2220–2224. :: (15) Runge, C. Uber die numerische Aufl::sung von Differentialgleichungen. Math. Ann. 1895, 46, 167. (16) Kutta, M. W. Ph. D. thesis, University of Munich, Munich, Germany, 1901. (17) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes: the Art of Scientific Computing; Cambridge University Press: New York, 1986; pp 547-577. (18) Norris, A. C. Computational Chemistry: An Introduction to Numerical Methods; John Wiley & Sons: New York, 1981. (19) Adomian, G. Systems of nonlinear partial differential equations. J. Math. Anal. Appl. 1986, 115, 235–238. (20) Adomian, G. Solution of coupled nonlinear partial differential equations by decomposition. Comput. Math. Appl. 1996, 31, 117–120. (21) Adomian, G. A review of the decomposition method in applied mathematics. J. Math. Anal. Appl. 1988, 135, 501–544. (22) Adomian, G. Solving Frontier Problems in Physics: the Decomposition Method; Kluwer Academic Publishers: Norwell, MA, 1994. (23) Olek, S. An accurate solution to the multispecies Lotka-Volterra equations. SIAM Rev. 1994, 36, 480–488. (24) Ghosh, S.; Roy, A.; Roy, D. An adaptation of adomian decomposition for numeric-analytic integration of strongly nonlinear and chaotic oscillators. Comput. Methods Appl. Mech. Eng. 2007, 196, 1133– 1153. (25) Kaya, D. A reliable method for the numerical solution of the kinetics problems. Appl. Math. Comput. 2004, 156, 261–270. (26) Abdelwahid, F. A mathematical model of Adomian polynomials. Appl. Math. Comput. 2003, 141, 447–453. (27) El-Wakil, S. A.; Abdou, M. New applications of Adomian decomposition method. Chaos Solitons Fractals 2007, 33, 513–522. (28) Gabet, L. The theoretical foundation of the Adomian method. Comput. Math. Appl. 1994, 27, 41–52. (29) Abbaoui, K.; Cherruault, Y.; Seng, V. Practical formulae for the calculus of multivariable Adomian polynomials. Math. Comput. Model. 1995, 22, 89–93. (30) The formulation for a single differential equation given by Abdelwahid is primarily found throughout the literature:26

ARTICLE

(36) Green, M. T.; Dawson, J. H.; Gray, H. B. Oxoiron(IV) in Chloroperoxidase is Basic: Implications of P450 Chemistry. Science 2004, 304, 1653–1656. (37) Tyreus, B. D.; Luyben, W. L.; Schiesser, W. E. Stiffness in Distillation Models and the Use of an Implicit Method to Reduce Computation Times. Ind. Eng. Chem. Des. Dev. 1975, 14, 427–433. (38) Pavlis, R. R. Kinetics without Steady State Approximations. J. Chem. Educ. 1997, 74, 1139–1140. (39) Hindmarsh, A. C.; Brown, P. N.; Grant, K. E.; Lee, S. L.; Serban, R.; Shumaker, D. E.; Woodward, C. S. SUNDIALS: Suite of Nonlinear and Differential/Algebraic Equation Solvers. ACM Trans. Math. Software 2005, 31, 363–396. (40) Nocedal, J.; Wright, S. Numerical Optimization; Springer: New York, 1999. (41) Fletcher, R.; Leyffer, S. Nonlinear programming without a penalty function. Math. Program., Ser. A 2002, 239–269. (42) W€achter, A.; Biegler, L. T. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Progr., Ser. A 2006, 106, 25–57. (43) Kocí, P.; Kubícek, M.; Marek, M. Modeling of Three-Way Catalyst Monolith Converters with Microkinetics and Diffusion in the Washcoat. Ind. Eng. Chem. Res. 2004, 43, 4503–4510. (44) Aranzabal, A.; Ayastuy-Arizti, J. L.; Gonzalez-Marcos, J. A.; Gonzalez-Velasco, J. R. Kinetics of the Catalytic Oxidation of Lean Trichloroethylene in Air over Pd/Alumina. Ind. Eng. Chem. Res. 2003, 42, 6007–6011. (45) Ayastuy, J. L.; Gutierrez-Ortiz, M. A.; Gonzalez-Marcos, J. A.; Aranzabal, A.; Gonzalez-Velasco, J. R. Kinetics of the Low-Temperature WGS Reaction over CuO/ZnO/Al2O3 Catalyst. Ind. Eng. Chem. Res. 2005, 44, 41–50. (46) Benaissa, M.; Carillo Le Roux, G.; Joulia, X.; Chaudhari, R. V.; Delmas, H. Kinetic Modeling of the Hydrogenation of 1,5,9-Cyclododecatriene on Pd/Al2O3 Catalyst Including Isomerization. Ind. Eng. Chem. Res. 1996, 35, 2091–2095. (47) Anikeev, V. I.; Yermakova, A.; Bogel-Lukasik, E.; Nunes da Ponte, M. Kinetics of Limonene Hydrogenation in High-Pressure CO2 at Variation of Hydrogen Pressure. Ind. Eng. Chem. Res. 2010, 49, 2084– 2090.

Aðn6¼ 0Þ ðyðtÞÞ n

¼

∑ v¼1

! nþ1-v 1 δðn, l1 þ l2 þ ::: þ lv Þ yðl1 Þ yðl2 Þ :::yðlv Þ g ðvÞ ðyð0ÞÞ v! l1 , l2 , :::, lv ¼ 1

∑

(31) Jiao, Y. C.; Yamamoto, Y.; Dang, C.; Hao, Y. An after treatment technique for improving the accuracy of Adomian’s decomposition method. Comput. Math. Appl. 2002, 43, 783–798. (32) Linz, P.; Wang, R. L. C. Exploring Numerical Methods: An Introduction to Scientific Computing Using MATLAB; Jones and Barlett Publishers: Sudbury, MA, 2003. (33) Hornbeck, R. W. Numerical Methods; Prentice-Hall: New York, 1975; pp 16-34. (34) Dunford, B. H. Heme Peroxidases; Wiley & Sons: New York, 1999. (35) Lambeir, A.-M.; Dunford, B.; Pickard, M. A. Kinetics of the oxidation of ascorbic acid, ferrocyanide and p-phenolsulfonic acid by chloroperoxidase compounds I and II. Eur. J. Biochem. 1987, 163, 123–127. 3109


Numerical Integration of the Chemical Rate Equations via a

Recommend Documents