# Matrix Exponentiation and the Frank-Kamenetskii Equation

Authors: E. Momoniat
Journal: Mathematical Problems in Engineering
 (2012)

---

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


---

## Body

## 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)∂u∂t=∂2u∂x2+kx∂u∂x+δeu,
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(1.2a)∂u∂x(0,t)=0,(1.2b)u(1,t)=0,and the initial condition
(1.3)u(x,0)=0.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
(1.4)c∂T∂t=kc∇2T+σQAexp(-ERT),
where c is the thermal capacity, kc 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
(1.5)θ=ERT02(T-T0),
where T0 is the ambient temperature into (1.4). The Laplacian ∇2 is given by
(1.6)∇2=∂2∂r2+N-1r∂∂r,
where the constant N is related to the geometry of the problem.We make a change of variables(1.7a)u=θ-θ0,(1.7b)x=r[δeθ0]1/2,where θ0 is a dimensionless temperature at the vessel walls and
(1.8)δ=[σQAkcERT02exp(-ERT0)].
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(1.9a)∂u∂x(0,t)=0,(1.9b)u(λ,t)=0,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]:
(1.10)∂u∂x(1,t)+αu(1,t)=0,
where α is a constant.Frank-Kamenetskii [1] has derived (1.1) by ignoring coefficients O(RT0/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 [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)dudt=ℒu+𝒩(u),
where ℒ is a linear operator and 𝒩 is a nonlinear function of u. 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𝒪(Δt2)+𝒪(Δx2) 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 𝒪(Δt4)+𝒪(Δx2). 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 Section2 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.


## 2. 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
(2.1)∂u∂t=∂2u∂x2+δ(1+u).
We make the assumption δ≪1 and investigate approximate solutions of the form
(2.2)u(x,t)≈u0(x,t)+δu1(x,t)+⋯.
Substituting (2.2) into (2.1) and separating to first order in δ we obtain the system(2.3a)∂u0∂t=∂2u0∂x2,(2.3b)∂u1∂t=∂2u1∂x2+(1+u0).The heat kernel solution to u0 is given by
(2.4)u0(x,t)=14πtexp(-x24t).
By inspection, we find that u1 admits the particular solution
(2.5)u1(x,t)=u0(x,t)(t+c0)+t,
where c0 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)uA=u(x,t)≈u0(x,t)+δ[u0(x,t)(t+c0)+t],uB=u(x,t)≈u0(x,t)+δ[u0(x,t)+u0(x,t)(t+c0)+t].It is immediately obvious fromuA and uB 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 u0. 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 c0 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.Plot comparing the approximate solutions (2.6).
(a)
(b)


## 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 atx=0 we approximate (1/x)∂u/∂x by ∂2u/∂x2 to write (1.1) as
(3.1)∂u∂t={(1+k)∂2u∂x2+δeu,x=0,∂2u∂x2+kx∂u∂x+δeu,x≠0.
We approximate the spatial derivatives by the central difference approximations
(3.2)∂u∂x|x=xi≈ui+1j-ui-1j2Δx+𝒪(Δx2),∂2u∂x2|x=xi≈ui+1j-2uij+ui-1jΔx2+𝒪(Δx2),
where uij=u(xi,tj). High-order central difference approximations to the derivatives are given by
(3.3)∂u∂x|x=xi≈-ui+2j+8ui+1j-8ui-1j+ui-2j12Δx+𝒪(Δx4),∂2u∂x2|x=xi≈-ui+2j+16ui+1j-30uij+16ui-1j-ui-2j12Δx2+𝒪(Δx4).
The derivative boundary condition (1.2a) gives the condition u-1j=u1j when i=0. Defining xi=iΔx and u=[u0,u1,…,um] where exp(u)=[eu0,eu1,…,eun] we write (1.1) in vector form as
(3.4)dudt=Au+δeu,
where
(3.5)A=1Δx2[-2(1+k)2(1+k)0000⋯000000α1-2β1000⋯⋮⋮⋮⋮⋮⋮a2b2-3012c2d20⋯⋮⋮⋮⋮⋮⋮0a3b3-3012c3d3⋯⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋯an-3bn-3-3012cn-3dn-30⋮⋮⋮⋮⋮⋮⋯0an-2bn-2-3012cn-2dn-2⋮⋮⋮⋮⋮⋮⋯000αn-1-2βn-1000000⋯000001],αi=(1-12i),βi=(1+12i),ai=-112(1-ki),bi=112(16-8ki),ci=112(16+8ki),di=-112(1+ki).
The matrix A 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 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
(3.6)ddt(exp(-At)·u)=δexp(-At)·exp(u),
where · is used to signify the multiplication of a matrix with a vector. Using a forward difference approximation to the time derivative given by
(3.7)dudt|t=tj≈uj+1-ujΔt+𝒪(Δt)
we obtain the Lawson-Euler scheme
(3.8)uj+1=exp(AΔt)·[uj+δΔtexp(uj)].
The truncation error of results obtained from the Lawson-Euler scheme (3.8) is 𝒪(Δt)+𝒪(Δx2).We implement the numerical scheme (3.8) in MATHEMATICA. We assume that a steady state has been reached when
(3.9)|u0j+1-u0j|<τ,
where τ is a tolerance. The quantity u0j=u(0,tj). When the tolerance condition (3.9) has been satisfied we set u0j=u(0,tj)=u(0), the temperature at the centre of the vessel. The boundary condition (1.2b) is imposed as unj=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)dudt|t=tj≈3uj+1-4uj+uj-12Δt+𝒪(Δt2).
We then obtain the scheme
(3.11)uj+1=13exp(AΔt)·[4uj+2δΔtexp(uj)]-13exp(A(2Δt))·uj-1.The numerical scheme (3.11) is implemented in MATHEMATICA. The scheme (3.11) requires two starting values u0 which we obtain from the initial value and u1 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 𝒪(Δt2)+𝒪(Δx2).We lastly consider a fourth-order Runge-Kutta scheme to solve (3.6). The resulting numerical scheme is given by
(3.12)uj+1=exp(AΔt)·uj+16(k1+2k2+2k3+k4)+𝒪(Δt4),
where
(3.13)k1=δΔtexp(AΔt)·exp(uj),k2=δΔtexp(AΔt2)·exp(uj+k12),k3=δΔtexp(AΔt2)·exp(uj+k22),k4=δΔtexp(uj+k3).Once again the scheme is implemented inMATHEMATICA. 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 𝒪(Δt4)+𝒪(Δx2).In Tables1 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 𝒪(10-6). For larger values of k and δ this difference increases to 𝒪(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, 𝒪(Δt4)+𝒪(Δx2) 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Δx=0.005, Δt=0.00001, and δ=1 for different values of k. The tolerance τ=10-6.
k
LE
BDF
RK4
0.5
0.449971
0.449957
0.449978
1.0
0.294548
0.294539
0.294551
1.5
0.221880
0.221871
0.221881
2.0
0.178816
0.178810
0.178817Table 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.
δ
LE
BDF
RK4
0.5
0.119555
0.119552
0.119555
1.0
0.294548
0.294539
0.294551
1.5
0.546525
0.546501
0.546541
2.0
1.165950
1.165710
1.166260To test the stability of each of the numerical schemes we consider the steady solution for the casek=1. For the case k=1 the model equation (1.1) must satisfy the steady solution u(x,t)=y(x) where y(x) satisfies the second-order ordinary differential equation
(3.14)d2ydx2+1xdydx+δey=0
and the boundary conditions
(3.15)dy(0)dx=0,y(1)=0.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
(3.16)y=y1=log[16ec1(2δ+ec1x2)2],y=y2=log[16ec2(1+2δec2x2)2],
where c1 and c2 are constants of integration given by
(3.17)c1=log[2(4-δ)±42(2-δ)],c2=log[4-δ±22(2-δ)2δ2].
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(3.18)y=y*=log[16ec*(1+2δec*x2)2],c*=log[4-δ-22(2-δ)2δ2],
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.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.
Stopping criteria
Time
u (0)
1
0
-
6
0.60297
0.294551
1
0
-
7
1.11334
0.314486
1
0
-
8
1.62445
0.316483
1
0
-
9
2.13563
0.316683
1
0
-
10
2.64682
0.316703
1
0
-
11
3.15802
0.316705
1
0
-
12
3.66920
0.316705In Figure2 we plot log|uij-yi*| 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.Plot comparing thelog|uij-yi*| 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).
(a)
(b)
(c)


## 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 Δ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 Table3 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 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)T=T0+RT02E(θ0+τ),
where T0 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 Figure3 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.Figure 3
Plot showing the decreasing value of the temperature at the centre of the vesselu(0) with the increasing value of the shape factor k.In Figure4 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.Figure 4
Plot showing the increasing value of the temperature at the centre of the vesselu(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(5.1)exp(AΔt)≈I+ΔtA+Δt22!A2+⋯,we find that the Lawson-Euler scheme (3.8) simplifies the FTCS scheme,(5.2)uj+1=(I+ΔtA)·uj+δΔtexp(uj)+𝒪(Δt2),
to 𝒪(Δ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 𝒪(Δt)+𝒪(Δx2). Approximating the time derivative with a backward differentiation formula (BDF) we obtain a numerical scheme with a truncation error 𝒪(Δt2)+𝒪(Δx2). We finally approximate the time derivative by a fourth-order Runge-Kutta approximation to obtain a numerical scheme with a truncation error 𝒪(Δt4)+𝒪(Δx2). For large values of the Frank-Kamenetskii parameter δ and the shape factor k there is a large discrepancy (numerical values differ by 𝒪(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 Figures3 and 4 are the same as those obtained by Harley [2] and Britz et al. [3].



---

*Source: 713798-2012-12-04.xml*