#### Abstract

Reaction-diffusion-advection equations provide precise interpretations for many important phenomena in complex interactions between natural and artificial systems. This paper studies second-order semi-discretizations for the numerical solution of reaction-diffusion-advection equations modeling quenching types of singularities occurring in numerous applications. Our investigations particularly focus at cases where nonuniform spatial grids are utilized. Detailed derivations and analysis are accomplished. Easy-to-use and highly effective second-order schemes are acquired. Computational experiments are presented to illustrate our results as well as to demonstrate the viability and capability of the new methods for solving singular quenching problems on arbitrary grid platforms.

#### 1. Introduction

Nonlinear reaction-diffusion-advection equations have been playing an important role in interactions between natural and artificial systems. The partial differential equations provide precisely mathematical interpretations for numerous natural phenomena, such as diffusion of heat and energy, burning or quenching of fuels, energy concentrations, and transformations in different environments. Nonlinear partial differential equations also deliver popular computational tools to cell biology and cancer treatment plans [1–4].

Consider an -dimensional heat engine, . Let be the temperature distribution inside of its combustion chamber and be the initial temperature distribution due to the sparks. Assume that combustion occurs at the unit temperature. Further, let , be the reciprocal chamber size index. The use of variable reflects the change of chamber size, which is typical when different chemical fuels are selected. Then, an idealized combustion model can be comprised as a nonlinear reaction-diffusion-advection initial-boundary value problem:where , is its boundary, is the -dimensional gradient vector, the coefficient is positive, , and is the stochastic indicator for [5, 6].

The nonlinear source function is strictly increasing with respect to and

A typical example of such reaction term iswhere is the internal combustion index utilized [3, 5–7].

The solution of (1)–(4) is said to *quench* if there exists a finite such that

Such a value is called the *quenching time* which is a possible combustion-reaction time [7, 8]. It has been shown that a necessary condition for quenching to occur is

The study of singular reaction-diffusion-advection equations of the form (1) and (5) can be traced back to Hale and Kawarada’s pioneering work [9, 10]. It was observed that when and , there exists a *critical value* such that the solution of (1)–(3) and (5) quenches if . Verifications of the existence and uniqueness to more general quenching models can be found in numerous recent publications including [1, 5, 9, 11].

To solve the nonlinear differential equation problem (1)–(3) numerically is an interesting, yet challenging, multitask since we often do not know if the solution of (1)–(3) will quench or not until a proper solution procedure is taking place. In other words, in addition to approximating the solution , a numerical method must be capable of evaluating simultaneously correct critical value , quenching time , quenching location , and extremely sensitive and possibly unbounded derivative function [5, 6, 11]. This can be difficult since values of remain bounded and well-behaved until time reaches a certain neighborhood of , if it exists. Furthermore, we may also observe from properties (6) and (7) that, while grows exponentially, or faster than exponentially, in the aforementioned neighborhoods, the solution itself continuously grows but stays bounded throughout the computation until a quenching suddenly erupts [7, 9, 10]. These coexisted and very distinct characters make it extremely hard to design an effective scheme.

Due to their extremely important applications in nature and societies, numerous investigations, in both theory and computations, have been carried out for nonlinear quenching problems including (1)–(3) in recent years. Most of the existing algorithms are constructed either based on the reduced problems, that is, the stationary problems by removing the variable , or by using fixed mesh steps (cf. [4, 7, 8, 12–14] and references therein). In these procedures, critical values, quenching times, solution , and the rate-of-change function are often approximated incorrectly. Computational procedures are less efficient and less reliable, in particular when multidimensional modeling equations are considered. Although adaptive strategies have been particularly in favor due to their great flexibility and geometric accuracy in capturing quenching singularities involved, orders of accuracies of the spatial discretization of existing adaptive methods have been limited to one due to the unpredictability of the nonuniform grids generated by adaptations [6, 8, 12].

On the other hand, effective new approaches in high-order approximations via finite difference, spectral, or finite element methods, have been introduced for solving other nonlinear partial differential equations [2, 15, 16]. Investigations have also been extended to fractional order partial differential equation problems [16, 17]. Since it is the high dimension that causes a tremendous increase in the computational cost, modern splitting techniques capable of converting higher dimensional problems into sets of single dimensional subproblems have reached a new height [13, 18]. Rigorous numerical analysis has also been given on split adaptations [18]. These have motivated our study of highly applicable and effective higher-order finite difference methods.

Note that problems (1)–(3) can be numerically treated through proper second-order stable dimensional splitting [18, 19]. The decomposition may effectively reduce the total number of operations from to , where is the number of dimensions and is the number of internal mesh points used [7, 18, 20]. Thus, in this paper, we shall focus on second-order semi-discretization for a one-dimensional modeling equation with and . In the circumstance, problems (1)–(3) can be simplified to

We consider an arbitrary partition of the spatial domain , that is, gridsfor which . Denote . Apparently, such sets are -dependent when they are used for solving (8)–(10). However, due to the feature of semi-discretization, we prefer dropping time location indicators for the simplicity of notations. We further assume that in general . At any spatial location , we adopt usual notations or without confusion.

#### 2. Semi-Discretization Approximations

At the time level , we setfor any available index . Without loss of generality, we drop the superindex for simplicity. We have expansions of the spatial derivatives

*Strategy 1. *From the above, we must haveNote that due to the homogeneous condition (9). Substituting the above into (8)–(10), we acquire thatwhereThe band matrix is large in size since is usually large. But it is relatively sparse and can be handled conveniently by many existing software packages such as “sparse” subroutines in MATLAB. Dropping the truncation error term in (19), we obtain a second-order semi-discretization scheme for solving (8)–(10):Therefore, the solution of (8)–(10) can be readily approximated numerically through the solution of the nonlinear system of ordinary differential equations (19) and (20).

*Strategy 2. *The coefficient matrix *A* in (19) can be extremely stiff particularly due to the use of one-sided finite difference formulas (13) and (15). To avoid using the formulas, we may assume first that, in addition to the boundary values , we also have values calculated via some other methods that are at least of second order in accuracy. In this way, we may only use (14) for our second derivative approximations. The fully centralized formulas may hopefully ease the stiffness [19]. This yields our semi-discretized new approximationwhereDropping the truncation error term from the differential system, we obtain another second-order semi-discretization scheme for solving (8)–(10):Therefore, the solution of (8)–(10) can be approximated numerically through solving the nonlinear system of ordinary differential equations (23) and (24), given that values of *u*_{−1}, *u*_{n+2} can indeed be calculated successfully.

#### 3. Parameter Determinations

Now, let us investigate the issue and see if values can be computed through some direct expansions. To this end, we have

We observe readily that

For the simplicity in calculations, we may set . Thus,

To use the above two formulations for a quadratic order approximation, we need to balance coefficients related to and at internal spatial mesh points. Recall (13) and (15) and set

Recall (25). To find the coefficients , we let

It follows therefore

It turns out that

From the last two equations,

Hence,

Consequently,

Similarly, to find coefficients , we consider

Balancing the coefficients, we have

It follows subsequently that

From the last two equations,

Hence,and furthermore,

Now, recall the expansion (25). To determine coefficients in (29), we have

From the above, we obtain the following linear system:

The system can be conveniently solved digitally via any existing solution package. By the same token, for evaluating coefficients , in (31), we find that

The above expansion leads to the following non-homogeneous linear system:

Now, recall (14). To determine coefficients , we have

Therefore, to guarantee overall *second-order approximations* of the required function values for , we must ask thatwhich can be compressed into a matrix formwhere

Furthermore, for (13), we observe that

Therefore, the result is as follows:

It can be rewritten aswhere

Now, for (15), we have

This gives us the linear systemwhere

Based on the above investigations, we acquire the following result.

Theorem 1. *If matrices , and are nonsingular based on the selection of , then the coefficient in (23) is uniquely determined.*

We note that properties, such as the spectrums, of matrices and still remain to be studied based on arbitrary sets of at each temporal level of . Bear in mind that nonuniform mesh steps in are determined through particular adaptive procedure employed, such as the arc-length monitoring functions investigated in [7, 13, 17, 19].

#### 4. Simulation Experiments

Let and . We consider two typical two-dimensional reaction-diffusion-advection modeling problems of (1)–(3) associated with (5).

*Example 1. *We are interested in natural quenching situations [21] with