Article pubs.acs.org/jchemeduc
Free Radical Addition Polymerization Kinetics without Steady-State Approximations: A Numerical Analysis for the Polymer, Physical, or Advanced Organic Chemistry Course H. Darrell Iler,*,† Amber Brown,† Amanda Landis,†,§ Greg Schimke,† and George Peters‡ †
Department of Chemistry and ‡Department of Mathematics, Greenville College, Greenville, Illinois 62246, United States S Supporting Information *
ABSTRACT: A numerical analysis of the free radical addition polymerization system is described that provides those teaching polymer, physical, or advanced organic chemistry courses the opportunity to introduce students to numerical methods in the context of a simple but mathematically stiff chemical kinetic system. Numerical analysis can lead students to an appreciation of the power of numerical methods applied to chemical systems while providing a better understanding of the steady-state approximation and giving a more complete picture and deeper understanding of reaction dynamics. The description and Mathcad and MATLAB tutorials provided in the Supporting Information offer a variety of levels of student engagement ranging from classroom presentations of numerical and steady-state comparisons to projects requiring students to set up, run, and interpret numerical analyses for a free radical addition polymerization system. KEYWORDS: Upper-Division Undergraduate, Polymer Chemistry, Physical Chemistry, Organic Chemistry, Computer-Based Learning, Kinetics
F
The justification for applying the SSA is based upon the vastly different rates of formation and destruction for the intermediate species as dictated by the system rate constants. For example, the rate constant associated with the destruction of chain radicals (k4) can be 8−10 orders of magnitude greater than the rate constant associated with the formation of primary radicals (k1). These conditions suggest that after a transient or induction period, the system should settle into a state characterized by a steady and nearly zero level of radical concentrations. For more than 50 years, investigators have been remarkably successful in using the SSA to help understand and predict properties of many free radical and other addition polymerization systems. However, there is a cost. The approximation eliminates both radicals (R• and M•) and a rate constant (k2) from the resulting steady-state equations, thereby giving an incomplete representation of the reaction dynamics. For example, the impact of the radicals on the development toward and establishment of steady state cannot be directly determined through this approximation. It should be noted that the primary radical (R•) undergoes other termination reactions with primary radicals and chain radicals that are not accounted for in eq 2. These reactions are diffusion controlled and very fast and, together with the k2[R•][M] term, result in R• rapidly achieving a steady-state concentration. Also, the k4 rate constant has a chain-length dependence that is not accounted for in eq 4. These factors are
ree radical addition polymerization (FRAP) is an important and widely used technique for producing polymers. The accepted mechanism for this process, as shown in Scheme 1, results in a set of four simultaneous rate equations.1−5 An analytical solution to this system does not exist. However, application of the steady-state approximation (SSA) allows one to set eqs 2 and 4 equal to zero resulting in constant concentration terms for the two radical species, R• and M•. These simplifications lead to equations for the time dependence of the monomer concentration as well as other important polymerization and polymer properties, such as the rate of polymerization, conversion (fraction of monomer converted to polymer), average radical lifetime, and average degree of polymerization. There are many published examples in this Journal of pedagogical applications of numerical simulations leading to validation of the SSA for relatively simple to complex kinetic systems.6−11 However, no such reports have been found with respect to the FRAP kinetic model. Furthermore, an inspection of popular physical chemistry and polymer chemistry textbooks3,5,12−15 reveals that their treatments of the free radical addition polymerization system are all limited to the SSA. In most cases, the approximation is applied and results are calculated without rigorous justification or verification. Therefore, there is precedent and need for a pedagogical report and student activities relating to the numerical analysis of the kinetics of this polymerization technique that is central to polymer education and polymer science in general. © 2014 American Chemical Society and Division of Chemical Education, Inc.
Published: January 27, 2014 374
dx.doi.org/10.1021/ed4005535 | J. Chem. Educ. 2014, 91, 374−379
Journal of Chemical Education
Article
Scheme 1. Free Radical Addition Polymerization Mechanism and Rate Equations
through the traditional steady-state approach. This affords a unique opportunity for advancing students’ knowledge and experience with the chemistry of free radical addition polymerization systems. Most undergraduate chemistry students will not have the background or desire to learn and understand the mathematical details of stiff solver algorithms. However, because the objective is for students to use numerical results to gain a deeper insight into the underlying chemistry, these details are not required.22 The typical polymer or physical chemistry student, with one or two semesters of calculus, is adequately equipped to follow an elementary description of the numerical methods and understand and interpret the results in terms of the FRAP system chemistry. In the polymer chemistry classroom, we have used these simulations to present data and plots comparing the numerical and steady-state results for a FRAP system and then discuss the validity of the SSA for this important and widely used polymerization technique. Other system aspects that are not accessible through a steady-state analysis are also addressed, such as the induction time, the dynamic role of radicals, and the k2 rate constant. For chemistry programs in which Mathcad or MATLAB is a part of the curriculum, students can easily use a set of instructions to set up, run, and interpret the numerical analysis of a FRAP system. We have used this as a physical chemistry project assignment in which students worked in groups of two or three and presented their results and interpretations to the class. There are several good reasons for introducing polymer and physical chemistry students to numerical analysis through the free radical addition polymerization model: • It reminds students that the steady-state analysis, like any approximation method, should be verified by independent means whenever possible. • It maintains the active presence and separate identity of the participating radicals, thereby giving a more complete picture of the reaction dynamics and promoting a deeper insight into the underlying mechanism.22,23 Basically, students use numerical analysis to learn chemistry.
beyond the scope of this paper. The authors have simplified the model by only considering one mode of consumption for the primary radicals (i.e., their reaction with monomer) and a rate constant for termination that is chain-length independent. The same kinetic features of the FRAP system, discussed above, that justify application of the SSA also present a special difficulty with respect to its analysis through numerical methods. That is, the broad range of rate constant values and a near zero and steady set of intermediate concentrations lead to a property of numerical systems known as stiffness instability.16−18 The stiffness of a system of equations is indicated by the complete failure of the computer program or the requirement of large amounts of computer time to converge to a solution when using standard numerical methods such as Runge−Kutta. Early investigators recognized that stiffness instability is a common feature for chemical kinetic systems and, therefore, had to rely on the steady-state approach for their analysis of those systems.19 However, the relatively recent development of powerful mathematical software programs that incorporate stiff solver algorithms can make generating numerical solutions to many complex kinetic problems a reasonable task, even for undergraduate chemistry students.6,20,21 This paper describes the Numerical Analysis of the full mechanistic form of the FRAP model using two widely available mathematical software products (Mathcad and MATLAB). Both products incorporate Runge−Kutta and stiff solver algorithms as part of their standard packages that provide essentially identical solutions for the problems described in this paper. Our quantitative and graphical comparisons of the numerical and steady-state results for the polymerization of polymethyl methacrylate (PMMA) validate the use of the SSA for this system. By fitting numerical solutions to experimental data for the polymerization of PMMA, k3 and k4 values and a k2 threshold value were found that give numerical solutions that are in excellent agreement with those from steady-state calculations and empirical data over a broad range of initial monomer concentrations. In addition, the numerical methods reveal a host of interesting mechanistic details not available 375
dx.doi.org/10.1021/ed4005535 | J. Chem. Educ. 2014, 91, 374−379
Journal of Chemical Education
■
Article
• It is simple. Many kinetic systems contain dozens or even hundreds of elementary reactions17 for which rate equations and constants would have to be determined and entered. This would be a daunting task for both students and instructors. The FRAP system, in its simplest form, has only four reactions whose rate constants have been experimentally determined. This results in a system whose numerical analysis is easy to set up, runs quickly (typically less than a minute), and lends itself to straightforward and simple interpretation. Thus, it provides a good undergraduate level introduction to the stiff systems that are so common in chemical kinetics. • It relates to familiar items. Students generally enjoy gaining a deeper technical understanding and knowledge of items and processes they encounter on a regular basis.24 Many familiar polymers, such as the polyethylene in plastic milk containers, the poly(vinyl chloride) in PVC pipe, and the PMMA in automobile headlight lenses, are produced using FRAP. Informing students of this and bringing example materials and products into the classroom can help make the concepts more relevant and engaging for them.
Table 1. Parameter Values Used for Numerical Simulations of PMMA Polymerization at 77 °C parameter
value
k1 k2 k3 k4 f
9.50 × 10−5/(1/s)a adjustable, but assumed to be ≥ k3 980 L/(mol s)b 2.80 × 107 L/(mol s)b 0.6c
a c
Value from ref 27. bValue from the best-fit results to empirical data. Value from ref 2.
that k3[M•] ≫ k2[R•] in eq 3. The latter assumption leads to −d[M]/dt = k3[M•][M], which with the [M•] substitution gives the steady-state rate of polymerization as seen in eq 8. Note that the k3[M•] ≫ k2[R•] assumption is not used in the numerical analysis. Also note that [M] and [M]o represent the time dependent and initial concentrations of monomer, respectively. Equations 11 and 12 give the time dependent [M] derived from steady-state assumptions and time dependent [I], respectively, both of which are used in eqs 8−10. Rpn = k 3[M][M•] + k 2[M][R•]
NUMERICAL ANALYSIS The goal of the analysis is to provide full numerical integration of the set of four simultaneous differential equations from Scheme 1 using algorithms that are available in commercial mathematical software products such as Mathcad and MATLAB. The polymerization of polymethyl methacrylate (PMMA) at 77 °C with the initiator azo-bis-isobutyrolnitrile (AIBN) was chosen as the test system for the simulations. However, other polymers, initiators, and temperatures work equally well (e.g., benzoyl peroxide-initiated polystyrene). Values for k1, k3, and k4 as a function of reaction temperature can be found in a variety of sources including the Polymer Handbook.25 The k1 constant is associated with the simple firstorder decay of the initiator and has been measured directly and independently of the polymerization reaction. The k3 and k4 constants were historically determined indirectly by combining steady-state assumptions with experimental measurements of rates of polymerization and chain radial lifetime.3 Due to the difficulty of measuring k3 and k4, the literature shows a relatively large range of values for each. Because the SSA equations do not incorporate k2, these values are not as broadly available. However, because the primary radical is usually more reactive than the chain radical toward monomer, it is assumed that k2 will have a value close to or greater than that of k3. Literature sources present k2 as being equal to or about 1 order of magnitude greater than k3.2,26 For this paper, k2 is used as an adjustable parameter in the simulations. The parameter values used for most simulations are listed in Table 1. The polymerization and polymer properties used for comparing the numerical and steady-state results in this paper and supplementary tutorials are the rate of polymerization (−d[M]/dt = Rp), conversion (X), and kinetic chain length (ν). The functional form for these properties, derived from numerical (Rpn, Xn, νn) and steady-state (Rps, Xs, vs) approaches, are given by eqs 5−10. Note that the numerical equations, eqs 5−7, retain the radical concentration terms that are absent in the steady-state equations. The steady-state equations, eqs 8−10, are derived from setting eqs 2 and 4 equal to zero to find expressions for [M•] and [R•] and by assuming
[M] [M]o
(6)
k 3[M][M•] k 2[R•][M]
(7)
Xn = 1 − vn =
⎛ fk [I] ⎞0.5 Rps = k 3[M]⎜ 1 ⎟ ⎝ k4 ⎠ Xs = 1 − e νs =
(5)
⎛ fk [I] ⎞ −k 3t ⎜ 1 ⎟ ⎝ k4 ⎠
(8)
0.5
(9)
k 3[M] 0.5
( )
2k 2
fk1[I] k4
(10) 0.5
( ) [M] = [M] fk1[I]o k4
o
e((k1t /2) − 1)
k1
[I] = [I]o e−k1t
(11) (12)
Runge−Kutta and Stiff Systems
The fourth-order Runge−Kutta method is widely used to provide numerical solutions for chemical kinetic systems involving ordinary differential equations7. Mathcad has two explicit versions of this method, RKfixed and RKadapt. These are explicit methods for which the solution value for each time step of the analysis is calculated in terms of the value of the previous time step. The RKfixed algorithm uses a fixed step size and fails to provide a reasonable solution for the PMMA FRAP system. At any k2 value above 10 L/(mol s), which is 2 orders of magnitude or more below its expected value, the system exhibits a rapidly increasing oscillatory behavior that will crash (fail) within the first few steps of the analysis. For example, within one step, the algorithm output for the primary radical concentration can go from of 10−7 to 1064 mol/L. This instability is characteristic of a stiff system when nonstiff 376
dx.doi.org/10.1021/ed4005535 | J. Chem. Educ. 2014, 91, 374−379
Journal of Chemical Education
Article
Stiff Solvers and Finding the Threshold Value for k2
methods are applied. For unrealistically small values of k2 (where k1 and k2 are more comparable in magnitude) the stiffness problem is apparently alleviated, and RKfixed provides a solution. However, under this condition, as can be seen in Figure 1, the numerical and steady-state solutions are not in agreement.
Stiff solver algorithms from Mathcad (the Bulirsch−Stoer method) and MATLAB (a modified Rosenbrock method) were applied to the FRAP problem. No mathematical details for these algorithms will be given here, but the Mathcad and MATLAB tutorial files in the Supporting Information provide examples and instructions for setting up the required matrices and equations for these algorithms. Both MATLAB and Mathcad stiff solvers give rapid and essentially identical solutions that agree with those from steady state as long as k2 exceeds a threshold value of about 1−2 L/(mol s). For k2 values above this threshold, the numerical analysis results are independent of k2. Figure 2 illustrates the quality of agreement through the plots of the numerical and steady-state results for the PMMA system with k2 set at 1000 L/(mol s). Quantitative comparisons show that after a short induction period of a few seconds the system properties (rate of polymerization, conversion, etc.) predicted numerically and by steady state are essentially identical across the range of reaction times examined. These results indicate that there exists a k2 threshold value below which the system will not exhibit steady-state behavior and above which the system becomes independent of k2. Because the rate of destruction of the primary radicals is proportional to k2 through reaction 2, it is reasonable to assume that k2 must be much larger than the k1 term associated with the formation of those same radicals via reaction 1. In the PMMA simulation, the k2 is required to be at least 4 orders of magnitude greater than k1. Figure 3 shows the difference of the numerical results compared to the steady-state results plotted versus k2. The difference reaches a minimum on the order of 10−8 L/(mol s) (compared to rates of polymerization on the order of 10−5 L/(mol s)) at a k2 value between 1 and 2 L/(mol s) and remains essentially constant thereafter. This k2 threshold value is independent of the initial concentrations of monomer and initiator used in the simulations.
Figure 1. Numerical and steady-state results for conversion plotted versus reaction time. The k2 value is set at 0.1 L/(mol s) for the numerical analysis using a fixed step-size Runge− Kutta algorithm.
RKadapt is a Runge−Kutta method with an adaptive step size of integration. The step size decreases when the solution changes quickly and increases when it changes slowly. This kind of step-size variation can reduce numerical instability but can also produce ridiculously small step sizes that greatly extend computation time.18 When applied to the FRAP problem, the RKadapt algorithm accommodates the larger k2 values, values for which RKfixed failed, but runs so slowly that it essentially does not allow one to obtain a numerical solution. The instability of the RKfixed and extensive run time of the RKadapt analyses indicate that the FRAP system is indeed stiff and will require a numerical method intended for stiff problems.
Induction Time
An obvious and important question that arises with the assumption of steady state is, “What is the time for the system to achieve this condition?” By assuming that steady state is reached when the rates of initiation and termination become equal, the development toward and achievement of steady state
Figure 2. Numerical and steady-state results for conversion (left) and rates of polymerization (right) plotted versus reaction time. The numerical results are generated using a stiff solver algorithm. 377
dx.doi.org/10.1021/ed4005535 | J. Chem. Educ. 2014, 91, 374−379
Journal of Chemical Education
Article
This time is independent of the initial concentrations of monomer and initiator used in the simulations. Two other methods for approximating the induction time for this simulation, based upon steady-state assumptions, give results that are on the same order of magnitude, 9 s28 and 3.1 s.3 The induction time is less than 0.03% of the monomer half-life. R i = k 2[M][R•]
(13)
R t = 2k4[M•]2
(14)
Fitting Numerical Results to Empirical Data
Because of the large number and underdetermined nature of the rate constants associated with many chemical kinetic systems, it is generally not recommended that numerical results be fit to experimental data for the purpose of determining those rate constants.17 However, the FRAP system is an exception to this rule. There are only four rate constants, and one of those (k1) is accurately known through direct measurements of initiator decomposition reactions. Also, as long as k2 is set above its threshold value, it will have no effect on the numerical results. That leaves only two rate constants (k3 and k4), which have typically been determined indirectly using steady-state assumptions. Numerical results were fit to PMMA initial rate of polymerization data run at 77 °C27 using a least-squares analysis code written in MATLAB. Six different sets of initial conditions of monomer and initiator concentrations from the literature were entered into the code with k1 and k2 set at 9.50 × 10−5 s−1 and 1000 L/(mol s), respectively. The k3 and k4 terms were used as fitting parameters with ranges for their values determined from literature.25 For each k3−k4 pair, the differential equations were solved for all six sets of initial conditions at a reaction time of 100 s and were considered to be the estimated initial rates. A least-squares analysis between the estimated and observed initial rates of polymerization was used to determine the minimal error. The k3 and k4 values that produced the best agreement with the observed data (minimal sum of the squares of the differences) are shown in Table 1. There is a line segment in the k3−k4 plane of solutions where the values produce nearly identically small errors. Table 2 lists the initial rate of polymerization values from the numerical, steady-state, and empirical studies for ten different sets of monomer and initiator concentrations.
Figure 3. Plots of the difference between rates of polymerization values for numerical and steady-state calculations versus k2 at reaction times of 1000 s and 10 000 s. The “threshold” refers to that value of k2 below which the system does not exhibit steady-state behavior.
can be graphically displayed using the numerical solutions. Figure 4 shows the plots of the numerical rates of initiation (Ri)
Figure 4. Numerical results for the rates of initiation and termination plotted against reaction time.
and termination (Rt). Both plots level off and coincide between 2 and 3 s into the reaction. This defines the time required for the system to reach steady-state behavior (the induction time).
Table 2. Initial Rates of Polymerization, Steady-State (SSA), and Numerical Analysis for Different Sets of Monomer [M] and Initiator [I] Concentrations
a
[M]/(mol/L)
[I] × 104/(mol/L)
Rp × 104/[(mol/L)/s] (empirical)a
Rp × 104/[(mol/L)/s] (steady state)
Rp × 104/[(mol/L)/s] (numerical)
2.07 3.26 4.17 4.22 4.75 4.96 6.13 7.19 8.63 9.04
2.11 2.45 5.81 2.30 1.92 3.13 2.28 2.55 2.06 2.35
0.415 0.715 1.30 0.857 9.37 1.22 1.29 1.65 1.70 1.93
0.422 0.715 1.41 0.897 9.23 1.23 1.30 1.61 1.74 1.94
0.422 0.716 1.41 0.898 9.23 1.23 1.30 1.61 1.74 1.94
Values from ref 27. 378
dx.doi.org/10.1021/ed4005535 | J. Chem. Educ. 2014, 91, 374−379
Journal of Chemical Education
■
Article
(8) Carresana García, J. A. ABC Kinetics. J. Chem. Educ. 2003, 80 (10), 1220. (9) Chemical Education Digital Library. www.chemeddl.org/ alfresco/service/org/chemeddl/symmath/apps?toc_id=13&guest= true (accessed Jan 2014). (10) Porter, M. D.; Skinner, G. B. The Steady-State Approximation in Free-Radical Calculations. A Numerical Example. J. Chem. Educ. 1976, 53, 366−367. (11) Bond, R. A.; Simoyl, R. H. The Quasi-Steady-State Approximation: Numerical Validation. J. Chem. Educ. 1998, 75, 1158−1165. (12) Odian, G. Principles of Polymerization; McGraw-Hill Companies, Inc.: Hoboken, NJ, 1970; p 161. (13) Levine, I. Physical Chemistry, 6th ed.; McGraw-Hill Companies, Inc.: New York, NY, 2009; p 554. (14) Cowie, J. Polymers: Chemistry and Physics of Modern Materials, 2nd ed.; Chapman and Hall: London, 1991; p 54. (15) Allcock, H.; Lampe, F. Contemporary Polymer Chemistry; Prentice-Hall Inc.: Englewood Cliffs, NJ, 1981; p 271. (16) Edelson, D. The New Look in Chemical Kinetics. J. Chem. Educ. 1975, 52, 642−644. (17) Hirst, D. A Computational Approach to Chemistry; Blackwell Scientific Publications: Oxford, United Kingdom, 1990; p 165. (18) Shampine, L. Numerical Solution of Oridnary Differential Equations; Chapman & Hall, Inc.: New York, NY, 1994; p 382. (19) Snow, R. A Chemical Kinetics Computer Program for Homogenous and Free-Radical Systems of Reactions. J. Phys. Chem. 1966, 70, 2780−2786. (20) Chesick, J. Interactive Program System for Integration for Reaction Rate Equations. J. Chem. Educ. 1988, 65, 599−602. (21) Goodman, J. How Do Approximations Affect the Solutions to Kinetic Equations? J. Chem. Educ. 1999, 76, 275−277. (22) Edelson, D. The Steady State Approximation: Fact or Fiction? J. Chem. Kinetics 1974, 6, 787. (23) Zielinski, T. Critical Thinking in Chemistry Using Symbolic Mathematics Documents. J. Chem. Educ. 2004, 81, 1553. (24) Iler, H.; Rutt, E.; Althoff, S. An Introduction to Polymer Processing, Morphology, and Property Relationships through Thermal Analysis of Plastic PET Bottles. Exercises Designed to Introduce Students to Polymer Physical Properties. J. Chem. Educ. 2006, 83, 439−442. (25) Polymer Handbook; Brandrup, J., Immergut, E. H., Grulke, E. A., Eds.; John Wiley & Sons, Inc.: Hoboken, NJ, 1999. (26) Beuermann, S.; Buback, M. Rate Coefficients of Free Radical Polymeriation Deduced from Pulsed Laser Experiments. Prog. Polym. Sci. 2002, 27, 191−254. (27) Arnett, L. Kinetics of the Polymerization of Methyl Methacrylate with Aliphatic Azobisnitriles as Initiators. J. Am. Chem. Soc. 1952, 74, 2027−2033. (28) Nozaki, K.; Bartlett, P. Rate Constants of the Steps in Addition Polymerization. I. The Induction Period in the Polymerization of Vinyl Acetate. J. Am. Chem. Soc. 1946, 68, 2377−2380.
CONCLUSIONS This work leads to the conclusion that the free radical addition polymerization mechanism results in a stiff system of coupled differential equations whose numerical analysis requires algorithms designed to accommodate stiffness. The commercially available mathematical software products, MATLAB and Mathcad, incorporate stiff solver algorithms that give essentially identical solutions, for this system, that are in excellent agreement with steady-state predictions. The numerical analysis leads to the determination of the induction time, the threshold value for the k2 rate constant, and a set of optimal k3 and k4 rate constants for the polymerization of PMMA at 77 °C. These parameters provide numerical rates of polymerization that are in excellent agreement with steady-state results and empirical data. This work has shown that presenting the numerical analysis of the free radical addition polymerization system to students makes good pedagogical sense. It offers the opportunity to introduce polymer chemistry or physical chemistry students to numerical methods in the context of a simple but very important chemical kinetic system. Numerical analysis can help students better understand free radical addition polymerization and give them a better understanding of the often valid and useful SSA. The paper, along with provided Mathcad and MATLAB tutorials, offer a variety of levels of student engagement ranging from classroom presentations of numerical and steady-state comparisons to projects requiring students to set up, run, and interpret numerical analysis for a free radical addition polymerization system.
■
ASSOCIATED CONTENT
S Supporting Information *
MATLAB and Mathcad files with basic coding instructions for setting up and running the numerical analysis of a FRAP system using stiff solvers. This material is available via the Internet at http://pubs.acs.org.
■
AUTHOR INFORMATION
Corresponding Author
*H. D. Iler. E-mail:
[email protected]. Present Address §
Department of Chemistry, University of Wyoming, Laramie, Wyoming 82071, United States. Notes
The authors declare no competing financial interest.
■
REFERENCES
(1) Flory, P. The Mechanism of Vinyl Polymerizations. J. Am. Chem. Soc. 1937, 59, 241−253. (2) Flory, P. Principles of Polymer Chemistry; Cornell University Press: Ithaca, NewYork, 1953; p 106. (3) Hiemenz, P. Polymer Chemistry: The Basic Concepts; Marcel Dekker, Inc.: New York, NY, 1984; p 345. (4) McGrath, J. Chain Reaction Polymerization. J. Chem. Educ. 1981, 58, 844−861. (5) Atkins, P. Physical Chemistry, 6th ed.; W. H. Freeman and Company: New York, 1998; p 803. (6) Pavlis, R. Kinetics Without Steady State Approximations. J. Chem. Educ. 1997, 74, 1139−1140. (7) Francl, M. M. Exploring Exotic Kinetics: An Introduction to the Use of Numerical Methods in Chemical Kinetics. J. Chem. Educ. 2004, 81 (10), 1535. 379
dx.doi.org/10.1021/ed4005535 | J. Chem. Educ. 2014, 91, 374−379
