Research Article  Open Access
Matrix Exponentiation and the FrankKamenetskii Equation
Abstract
Long time solutions to the FrankKamenetskii partial differential equation modelling a thermal explosion in a vessel are obtained using matrix exponentiation. Spatial derivatives are approximated by highorder finite difference approximations. A forward difference approximation to the time derivative leads to a LawsonEuler scheme. Computations performed with a BDF approximation to the time derivative and a fourthorder RungeKutta approximation to the time derivative are compared to results obtained with the LawsonEuler scheme. Variation in the central temperature of the vessel corresponding to changes in the shape parameter and FrankKamenetskii parameter are computed and discussed.
1. Introduction
In this paper we use matrix exponentiation to determine numerical solutions of the FrankKamenetskii partial differential equation (FKPDE) modelling a thermal explosion in a vessel. The FKPDE is given by the nonlinear partial differential equation: where is the nondimensional temperature in a vessel, is a coefficient that represents the shape of the vessel, and is the FrankKamenetskii parameter [1]. The FKPDE (1.1) is solved subject to the boundary conditionsand the initial condition
In this paper the spatial derivatives of the FKPDE (1.1) are approximated by highorder finite difference schemes. The resulting system of equations is approximated using the method of lines and is then integrated using matrix exponentiation. The resulting matrix approximation is pentadiagonal because we have approximated the spatial derivatives using highorder finite difference schemes. Comparisons are made between two highorder approximations to the time derivative by comparing a BDF approximation and a RungeKutta approximation to the time derivative. The results obtained compare favorably with the results obtained by Harley [2] and Britz et al. [3].
The FKPDE (1.1) is derived from a heat balance equation given by where is the thermal capacity, is the thermal conductivity, the density, is the heat of reaction, is the frequency factor, is the energy of activation of the chemical reaction, is the universal gas constant, and is the gas temperature. The heat balance equation (1.4) is made nondimensional by substituting where is the ambient temperature into (1.4). The Laplacian is given by where the constant is related to the geometry of the problem.
We make a change of variableswhere is a dimensionless temperature at the vessel walls and The heat balance equation (1.4) reduces to the FKPDE equation (1.1) where boundary conditions for a fixed constant temperature at the cylinder wall are given bywhere is the boundary of vessel and . We choose to obtain the boundary conditions (1.2a) and (1.2b). For welldefined geometries we have for a rectangular slab, for an infinite circular cylinder, and for a sphere. Balakrishnan et al. [4] have solved the steady problem for nonstandard geometries, that is, is not an integer. The problem (1.1) has been solved using the Robin boundary conditions at the boundary given by [5–8]: where is a constant.
FrankKamenetskii [1] has derived (1.1) by ignoring coefficients from the heat balance equation (1.4). The motivation for deriving (1.1) comes from the study of heat generated in chemical reactions that follow an Arrhenius law. When the heat generated by a chemical reaction is far greater than the heat lost to the walls of the vessel in which the reaction is taking place, a thermal explosion occurs. Experimental work on thermal explosions has been conducted by Rice [9] and Rice and Campbell [10]. Steggerda [11] has extended the work done by FrankKamenetskii [1] on the criteria under which a thermal explosion can take place. In the models presented by Zhang et al. [8, 12] and Abd ElSalam and Shehata [13] the fuel used in the chemical reaction is consumed during the reaction.
Harley [2] has compared hopscotch, method of lines, and CrankNicolson schemes to solve the nonlinear partial differential equation (1.1). Britz et al. [3] improve on the work of Harley [2] by treating the nonlinear source term in (1.1) in both an implicit and explicit form. Britz et al. [3] include the effects of the critical activation parameter which we ignore. In this paper we use matrix exponentiation to improve on the work of Harley [2] and Britz et al. [3]. The nonlinear source term in (1.1) is treated explicitly in a natural way. We derive numerical schemes that are accurate in both time and space by considering a highorder finite difference approximations to the spatial derivatives coupled with BDF and RungeKutta approximations to the time derivative.
We reduce (1.1) to the semilinear system of equations where is a linear operator and is a nonlinear function of . The firstorder semilinear equation (1.11) can be solved by multiplication with an integrating factor. The resulting solution contains a matrix exponential that needs to be evaluated. The interested reader is referred to the reviews of Vanden Bosch et al. [14], Moler and Van Loan [15], and Hochbruck and Ostermann [16] on the construction of matrix exponentials for semilinear systems of the form (1.11). Recent developments have focused on constructing matrix exponentials using Krylov subspace approximations [17] and polar decompositions [18]. In this paper we use MATHEMATICA to determine the relevant matrix exponentials. Alternate computer packages to determine matrix exponentials have been developed by Sidje [19] and Berland et al. [20].
A backward differentiation formula (BDF) approximation to the time derivative leads to a numerical scheme with a truncation error of where is the spatial step and is the time step. A RungeKutta approximation to the time gives a scheme with a truncation error . A BDF approximation to the time derivative leads to Astable multistep methods [21]. A fourthorder RungeKutta approximation to the time derivative is computationally expensive because it requires four additional computations per time step. The interested reader is referred to Gottlieb et al. [22] and Kassam and Trefethen [23] for discussions about highorder approximations to the time derivate and their respective properties.
The paper is divided up as follows. In Section 2 we consider approximate solutions for the case in terms of the heat kernel to investigate the influence of the parameter . In Section 3 we derive three numerical schemes based on matrix exponentiation to determine numerical solutions admitted by (1.1). Concluding remarks are made in Section 4.
2. Analytical Solution
In this section we consider an approximate analytical solution admitted by (1.1) for the case in terms of the heat kernel. The heat kernel solution is a good approximation to the behaviour of the solutions of (1.1) as the solution satisfies the boundary conditions (1.9a) and (1.9b) for . It will be useful to use the heat kernel solution to examine the effects of the FrankKamenetskii parameter on the temperature at the centre of the vessel. Momoniat [24] has used the Lie symmetry approach to determine group invariant solutions admitted by (1.1) for . The solutions obtained by Momoniat [24] are valid after blowup.
In this section we determine approximate solutions admitted by (1.1) by using the approximation . As discussed in Harley [2] this approximation is valid close to the centre of the cylindrical vessel and yields results that give an error when compared to the numerical solution at the steady state. The model equation simplifies to We make the assumption and investigate approximate solutions of the form Substituting (2.2) into (2.1) and separating to first order in we obtain the systemThe heat kernel solution to is given by By inspection, we find that admits the particular solution where is a constant. The heat kernel solution (2.4) is also a particular solution as is any solution of the linear phenomenological diffusion equation (2.3a). The linear nature of (2.3b) implies that a linear combination of these solutions is also a solution of (2.3b). We can write these solutions as
It is immediately obvious from and that the main effect of the FrankKamenetskii parameter is to increase the temperature at the center of the vessel when compared to the temperature for the phenomenological diffusion equation . We can thus expect that an increase in the critical parameter will yield an increase in the central temperature irrespective of the shape of the vessel. We plot the solutions (2.6) in Figure 1. The curves in Figure 1 confirm the effect of to increase the central temperature of the vessel. Positive values of the constant also cause a corresponding increase in the core temperature of the vessel. In the next section we consider numerical solutions admitted by (1.1) obtained from matrix exponentiation.
(a)
(b)
3. Derivation and Implementation of the Numerical Scheme
In this section we compare three numerical schemes for solving the FrankKamenetskii partial differential equations (1.1) based on matrix exponentiation. In the first instance we approximate the time derivative by a forward difference approximation. We then approximate the time derivative by a backward difference formula. We finally use a RungeKutta approximation to the time derivative.
Using L’Hospitals rule at we approximate by to write (1.1) as We approximate the spatial derivatives by the central difference approximations where . Highorder central difference approximations to the derivatives are given by The derivative boundary condition (1.2a) gives the condition when . Defining and where we write (1.1) in vector form as where The matrix is different to the matrix obtained by Harley [2] in that we have used highorder finite difference approximations to the spatial derivatives to improve the accuracy of the numerical scheme. The matrix A is pentadiagonal and not tridiagonal as obtained by Harley [2]. Loworder approximations are used in row positions and because no information is known about the temperature at the grid points that occur at two steps on either side of the boundaries. The entries in row positions and are determined from the boundary conditions (1.2a) and (1.2b).
We multiply (3.4) by the integrating factor to obtain where is used to signify the multiplication of a matrix with a vector. Using a forward difference approximation to the time derivative given by we obtain the LawsonEuler scheme The truncation error of results obtained from the LawsonEuler scheme (3.8) is .
We implement the numerical scheme (3.8) in MATHEMATICA. We assume that a steady state has been reached when where is a tolerance. The quantity . When the tolerance condition (3.9) has been satisfied we set , the temperature at the centre of the vessel. The boundary condition (1.2b) is imposed as at each iteration.
An improvement on the accuracy of the LawsonEuler scheme (3.8) can be made by considering a backward difference formula (BDF) to the time derivative given by [21] We then obtain the scheme
The numerical scheme (3.11) is implemented in MATHEMATICA. The scheme (3.11) requires two starting values which we obtain from the initial value and which is obtained from the LawsonEuler scheme. The scheme terminates when condition (3.9) is satisfied. The truncation error of results obtained from the BDF formulation of the time derivative is .
We lastly consider a fourthorder RungeKutta scheme to solve (3.6). The resulting numerical scheme is given by where
Once again the scheme is implemented in MATHEMATICA. Unlike the scheme (3.11) the numerical scheme (3.12) does not require two starting values. The scheme does however require four evaluations at each time step. The truncation error of results obtained using the numerical scheme (3.12) is .
In Tables 1 and 2 we compare output from numerical schemes (3.8), (3.11), and (3.12) for different values of and . In Table 1 we fix the value of and vary , while in Table 2 we fix the value of and vary . We note from the results in Tables 1 and 2 that for small values of and the difference in the output from the numerical schemes is . For larger values of and this difference increases to . Due to this large discrepancy that occurs in the outputs from the different numerical schemes for large values of and we implement the fourthorder RungeKutta numerical scheme (3.12) in the rest of the paper. The numerical scheme (3.12) based on a fourthorder RungeKutta scheme that gives results with the smallest truncation error, that is, and are hence the most accurate.


To test the stability of each of the numerical schemes we consider the steady solution for the case . For the case the model equation (1.1) must satisfy the steady solution where satisfies the secondorder ordinary differential equation and the boundary conditions
FrankKamenetskii [1] (see also Chambré [25], Dresner [26], and Momoniat and Harley [27]) have shown that (3.14) solved subject to (3.15) admits the solutions where and are constants of integration given by These multiple solutions lead to the wellknown bifurcation that occurs when plotting the temperature at the centre of the vessel against the critical parameter .
We consider the solution that has a maximum . This corresponds to the results given in Table 3 for the central temperature at a tolerance of . We label the numerical scheme (3.8) corresponding to the LawsonEuler scheme as the LE scheme, the numerical scheme (3.11) corresponding to the BDF approximation to the time derivative as the BDF scheme, and the numerical scheme (3.12) corresponding to a RungeKutta approximation to the time integration as the RK4 scheme. We use (3.18) as the initial value for the numerical schemes and consider 10, 100, and 1000 iterations, respectively.

In Figure 2 we plot for , and 1000, respectively, where we have chosen , , , and . We note from the results plotted in Figure 2 that the LawsonEuler scheme (3.8) and the RungeKutta scheme (3.12) produce the most stable results. The results corresponding to the BDF scheme (3.11) are not as stable. We note that all three schemes produce results that are accurate to more than 12 decimal places after 1000 iterations.
(a)
(b)
(c)
4. Results
The results in this section have been determined using the numerical scheme (3.12) corresponding to the fourthorder RungeKutta approximation for the time integration. As indicated in the previous section the RK4 scheme produces results that are stable and very accurate because it is a highorder scheme. An important aspect of the problem is the effect of the tolerance on the time to reach the steady state and the corresponding temperature at the centre of the vessel. Choosing and for and we evaluate the numerical scheme (3.12) for different tolerances. In Table 3 we note the stopping time and corresponding temperature at the centre of the vessel, , for each value of the tolerance.
From the results in Table 3 we note that the temperature converges to but the stopping time increases from to . This big increase in the stopping time reflects the large number of iterations required for the stopping criteria to be satisfied. What is important in interpreting these results is the scaling of the temperature. The tolerance, , is nondimensional, but it has the same order of behaviour as the temperature. To understand the effect of the tolerance on the temperature we transform back to the original variables using (1.7a), (1.7b), and (1.5) to obtain where is the ambient temperature, is the dimensionless temperature at the vessel wall, is the universal gas constant, and is the energy of the activation of the chemical reaction. Therefore the tolerance can be physically interpreted as a change in the temperature at the wall of the vessel and the temperature at the centre of the vessel corresponding to the wall temperature . As the temperature at the centre of the vessel stops having to adjust; hence we get convergence as is indicated in Table 3.
The large variation in time shows that using (1.1) to determine the time of a thermal explosion is not appropriate. It is more practical to monitor the temperature at the centre of the vessel or at the wall of the vessel. Once the critical temperature is reached a thermal explosion takes place. In the results for the rest of the paper we have used the tolerance for the steady state.
In Figure 3 we plot the temperature at the centre of the vessel against the shape factor for fixed . We note that for increasing the value of the shape factor the steady temperature decreases for the tolerance .
In Figure 4 we plot the temperature at the centre of the vessel against the increasing value of the FrankKamenetskii parameter for fixed . We note that for the increasing value of the steady temperature increases for the tolerance . These results match up with the approximate analytical results predicted by the study of the approximate heat kernel solution in Section 2.
5. Concluding Remarks
In this paper we have shown how matrix exponentiation can be used to solve a nonlinear partial differential equation modelling a thermal explosion in a vessel. Matrix exponentiation has the advantage over a standard forwardtime central space (FTCS) scheme in that the scheme inherently contains higherorder timestep terms. This can be seen if we consider the approximation we find that the LawsonEuler scheme (3.8) simplifies the FTCS scheme, to . The advantage of the approach taken in this paper is that we have used matrix exponentiation in combination with highorder spatial and highorder time discretizations. Future work will investigate the use of multistep exponential integrators of the type introduced by Calvo and Palencia [28] and Ostermann et al. [29] to improve the accuracy of our results.
We have reduced the nonlinear partial differential equation (1.1) to a semilinear form. The resulting semilinear equation is reduced to an integrable form by the multiplication of an integration factor. We approximate the time derivative by a forward difference approximation to obtain a LawsonEuler numerical scheme with a truncation error . Approximating the time derivative with a backward differentiation formula (BDF) we obtain a numerical scheme with a truncation error . We finally approximate the time derivative by a fourthorder RungeKutta approximation to obtain a numerical scheme with a truncation error . For large values of the FrankKamenetskii parameter and the shape factor there is a large discrepancy (numerical values differ by ) in the values obtained for the temperature at the centre of the vessel, , when the tolerance is set at between the three numerical schemes. To ensure that we obtain numerical values that are accurate we implement the numerical scheme obtained from approximating the time derivative by a fourthorder RungeKutta approximation.
We have shown that there is a strong dependence on the value obtained for the temperature at the centre of the vessel and the tolerance. A consequence of this strong dependence is that it is not practical to use the time as a measure of when a thermal explosion will take place. It is more practical to measure the temperature at the centre of the vessel or at the vessel wall. We have also shown that by changing the shape from a rectangular slab to a sphere decreases the temperature at which a thermal explosion will take place. Similarly, we have shown that increasing the FrankKamenetskii parameter will increase the temperature at which a thermal explosion will take place. The results indicted in Figures 3 and 4 are the same as those obtained by Harley [2] and Britz et al. [3].
Acknowledgment
E. Momoniat acknowledges support from the National Research Foundation of South Africa.
References
 D. A. FrankKamenetskii, Diffusion and Heat Transfer in Chemical Kinetics, Plenum Press, New York, NY, USA, 1969.
 C. Harley, “Hopscotch method: the numerical solution of the FrankKamenetskii partial differential equation,” Applied Mathematics and Computation, vol. 217, no. 8, pp. 4065–4075, 2010. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 D. Britz, J. Strutwolf, and O. Østerby, “Digital simulation of thermal reactions,” Applied Mathematics and Computation, vol. 218, no. 4, pp. 1280–1290, 2011. View at: Publisher Site  Google Scholar
 E. Balakrishnan, A. Swift, and G. C. Wake, “Critical values for some nonclass A geometries in thermal ignition theory,” Mathematical and Computer Modelling, vol. 24, no. 8, pp. 1–10, 1996. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 N. W. Bazley and G. C. Wake, “The dependence of criticality on activation energy when reactant consumption is neglected,” Combustion and Flame, vol. 33, no. 2, pp. 161–168, 1978. View at: Publisher Site  Google Scholar
 C. S. Chen, “The method of fundamental solutions for nonlinear thermal explosions,” Communications in Numerical Methods in Engineering, vol. 11, no. 8, pp. 675–681, 1995. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. Kubíček and V. Hlaváček, “Direct evaluation of branching points for equations arising in the theory of explosions of solid explosives,” Journal of Computational Physics, vol. 17, pp. 79–86, 1975. View at: Google Scholar
 G. Zhang, J. H. Merkin, and S. K. Scott, “Reactiondiffusion model for combustion with fuel consumption. II. Robin boundary conditions,” IMA Journal of Applied Mathematics, vol. 51, no. 1, pp. 69–93, 1993. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 O. K. Rice, “The role of heat conduction in thermal gaseous explosions,” The Journal of Chemical Physics, vol. 8, no. 9, pp. 727–733, 1940. View at: Publisher Site  Google Scholar
 O. K. Rice and H. C. Campbell, “The explosion of ethyl azide in the presence of diethyl ether,” The Journal of Chemical Physics, vol. 7, no. 8, pp. 700–709, 1939. View at: Publisher Site  Google Scholar
 J. J. Steggerda, “Thermal stability: an extension of FrankKamenetsky's theory,” The Journal of Chemical Physics, vol. 43, no. 12, pp. 4446–4448, 1965. View at: Publisher Site  Google Scholar
 G. Zhang, J. H. Merkin, and S. K. Scott, “Reactiondiffusion model for combustion with fuel consumption. I. Dirichlet boundary conditions,” IMA Journal of Applied Mathematics, vol. 47, no. 1, pp. 33–60, 1991. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. R. Abd ElSalam and M. H. Shehata, “The numerical solution for reactiondiffusion combustion with fuel consumption,” Applied Mathematics and Computation, vol. 160, no. 2, pp. 423–435, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 P. M. Vanden Bosch, D. C. Dietz, and E. A. Pohl, “Choosing the best approach to matrix exponentiation,” Computers & Operations Research, vol. 26, no. 9, pp. 871–882, 1999. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 C. Moler and C. Van Loan, “Nineteen dubious ways to compute the exponential of a matrix, twentyfive years later,” SIAM Review, vol. 45, no. 1, pp. 3–49, 2003. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. Hochbruck and A. Ostermann, “Exponential integrators,” Acta Numerica, vol. 19, pp. 209–286, 2010. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 Y. Saad, “Analysis of some Krylov subspace approximations to the matrix exponential operator,” SIAM Journal on Numerical Analysis, vol. 29, no. 1, pp. 209–228, 1992. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 A. Iserles and A. Zanna, “Efficient computation of the matrix exponential by generalized polar decompositions,” SIAM Journal on Numerical Analysis, vol. 42, no. 5, pp. 2218–2256, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 R. B. Sidje, “Expokit: a software package for computing Matrix Exponentials,” ACM Transactions on Mathematical Software, vol. 24, no. 1, pp. 130–156, 1998. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 H. Berland, B. Skaflestad, and W. M. Wright, “EXPINTA MATLAB package for exponential integrators,” ACM Transactions on Mathematical Software, vol. 33, pp. 1–17, 2007. View at: Google Scholar
 A. Iserles, A first Course in the Numerical Analysis of Differential Equations, Cambridge University Press, Cambridge, UK, 1996.
 S. Gottlieb, C.W. Shu, and E. Tadmor, “Strong stabilitypreserving highorder time discretization methods,” SIAM Review, vol. 43, no. 1, pp. 89–112, 2001. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 A.K. Kassam and L. N. Trefethen, “Fourthorder timestepping for stiff PDEs,” SIAM Journal on Scientific Computing, vol. 26, no. 4, pp. 1214–1233, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 E. Momoniat, “A thermal explosion in a cylindrical vessel: a nonclassical symmetry approach,” International Journal of Modern Physics B, vol. 23, no. 14, pp. 3089–3099, 2009. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 P. L. Chambré, “On the solution of the poissonboltzmann equation with application to the theory of thermal explosions,” The Journal of Chemical Physics, vol. 20, no. 11, pp. 1795–1797, 1952. View at: Publisher Site  Google Scholar
 L. Dresner, “Phaseplane analysis of nonlinear, secondorder, ordinary differential equations,” Journal of Mathematical Physics, vol. 12, no. 7, pp. 1339–1348, 1971. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 E. Momoniat and C. Harley, “An implicit series solution for a boundary value problem modelling a thermal explosion,” Mathematical and Computer Modelling, vol. 53, no. 12, pp. 249–260, 2011. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 M. P. Calvo and C. Palencia, “A class of explicit multistep exponential integrators for semilinear problems,” Numerische Mathematik, vol. 102, no. 3, pp. 367–381, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 A. Ostermann, M. Thalhammer, and W. M. Wright, “A class of explicit exponential general linear methods,” BIT, vol. 46, no. 2, pp. 409–431, 2006. View at: Publisher Site  Google Scholar  Zentralblatt MATH
Copyright
Copyright © 2012 E. Momoniat. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.