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 obtainedwith 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.


Introduction
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 u x, t is the nondimensional temperature in a vessel, k is a coefficient that represents the shape of the vessel, and δ is the Frank-Kamenetskii parameter 1 .The FKPDE 1.1 is solved subject to the boundary conditions ∂u ∂x 0, t 0, 1.2a u 1, t 0, 1.2b and the initial condition u x, 0 0.

1.3
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 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 2 and Britz et al. 3 .
The FKPDE 1.1 is derived from a heat balance equation given by where c is the thermal capacity, k c is the thermal conductivity, σ the density, Q is the heat of reaction, A is the frequency factor, E is the energy of activation of the chemical reaction, R is the universal gas constant, and T is the gas temperature.The heat balance equation 1.4 is made nondimensional by substituting where T 0 is the ambient temperature into 1.4 .The Laplacian ∇ 2 is given by N − 1 r ∂ ∂r , 1.6 where the constant N is related to the geometry of the problem.We make a change of variables u θ − θ 0 , 1.7a x r δe θ 0 1/2 , 1.7b where θ 0 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 by where λ is the boundary of vessel and N k 1.We choose λ 1 to obtain the boundary conditions 1.2a and 1.2b .For well-defined geometries we have k 0 for a rectangular slab, k 1 for an infinite circular cylinder, and k 2 for a sphere.Balakrishnan et al. 4 have solved the steady problem for nonstandard geometries, that is, k is not an integer.The problem 1.1 has been solved using the Robin boundary conditions at the boundary x 1 given by 5-8 : where α is a constant.Frank-Kamenetskii 1 has derived 1.1 by ignoring coefficients O RT 0 /E 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  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.The paper is divided up as follows.In Section 2 we consider approximate solutions for the case k 0 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.

Analytical Solution
In this section we consider an approximate analytical solution admitted by 1.1 for the case k 0 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 24 has used the Lie symmetry approach to determine group invariant solutions admitted by 1.1 for k 1.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 exp u ≈ 1 u.As discussed in Harley 2 this approximation is valid close to the centre of the cylindrical vessel x ≈ 0 and yields results that give an error O 10 −2 when compared to the numerical solution at the steady state.The model equation simplifies to We make the assumption δ 1 and investigate approximate solutions of the form Substituting 2.2 into 2.1 and separating to first order in δ we obtain the system The heat kernel solution to u 0 is given by By inspection, we find that u 1 admits the particular solution where c 0 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

2.6
It is immediately obvious from u A and u B 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 u 0 .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 c 0 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.

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 x 0 we approximate 1/x ∂u/∂x by ∂ 2 u/∂x 2 to write 1.1 as

3.1
We approximate the spatial derivatives by the central difference approximations where u j i u x i , t j .High-order central difference approximations to the derivatives are given by

3.3
The derivative boundary condition 1.2a gives the condition u j −1 u j 1 when i 0. Defining x i iΔx and u u 0 , u 1 , . . ., u m where exp u e u 0 , e u 1 , . . ., e u n we write 1.1 in vector form as where

3.5
The matrix A 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 .Low-order approximations are used in row positions i 2 and i n − 1 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 i 0 and i n are determined from the boundary conditions 1.2a and 1.2b .We multiply 3.4 by the integrating factor exp −At 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 u j 1 exp AΔt • u j δΔt exp u j .

3.8
The truncation error of results obtained from the Lawson-Euler scheme 3.

Mathematical Problems in Engineering
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 u j 0 u 0, t j .When the tolerance condition 3.9 has been satisfied we set u j 0 u 0, t j u 0 , the temperature at the centre of the vessel.The boundary condition 1.2b is imposed as u j n 0 at each iteration.An improvement on the accuracy of the Lawson-Euler scheme 3.8 can be made by considering a backward difference formula BDF to the time derivative given by 21

3.10
We then obtain the scheme

3.11
The numerical scheme 3.11 is implemented in MATHEMATICA.The scheme 3.11 requires two starting values u 0 which we obtain from the initial value and u 1 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 O Δt 2 O Δx 2 .We lastly consider a fourth-order Runge-Kutta scheme to solve 3.6 .The resulting numerical scheme is given by where

3.13
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 O Δt 4 O Δx 2 .
In Tables 1 and 2 we compare output from numerical schemes 3.8 , 3.11 , and 3.12 for different values of δ and k.In Table 1 we fix the value of δ and vary k, while in Table 2 we fix the value of k and vary δ.We note from the results in Tables 1 and 2 that for small values of k and δ the difference in the output from the numerical schemes is O 10 −6 .For where c 1 and c 2 are constants of integration given by 3.17 These multiple solutions lead to the well-known bifurcation that occurs when plotting the temperature at the centre of the vessel y 0 against the critical parameter δ.
We consider the solution that has a maximum y 0 0.316694.This corresponds to the results given in Table 3 for the central temperature at a tolerance of τ 10 −12 .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 log |u j i − y * i | for j 10, 100, and 1000, respectively, where we have chosen δ 1, Δx 0.005, Δt 0.00001, and k 1.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.

Results
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 Δx 0.005 and Δt 0.00001 for k 1 and δ 1 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, u 0 , for each value of the tolerance.
From the results in Table 3 we note that the temperature converges to u 0 0.316705 but the stopping time increases from t 0.60297 to t 3.66920.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 We have used 3.18 as the initial value of the numerical schemes where LE scheme corresponds to the Lawson-Euler scheme 3.8 , BDF scheme corresponds to the numerical scheme 3.11 , and RK4 scheme corresponds to the numerical scheme 3.12 .
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 T 0 is the ambient temperature, θ 0 is the dimensionless temperature at the vessel wall, R is the universal gas constant, and E 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 u 0 at the centre of the vessel corresponding to the wall temperature θ 0 τ.As τ → 0 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 τ 10 −6 for the steady state.In Figure 3 we plot the temperature at the centre of the vessel u 0 against the shape factor k for fixed δ 1.We note that for increasing the value of the shape factor k the steady temperature u 0 decreases for the tolerance τ 10 −6 .
In Figure 4 we plot the temperature at the centre of the vessel against the increasing value of the Frank-Kamenetskii parameter δ for fixed k 1.We note that for the increasing value of δ the steady temperature u 0 increases for the tolerance τ 10 −6 .These results match up with the approximate analytical results predicted by the study of the approximate heat kernel solution in Section 2.

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 O Δt .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 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 Lawson-Euler numerical scheme with a truncation error O Δt O Δx 2 .Approximating the time derivative with a backward differentiation formula BDF we obtain a numerical scheme with a truncation error O Δt 2 O Δx 2 .We finally approximate the time derivative by a fourth-order Runge-Kutta approximation to obtain a numerical scheme with a truncation error O Δt 4 O Δx 2 .For large values of the Frank-Kamenetskii parameter δ and the shape factor k there is a large discrepancy numerical values differ by O 10 −3 in the values obtained for the temperature at the centre of the vessel, u 0 , when the tolerance is set at τ 10 −6 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 9 and Rice and Campbell 10 .Steggerda 11 has extended the work done by Frank-Kamenetskii 1 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 13 the fuel used in the chemical reaction is consumed during the reaction.Harley 2 has compared hopscotch, method of lines, and Crank-Nicolson 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.

Figure 2 :
Figure 2: Plot comparing the log |u j i − y * i | for j 10, 100, and 1000, respectively, where we have chosen δ 1, Δx 0.005, Δt 0.00001, and k 1.We have used 3.18 as the initial value of the numerical schemes where LE scheme corresponds to the Lawson-Euler scheme 3.8 , BDF scheme corresponds to the numerical scheme 3.11 , and RK4 scheme corresponds to the numerical scheme 3.12 .

Figure 3 :Figure 4 :
Figure 3: Plot showing the decreasing value of the temperature at the centre of the vessel u 0 with the increasing value of the shape factor k.
are the same as those obtained by Harley 2 and Britz et al. 3 .
Δt 2 O Δx 2 where Δx is the spatial step and Δt is the time step.A Runge-Kutta approximation to the time gives a scheme with a truncation error O Δt 4 O Δx 2 .A BDF approximation to the time derivative leads to A-stable multistep methods 21 .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. 22 and Kassam and Trefethen 23 for discussions about high-order approximations to the time derivate and their respective properties.

Table 1 :
Table comparing values obtained from the Lawson-Euler LE scheme, BDF scheme, and the fourth-order Runge-Kutta RK4 scheme where Δx 0.005, Δt 0.00001, and δ 1 for different values of k.The tolerance τ 10 −6 .

Table 2 :
Table comparing values obtained from the Lawson-Euler LE scheme, BDF scheme, and the fourth-order Runge-Kutta RK4 scheme where Δx 0.005, Δt 0.00001, and k 1 for different values of δ.The tolerance τ 10 −6 .Δt 4O Δx 2 and are hence the most accurate.To test the stability of each of the numerical schemes we consider the steady solution for the case k 1.For the case k 1 the model equation 1.1 must satisfy the steady solution larger values of k and δ this difference increases to O 10 −3 .Due to this large discrepancy that occurs in the outputs from the different numerical schemes for large values of k 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, O

Table 3 :
Table comparing the temperature at the centre of the vessel for different stopping criteria and the corresponding stopping times where Δx 0.005, Δt 0.00001, k 1, and δ 1.