Hopscotch method: The numerical solution of the Frank-Kamenetskii partial differential equation

Applied Mathematics and Computation 217 (2010) 4065–4075

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Hopscotch method: The numerical solution of the Frank-Kamenetskii partial differential equation C. Harley Centre for Differential Equations, Continuum Mechanics and Applications, School of Computational and Applied Mathematics, University of the Witwatersrand, Johannesburg, Private Bag 3, Wits 2050, South Africa

a r t i c l e

i n f o

Keywords: Hopscotch scheme Thermal explosion Nonlinear source term Linear stability analysis

a b s t r a c t Numerical solutions to the Frank-Kamenetskii partial differential equation modelling a thermal explosion in a cylindrical vessel are obtained using the hopscotch scheme. We observe that a nonlinear source term in the equation leads to numerical difficulty and hence adjust the scheme to accommodate such a term. Numerical solutions obtained via MATLAB, MATHEMATICA and the Crank–Nicolson implicit scheme are employed as a means of comparison. To gain insight into the accuracy of the hopscotch scheme the solution is compared to a power series solution obtained via the Lie group method. The numerical solution is also observed to converge to a well-known steady state solution. A linear stability analysis is performed to validate the stability of the results obtained. Ó 2010 Elsevier Inc. All rights reserved.

1. Introduction In this paper we consider numerical solutions to the Frank-Kamenetskii [12] partial differential equation which is capable of modelling a thermal explosion in a cylindrical vessel. The heat balance equation, neglecting reactant consumption, is given by

c

  @T E : ¼ jr2 T þ rQA exp  @t RT

ð1Þ

In this formulation, c is the thermal capacity, j the thermal conductivity, r the density, Q the heat of the reaction, A the frequency factor, E the energy of activation of the chemical reaction, R the universal gas constant and T the gas temperature. For this kind of chemical process it must be noted that c, j, r and Q would likely evolve as the heat increases. However it has been noted in the work done by Frank-Kamenetskii [12] that Eq. (1) is a reasonable model for the physical process in question and able to provide us with reasonable insight into the temperature distribution. A thermal explosion occurs 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. The increase in the heat generated by the reaction increases the rate of the chemical reaction exponentially following the Arrhenius equation until a thermal explosion occurs. Eq. (1) was derived by Frank-Kamenetskii [12]. Further work was done by Steggerda [33] on Frank-Kamenetskii’s original criterion for a thermal explosion showing that a more detailed consideration of the situation is possible. Rice [31] has applied the theory developed by Frank-Kamenetskii for thermal explosions in which heat is removed by conduction only to the azomethane, ethyl azide and methyl nitrate explosions. As a precursor to this work it was concluded by Rice and Campbell [32] that the methyl nitrate explosion appears to not be a thermal explosion. Certain criticisms directed by Frank-Kamenetskii against the determination of the heat of reaction from the induction period, as carried out by Rice and Campbell [32], were considered and refuted. E-mail addresses: [email protected], [email protected] 0096-3003/$ - see front matter Ó 2010 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2010.10.020

4066

C. Harley / Applied Mathematics and Computation 217 (2010) 4065–4075

Eq. (1) is made non-dimensional by introducing the following dimensionless temperature difference

h¼

E RT 20

ðT  T 0 Þ;

ð2Þ

where T0 is the ambient temperature. The model Eq. (1) reduces to

  c @h h ; ¼ r2 h þ d exp k @t 1 þ h

ð3Þ

where

"

 # E exp  RT 0 RT 20

rQA E

d¼

k

ð4Þ

is known as the Frank-Kamenetskii [12] parameter and

¼

RT 0 : E

ð5Þ

Following the work of Frank-Kamenetskii [12] it is assumed that

r2 ¼

@2 N  1 @ ; þ r @r @r 2

  1. The Laplacian operator r2 is given by ð6Þ

where N = 1 for an infinite slab, N = 2 for an infinite circular cylinder and N = 3 for a sphere. By making a further change of variables the explicit dependence of the model Eq. (3) on the Frank-Kamenetskii parameter d can be removed

u ¼ h  h0 ; ðaÞ x ¼ r½deh0 1=2 ; ðbÞ;

ð7Þ

where h0 is a dimensionless temperature at the vessel walls. Then by rescaling time t, ignoring coefficients O() and taking N = 2 (3) is written as

  @u 1 @ @u ¼ x þ eu : @t x @x @x

ð8Þ

The model Eq. (8) is solved subject to

uðx; 0Þ ¼ constant ¼ h0 ;

ð9Þ

@u ð0; tÞ ¼ 0 uðR; tÞ ¼ h0 ; @x

ð10Þ

R ¼ ½deh0 1=2 :

ð11Þ

where

The boundary condition (11) is now dependent on d. The Eq. (8) describes a thermal explosion in a cylindrical vessel and the boundary conditions (9) and (10) imply that the temperature at the vessel walls is kept fixed and the solution is symmetric about the origin. Appropriate boundary conditions for a fixed constant temperature at the cylinder wall and symmetry about the centre of the cylinder are given by

uðx; 0Þ ¼ 0; @u ð0; tÞ ¼ 0; @x

ð12Þ ðaÞ uð1; tÞ ¼ 0;

ðbÞ;

ð13Þ

where R = 1 is the radius of the cylinder. In Fig. 1 a numerical solution to (8) solved subject to (12) and (13) obtained using pdepe in MATLAB is plotted. We note that as t increases the solution becomes singular at x = 0 tending toward blow-up. The numerical solution does not exhibit what happens after blow-up. Eq. (8) is a nonlinear diffusion equation with an exponential source term. In the case of a thermal explosion the Arrhenius law is maintained by the introduction of the exponential term which acts as a source for the heat generated by the chemical reaction. Frank-Kamenetskii [12] obtained a steady state solution to this problem. Zeldovich et al. [37] considered similarity solutions admitted by (8) that exhibit blow-up in finite time. These kinds of solutions, while noteworthy, have limited significance due to the restricted form of the initial profiles compatible with the similarity solutions. These solutions correspond to very special initial conditions for the temperature evolution profile, limiting the degree to which results obtained in this manner are applicable. This disadvantage has been noted by Anderson et al. [1] while analytically investigating the time evolution of the one-dimensional temperature profile in a fusion reactor plasma. A solution which also models blow-up in finite time has been obtained by Harley and Momoniat [20] via nonlocal symmetries of the steady-state equation.

4067

C. Harley / Applied Mathematics and Computation 217 (2010) 4065–4075

0.25 T = 0.3 0.2

T = 0.2

u

0.15 T = 0.1 0.1

0.05

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x Fig. 1. Numerical solution of (8) solved subject to (12) and (13).

In this paper a well-known numerical algorithm is introduced for solving the Frank-Kamenestkii [12] partial differential equation. The efficiency and stability of the method is investigated and compared to the pdepe function in MATLAB, the method of lines implemented in MATHEMATICA and the Crank–Nicolson scheme. The hopscotch scheme implemented in this paper is the explicit–implicit member of a class of hopscotch algorithms as defined by Gourlay [15], but will simply be referred to as the hopscotch scheme. It is shown that, even though the algorithm is modified due to the presence of a nonlinear term in the equation, the hopscotch scheme is applicable to the relevant equation and still able to perform well in comparison to the other numerical schemes employed. The hopscotch method has been discussed in depth by Gourlay [16] and the findings of this paper support his view that the method is a fast general-purpose algorithm, capable of producing answers with minimum machine time and effort. The Crank–Nicolson and hopscotch scheme indicate a similar level of efficiency in terms of running time, whereas in comparison with the pdepe function in MATLAB and the method of lines implemented in MATHEMATICA the hopscotch scheme computes much faster. As pointed out by Gourlay [16] another merit of the hopscotch scheme is that it is a scheme constructed for the solution of a general equation due to the fact that it appeals to properties of the mesh instead of properties of the differential equation. In contrast however the Crank–Nicolson scheme is very much problem-dependent. This paper shows that the hopscotch scheme performs as well as the Crank–Nicolson scheme when implemented for the Frank-Kamenestkii partial differential equation in most respects. However the numerical results obtained correlate with those presented by Feldberg [11], indicating that large values of Dt/Dx2 produce the problem of propagational inadequacy. As the value of Dt/Dx2 increases we find that the hopscotch scheme is outperformed by the Crank–Nicolson scheme indicating that the scheme is possibly not as accurate. This paper is set out as follows. In Section 2 a numerical solution is obtained for the model equation using the hopscotch scheme. This solution is compared with numerical solutions obtained via the pdepe function in MATLAB, the method of lines implemented in MATHEMATICA and the Crank–Nicolson scheme. In Section 3 a comparison is made between the numerical solution obtained through the hopscotch scheme and two analytical solutions: a power series solution obtained via the Lie group method and the steady state solution obtained by Frank-Kamenetskii [12]. In Section 4 a linear stability analysis is performed with concluding remarks made in Section 5. 2. Comparison of numerical solutions 2.1. Method of lines In this section we solve Eq. (8) numerically using the method of lines. The domain [0, 1] is sub-divided into N intervals of equal length, Dx, i.e. 0 = x0 < x1 < x2<  
ð14Þ

4068

C. Harley / Applied Mathematics and Computation 217 (2010) 4065–4075

 um @u um n1  nþ1 : @x 2Mx

ð15Þ

The time derivative is approximated by the forward-difference approximation

@u umþ1  um n  n ; @t Mt

ð16Þ

where um n ¼ uðxn ; t m Þ with xn = x0 + nDx and tm = t0 + mDt for n = 0, 1, 2, . . ., N and m = 0, 1, 2, . . ., M. The constant N refers to the number of intervals into which the domain is divided while M is the number of iterations. To impose zero-shear boundary conditions at the edges we approximate the spatial first derivative by the central-difference approximation (16) which leads to the following condition: m um 1 ¼ u1 :

ð17Þ

Additionally the boundary condition u(1, t) = 0 leads to uN(t) = 0. The fundamental idea behind the method of lines is the replacement of the spatial derivatives with approximations such as (14)–(16). This means that the spatial derivatives are no longer explicitly stated in terms of the spatial independent variables. This leaves us with a system of initial value ordinary differential equations to solve through numerical integration which should approximate the original partial differential equation. Consequently, considering Eq. (8), according to the method of lines we need to solve the following system of equations [25,26]:

u0 ðtÞ ¼ AuðtÞ þ euðtÞ

ð18Þ

where

2 6 6 6 6 6 6 A¼6 6 6 6 6 4

4j

4j

0

0

...

0

0

0

2j

4j

2j

0

...

0

0

0

0 .. .

2j .. .

4j 2j . . . .. .. . . ...

0 .. .

0 .. .

0 .. .

0 0 0

0 0 0

0 0 0

0 0 0

. . . 2j 4j 2j ... 0 2j 4j ... 0 0 0

0

3

0 7 7 7 0 7 7 .. 7 . 7 7 0 7 7 7 2j 5 1

with j ¼ Mx1 2 . If the system had been homogenous, solving (18) would simply be a matter of solving the matrix system which is generally structured as Cu0 (t) = Du(t) through explicit or implicit integration. Since we find that (18) is nonhomogeneous, each equation will be solved individually [26]. To implement the method of lines we use the in-built function NDSolve in MATHEMATICA. The system (18) is solved for xN = 1 where the minimum number of points for each dimension of the grid is chosen as N = 2000. As is usually the case the step-size, Dx, is related to the number of grid-points chosen, i.e. Dx = (xN  x0)/N. It is important to note that the interval over which we solve Eq. (8) leads to a singularity. In order to overcome this we structure 2 the code to account for two instances: x = 0 and x – 0. Using L’Hopitals rule we make the simplification ð1=xÞ @u  @@x2u at x = 0 @x for Eq. (8). Doing this we attain the following approximation:

@u @2u ¼ 2 2 þ eu @t @x

ð19Þ

to Eq. (8) at x = 0. This approximation has been taken into account in the matrix A of the system of equations represented by (18). This approximation is used in many numerical algorithms and has been implemented by Crank and Furzeland [8]. In [8] they presented a modified finite-difference method which eliminates inaccuracies that occur in the standard numerical solution near singularities. The approximation has also been used by Harley and Momoniat [21] to generate a consistency criteria for initial values at x = 0 for a Lane–Emden equation of the second-kind. In this instance it allows the NDSolve program in MATHEMATICA to solve the relevant equation while avoiding the singularity. The solution obtained by MATHEMATICA is an interpolating function which represents an approximate function whose values are found by interpolation within the relevant domains for x and t respectively. The solution found in this manner is evaluated on the interval [0, 1] at tM = T = 0.3 and will be compared to other numerical solutions which are yet to be discussed. 2.2. Crank–Nicolson implicit scheme We now consider an implicit Crank–Nicolson scheme for Eq. (8) subject to (12) and (13) to serve as a means of comparison against the hopscotch scheme. The nonlinear term remains explicit in the formulation of the scheme in this section. Choosing to maintain the explicit nature of the nonlinear term is a natural way of employing the Crank–Nicolson scheme when one allows for only the averaging of the derivatives over time. The reasoning behind this choice will be made clear in the next section when the hopscotch scheme is constructed.

4069

C. Harley / Applied Mathematics and Computation 217 (2010) 4065–4075

As clarified in Section 2.1 central-difference approximations are used for the spatial derivatives while a forwarddifference approximation is used for the time derivative. We implement a Crank–Nicolson scheme by approximating the second-derivative on the right-hand side of (8) by the implicit Crank–Nicolson [7] approximation m  2umþ1 þ umþ1 um  2um @ 2 u umþ1 n n þ un1 n1  nþ1 þ nþ1 : 2 2 2 @x 2Mx 2Mx

ð20Þ

In a similar fashion the first-derivative on the right-hand side becomes

 umþ1 um  um @u umþ1 n1 n1  nþ1 þ nþ1 : @x 4Mx 4Mx

ð21Þ

To impose zero-shear boundary conditions at the edges we approximate the spatial first-derivative by the central-difference approximation (16) which leads to the following condition: m um 1 ¼ u1 :

ð22Þ 2

Due to the fact that the x0 = 0 would lead to a singularity in Eq. (8) we use the fact that ð1=xÞ @u  @@x2u at x = 0 to obtain the @x approximation (19) to Eq. (8) as discussed previously. From this equation an initial condition for u(x, t) is obtained at x0 = 0 giving the following: m m1 ð1 þ 2bÞum þ 2bum1 ; 0  2bu1 ¼ ð1 þ 2bÞu0 1

ð23Þ

Mt , Mx2

where b ¼ as the initial difference scheme. Implementing the difference approximations discussed above we obtain the general numerical scheme



kn m cn m kn c þ ð1 þ bÞum ¼ um1 þ ð1  bÞum1 þ n um1 ; u u n  n 2 n1 2 nþ1 2 n1 2 nþ1

ð24Þ

    Mt 1 1 where xn = nDx and b ¼ Mx and kn ¼ b 1 þ 2n . This difference scheme (24), including the initial dif2 such that cn ¼ b 1  2n ference condition (23), form a system of equations which are to be solved iteratively. As indicated by the boundary conditions (13) we consider the problem for x 2 [0, 1] and t 2 [0, T]. The domain [0, 1] is subdivided into N equidistant intervals termed Dx, i.e. 0 = x0 < x1 < x2<  
Aum ¼ Bum1 þ Mteu

ð25Þ

and is solved as follows:

  m1 um ¼ ðAÞ1 Bum1 þ Mteu :

ð26Þ

The inverse of A is calculated using the inv function in MATLAB. The results produced are extremely similar to those obtained by the in-built function pdepe in MATLAB. This will be discussed further in a subsequent section. 2.3. Hopscotch scheme The hopscotch scheme is fashioned in such a way that the unknown grid points are obtained at two stages. The first stage involves solving (8) for all points on time level m such that (m + n) are even values using an explicit scheme. At odd values of (m + n) the points are calculated at time level m using an ‘effectively’ implicit scheme. By using the values newly calculated from the first stage we observe that the same explicit formulation used in stage two becomes effectively an implicit scheme. Diagrammatically the hopscotch scheme attains the solutions to the unknown points as pointed out by Fig. 2 where the crosses denote the points that were attained using the explicit scheme and the squares indicate points attained using the

X

X

X

Fig. 2. Diagrammatical representation of the hopscotch scheme.

4070

C. Harley / Applied Mathematics and Computation 217 (2010) 4065–4075

implicit scheme. In this manner no two points lying next to each other, in either direction of the axis, are determined only explicitly or implicitly. For more information regarding this method refer to [10,15–19,28,29]. To structure the hopscotch scheme as required we implement the difference approximations indicated by (14)–(16). As discussed in Section 2.2 we impose zero-shear boundary conditions (22) and use the approximate Eq. (19) to overcome the singularity of the equation. In this manner we obtain an initial difference scheme at x0 = 0 for the explicit and implicit cases respectively m1

m1 um þ 4bum1 þ Mteu0 ; 0 ¼ ð1  4bÞu0 1 m1

m m1 ð1 þ 4bÞum þ Mteu0 : 0  4bu1 ¼ u0

ð27Þ

Mt where b ¼ Mx 2 . These two initial difference Eq. (27) will be incorporated into a system of equations which will be solved iteratively. The system generated through the implementation of the hopscotch scheme on Eq. (8) subject to (12) and (13) is predominantly structured in two steps

(a) an explicit scheme m1 m1 um1  2um1 þ um1 um 1 um1 m1 n n1 nþ1  un1 n  un þ ¼ nþ1 þ eun xn Mt Mx2 2Mx

ð28Þ

for even values of (m + n), i.e. m and n must either both be even or both be odd, which leads to the difference scheme m1

m1 m1 un um þ kn um1 ; n ¼ cn un1 þ ð1  2bÞun nþ1 þ Mte

ð29Þ

    1 1 where xn = nDx and b ¼ such that cn ¼ b 1  2n and kn ¼ b 1 þ 2n (b) and an implicit scheme Mt Mx2

m m m1 um  2um um 1 um m1 n þ un1 nþ1  un1 n  un þ ¼ nþ1 þ eun 2 xn Mt Mx 2Mx

ð30Þ

for odd values of (m + n), i.e. when m is even then n must be odd and vice versa, which produces the difference scheme m1

m m m1 cn um þ Mteun n1 þ ð1 þ 2bÞun  kn unþ1 ¼ un

ð31Þ

with b, cn and kn are defined as before. Implementing the hopscotch scheme in the usual fashion would require the exponential term to be evaluated at time level m which causes numerical difficulty. To implement the implicit scheme we have m removed the complication by evaluating the nonlinear term eun at time level m  1 in Eqs. (30) and (31) instead of at time level m. We force the nature of the nonlinear term to remain explicit while performing the hopscotch scheme but not when implementing the Crank–Nicolson scheme. Maintaining an explicit nonlinear term is a natural way of implementing the Crank–Nicolson scheme as discussed in Section 2.2. In this way we are able to compare two schemes which both iterate linearly. Hence while a degradation of the performance for the hopscotch method may come about this is not the case for the Crank–Nicolson scheme. In fact in this manner the comparison between the two schemes can be said to be of some significance since both are the same ‘type’ of iterative scheme. A similar methodology with regard to the nonlinear term was employed by Gaidamauskaite_ and Baronas [13] and discussed further in Britz et al. [5] where the nonlinearities were in fact maintained while considering the Crank–Nicolson scheme through the use of a Newton iteration for the nonlinear term. This issue will not be considered in more depth in this paper due to the complication rendering the nonlinear term implicit presents. In not following the work of [5] we are in fact able to show that even though the hopscotch scheme has been modified, which may lead to poorer performance, the scheme can still be effectively employed for Eq. (8). For consistency we evaluate the system of equations obtained via the hopscotch scheme with Eqs. (29) and (31) using step-sizes Dx = 0.1, Dt = 0.0001 where x 2 [0, 1] and t 2 [0, T]. The explicit system given by (29) has the following structure: m1

um ¼ Bum1 þ Mteu

ð32Þ

and the implicit scheme (31) becomes

  m1 : um ¼ A1 Bum1 þ Mteu

ð33Þ

The inverse of A is calculated using the inv function in MATLAB. The system will iterate until tm + Dt = T which means that the system will iterate for M steps. At each iteration a column of values are determined in the grid using both the explicit and implicit scheme. If the mth column is odd then for odd n, i.e. odd valued rows, the explicit method (32) is employed and for even n, i.e. even valued rows, we use the implicit method (33). For an even-valued column, m, the inverse occurs – for n odd the implicit method is employed and n even the explicit method is implemented.

C. Harley / Applied Mathematics and Computation 217 (2010) 4065–4075

4071

2.4. Remarks The solutions to Eq. (8) subject to (12) and (13) obtained using the method of lines, Crank–Nicolson implicit scheme and hopscotch scheme show extreme similarity to Fig. 1 when represented graphically for T = 0.1, 0.2, 0.3. As a consequence we consider an error profile of the four schemes where we compare the pdepe function in MATLAB, the method of lines and the Crank–Nicolson scheme with the hopscotch scheme at T = 0.3. The differences between the solutions obtained from the various methods and the hopscotch scheme quite consistently indicate a decreased error profile as per Table 1. This may be a function of the boundary condition (13b) which is enforced by the numerical scheme. It was noticed that as T increases so does the overall error profile across the values of x for each difference tabulated. It is also clear that the solutions obtained by the Crank–Nicolson and hopscotch scheme are a close match whereas those obtained by the pdepe function in MATLAB and the method of lines implemented in MATHEMATICA are noticeably different from the solution obtained via the hopscotch scheme. The error profiles for the method of lines or Crank–Nicolson scheme and the hopscotch scheme indicate a consistency in the error components across the x values in terms of magnitude – for the method of lines the error predominantly remains in the fourth decimal place whereas for the Crank–Nicolson scheme it lies in the eighth and eventually ninth decimal place. The pdepe function in MATLAB however, indicates a larger error at the origin with an eventual decline as the x values increase. Due to the fact that the Crank–Nicolson and hopscotch scheme are both difference schemes it is difficult to ascertain whether the small error difference observed is indicative of the accuracy of the methods. It is beyond the scope of this paper to specify which of the numerical methods implemented have produced the most accurate results. In terms of performance time the hopscotch scheme is able to outperform the pdepe function in MATLAB and the method of lines implemented in MATHEMATICA and perform almost as well as the Crank–Nicolson scheme. To obtain average running times for the relevant methods we compute the running time using in-built functions in MATLAB and MATHEMATICA for each method for up to 30 computing instances at T = 0.3. This has indicated a consistent running time for each method over subsequent iterations. It is noticed that the pdepe function in MATLAB takes on average 0.6560 s at a relative tolerance of 103. This tolerance is a measure of the error relative to the size of each solution component of the vector u(x, t). In turn the absolute error tolerance, which determines the accuracy when the solution approaches zero, is chosen as 106. These tolerances influence the computing time of MATLAB – an increase in either of the tolerances leads to an increase in the time required to compute the solution. The method of lines employed in MATHEMATICA has a running time which is much higher – approximately 3.8 s. This is most likely influenced by the minimum amount of grid points for the spatial discretisation which was specified as 2000, however this alone cannot account for the computing time observed. The Crank–Nicolson scheme has indicated a shorter computing time than the hopscotch scheme. While the hopscotch scheme takes on average 0.2350 s to compute the Crank–Nicolson scheme takes an average of 0.2190 s. Both of these methods indicate a noticeable improvement on the other methods used to solve Eq. (8) in terms of machine time. Since the differences in running time and solutions of the Crank–Nicolson and hopscotch scheme are small it seems that the hopscotch method is able to perform at least as well as the Crank–Nicolson scheme as a means of solving Eq. (8).

3. Comparison with analytical solutions In this section a solution obtained to Eq. (8) using the Lie group approach will be compared to the numerical solution obtained above. The method employed will not be discussed in this section, however for further information regarding its implementation the reader may refer to [2–4,22,24,30,34]. Using the classical symmetry method Momoniat [27] obtained the following Lie point symmetries admitted by the model Eq. (8)

X1 ¼ @t ;

X 2 ¼ 2t@ t þ x@ x  2@ u ;

ð34Þ

Table 1 Absolute errors between the pdepe function in MATLAB, the method of lines implemented in MATHEMATICA (MATH), the Crank–Nicolson scheme (CNS) and the hopscotch scheme (HS) at T = 0.3. x

pdepe vs. HS:

MATH vs. HS:

CNS vs. HS:

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.002927 0.002952 0.001721 0.0010670 0.000661 0.000396 0.000222 0.000110 4.382801e05 1.018668e05 0

0.000514 0.000506 0.000481 0.000441 0.000388 0.000326 0.000257 0.000187 0.000118 5.504498e005 0

2.386976e08 1.163558e08 2.208926e08 1.020220e08 1.822757e08 7.809201e09 1.262279e08 4.755427e09 6.189647e09 1.520620 e09 0

4072

C. Harley / Applied Mathematics and Computation 217 (2010) 4065–4075

where @ t = @/@ t, @ x = @/@ x and @ u = @/@ u. By taking a linear combination of the symmetries a group invariant solution indicating the ‘‘blow-up” time T was found. Substituting this into the model equation produced a second-order ordinary differential equation from which a power series solution which converges on a small interval around the central axis of the cylinder was obtained. As a consequence it was shown that the temperature at the centre of the vessel is [27]

uð0; tÞ ¼ log



 T : T t

ð35Þ

In Fig. 3 we can see that the numerical solution obtained in the previous section matches Eq. (35) over a very small interval. It is specified by Momoniat [27] that the power series solution obtained satisfies the boundary condition (12) approximately on a small interval close to the origin. The boundary condition (13b) which has been implemented for the numerical scheme used in this paper was not imposed by Momoniat thus the restriction, R  1, on the radius of the cylinder was removed. These dissimilar boundary conditions are able to explain the divergence of the solutions. The numerical solution adheres to the condition at t = 0, (12), and on the boundary, (13b), unlike the power series solution which produces (35). This validates the numerical solution obtained in Section 2 since the solutions are similar close to the origin – the only region which is relevant due to the comments made by Momoniat [27] concerning the region on which the boundary conditions are satisfied. The error remains of order 104 until t = 0.002 at which stage the error starts to increase leading to a drastic divergence as noted in Fig. 3. Analytical solutions for Eq. (8) have been obtained by Frank-Kamenetskii for the steady state case of the temperature distribution. Both of the solutions he obtained are able to satisfy the condition y0 (0) = 0 [12,6,9]

" yðxÞ ¼ log

16ec1

ð2d þ ec1 x2 Þ2

# ;

" yðxÞ ¼ log

#

16ec2 ð1 þ 2dec2 x2 Þ2

ð36Þ

:

The constants c1 and c2 can be determined by imposing the boundary condition y(1) = 0 on (36)

h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii c1 ¼ log 2ð4  dÞ  4 2ð2  dÞ ; " pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi# 4  d  2 2ð2  dÞ : c2 ¼ log 2d2

ð37Þ

To attempt a comparison we consider the numerical solution of Eq. (8) obtained via the hopscotch scheme for a range of time intervals such that Te[0.1, 2], which allows us to observe the convergence of the numerical solution. We find that the system of solutions converges to a steady state solution as indicated in Fig. 4. In fact two of the analytical solutions, where ± is chosen as – for c1 and c2, obtained by Frank-Kamenetskii [12] match the solution observed through convergence closely with errors of order 104. This is also an indication that the scheme implemented for Eq. (8) is reasonably stable for Dx = 0.1 and Dt = 0.0001.

12

10

u

8

6

4

2

0

0

0.2

0.4

0.6

t

0.8

1

1.2

1.4

Fig. 3. Numerical solution of (8) at the centre of the vessel (———) compared to the temperature at the centre of the vessel given by (35) ((* * * *).

4073

C. Harley / Applied Mathematics and Computation 217 (2010) 4065–4075

0.35 0.3 0.25

u

0.2 0.15 0.1 0.05 0

0

0.1

0.2

0.3

0.4

0.5

x

0.6

0.7

0.8

0.9

1

Fig. 4. Convergence of numerical solutions of (8) obtained using the hopscotch numerical scheme for te[0.1, 1] and the solution (36) obtained by FrankKamenetskii (e  e  e  e).

4. Linear stability analysis To elaborate on the stability of this particular scheme we consider a linearised version of Eq. (8)

  @u 1 @ @u þ u þ 1; ¼ x @t x @x @x

ð38Þ

where we have used the approximation eu = 1 + u as per the Taylor series expansion since the stability analysis we are attempting to perform is only relevant for linear equations. Though this may seem to be a poor approximation, the magnitude of the numbers and the numerical range illustrated by the numerical solutions obtained above is an indication that we can reasonably assume that in this context the nonlinear exponential term is bounded appropriately, particularly close to the origin. The deviation of 1 + u from eu, where u is the numerical solution obtained for Eq. (8) using the hopscotch scheme, has been noted to grow with time and as a consequence the value of T impacts on the maximum error observed. For T = 0.1 the maximum absolute error obtained is jeu(x,t)  (1 + u(x, t))jT=0.1 = 0.0053 at x = 0, whereas for T = 0.3 the error increases to jeu(x,t)  (1 + u(x, t))jT=0.3 = 0.0283 at x = 0. Hence though not an error-free approximation we still deem it reasonable enough to employ. We consider the work by Hundsdorfer and Verwer [23] on the linear stability of the hopscotch scheme in the hope of determining whether the scheme has a bounded error propagation. As per [23] this will be true if and only if

sup kT n k 6 c

ð39Þ

nP0

for some c  (0, 1). The close relation this condition holds with the spectral radius q() allows us to use the inequality q() < 1 as sufficient for (39) [23]. To be able to structure the matrix T we consider Eq. (8) in the following form:

@u @u @2u ¼ a þ b 2 þ u þ 1; @t @x @x

ð40Þ

where the parameters a and b are defined as 1/x and 1 respectively. Eq. (40) can be rewritten as

d uðtÞ ¼ AuðtÞ; dt

t P 0;

uð0Þ ¼ u0 ;

ð41Þ

where the matrix A is given by

A¼

a b ðE  ET Þ þ 2 ðE  2I þ ET Þ þ I 2Mx Mx

ð42Þ

for a central-difference approximation with I = diag(1, 1, 1, . . ., 1, 1, 1) and E defined as the shift operator

0

0 1 0 ... 0 0

0

1

B0 0 1 ... 0 0 0C C B C B. . . . .. .. .. . . ... ... ... C: E¼B C B C B @0 0 0 ... 0 0 1A 0 0

0 ... 0 0

0

ð43Þ

4074

C. Harley / Applied Mathematics and Computation 217 (2010) 4065–4075

Table 2 Comparison of the spectral radii for different values of Dt.

Dt:

0.0001

0.00015

0.001

0.0015

0.01

0.015

Spectral radius:

0.9978

0.9967

0.9780

0.9671

0.7756

0.5643

Using the hopscotch scheme as discussed in [23] with I1 and I2 defined as two m  m diagonal matrices, whose diagonal elements are either 0 or 1 such that I1 + I2 = I (with I the unit matrix), we obtain the following:

T ¼ ðI  A1 Þ1 ðI þ A2 ÞðI  A2 Þ1 ðI þ A1 Þ:

ð44Þ

The conditions for stability of the hopscotch scheme proposed in this section are governed by the matrix T, indicated in Eq. (44), and hence using Eq. (39) and the consequent spectral radius condition we are able to comment on the stability of said scheme. Table 2 shows the values for the spectral radius over a range of values for Dt. We see that the spectral condition is maintained over the range of values for Dt with Dx = 0.1. From these results it is reasonable to maintain that the hopscotch scheme employed in this paper, with specific values assigned to Dx and Dt, has a bounded error propagation. Mt In the numerical work performed thus far the value for b, where b ¼ Mx 2 , has remained small – more specifically it was chosen to be 0.01 for Dt = 0.0001 and Dx = 0.1. We now implement the scheme for larger values of b as a check for the performance. It has been pointed out by Feldberg [11] that large values of b produce the problem of propagational inadequacy which leads to inaccuracies. This is a disturbing result due to the fact that an advantage of the scheme was that it was thought to be able to accommodate large values of b, more particularly large time-intervals, while maintaining stability. Within the context of this work however it was observed that for b-values of 0.2, 0.4 and 0.6 the system produces accurate results, attaining the same height and similar shape – dependent on the value of Dx – to those solutions obtained for a small value of b. However for b = 0.8 the height obtained is higher than that obtained for smaller values of b, indicating – as would be expected – that if b were to increase drastically it would lead to poor performance. This however was not observed for the Crank–Nicolson scheme which does not develop the problem of propagational inadequacy and hence presented accurate results even as b increased. 5. Concluding remarks A numerical solution was obtained for the Frank-Kamenetskii partial differential equation using the hopscotch scheme. The solution closely matched solutions obtained via the in-built MATLAB function pdepe, the method of lines implemented in MATHEMATICA and in particular the Crank–Nicolson implicit difference scheme. We have shown, even though a modification of the hopscotch scheme is employed by evaluating the nonlinear term explicitly instead of implicitly as is usually required, that the hopscotch scheme can be implemented with confidence for the solution of Eq. (8). It is found that the hopscotch scheme performs as well as the Crank–Nicolson scheme when the value of b remains below 0.6 and outperforms the pdepe function in MATLAB and method of lines in terms of running time. The numerical solution was also compared to a power series solution obtained via Lie group theory which indicated a close match near the origin. The steady state solution obtained by Frank-Kamenetskii [12] was then compared to the convergent numerical solution obtained by considering a range of T-values. Both of these comparisons indicated that the numerical solution obtained in this paper is relevant and capable of producing a reasonable model for the process of temperature distribution. Analysing the hopscotch scheme is complicated due to the nature of the scheme and the addition of a nonlinear source term in Eq. (8) only compounds the problem. As discussed in [23] unfortunately desirable stability inequalities are very hard to obtain for the hopscotch scheme even in relatively simple cases. Work has been done by [10,35,23] with particular emphasis on the linear stability of the scheme. Further discussions can be found in [14,15,29,36]. In an attempt to comment on the stability of the hopscotch numerical scheme used in this paper we considered a linearised version of the equation and showed that it maintained the condition that the spectral radius of the coefficient matrix be bounded by one indicating bounded error propagation. It is important to note that the stability results obtained in this paper were for specific values of Dx and Dt. It was shown that for reasonable values of b the hopscotch scheme behaves well, however at some point as the value of b increases the scheme would start producing inaccurate results. A thorough analysis needs to be done to ascertain the domain of b which maintains a stable scheme. This will be a topic of further research. References [1] D. Anderson, H. Hamnén, M. Lisak, T. Elevant, H. Persson, Transition to thermonuclear burn in fusion plasmas, Plasma Phys. Contr. Fusion 33 (10) (1991) 1145–1159. [2] G. Baumann, Symmetry Analysis of Differential Equations with Mathematica, Springer-Verlag, Berlin, 2000. [3] G.W. Bluman, S. Kumei, Symmetries and Differential Equations, Springer, Berlin, 1989. [4] G.W. Bluman, J.D. Cole, J. Math. Mech. 18 (1969) 1025. [5] D. Britz, R. Baronas, E. Gaidamauskaite˙, F. Ivanauskas, Further comparisons of finite difference schemes for computational modelling of biosensors, Nonlinear Anal.: Model. Contr. 14 (4) (2009) 419–433.

C. Harley / Applied Mathematics and Computation 217 (2010) 4065–4075

4075

[6] P.L. Chambré, On the solution of the Poisson–Boltzmann equation with application to the theory of thermal explosions, J. Chem. Phys. 20 (1952) 1795– 1797. [7] J. Crank, E. Nicolson, A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type, Proc. Camb. Phil. Soc. 43 (1947) 50–67. [8] J. Crank, R.M. Furzeland, The treatment of boundary singularities in axially symmetric problems containing discs, J. Inst. Math. Appl. 20 (3) (1977) 355– 370. [9] L. Dresner, Phase-plane analysis of nonlinear, second-order, ordinary differential equations, J. Math. Phys. 12 (1971) 1339–1348. [10] D.J. Evans, A. Danaee, A new group Hopscotch method for the numerical solution of partial differential equations, SIAM J. Numer. Anal. 19 (3) (1982). [11] S.W. Feldberg, J. Electroanal. Chem. 222 (1987) 101. [12] D.A. Frank-Kamenetskii, Diffusion and Heat Transfer in Chemical Kinetics, Plenum Press, New York, 1969. [13] E. Gaidamauskaite˙, R. Baronas, A comparison of finite difference schemes for computational modelling of biosensors, Nonlinear Anal.: Model. Contr. 12 (3) (2007) 359–369. [14] P. Gordon, Nonsymmetric difference equations, SIAM J. Appl. Math. 13 (1965) 667–673. [15] A.R. Gourlay, Hopscotch: a fast second order partial differential equation solver, J. Inst. Math. Appl. 6 (1970) 375–390. [16] A.R. Gourlay, Some recent methods for the numerical solution of time-dependent partial differential equations, Proc. Royal Soc. London A 323 (1971) 219–235. [17] A.R. Gourlay, Splitting methods for time dependent partial differential equations, in: D. Jacobs (Ed.), The State of the Art in Numerical Analysis, Academic Press, New York, 1977, pp. 757–791. [18] A.R. Gourlay, J.LI. Morris, Hopscotch difference methods for nonlinear hyperbolic systems, IBM J. Res. Develop. 16 (1972) 349–353. [19] A.R. Gourlay, Recent developments of the hopscotch method, in: R. Vichnevetsky (Ed.), Advances in Computer Methods for Partial Differential Equations II, IMACS, New Brunswick, NJ, 1977, pp. 1–6. [20] C. Harley, E. Momoniat, Steady state solutions for a thermal explosion in a cylindrical vessel, Mod. Phys. Lett. B (MPLB) 21 (14) (2007) 831–841. [21] C. Harley, E. Momoniat, Instability of invariant boundary conditions of a generalized Lane–Emden equation of the second-kind, Appl. Math. Comput. 198 (2008) 621–633. [22] A.K. Head, Comput. Phys. Commun. 77 (1993) 241. [23] W.H. Hundsdorfer, J.G. Verwer, Linear Stability of the Hopscotch scheme, Appl. Numer. Math. 5 (1989) 423–433. [24] N.H. Ibragimov, CRC Handbook of Lie Group Analysis of Differential Equations, vol. 1, CRC Press Inc., Boca Raton, Florida, 1994. [25] M.N. Mikhail, On the validity and stability of the method of lines for the solution of partial differential equations, Appl. Math. Comput. 22 (1987) 89– 98. [26] E. Momoniat, R. McIntyre, R. Ravindran, Numerical inversion of a Laplace transform solution of a diffusion equation with a mixed derivative term, Appl. Math. Comput. 209 (2009) 222–229. [27] E. Momoniat, A thermal explosion in a cylindrical vessel: a non-classical symmetry approach, Int. J. Mod. Phys. B 23 (14) (2009) 3089–3099. [28] G.R. McGuire, M.Sc. Thesis, University of Dundee, 1970. [29] G.R. McGuire, J. LI. Morris, Explicit–implicit schemes for the numerical solution of nonlinear hyperbolic systems, Math. Comput. 29 (130) (1975) 407– 424. [30] L.V. Ovsiannikov, Group Analysis of Differential Equations, Academic Press, New York, 1982. [31] O.K. Rice, The role of heat conduction in thermal gaseous explosions, J. Chem. Phys. 8 (1940) 727–733. [32] O.K. Rice, H.C. Campbell, The explosion of ethyl azide in the presence of diethyl ether, J. Chem. Phys. 7 (1939) 700–709. [33] J.J. Steggerda, Thermal stability: an extension of Frank-Kamenetskii’s theory, J. Chem. Phys. 43 (1965) 4446–4448. [34] J. Sherring, A.K. Head, G.E. Prince, Dimsym and LIE: symmetry determining packages, Math. Comput. Mod. 25 (1997) 153–164. [35] G.L.G. Sleijpen, Strong Stability Results for the Hopscotch Method with Applications to Bending Beam Equations, Utrecht. [36] T.R. Taha, M.J. Ablowitz, Analytical and numerical aspects of certain nonlinear evolution equations II, Numer. Nonlinear Schrodinger Equat. Reprinted J. Comput. Phys., vol. 55, Academic Press, New York and London, 1984. 2. [37] Y.B. Zeldovich, G.I. Barenblatt, V.B. Librovich, G.M. Makhviladze, The Mathematical Theory of Combustion and Explosions, Consultants Bureau, New York, 1985.