Second-Order Semi-Discretized Schemes for Solving Stochastic Quenching Models on Arbitrary Spatial Grids

Department of Mathematics and Center for Astrophysics, Space Physics and Engineering Research, Baylor University, Waco, TX 76798-7328, USA Department of Mathematics, School of Digital Technologies, Tallinn University, Narva Rd. 25, 10120 Tallinn, Estonia Departamento de Matemáticas y Fı́sica, Universidad Autónoma de Aguascalientes, Avenida Universidad 940, Ciudad Universitaria, Aguascalientes, Ags. 20131, Mexico


Introduction
Nonlinear reaction-diffusion-advection equations have been playing an important role in interactions between natural and artificial systems. e 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][2][3][4].
Consider an m-dimensional heat engine, m ≥ 1. Let u be the temperature distribution inside of its combustion chamber and u 0 be the initial temperature distribution due to the sparks. Assume that combustion occurs at the unit temperature. Further, let 0 < σ ≤ 1, t > 0, be the reciprocal chamber size index. e use of variable σ reflects the change of chamber size, which is typical when different chemical fuels are selected. en, an idealized combustion model can be comprised as a nonlinear reaction-diffusion-advection initial-boundary value problem: where D ⊂ R m , zD is its boundary, ∇ is the m-dimensional gradient vector, the coefficient a � a(x) ∈ C 1 (D) is positive, 0 ≤ u 0 ≪ 1, and X(u) is the stochastic indicator for u < 1 [5,6]. e nonlinear source function f(u) > 0 is strictly increasing with respect to 0 ≤ u < 1 and A typical example of such reaction term is where κ > 0 is the internal combustion index utilized [3,[5][6][7]. e solution u of (1)-(4) is said to quench if there exists a finite T σ such that Such a value T σ 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 e 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 m � κ � 1 and a ≡ 1, there exists a critical value σ * ≈ 0.65340 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 u, a numerical method must be capable of evaluating simultaneously correct critical value σ * , quenching time T * σ , quenching location x * , and extremely sensitive and possibly unbounded derivative function u t [5,6,11]. is can be difficult since values of max x∈D u t remain bounded and well-behaved until time t reaches a certain neighborhood of T * , if it exists. Furthermore, we may also observe from properties (6) and (7) that, while max x∈D u t grows exponentially, or faster than exponentially, in the aforementioned neighborhoods, the solution u itself continuously grows but stays bounded throughout the computation until a quenching suddenly erupts [7,9,10]. ese 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 t, or by using fixed mesh steps (cf. [4,7,8,[12][13][14] and references therein). In these procedures, critical values, quenching times, solution u, and the rate-ofchange function u t 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 highorder 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]. ese 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]. e decomposition may effectively reduce the total number of operations from O(n 3m ) to O(n 3 ), where m is the number of dimensions and n is the number of internal mesh points used [7,18,20]. us, in this paper, we shall focus on second-order semi-discretization for a one-dimensional modeling equation with m � 1 and a ≡ 1. In the circumstance, problems (1)-(3) can be simplified to We consider an arbitrary partition of the spatial domain [0, 1], that is, grids Apparently, such sets Ω n , H n are t-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 h k+1 ≠ h k , k � 0, 1, . . . , n. At any spatial location x k , we adopt usual notations ϕ(x k , t) � ϕ k (t) or ϕ(x k , t) � ϕ k without confusion.

Semi-Discretization Approximations
At the time level t > 0, we set for any available index k. Without loss of generality, we drop the superindex for simplicity. We have expansions of the spatial derivatives 2 Discrete Dynamics in Nature and Society Strategy 1. From the above, we must have Note that u 0 � u n+1 � 0 due to the homogeneous condition (9). Substituting the above into (8)-(10), we acquire that where ⋱ ⋱ ⋱ ⋱ ⋱ a n−2,−2 a n−2,−1 a n−2,0 a n−2,1 a n−2,2 a n−1,−2 a n−1,−1 a n−1,0 a n−1,1 c n,−3 c n,−2 c n,−1 c n,0 e band matrix A is large in size since n 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): Discrete Dynamics in Nature and Society erefore, 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.
e 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 u 0 � u n+1 � 0, we also have values u −1 , u n+2 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. e fully centralized formulas may hopefully ease the stiffness [19]. is yields our semidiscretized new approximation a n−2,−2 a n−2,−1 a n−2,0 a n−2,1 a n−2,2 a n−1,−2 a n−1,−1 a n−1,0 a n−1,1 a n,−2 a n,−1 a n,0 Dropping the truncation error term from the differential system, we obtain another second-order semi-discretization scheme for solving (8)-(10): erefore, 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.

Parameter Determinations
Now, let us investigate the issue and see if values u −1 , u n+2 can be computed through some direct expansions. To this end, we have 4 Discrete Dynamics in Nature and Society We observe readily that For the simplicity in calculations, we may set To use the above two formulations for a quadratic order approximation, we need to balance coefficients related to u 0 ′ , u 0 ″ and u n ′ , u n ″ at internal spatial mesh points. Recall (13) and (15) and set Recall (25). To find the coefficients β 1,0 , β 1,1 , β 1,2 , we let Discrete Dynamics in Nature and Society It follows therefore It turns out that From the last two equations, Hence, Consequently, Similarly, to find coefficients c 1,n−2 , c 1,n−1 , c 1,n , we consider Balancing the coefficients, we have It follows subsequently that From the last two equations, Hence, and furthermore, 6 Discrete Dynamics in Nature and Society Now, recall the expansion (25). To determine coefficients β 2,0 , β 2,1 , β 2,2 , β 2,3 , β 2,4 in (29), we have From the above, we obtain the following linear system: e system can be conveniently solved digitally via any existing solution package. By the same token, for evaluating coefficients c 2,n−4 , c 2,n−3 , c 2,n−2 , c 2,n−1 , c 2,n , in (31), we find that Discrete Dynamics in Nature and Society u n ″ � +c 2,n−4 u n − h n−4 + h n−3 + h n−2 + h n−1 u n ′ + e above expansion leads to the following non-homogeneous linear system: c 2,n−4 + c 2,n−3 + c 2,n−2 + c 2,n−1 + c 2,n � 0, c 2,n−4 h n−4 + h n−3 + h n−2 + h n−1 + c 2,n−3 h n−3 + h n−2 + h n−1 + c 2,n−2 h n−2 + h n−1 + c 2,n−1 h n−1 � 0, Now, recall (14). To determine coefficients a k−2 , a k−1 , a k , a k+1 , a k+2 , we have 8 Discrete Dynamics in Nature and Society u k ″ � a k−2 + a k−1 + a k + a k+1 + a k+2 u k erefore, to guarantee overall second-order approximations of the required function values for k � 1, 2, . . . , n, we must ask that a k−2 + a k−1 + a k + a k+1 + a k+2 � 0, Furthermore, for (13), we observe that (52)

Theorem 1. If matrices M k , M l , and M r are nonsingular based on the selection of H n , then the coefficient B in (23) is uniquely determined.
We note that properties, such as the spectrums, of matrices M k , M l , M r and A, B still remain to be studied based on arbitrary sets of H n at each temporal level of t. Bear in mind that nonuniform mesh steps in H n are determined through particular adaptive procedure employed, such as the arc-length monitoring functions investigated in [7,13,17,19].

Example 1.
We are interested in natural quenching situations [21] with σ � � 5 √ and We further set X(u) ≡ 1 for the simplicity in illustrations. is problem can be viewed as a simplification of the internal combustion when dual ignition sparks are utilized [2,3]. e partial differential equation problem is first split to one-dimensional subproblems via exponential splitting and then Strategies 1 and 2 are applied, respectively [8,14,18]. e semi-discretization mesh density controller is chosen to be n � 200, while the temporal discretization parameter τ � λh is fixed, where λ � (1/10) > 0 is the Courant-Friedrichs-Lewy constant, and h � max k h k is taken in an appropriate temporal level, since temporal adaptations are not a goal of the current study. Standard exponentially evolving grid adaption is adopted [6,12]. Cubic splines are used for passing data on each temporal level since spatial grids are moved. To see more clearly solution profiles, surface plots are only shown on internal mesh points without the fixed homogeneous boundary data (2).
It is observed that, due to the strong smoothness effect of the diffusion operator, the initial temperature distribution u 0 due to the dual ignition quickly converges to a smooth surface as time increases. e surface further forms a symmetric distribution with respect to the center of the spatial domain. Effects of the two sparks become hardly noticeable. e phenomenon is demonstrated in Figure 1 at the second, 10,000th, 20,000th, and 30,000th temporal levels, respectively. It has also been found that the numerical solution is monotonically increasing pointwise as t increases.
Another extremely interesting physical function to monitor in quenching-combustion situations is the temporal derivative, or rate of change function, u t [7,21]. e evolution of u t accelerates rapidly as time t is approaching the quenching time T * , which is approximately 69727τ ≈ 0.47624479. To see more details, we show the rate of change function u t in four temporal levels: T � 69600τ, 69700τ, 69720τ, and 69726τ under the same scale in Figure 2. A rapid increase of the peak values can be observed. ese consistently match the latest observations reported in [7,21].
To see more precisely profiles of u and u t , we plot evolutional trajectories of the maximal values and mean values of them for T ∈ (0, T * ) in Figure 3. It can be seen that the maximum of the temperature field derivative u t increases exponentially as t ⟶ T * , while the solution quenches peacefully. e rapid increase of averaging temperature in the combustion chamber is particularly important to engineering applications, since it indicates the success of the quick release of the chemical energy from fuel [3]. A y-semilogarithmic scale is used to show greater details of the phenomenon, especially for the derivative function max 0≤x,y≤1 u t (x, y, T).
In Figure 4, three-dimensional surfaces of u and u t , together with their x-(u/u t ) plane projections, immediately before the quench are given. It can be seen that while the peak temperature, which is at the center of the combustion chamber, gradually approaches the fuel ignition, distribution of u tends to occupy the entire space of the chamber. Original influence of the dual sparks disappears. On the other hand, however, rates of the temperature change, u t , increase exponentially to the infinity at the center. e semidiscretized scheme is highly stable and captures the phenomenon satisfactorily in about seventy thousand time steps.
Example 2. As a highly realistic extension, we consider another favorable natural quenching case with σ � � 5 √ and a much simplified four ignition spark initial condition [2,3]: We continue setting X(u) ≡ 1. Again, we consider corresponding one-dimensional subproblems via exponential splitting and then adopting Strategies 1 and 2, respectively. e same spatial mesh density controller and Courant-Friedrichs-Lewy constant are used for most of the computations until t is extremely close to the quenching time. Standard exponentially evolving grid adaption is adopted in the final stages of solution calculations [12]. Cubic splines are again used for passing data on each temporal level since spatial grids are moved. We only plot solutions on internal mesh points without showing the homogeneous boundary data (2). Figure 5 is devoted to four key moments of the solution procedures. Distributions of the temperature u and its temporal derivative u t inside the rectangular combustion chamber are simulated. e evidence of initial multiple sparks is obvious at T ≈ 0.00683013. However, the influence apparently becomes fuzzy at T ≈ 0.01366026 and fuzzier at T ≈ 0.48121712, while the maximal temperature, Discrete Dynamics in Nature and Society max 0≤x,y≤1 u(x, y, T), increases monotonically from 0.13965273 to 0.20441452 and then 0.97828162. We also notice that the maximal value locations move from the initial four sparks to the central of the spatial domain. e phenomenon can be seen more clearly in the last figure of u, for which T ≈ 0.48142203. is is further supported by the corresponding set of u t surfaces in Figure 5. e facts indicate a unique quenching/blow-up position in the combustion chamber which is physically correct [4,5,11].
In the final phase at T ≈ 0.48142203, values of u ≈ 1, which is the singular point of the source function f(u), in almost the entire combustion domain. On the other hand, we also observe that, while the distribution profile of u t is similar to that of u, the amplitude of rate function u t jumps to approximately 1.5 × 10 5 . is is an indication that, during the final ignition, the fuel combustion speed must be tremendously high, while the dimensionless temperature stays under the unity.
We show both the maximal and mean values of the numerical solution u and its temporal derivative u t , (x, y) ∈ D, in Figure 6. It is interesting that, while both mean values increase monotonically, the maximal derivative function decreases slightly in the beginning but quickly picks up the momentum and increases rapidly and becomes unbounded as t ⟶ T * . However, both maximal values remain positive throughout computations. e simulation features well agree with the monotone solution property predicted in the theory of thermal combustion [2,3].    Figure 4: Numerical solution u(x, y, T) (a, b) and its temporal rate of change function u t (x, t) (c, d), together with their respective x − u/u t plane projections, are given at the reference T ≈ 69727τ ≈ 0.47624479 immediately before quenching. We have max 0≤x,y≤1 u(x, y, T) ≈ 98952207, max 0≤x,y≤1 u t (x, y, T) ≈ 87.32452308. Again, plots are over internal mesh points without showing the boundary data (2). Finally, in Figure 7, profiles of the solution u and its temporal derivative u t are given at y � 0.5, T � 70483τ ≈ 0.48140837 immediately before quenching. It can be seen that while the temperature inside the chamber increases uniformly towards the unity, the rate-of-change function u t , in other words, the velocity of the temperature field evolution, grows rapidly to the infinite at the quenching location (x * , y * ) � (0.5, 0.5) ∈ D.  )

Conclusions and Future Endeavours
In this article, we have built and investigated two interconnected second-order finite difference schemes for the numerical solution of the stochastic quenching-combustion partial differential equation (1) equipped with Dirichlet boundary conditions and initial condition (2) and (3). is is for the first time, to the authors' best knowledge, to use a second-order finite difference scheme for effective approximations of quenching solutions on arbitrary spatial meshes. Numerical results obtained are not only consistent with existing results, but also offering more structure related details. Parameter determination procedures are discussed and analyzed.
Computer simulation experiments are focused on two examples with solid combustion-reaction backgrounds [1,5]. Evenly distributed multiple sparks are deployed in combustion chambers. It is found that both the temperature and its temporal derivative distributions increase monotonically inside the chamber once the ignition starts. e evolution of the temperature field remains smooth and bounded during computations. However, on the other hand, the derivative distribution in the chamber increases rapidly and becomes unbounded as t reaches the neighborhood of the quenching time T * . is is an indication of a physical combustion reaction of the fuel utilized.
e new second-order methods developed on arbitrary spatial grids are extremely straightforward and easy to use in programming implementations. ey capture successfully the unique quenching location of the solution of (1)-(3). e numerical results acquired correctly reflect the fact that the final quenching-combustion location is unique and in the center of the camber domain, rather than the multiple initial ignition spark positions suggested by initial functions.
Future endeavors in this direction involve continuing numerical analysis and the implementation of numerical algorithms for more rigorous stochastic structure coalitions and conditions. e authors will also pay special attention to the degeneracies, since these are indications of engine metal fatigue [12,14]. ey are extremely important in applications. Higher-order continuous and discontinuous Galerkin methods [15,16] may also be explored for solving singular quenching problems including (1)- (3). It is also the authors' intention to inspire more collaborations in this very meaningful and promising territory of advanced theory and numerical methods for the natural world.
Data Availability e data are available from the corresponding author upon any reasonable request.

Conflicts of Interest
e authors declare no conflicts of interest.