Numerical Convergence Analysis of the Frank–Kamenetskii Equation

This work investigates the convergence dynamics of a numerical scheme employed for the approximation and solution of the Frank–Kamenetskii partial differential equation. A framework for computing the critical Frank–Kamenetskii parameter to arbitrary accuracy is presented and used in the subsequent numerical simulations. The numerical method employed is a Crank–Nicolson type implicit scheme coupled with a fourth order spatial discretisation as well as a Newton–Raphson update step which allows for the nonlinear source term to be treated implicitly. This numerical implementation allows for the analysis of the convergence of the transient solution toward the steady-state solution. The choice of termination criteria, numerically dictating this convergence, is interrogated and it is found that the traditional choice for termination is insufficient in the case of the Frank–Kamenetskii partial differential equation which exhibits slow transience as the solution approaches the steady-state. Four measures of convergence are proposed, compared and discussed herein.


Introduction
Thermal ignition has been of a topic of interest for many years, and continues to be an area of scientific relevance. The development of such explosions and solutions dependent upon various geometries present numerous challenges. Exothermic reactions are rife within the real world with many pertaining to sudden increases in pressure confined to physical boundaries. The impact of such constraints forces thermal reactions, heating systems to the point of explosion or steady-state.
This initial steady-state was originally described by Frank-Kamenetskii [1], extended further in 1969 to the time development of a thermal reaction [2]. The governing equation for this process under alternative notation is [2,3], where X the spatial variable, T is the temperature, τ the time, K the thermal diffusivity (m 2 s −1 ), c the specific heat capacity (J kg −1 K −1 ), Q the heat of the reaction (J kg −1 ), z the pre-exponential Arrhenius parameter (s −1 ), E the activation energy of the reaction and R the gas constant. The spatial variable describes the distance from the centre of the body for either (D = 1) a slab with thickness 2a and large length, (D = 2) a cylinder with radius a and large length or (D = 3) a sphere of radius a. Thus D specifies the geometry in which the reaction is taking place. Boundary conditions are provided as ∂T(0, τ) Other boundary conditions have been employed (see [3] and references therein) however we will employ those provided by ( (3) and (4)).
Normalisations are as follows [3] The Frank-Kamenetskii parameter δ is defined as [1,2,4], where expresses the rate of the reaction at ambient temperature T 0 and is often 1. This results in the dimensionless partial differential equation (PDE) ∂u ∂t where k = D − 1. The Frank-Kamenetskii parameter is of considerable importance given that at values below the critical value thereof, δ cr , a steady-state solution is reached for respective geometry and boundary conditions, while a thermal explosion presents itself for values above it. The Frank-Kamenetskii partial differential equation (FKPDE) has received a great deal of attention in recent years both in terms of computational and analytical work. Momoniat [5] obtained new solutions to the FKPDE via a non-classical symmetry approach. Numerical solutions have been obtained by Britz et al. [3], Harley [6], Momoniat [7] and Soliman [8] via various numerical techniques. Motsa and Sibanda [9] present a quasi-linearisation technique for nonlinear differential equations with higher-order convergence which is applied to the FKPDE. The presented method verifies the critical δ value of 2 obtained by Britz et al. [3] for the cylindrical geometry and activation parameter = 0.
It has been observed that numerical solutions to the steady equation [3] and those obtained from the transient equation [7] exhibited some discrepancies. In this work, we consider the convergence performance of four terminating criteria with high-order spatial discretisations for a Crank-Nicolson time-marching scheme, used for direct calculation of the steady-state solution. The motivation behind this work comes from uncertainty regarding the results obtained in [7] which centred around the use of higher-order discrete approximations. In particular, although central-difference five-point approximations to spatial derivatives were applied in the bulk of the propagation matrix, lower-order approximations were used at the boundaries. We now propose the use of asymmetric higher-order approximations at these points as a means of maintaining the order accuracy of the scheme. This is an important consideration given that due to symmetry, the temperature profile depicts the behaviour of the temperature at the centre of the vessel on one of the boundaries. Furthermore, we seek to clarify some issues surrounding the steady-state solution of the FKPDE and the termination criterion employed in the numerical schemes considered. We note here that the term 'convergence' is traditionally used in numerical analysis to describe the convergence of a numerical approximation to the true underlying exact solution in the limit of some discretisation parameters. In this work we interrogate the convergence of the numerical solution toward the steady-state solution for long time simulations and specifically the termination criteria employed dictating this convergence to some tolerance.
The FKPDE exhibits a vast range of dynamics given the parameters, however, the focus of this research is a comparative analysis of termination criteria which is conducted across three geometries (k = 0 for the large slab geometry, k = 1 for the cylindrical geometry and k = 2 for the spherical geometry) while taking the activation parameter to be 0. Then Equation (8) in which u is the normalised temperature, t the normalised time and x the spatial variable where 0 ≤ x ≤ 1 and δ is defined as the Frank-Kamenetskii parameter. We impose the boundary conditions and initial condition u(x, 0) = 0 .
The remainder of the paper is structured as follows. Section 2 discusses the candidacy of various methodologies for the solution of the Frank-Kamenetskii equation. Section 3 provides an analysis of the critical value of the Frank-Kamenetskii equation. Section 4 introduces the Crank-Nicolson scheme used in the structuring of a finite difference scheme for the solution of the FKPDE. Section 5 proposes several terminating criteria used for the numerical experiments conducted, with the results of these criteria discussed in Section 6. Finally, conclusions are drawn in Section 7.

Candidate Methodologies
This section presents some additional details of potential numerical methodologies amenable to the numerical solution of the Frank-Kamenetskii equation as well as supporting reasoning behind the selection of the implemented numerical method.
The Lawson-Euler [10] methodology implemented by Momoniat in [7] essentially treats the linear part of the Frank-Kamenetskii equation exactly by way of an integrating factor. The truncation error of the Lawson-Euler scheme is O( t) + O( x 4 ).
By replacing the time integration procedure of the Lawson-Euler scheme with a fourth-order Runge-Kutta scheme one may obtain (details of the implementation may be found in [5]) a scheme with truncation error of O( t 4 ) + O( x 4 ). This scheme enjoys a higher accuracy in the temporal discretisation at the cost of four additional evaluations at each time step.
One possible solution methodology distinct from finite difference schemes are pseudo-spectral methods. This class of methodologies is attractive due to their geometric convergence and hence it is generally acceptable to use fewer gridpoints while still obtaining accurate solutions. The applicability of a Chebyshev collocation spatial discretization was explored and it was found that due to the dense differentiation matrix the method can become computationally expensive. Moreover, as the number of collocation points increases the truncation error decreases but round-off errors accumulate and dominate the error [11,12]. In addition to the accumulation of these errors, instabilities in the scheme require an appropriately small choice of ∆t which hinders the convergence of the method.
The above methodologies produce time dynamic numerical simulations to the Frank-Kamenetskii equation. Given the nature of the equation and the potential interest in the steady-state solution to the equation, an accurate numerical solution to the steady solution may be readily obtained using the shooting method with Richardson's extrapolation, as conducted by Britz et al. [3]. Our interest is in the behaviour of time-marching schemes, and hence this method is not of interest in this work.
The discussed methodologies were implemented and compared for the Crank-Nicolson scheme (the details of which are discussed in Section 4.2). It was found that unconditional stability afforded by the Crank-Nicolson scheme rendered it most amenable to the long-time simulation of the Frank-Kamenetskii equation. Moreover the Crank-Nicolson implementation presented in this work enjoys O( t 2 ) + O( x 4 ) accuracy.

Frank-Kamenetskii Parameter δ
In this section we will begin by establishing a robust approach for computing the critical value of the Frank-Kamenetskii parameter δ, below which steady solutions are obtained, with arbitrary precision, to be used for the numerical simulations to follow in subsequent Sections. Beginning with the steady-state solution derived by Frank-Kamenetskii [2,4] for the large slab geometry, where a is the solution to the auxiliary equation we define so that (a, δ) pairs that satisfy f (a, δ) = 0 are of particular interest. Figure 1 visualises this function close to the critical value of δ, with a zero plane to illustrate the parabolic nature of δ near the criticality. Rearranging we find that satisfies f (a, δ) = 0, parameterising the parabolic shape seen in Figure 1. The value of a which maximises δ is defined in terms of the roots of a complicated polynomial, which is unfortunately intractable analytically. A more pragmatic approach is to maximise Equation (16) by solving for a when dδ da = 0, given explicitly by Finding the root of this function can be done numerically to arbitrary accuracy. Applying the bisection method with a tolerance of 10 −16 (machine precision) initially on the domain a ∈ [3, 4] we obtain in 54 iterations, corresponding to a = 3.2767175312280727.
Substituting the above values for a and δ into Equation (13) we can find the corresponding left-hand boundary value to be We note here that the critical value for δ was investigated by Britz et al. [3] and was accurately found to eight significant figures utilising the shooting method on the steady equation which is widely agreed upon in the literature. The result obtained in (18) is accurate to seventeen significant figures and the methodology does not require the discretisation of the differential equation, although it does require the exact steady solution. The exact critical value for δ in the cylindrical case is readily found to be 2 from the exact solution presented by Harley [6], and in the spherical case the critical value of δ = 3.3219921, as reported in [3], is used.

Crank-Nicolson
The Crank-Nicolson scheme is an average of the explicit and implicit finite difference spatial discretisation. The numerical scheme enjoys the stability properties of the implicit finite difference discretisation but also exhibits a higher order temporal accuracy than either the explicit or implicit finite difference methods. The numerical implementation selected for the experiments conducted in Section 6 is a Crank-Nicolson type scheme with fourth order spatial discretisations which ensures the scheme exhibits sufficient accuracy so as not to obfuscate the convergence behaviour.

Spatial Discretisation
In the following sections, our space domain is divided into N intervals of size 1/N and we make use of fourth-order discretisations for the spatial derivatives. The central difference discretisation is given as and Given the symmetric nature of the geometry under consideration we have employed a Neumann boundary condition at x = 0, which implies The left hand boundary approximations thus become and and at and We note here that the Neumann boundary condition discretisations (23) are second-order accurate hence the approximations at the boundaries also have second-order accuracy. At x N−1 = 1 − x we cannot invoke any physical symmetry, nor can we apply (22) directly since we will require values outside of the grid space. Here we employ a fifth-order asymmetric six-point approximation as follows, and which provides better accuracy, which a fourth-order five-point approximation would not [13]. The coefficients were obtained using the Fornberg algorithm [14]. Finally at x N = 1 we impose the condition u N = 0, noting that no derivative is needed there for the time march or direct solution.
We discretise Equation (9), using Equations (21), (22) and (24)-(29) to obtain, in vector form, and The value u N is set to zero because of the boundary conditions.

Implementation
Implementing the Crank-Nicolson scheme on Equation (30) produces where we note that u j+1 denotes the next, unknown, value. Incorporating the relevant boundary conditions as discussed in Section 4.1 and given u N = 0, we can rearrange Equation (33) to construct the following matrix-vector equation where with I the identity matrix and In this manner, we can implicitly treat the nonlinear source term resulting in a nonlinear system to be solved by Newton's method for each iterate. To achieve this we iteratively solve the linearised system given by where S is initialised to u j for each application of Newton's method and J is the Jacobian given by with The Newton iteration is stopped when the norm of ∆S < 10 −8 , generally, 3-4 iterations is sufficient due to the quadratic convergence of Newton's method. The truncation error associated with this Crank-Nicolson scheme is then O( t 2 ) + O( x 4 ).

Convergence Criteria
In this section, we investigate discrepancies in convergence criteria and critically analyse their applicability to the numerical solution of the FKPDE for which dependency on convergence is critical in defining the barrier between stable solutions and physical blow up. A common practice within the literature is to deem that a numerical simulation has "converged" when where tol refers to the chosen tolerance, is satisfied. Indeed, this is the criterion utilised by Momoniat [7]. The problem with Equation (40) is that if the order of the method used to determine u is O(∆t n ). Then Equation (40) approximately yields, Explicitly, if tol is 10 −6 then n > 5. If the method used doesn't guarantee this the termination criterion may yield early termination.
We propose that this inconsistency is attributed to the convergence condition and not the order of discretization. Rather, equations such as the Frank-Kamenetskii equation that exhibit slow temporal dynamics can trigger premature satisfaction of the above criterion. The convergence condition quoted in Equation (40) is derived from an O(∆t) approximation to du dt = 0 and is insufficient in obtaining the true steady solution. In this section, we introduce three new alternative convergence conditions in an attempt to address this early termination.
The first new criterion proposed is derived from a fourth-order backward difference approximation to the temporal derivative and is given by yielding the criterion Taking into account the slow temporal dynamics exhibited by the Frank-Kamenetskii equation, the second convergence criterion proposed is based on the spatial solution of the equation. This measure assesses how closely the current time step satisfies the steady equation and is given by Our final proposition for convergence is comparing the numerical solution at the current time step to the exact steady-state solution. While this condition is reliant on the known exact solution to the steady-state we use this in the current study as a benchmark against which the other convergence criteria may be assessed. This condition is given by where U(x) is the exact steady solution given by Equation (13) with x ∈ [0, 1] and discretised as previously stated.

Results
All results were obtained on a AMD Ryzen TM Threadripper TM 1950x CPU @ 3.4GHz, with 128 GB RAM @ 2400 MHz on Ubuntu 19.10 implemented in MATHEMATICA 12. Figure 2 serves as an aid for the reader to visualise the equation dynamics, this figure illustrates the convergent dynamics as time tends to infinity. Particularly, the dynamics of the left-hand boundary value which is indicative of the overall convergence of the solution. The Crank-Nicolson scheme derived in Section 4 was implemented for the large slab geometry and each of the four terminating criteria discussed in Section 5. The obtained results for each termination criterion are given in Tables 1-4. For consistency all experimental parameters were fixed, ∆x = ∆t = 1/200, δ = 0.8784576797812904 as described and calculated in Section 3, with only the termination criterion changing for various experiments. The tables presented below illustrate, column-wise, the tolerance for which the termination criterion is satisfied, the resulting left-hand boundary value, the given measure of convergence, number of iterations taken to satisfy that measure of convergence, absolute CPU time elapsed and finally the absolute distance between the simulated boundary and the exact steady boundary value. Table 1 illustrates the early termination exhibited by the traditionally [7] used measure of convergence. This measure reports 'convergence' to O(10 −8 ) when the solution is only accurate up to O(10 −3 ) when compared to the steady-state value at the left-hand boundary.
The second convergence criterion assessed improves on the traditional selection by measuring the rate of change to O(∆ 4 ). This choice again leads to premature termination of the numerical simulation as can be seen in the last row of the table, with the measure reporting convergence to O(10 −8 ) when the solution is only O(10 −4 ) accurate when compared to the steady-state value at the left-hand boundary.
The dynamics of the Frank-Kamenetskii equation give rise to extremely slowly evolving solutions. As such, the third measure of convergence is not based on measuring the rate of evolution but rather on 'how well' the approximate solution satisfies the steady equation. This yields a more accurate convergence criterion than the previous experiments, reporting convergence to O(10 −5 ) from the steady-state value at the left-hand boundary when the criterion indicates convergence to O(10 −8 ). We note here that in the event that the steady-state solution is unknown this convergence criterion may still be employed yielding a more robust measure of convergence when numerical solutions are evolving slowly.
Results for the final experiment are presented in Table 4. We note that columns three and six are commensurate, indicating this measure is absolute in measuring the distance of the solution from the steady-state. It is also worth noting that the infinity norm measures maximum absolute distance which, in this case, always occurs at the left-hand boundary supporting the exact agreement between columns three and six.
A visual comparison of the early termination is presented in Figure 3. This figure illustrates the left hand boundary value at the time of convergence as dictated by each of the four discussed metrics for convergence for varying tolerances. The results presented in this figure support the findings reported in Tables 1-4. Figure 4 presents the computational time and number of iterations taken to reach a prescribed tolerance for each of the convergence metrics. As expected, the computational time and number of iterations are highly correlated, with small discrepancies occurring at very low tolerances due to algorithmic overhead. Interestingly the slope of the relationship between the tolerance and number of iterations is the same for metrics one through three, whereas metric four has a steeper gradient (implying the need for proportionally more iterations to satisfy the next tolerance order) for lower tolerances and a much flatter gradient for higher tolerances. Again, this is evident in Table 4 where satisfying the higher tolerances occurs after fewer iterations when compared with satisfying lower tolerances (specifically requiring 10 million iterations moving from O(10 −4 ) to O(10 −5 ) and to O(10 −6 ) but only 1.2 million to then move to O(10 −7 ) and 100 thousand iterations to get to O(10 −8 )).  (43)) with k = 0.   (45)) with k = 0.  6.2. Cylindrical Geometry (k = 1) We now assess the convergence of the Crank-Nicolson scheme for the cylindrical geometry. Again, for consistency all experimental parameters were fixed, ∆x = ∆t = 1/1000, δ = 2, with only the termination criterion changing for various experiments. Tables 5-8 support the same findings as in the case of the large slab, specifically that a poor choice of convergence criterion leads to early termination of the numerical simulation. Table 5. Results for Metric 1 (Equation (40)) with k = 1.  Table 6. Results for Metric 2 (Equation (43)) with k = 1.  Table 7. Results for Metric 3 (Equation (44)) with k = 1.  (45)) with k = 1. The cylindrical geometry gives rise to interesting dynamics since the critical value of δ is known exactly. This value is asymptotic, below which stable solutions exist and above which blow up occurs. The δ value used for the other geometries is approximated numerically (to high precision) and tends to the true critical value from the left which ensures stable simulations. In the cylindrical geometry, since we are exactly at the critical value, numerical simulations will cross the threshold of stable computation after about 400,000 iterations, in the current experimental setup. This necessitates the use of a lower tolerance in Figure 5 for a meaningful comparison between the candidate metrics. As before we observe metrics 1 and 2 result in premature termination of the numerical simulation, while metric 3 most accurately predicts convergence to steady-state.   Figure 5. The convergence of terminating criteria over varying tolerance.

Spherical Geometry (k = 2)
Finally, we consider the convergence dynamics for the spherical geometry. Once more, all experimental parameters were fixed, ∆x = ∆t = 1/1000, δ = 3.3219921, with only the termination criterion changing for various experiments. Tables 9-11 illustrate the convergence of the Crank-Nicolson scheme under the candidate convergence criteria. To the best of the authors' knowledge no exact solution exists for the spherical geometry steady equation, as such we are only able to present results for the first three metrics. The final column in Tables 9-11 present the absolute distance between the transient solution and the steady solution reported by Britz et al. [3]. Figure 6 illustrates the critical temperature at convergence as dictated by the various metrics for varying tolerances.  Figure 6. The convergence of terminating criteria over varying tolerance.

Discussion
The results presented in this section highlight the 'subjectivity of convergence' introduced by the choice of criterion. Since all experiments were conducted with the same discretisation parameters, the number of iterations taken to reach convergence yield an absolute measure of premature termination. When comparing the first three metrics to the final metric, it is evident that the selected measure terminates prematurely by some margin (two million iterations in case three versus 25 million iterations in the final case). This discrepancy and sensitivity are important considerations when assessing the requirements of a numerical method for the solution of an equation which exhibits slowing temporal evolutions.
Moreover, the importance of such analysis is emphasized in the cases where no exact solution is known. In the absence of an exact solution, a naive choice of convergence criteria, in spite of a highly accurate, stable numerical method, can lead to erroneous results with no problematic behaviour in solution dynamics.

Conclusions
In this work, we have provided a semi-analytical framework for computing the critical value of δ in slab geometry to arbitrary precision of the FKPDE which is a key parameter in understanding the modelling process of the equation. This framework circumvents the need to simulate the differential equation directly but does rely on the known steady-state solution. This analysis aids the subsequent numerical experimentation around the "subjectivity" of a convergence criteria.
Moreover, this work explores the sensitivity of numerical simulations of the Frank-Kamenetskii equation in all three geometries and has shown that a poor choice of termination criteria in the time-dependent simulation can lead to erroneous results which have given rise to discrepancies in the results obtained in previous studies. The numerical solutions obtained are based on a fourth-order spatial discretisation Crank-Nicolson methodology which was selected from a variety of methodologies (discussed in Section 2) due to its stability properties as well as its purported accuracy of O(∆t 2 + ∆x 4 ).
While this work is concerned with the dynamics of the Frank-Kamenetskii equation, the analysis is applicable across a broad range of models which employ balanced forces in the limit and sensitivity to critical parameter choices. Many problems, for which an exact steady solution is intractable, require careful consideration of whether the chosen termination criterion of the numerical method is appropriate. Such is the purpose of this work: to highlight the need for this analysis which has been overlooked and can be detrimental to the obtained numerical solution with no obvious indication of erroneous behaviour.