Research Article | Open Access
Matrix Exponentiation and the Frank-Kamenetskii Equation
Long time solutions to the Frank-Kamenetskii partial differential equation modelling a thermal explosion in a vessel are obtained using matrix exponentiation. Spatial derivatives are approximated by high-order finite difference approximations. A forward difference approximation to the time derivative leads to a Lawson-Euler scheme. Computations performed with a BDF approximation to the time derivative and a fourth-order Runge-Kutta approximation to the time derivative are compared to results obtained with the Lawson-Euler scheme. Variation in the central temperature of the vessel corresponding to changes in the shape parameter and Frank-Kamenetskii parameter are computed and discussed.
In this paper we use matrix exponentiation to determine numerical solutions of the Frank-Kamenetskii 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 Frank-Kamenetskii parameter . 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 high-order 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 high-order finite difference schemes. Comparisons are made between two high-order approximations to the time derivative by comparing a BDF approximation and a Runge-Kutta approximation to the time derivative. The results obtained compare favorably with the results obtained by Harley  and Britz et al. .
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 well-defined geometries we have for a rectangular slab, for an infinite circular cylinder, and for a sphere. Balakrishnan et al.  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.
Frank-Kamenetskii  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  and Rice and Campbell . Steggerda  has extended the work done by Frank-Kamenetskii  on the criteria under which a thermal explosion can take place. In the models presented by Zhang et al. [8, 12] and Abd El-Salam and Shehata  the fuel used in the chemical reaction is consumed during the reaction.
Harley  has compared hopscotch, method of lines, and Crank-Nicolson schemes to solve the nonlinear partial differential equation (1.1). Britz et al.  improve on the work of Harley  by treating the nonlinear source term in (1.1) in both an implicit and explicit form. Britz et al.  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  and Britz et al. . 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 high-order finite difference approximations to the spatial derivatives coupled with BDF and Runge-Kutta 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 first-order 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. , Moler and Van Loan , and Hochbruck and Ostermann  on the construction of matrix exponentials for semi-linear systems of the form (1.11). Recent developments have focused on constructing matrix exponentials using Krylov subspace approximations  and polar decompositions . In this paper we use MATHEMATICA to determine the relevant matrix exponentials. Alternate computer packages to determine matrix exponentials have been developed by Sidje  and Berland et al. .
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 Runge-Kutta approximation to the time gives a scheme with a truncation error . A BDF approximation to the time derivative leads to A-stable multistep methods . A fourth-order Runge-Kutta 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.  and Kassam and Trefethen  for discussions about high-order 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 Frank-Kamenetskii parameter on the temperature at the centre of the vessel. Momoniat  has used the Lie symmetry approach to determine group invariant solutions admitted by (1.1) for . The solutions obtained by Momoniat  are valid after blowup.
In this section we determine approximate solutions admitted by (1.1) by using the approximation . As discussed in Harley  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 Frank-Kamenetskii 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.
3. Derivation and Implementation of the Numerical Scheme
In this section we compare three numerical schemes for solving the Frank-Kamenetskii 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 Runge-Kutta 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 . High-order 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  in that we have used high-order 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 . Low-order 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 Lawson-Euler scheme The truncation error of results obtained from the Lawson-Euler 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.
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 Lawson-Euler 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 fourth-order Runge-Kutta 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 fourth-order Runge-Kutta numerical scheme (3.12) in the rest of the paper. The numerical scheme (3.12) based on a fourth-order Runge-Kutta 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 second-order ordinary differential equation and the boundary conditions
Frank-Kamenetskii  (see also Chambré , Dresner , and Momoniat and Harley ) 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 well-known 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 Lawson-Euler 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 Runge-Kutta 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 Lawson-Euler scheme (3.8) and the Runge-Kutta 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.
The results in this section have been determined using the numerical scheme (3.12) corresponding to the fourth-order Runge-Kutta 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 high-order 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 Frank-Kamenetskii 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 forward-time central space (FTCS) scheme in that the scheme inherently contains higher-order time-step terms. This can be seen if we consider the approximation we find that the Lawson-Euler 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 high-order spatial and high-order time discretizations. Future work will investigate the use of multistep exponential integrators of the type introduced by Calvo and Palencia  and Ostermann et al.  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 Lawson-Euler 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 fourth-order Runge-Kutta approximation to obtain a numerical scheme with a truncation error . For large values of the Frank-Kamenetskii 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 fourth-order Runge-Kutta 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 Frank-Kamenetskii 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  and Britz et al. .
E. Momoniat acknowledges support from the National Research Foundation of South Africa.
- D. A. Frank-Kamenetskii, Diffusion and Heat Transfer in Chemical Kinetics, Plenum Press, New York, NY, USA, 1969.
- C. Harley, “Hopscotch method: the numerical solution of the Frank-Kamenetskii partial differential equation,” Applied Mathematics and Computation, vol. 217, no. 8, pp. 4065–4075, 2010.
- D. Britz, J. Strutwolf, and O. Østerby, “Digital simulation of thermal reactions,” Applied Mathematics and Computation, vol. 218, no. 4, pp. 1280–1290, 2011.
- E. Balakrishnan, A. Swift, and G. C. Wake, “Critical values for some non-class A geometries in thermal ignition theory,” Mathematical and Computer Modelling, vol. 24, no. 8, pp. 1–10, 1996.
- 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.
- C. S. Chen, “The method of fundamental solutions for non-linear thermal explosions,” Communications in Numerical Methods in Engineering, vol. 11, no. 8, pp. 675–681, 1995.
- 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.
- G. Zhang, J. H. Merkin, and S. K. Scott, “Reaction-diffusion model for combustion with fuel consumption. II. Robin boundary conditions,” IMA Journal of Applied Mathematics, vol. 51, no. 1, pp. 69–93, 1993.
- 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.
- 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.
- J. J. Steggerda, “Thermal stability: an extension of Frank-Kamenetsky's theory,” The Journal of Chemical Physics, vol. 43, no. 12, pp. 4446–4448, 1965.
- G. Zhang, J. H. Merkin, and S. K. Scott, “Reaction-diffusion model for combustion with fuel consumption. I. Dirichlet boundary conditions,” IMA Journal of Applied Mathematics, vol. 47, no. 1, pp. 33–60, 1991.
- M. R. Abd El-Salam and M. H. Shehata, “The numerical solution for reaction-diffusion combustion with fuel consumption,” Applied Mathematics and Computation, vol. 160, no. 2, pp. 423–435, 2005.
- 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.
- C. Moler and C. Van Loan, “Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later,” SIAM Review, vol. 45, no. 1, pp. 3–49, 2003.
- M. Hochbruck and A. Ostermann, “Exponential integrators,” Acta Numerica, vol. 19, pp. 209–286, 2010.
- 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.
- 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.
- R. B. Sidje, “Expokit: a software package for computing Matrix Exponentials,” ACM Transactions on Mathematical Software, vol. 24, no. 1, pp. 130–156, 1998.
- H. Berland, B. Skaflestad, and W. M. Wright, “EXPINT-A MATLAB package for exponential integrators,” ACM Transactions on Mathematical Software, vol. 33, pp. 1–17, 2007.
- 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 stability-preserving high-order time discretization methods,” SIAM Review, vol. 43, no. 1, pp. 89–112, 2001.
- A.-K. Kassam and L. N. Trefethen, “Fourth-order time-stepping for stiff PDEs,” SIAM Journal on Scientific Computing, vol. 26, no. 4, pp. 1214–1233, 2005.
- E. Momoniat, “A thermal explosion in a cylindrical vessel: a non-classical symmetry approach,” International Journal of Modern Physics B, vol. 23, no. 14, pp. 3089–3099, 2009.
- P. L. Chambré, “On the solution of the poisson-boltzmann equation with application to the theory of thermal explosions,” The Journal of Chemical Physics, vol. 20, no. 11, pp. 1795–1797, 1952.
- L. Dresner, “Phase-plane analysis of nonlinear, second-order, ordinary differential equations,” Journal of Mathematical Physics, vol. 12, no. 7, pp. 1339–1348, 1971.
- 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. 1-2, pp. 249–260, 2011.
- 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.
- 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.
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.