Matrix Exponentiation and the Frank-Kamenetskii Equation

Abstract

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.

1. 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: 𝜕𝑢=𝜕𝜕𝑡²𝑢𝜕𝑥²+𝑘𝑥𝜕𝑢𝜕𝑥+𝛿𝑒^𝑢, (1.1) where 𝑢(𝑥,𝑡) is the nondimensional temperature in a vessel, 𝑘 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𝜕𝑢()𝑢()𝜕𝑥0,𝑡=0,1,𝑡=0, (1.2a)(1.2b)and the initial condition 𝑢(𝑥,0)=0. (1.3)

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 [2] and Britz et al. [3].

The FKPDE (1.1) is derived from a heat balance equation given by 𝑐𝜕𝑇𝜕𝑡=𝑘_𝑐∇²−𝐸𝑇+𝜎𝑄𝐴exp,𝑅𝑇 (1.4) 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 𝐸𝜃=𝑅𝑇²₀𝑇−𝑇₀, (1.5) where 𝑇₀ is the ambient temperature into (1.4). The Laplacian ∇² is given by ∇²=𝜕²𝜕𝑟²+𝑁−1𝑟𝜕,𝜕𝑟 (1.6) where the constant 𝑁 is related to the geometry of the problem.

We make a change of variables𝑢=𝜃−𝜃₀,𝑥=𝑟𝛿𝑒^𝜃01/2, (1.7a)(1.7b)where 𝜃₀ is a dimensionless temperature at the vessel walls and 𝛿=𝜎𝑄𝐴𝑘_𝑐𝐸𝑅𝑇²₀−𝐸exp𝑅𝑇₀. (1.8) 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𝜕𝑢()𝑢()𝜕𝑥0,𝑡=0,𝜆,𝑡=0, (1.9a)(1.9b)where 𝜆 is the boundary of vessel and 𝑁=𝑘+1. We choose 𝜆=1 to obtain the boundary conditions (1.2a) and (1.2b). For well-defined geometries we have 𝑘=0 for a rectangular slab, 𝑘=1 for an infinite circular cylinder, and 𝑘=2 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 𝑥=1 given by [5–8]: 𝜕𝑢()()𝜕𝑥1,𝑡+𝛼𝑢1,𝑡=0, (1.10) where 𝛼 is a constant.

Frank-Kamenetskii [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 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.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 𝑑𝐮(𝐮),𝑑𝑡=ℒ𝐮+𝒩 (1.11) 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. [14], Moler and Van Loan [15], and Hochbruck and Ostermann [16] 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 [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 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 [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.

The paper is divided up as follows. In Section 2 we consider approximate solutions for the case 𝑘=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.

2. Analytical Solution

In this section we consider an approximate analytical solution admitted by (1.1) for the case 𝑘=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 𝑘=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(𝑢)≈1+𝑢. As discussed in Harley [2] this approximation is valid close to the centre of the cylindrical vessel 𝑥≈0 and yields results that give an error 𝑂(10⁻²) when compared to the numerical solution at the steady state. The model equation simplifies to 𝜕𝑢=𝜕𝜕𝑡²𝑢𝜕𝑥²+𝛿(1+𝑢). (2.1) We make the assumption 𝛿≪1 and investigate approximate solutions of the form 𝑢(𝑥,𝑡)≈𝑢₀(𝑥,𝑡)+𝛿𝑢₁(𝑥,𝑡)+⋯. (2.2) Substituting (2.2) into (2.1) and separating to first order in 𝛿 we obtain the system𝜕𝑢₀=𝜕𝜕𝑡²𝑢₀𝜕𝑥²,𝜕𝑢₁=𝜕𝜕𝑡²𝑢₁𝜕𝑥²+1+𝑢₀. (2.3a)(2.3b)The heat kernel solution to 𝑢₀ is given by 𝑢₀()=1𝑥,𝑡√−𝑥4𝜋𝑡exp².4𝑡 (2.4) By inspection, we find that 𝑢₁ admits the particular solution 𝑢₁(𝑥,𝑡)=𝑢₀(𝑥,𝑡)𝑡+𝑐₀+𝑡, (2.5) 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 𝑢_𝐴=𝑢(𝑥,𝑡)≈𝑢₀(𝑥,𝑡)+𝛿𝑢₀(𝑥,𝑡)𝑡+𝑐₀,𝑢+𝑡_𝐵=𝑢(𝑥,𝑡)≈𝑢₀(𝑥,𝑡)+𝛿𝑢₀(𝑥,𝑡)+𝑢₀(𝑥,𝑡)𝑡+𝑐₀+𝑡. (2.6)

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.

(a)

(b)

Figure 1 

Plot comparing the approximate solutions (2.6).

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 𝑥=0 we approximate (1/𝑥)𝜕𝑢/𝜕𝑥 by 𝜕²𝑢/𝜕𝑥² to write (1.1) as 𝜕𝑢=()𝜕𝜕𝑡1+𝑘²𝑢𝜕𝑥²+𝛿𝑒^𝑢𝜕,𝑥=0,₂𝑢𝜕𝑥²+𝑘𝑥𝜕𝑢𝜕𝑥+𝛿𝑒^𝑢,𝑥≠0. (3.1) We approximate the spatial derivatives by the central difference approximations 𝜕𝑢𝜕𝑥_{𝑥=𝑥𝑖}≈𝑢^𝑗_𝑖+1−𝑢^𝑗_𝑖−12Δ𝑥+𝒪Δ𝑥²,𝜕²𝑢𝜕𝑥²_{𝑥=𝑥𝑖}≈𝑢^𝑗_𝑖+1−2𝑢^𝑗_𝑖+𝑢^𝑗_𝑖−1Δ𝑥²+𝒪Δ𝑥², (3.2) where 𝑢^𝑗_𝑖=𝑢(𝑥_𝑖,𝑡_𝑗). High-order central difference approximations to the derivatives are given by 𝜕𝑢𝜕𝑥_{𝑥=𝑥𝑖}≈−𝑢^𝑗_𝑖+2+8𝑢^𝑗_𝑖+1−8𝑢^𝑗_𝑖−1+𝑢^𝑗_𝑖−212Δ𝑥+𝒪Δ𝑥⁴,𝜕²𝑢𝜕𝑥²_{𝑥=𝑥𝑖}≈−𝑢^𝑗_𝑖+2+16𝑢^𝑗_𝑖+1−30𝑢^𝑗_𝑖+16𝑢^𝑗_𝑖−1−𝑢^𝑗_𝑖−212Δ𝑥²+𝒪Δ𝑥⁴. (3.3) The derivative boundary condition (1.2a) gives the condition 𝑢^𝑗₋₁=𝑢^𝑗₁ when 𝑖=0. Defining 𝑥_𝑖=𝑖Δ𝑥 and 𝐮=[𝑢₀,𝑢₁,…,𝑢_𝑚] where exp(𝐮)=[𝑒^𝑢0,𝑒^𝑢1,…,𝑒^𝑢𝑛] we write (1.1) in vector form as 𝑑𝐮𝑑𝑡=𝐀𝐮+𝛿𝑒^𝐮, (3.4) where 1𝐀=Δ𝑥²()2()𝛼−21+𝑘1+𝑘0000⋯000000₁−2𝛽₁𝑎000⋯⋮⋮⋮⋮⋮⋮₂𝑏₂−30𝑐12²𝑑₂0⋯⋮⋮⋮⋮⋮⋮0𝑎₃𝑏₃−30𝑐12³𝑑₃⋯⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋯𝑎_𝑛−3𝑏_𝑛−3−30𝑐12^𝑛−3𝑑_𝑛−30⋮⋮⋮⋮⋮⋮⋯0𝑎_𝑛−2𝑏_𝑛−2−30𝑐12^𝑛−2𝑑_𝑛−2⋮⋮⋮⋮⋮⋮⋯000𝛼_𝑛−1−2𝛽_𝑛−1,𝛼000000⋯000001_𝑖=11−2𝑖,𝛽_𝑖=11+,𝑎2𝑖_𝑖1=−𝑘121−𝑖,𝑏_𝑖=11216−8𝑘𝑖,𝑐_𝑖=11216+8𝑘𝑖,𝑑_𝑖1=−𝑘121+𝑖. (3.5) The matrix 𝐀 is different to the matrix obtained by Harley [2] 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 [2]. Low-order approximations are used in row positions 𝑖=2 and 𝑖=𝑛−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 𝑖=0 and 𝑖=𝑛 are determined from the boundary conditions (1.2a) and (1.2b).

We multiply (3.4) by the integrating factor exp(−𝐀𝑡) to obtain 𝑑()()(𝐮),𝑑𝑡exp−𝐀𝑡⋅𝐮=𝛿exp−𝐀𝑡⋅exp (3.6) where ⋅ is used to signify the multiplication of a matrix with a vector. Using a forward difference approximation to the time derivative given by 𝑑𝐮𝑑𝑡_{𝑡=𝑡𝑗}≈𝐮^𝑗+1−𝐮^𝑗()Δ𝑡+𝒪Δ𝑡 (3.7) we obtain the Lawson-Euler scheme 𝐮^𝑗+1=exp(𝐀Δ𝑡)⋅𝐮^𝑗+𝛿Δ𝑡exp𝐮^𝑗. (3.8) 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 𝑢₀^𝑗+1−𝑢^𝑗₀<𝜏, (3.9) where 𝜏 is a tolerance. The quantity 𝑢^𝑗₀=𝑢(0,𝑡_𝑗). When the tolerance condition (3.9) has been satisfied we set 𝑢^𝑗₀=𝑢(0,𝑡_𝑗)=𝑢(0), the temperature at the centre of the vessel. The boundary condition (1.2b) is imposed as 𝑢^𝑗_𝑛=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𝐮^𝑗+1−4𝐮^𝑗+𝐮^𝑗−12Δ𝑡+𝒪Δ𝑡². (3.10) We then obtain the scheme 𝐮^𝑗+1=13exp(𝐀Δ𝑡)⋅4𝐮^𝑗+2𝛿Δ𝑡exp𝐮^𝑗−13exp(𝐀(2Δ𝑡))⋅𝐮^𝑗−1. (3.11)

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 𝐮^𝑗+1=exp(𝐀Δ𝑡)⋅𝐮^𝑗+16𝑘₁+2𝑘₂+2𝑘₃+𝑘₄+𝒪Δ𝑡⁴, (3.12) where 𝑘₁=𝛿Δ𝑡exp(𝐀Δ𝑡)⋅exp𝐮^𝑗,𝑘₂=𝛿Δ𝑡exp𝐀Δ𝑡2⋅exp𝐮^𝑗+𝑘₁2,𝑘₃=𝛿Δ𝑡exp𝐀Δ𝑡2⋅exp𝐮^𝑗+𝑘₂2,𝑘₄=𝛿Δ𝑡exp𝐮^𝑗+𝑘₃. (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 𝒪(Δ𝑡⁴)+𝒪(Δ𝑥²).

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 𝒪(10⁻⁶). For larger values of 𝑘 and 𝛿 this difference increases to 𝒪(10⁻³). 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.

Table 1 

Table comparing values obtained from the Lawson-Euler (LE) scheme, BDF scheme, and the fourth-order Runge-Kutta (RK4) scheme where Δ𝑥=0.005, Δ𝑡=0.00001, and 𝛿=1 for different values of 𝑘. The tolerance 𝜏=10⁻⁶.

Table 2 

Table comparing values obtained from the Lawson-Euler (LE) scheme, BDF scheme, and the fourth-order Runge-Kutta (RK4) scheme where Δ𝑥=0.005, Δ𝑡=0.00001, and 𝑘=1 for different values of 𝛿. The tolerance 𝜏=10⁻⁶.

To test the stability of each of the numerical schemes we consider the steady solution for the case 𝑘=1. For the case 𝑘=1 the model equation (1.1) must satisfy the steady solution 𝑢(𝑥,𝑡)=𝑦(𝑥) where 𝑦(𝑥) satisfies the second-order ordinary differential equation 𝑑²𝑦𝑑𝑥²+1𝑥𝑑𝑦𝑑𝑥+𝛿𝑒^𝑦=0 (3.14) and the boundary conditions 𝑑𝑦(0)(1)𝑑𝑥=0,𝑦=0. (3.15)

Frank-Kamenetskii [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 𝑦=𝑦₁=log16𝑒^𝑐12𝛿+𝑒^𝑐1𝑥²²,𝑦=𝑦₂=log16𝑒^𝑐21+2𝛿𝑒^𝑐2𝑥²², (3.16) where 𝑐₁ and 𝑐₂ are constants of integration given by 𝑐₁=log2(4−𝛿)±42(2−𝛿),𝑐₂=log4−𝛿±22(2−𝛿)2𝛿². (3.17) These multiple solutions lead to the well-known bifurcation that occurs when plotting the temperature at the centre of the vessel 𝑦(0) against the critical parameter 𝛿.

We consider the solution 𝑦=𝑦^∗=log16𝑒^𝑐∗1+2𝛿𝑒^𝑐∗𝑥²²,𝑐^∗=log4−𝛿−22(2−𝛿)2𝛿², (3.18) that has a maximum 𝑦(0)=0.316694. This corresponds to the results given in Table 3 for the central temperature at a tolerance of 𝜏=10⁻¹². 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.

Table 3 

Table comparing the temperature at the centre of the vessel for different stopping criteria and the corresponding stopping times where Δ𝑥=0.005, Δ𝑡=0.00001, 𝑘=1, and 𝛿=1.

In Figure 2 we plot log|𝑢^𝑗_𝑖−𝑦^∗_𝑖| for 𝑗=10,100, and 1000, respectively, where we have chosen 𝛿=1, Δ𝑥=0.005, Δ𝑡=0.00001, and 𝑘=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.

(a)

(b)

(c)

Figure 2 

Plot comparing the log|𝑢𝑗𝑖−𝑦∗𝑖| for 𝑗=10,100, and 1000, respectively, where we have chosen 𝛿=1, Δ𝑥=0.005, Δ𝑡=0.00001, and 𝑘=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).

4. 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 Δ𝑥=0.005 and Δ𝑡=0.00001 for 𝑘=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, 𝑢(0), for each value of the tolerance.

From the results in Table 3 we note that the temperature converges to 𝑢(0)=0.316705 but the stopping time increases from 𝑡=0.60297 to 𝑡=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 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 𝑇=𝑇₀+𝑅𝑇²₀𝐸𝜃₀+𝜏, (4.1) 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 𝑢(0) at the centre of the vessel corresponding to the wall temperature 𝜃₀+𝜏. 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⁻⁶ for the steady state.

In Figure 3 we plot the temperature at the centre of the vessel 𝑢(0) against the shape factor 𝑘 for fixed 𝛿=1. We note that for increasing the value of the shape factor 𝑘 the steady temperature 𝑢(0) decreases for the tolerance 𝜏=10⁻⁶.

Figure 3 

Plot showing the decreasing value of the temperature at the centre of the vessel 𝑢(0) with the increasing value of the shape factor 𝑘.

In Figure 4 we plot the temperature at the centre of the vessel against the increasing value of the Frank-Kamenetskii parameter 𝛿 for fixed 𝑘=1. We note that for the increasing value of 𝛿 the steady temperature 𝑢(0) increases for the tolerance 𝜏=10⁻⁶. These results match up with the approximate analytical results predicted by the study of the approximate heat kernel solution in Section 2.

Figure 4 

Plot showing the increasing value of the temperature at the centre of the vessel 𝑢(0) with the increasing value of the critical parameter 𝛿.

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 exp(𝐀Δ𝑡)≈𝐈+Δ𝑡𝐀+Δ𝑡²𝐀2!²+⋯, (5.1)we find that the Lawson-Euler scheme (3.8) simplifies the FTCS scheme,𝐮^𝑗+1=(𝐈+Δ𝑡𝐀)⋅𝐮^𝑗+𝛿Δ𝑡exp𝐮^𝑗+𝒪Δ𝑡², (5.2) 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 [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 𝒪(Δ𝑡)+𝒪(Δ𝑥²). 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 𝒪(10⁻³)) in the values obtained for the temperature at the centre of the vessel, 𝑢(0), when the tolerance is set at 𝜏=10⁻⁶ 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 [2] and Britz et al. [3].

Acknowledgment

E. Momoniat acknowledges support from the National Research Foundation of South Africa.