Numerical solution of degenerate stochastic Kawarada equations via a semi-discretized approach

TL;DR

This paper develops a semi-discretized numerical method with exponential splitting for solving complex degenerate stochastic Kawarada equations, ensuring stability and key solution features on adaptive grids.

Contribution

It introduces a novel semi-discretized approach combined with exponential splitting for degenerate stochastic PDEs, preserving positivity and stability.

Findings

Key solution features are preserved under modest restrictions.

The splitting method maintains stability without extra restrictions.

Numerical experiments demonstrate convergence and effectiveness.

Abstract

The numerical solution of a highly nonlinear two-dimensional degenerate stochastic Kawarada equation is investigated. A semi-discretized approximation in space is comprised on arbitrary nonuniform grids. Exponential splitting strategies are then applied to advance solutions of the semi-discretized scheme over adaptive grids in time. It is shown that key quenching solution features including the positivity and monotonicity are well preserved under modest restrictions. The numerical stability of the underlying splitting method is also maintained without any additional restriction. Computational experiments are provided to not only illustrate our results, but also provide further insights into the global nonlinear convergence of the numerical solution.

Figures

Equations

$$\mathcal{L}u=f(\epsilon ,u),(x,y,t)\in \Omega ,$$

$$u(x,y,t)=0,(x,y,t)\in \text{SS},$$

$$u(x,y,t_0)=u_0(x,y)\ll c,(x,y)\in \mathcal{D},$$

$$\mathcal{L}:=\sigma (x,y){\partial }_t-\Delta ,$$

$$f(\epsilon ,0)=f_0(\epsilon )>0,f_u(\epsilon ,u)>0,\underset{u\to c^{-}}{\lim }f(\epsilon ,u)=\infty ,{\int }_0^{c}f(\epsilon ,u)du=\infty$$

$$0<f_{\min }^{\epsilon }\le \underset{\epsilon }{\min }\{f(\epsilon ,u)\}\text{mbox}and\underset{\epsilon }{\max }\{f(\epsilon ,u)\}\le f_{\max }^{\epsilon }\ll \infty$$

$$\sup \left\{u_t(x,y,t):(x,y)\in \mathcal{D}\right\}\to \infty \text{mbox}ast\to T_q^{-}.$$

$$\max \left\{u(x,y,t):(x,y)\in \bar{\mathcal{D}}\right\}\to c^{-}\text{mbox}ast\to T_q^{-}.$$

$$u_t=\frac{1}{a^2\sigma }{u}_{xx}+\frac{1}{b^2\sigma }{u}_{yy}+g(\epsilon ,u),(x,y,t)\in \Omega ,$$

$$u(x,y,t)=0,(x,y,t)\in \text{SS},$$

$$u(x,y,t_0)=u_0(x,y)\ll 1,(x,y)\in \mathcal{D},$$

$${\mathcal{D}}_h$$

$${\mathcal{D}}_h^{\circ }$$

$${\frac{{\partial }^{2}u}{\partial x^2}|}_{(x_i,y_j,t)}$$

$${\frac{{\partial }^{2}u}{\partial y^2}|}_{(x_i,y_j,t)}$$

$$v^{\prime }(t)$$

$$v(t_0)$$

$$M_x=\frac{1}{a^2}B(I_{N_y}\otimes T_x),M_y=\frac{1}{b^2}B(T_y\otimes I_{N_x}),M=M_x+M_y,$$

$$B$$

$$T_{\alpha }$$

$$l_{\alpha ,j}$$

$$m_{\alpha ,j}$$

$$v(t)=E(tM)v_0+{\int }_{t_0}^{t}E((t-s)M)g(\epsilon ,v(s))ds,t_0<t<T_q,$$

$$v_{k+1}=E({\tau }_{k}M)v_k+{\int }_{t_k}^{t_{k+1}}E((t_{k+1}-s)M)g(\epsilon ,v(s))ds,k=0,1,\dots$$

$${\int }_{t_k}^{t_{k+1}}E((t_{k+1}-s)M)g(\epsilon ,v(s))ds-M^{-1}[E({\tau }_{k}M)-I][\theta g(\epsilon ,v_k)+(1-\theta )g(\epsilon ,v_{k+1})]=\mathcal{O}\left({\tau }_k^p\right),$$

$${\eta }_k={\int }_{t_k}^{t_{k+1}}E((t_{k+1}-s)M)g(\epsilon ,v(s))ds-{\int }_{t_k}^{t_{k+1}}E((t_{k+1}-s)M)[\theta g(\epsilon ,v_k)+(1-\theta )g(\epsilon ,v_{k+1})]ds.$$

$${\eta }_k=\frac{{\tau }_k}{2}\left[g(\epsilon ,v_{k+1})+E({\tau }_{k}M)g(\epsilon ,v_k)\right]-M^{-1}[E({\tau }_{k}M)-I][\theta g(\epsilon ,v_k)+(1-\theta )g(\epsilon ,v_{k+1})]+\mathcal{O}\left({\tau }_k^3\right),$$

$${\eta }_k=\frac{{\tau }_k^2}{2}{g}_{v_k}(\epsilon ,v_k)-(1-\theta ){\tau }_k^2{g}_{v_k}(\epsilon ,v_k)+\mathcal{O}\left({\tau }_k^3\right),{\tau }_k\to 0^{+},$$

$$v_{k+1}$$

$$x_i=g^{(x)}({\omega }_i^{(x)}),i=0,\dots ,N+1$$

$${\psi }_{i,j}^{(x)}(t)$$

$$(x_{i}-x_{i-1})-(x_{i+1}-x_{i})=-k_{1}^{2}\frac{x_{i-1}-2x_{i}+x_{i+1}}{k_{1}^{2}}=-k_{1}^{2}\frac{d^{2}g^{(x)}}{d(\omega^{(x)})^{2}}\bigg{|}_{\omega^{(x)}_{i}}+{\cal{O}}(k_{1}^{2}).$$

$$\displaystyle h_{x,i-1}=x_{i}-x_{i-1}=g^{(x)}(\omega_{i})-g^{(x)}(\omega_{i-1})=\frac{\partial g^{(x)}}{\partial\omega}\bigg{|}_{\xi_{i}}k_{1}$$

$$E({\tau }_{k}M)=E({\tau }_k(N-\beta I))=E(-{\tau }_k\beta I)E({\tau }_{k}N).$$

$$\mu (A)=\underset{h\to 0^{+}}{\lim }\frac{\Vert I_K+hA\Vert -1}{h}.$$

Full text

††journal: Applied Mathematics and Computation

1 Introduction

Kawarada partial differential equations arise in the modeling of highly sophisticated, yet important, natural phenomena where singularities may develop as the solution evolves in time. Such singularities often result from the energy of a system concentrating and approaching its activation criterion [1, 2].

Consider a typical solid-fuel ignition process. If the combustion chamber with fuel and air are appropriately mixed, the temperature in the chamber may increase monotonically until a certain critical value is reached. However, rates of such temperature changes can occur in a nonlinear manner throughout the media. This nonuniform distribution of heat may result in high temperatures being extremely localized within the combustion chamber and lead to an ignition once the peak temperature reaches a certain threshold [1]. This phenomenon is carefully characterized by a Kawarada model in which temporal derivatives of the solution may grow at an explosive rate, while the solution itself remains bounded [3, 4, 5, 6]. This strong nonlinear singularity, which is referred to as a quenching singularity, has been the backbone of the modeling equations. Applications of Kawarada equations can be found in fuel combustion simulations, enhanced thermionic emission optimization processes, electric current transients in polarized ionic chambers, and enzyme kinematics.

Let $\displaystyle{\cal{D}}=(-a,a)\times(-b,b)\subset{\mathbb{R}}^{2}$ be an idealized combustion chamber, $\displaystyle\partial{\cal{D}}=\bar{{\cal{D}}}\setminus{\cal{D}}$ be its boundary, and $\displaystyle\Omega={\cal{D}}\times(t_{0},T_{q}),~{}\SS=\partial{\cal{D}}\times(t_{0},T_{q}),$ where $\displaystyle a,b>0,~{}0\leq t_{0}<T_{q}<\infty$ with $\displaystyle T_{q}$ being the terminal time. Further, let $\displaystyle c>0$ be the ignition temperature. If $\displaystyle u(x,y,t)$ denotes the temperature distribution within the chamber, we have the following 2D model generalized from the original Kawarada configuration [4]:

[TABLE]

where

[TABLE]

and $\displaystyle\Delta$ is the standard 2D Laplacian. The degeneracy function $\displaystyle\sigma(x,y)\geq 0$ for $\displaystyle x\in\bar{{\cal{D}}}$ with equality occurring only for $\displaystyle(x,y)\in\partial{\cal{D}}_{0}\subseteq\partial{\cal{D}}.$ The nonlinear reaction term, $\displaystyle f(\varepsilon,u)$ with a bounded stochastic influence $\displaystyle\varepsilon(x,y),$ satisfies the following properties

[TABLE]

for $\displaystyle u\in[0,c).$ The stochastic influence is characterized by

[TABLE]

for each $\displaystyle u\in[0,c).$ The solution $\displaystyle u$ of $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ is said to quench if there exists $\displaystyle T_{q}<\infty$ such that

[TABLE]

$\displaystyle T_{q}$ is then referred to as the quenching time [3, 7, 8]. It has been shown that a necessary condition for quenching to occur with the above conditions placed on the reaction term is

[TABLE]

We may note that quenching is a mathematical interpretation of the fuel ignition/combustion process in physics. Since $\displaystyle T_{q}<\infty$ only when certain spatial references, such as the size and shape of $\displaystyle{\cal{D}},$ reach their critical limits, a domain $\displaystyle\Omega_{*}$ is called a critical domain of $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ if the solution of $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ exists globally as $\displaystyle\Omega\subseteq\Omega_{*}$ and $\displaystyle(\ref{a2})$ occurs as $\displaystyle\Omega_{*}\subseteq\Omega.$ We note that such critical domains are not unique in multi-dimensional circumstances, since both the size and shape of $\displaystyle\Omega$ affect combustion [9, 10, 11, 12].

Highly effective and efficient algorithms for solving Kawarada equations such as $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ have been difficult to obtain due to the highly nonlinear mechanism. There has been much effort to develop algorithms that may accurately predict the ignition location and time, solution profiles, and critical domain characteristics, while reasonably conserving the solution positivity and monotonicity [12, 13, 14, 15, 16, 17]. While explorations have been carried out on both uniform and nonuniform grids, as well as via adaptations, rigorous analysis of different schemes involving nonuniform grids has been incomplete. To improve computational efficiency, splitting strategies are introduced and incorporated with adaptive mesh structures [13, 18, 11]. The success of such combinations has also become a key motivation for the study to be presented in this paper.

Rigorous mathematical analysis is implemented for a highly vibrant semi-discretization oriented method for solving $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ in our study. In the next section, we verify carefully the order of accuracy of the scheme proposed. In Sections 3 and 4, we prove that the solution of the semi-discretized system, as well as the corresponding scheme, preserve the expected solution positivity and monotonicity under reasonable restrictions. Section 5 is devoted to standard and reinforced stability analysis of the numerical method. Experimental results to validate our analysis and explore global nonlinear convergence of the stochastic singular Kawarada solutions are presented in Section 6. Finally, Section 7 is tailored for concluding our remarks, concerns, and forthcoming endeavors.

2 Derivation of Scheme

Without loss of generality, we set $\displaystyle c=1.$ Utilizing transformations $\displaystyle\tilde{x}=x/a,~{}\tilde{y}=y/b,$ and reusing the original variables for notation simplicity, we may reformulate $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ as

[TABLE]

where $\displaystyle{\cal{D}}=(-1,1)\times(-1,1)\subset{\mathbb{R}}^{2}$ and $\displaystyle g(\varepsilon,u)=f(\varepsilon,u)/\sigma(x,y).$ For fixed $\displaystyle N_{x},N_{y}\gg 1$ we define

[TABLE]

where $\displaystyle x_{0}=y_{0}=-1,~{}x_{i+1}-x_{i}=h_{x,i}>0,~{}y_{j+1}-y_{j}=h_{y,j}>0$ for $\displaystyle 0\leq i\leq N_{x},~{}0\leq j\leq N_{y},$ and $\displaystyle x_{N_{x}+1}=y_{N_{y}+1}=1.$ Let $\displaystyle u_{i,j}(t)$ denote an approximation of the solution of $\displaystyle(\ref{c1})$-$\displaystyle(\ref{c2b})$ at $\displaystyle(x_{i},y_{j},t)\in{\cal{D}}_{h}\times[t_{0},T_{q}).$ We consider the following nonuniform finite difference formulas [13],

[TABLE]

for $\displaystyle(x_{i},y_{j})\in{\cal{D}}_{h}^{\circ}.$ Further, for a standard lexicographical ordering, we denote $\displaystyle v(t)=(u_{1,1},u_{2,1},\dots,$ $\displaystyle u_{N_{x},1},u_{1,2},\dots,u_{N_{x},N_{y}})^{\intercal}\in{\mathbb{R}}^{N_{x}N_{y}},$ and let $\displaystyle g(\varepsilon,v)$ be a suitable discretization of the nonlinear stochastic reaction term. We derive from from $\displaystyle(\ref{c1})$-$\displaystyle(\ref{c2b})$ our semi-discretized system

[TABLE]

Denote $\displaystyle h_{\min}=\textstyle\min_{1\leq i\leq N_{x};1\leq j\leq N_{y}}\{h_{x,i-1}h_{x,i},h_{y,j-1}h_{y,j}\}$ and

[TABLE]

where $\displaystyle\otimes$ stands for the Kronecker product, $\displaystyle I_{N_{\alpha}}\in{\mathbb{R}}^{N_{\alpha}\times N_{\alpha}},~{}\alpha=x,y,$ are identity matrices. Further, let

[TABLE]

and for the above,

[TABLE]

The formal solution of $\displaystyle(\ref{c8})$, $\displaystyle(\ref{c8b})$ is

[TABLE]

where $\displaystyle E(\cdot)$ is the corresponding matrix exponential [13, 19, 20]. An aim of this paper is to investigate a highly efficient numerical procedure for the problem $\displaystyle(\ref{c8})$-$\displaystyle(\ref{rs1})$.

To this end, we consider the following realization of $\displaystyle(\ref{rs1})$ on variable temporal grids:

[TABLE]

where $\displaystyle v_{k}$ and $\displaystyle v_{k+1}$ are approximations to $\displaystyle v(t_{k})$ and $\displaystyle v(t_{k+1}),$ respectively, $\displaystyle v_{0}$ is the initial vector, $\displaystyle t_{k}=t_{0}+\textstyle\sum_{j=0}^{k-1}\tau_{j},$ $\displaystyle 0<\tau_{k}\ll 1,$ $\displaystyle k=0,1,2,\dots,$ and $\displaystyle\{\tau_{k}\}_{k\geq 0}$ is a set of adaptive temporal steps. In fact, a rapid temporal step refinement mechanism must be triggered robotically as $\displaystyle t_{k}\to T_{q}^{-},$ since $\displaystyle u_{t}$ may grow exponentially and affect computational accuracy near quenching points [1, 3, 21]. On the other hand, an overly refined mesh should also be avoided [22, 23]. In light of these concerns, a standard arc-length monitoring function may be implemented to handle the strong nonlinearity at quenching [13, 14, 18, 22, 12, 24].

Further, to avoid the need for nonlinear iterations, in computations we replace $\displaystyle g(\varepsilon,v_{k+1})$ by $\displaystyle g(\varepsilon,w_{k}),$ where $\displaystyle w_{k}$ is a suitable approximation to $\displaystyle v_{k+1}$ with a proper order of accuracy [10, 5, 6]. However, no such approximation is utilized in the following analysis.

We have the following weighted approximation.

Lemma 2.1**.**

Assume that $\displaystyle f(\varepsilon,u)\in C_{c}^{2}[0,1)$ and $\displaystyle\tau_{k}\rightarrow 0^{+}.$ Then for $\displaystyle\theta\in[0,1]$ we have

[TABLE]

where $\displaystyle p=2$ when $\displaystyle\theta=1/2,$ and $\displaystyle p=1$ otherwise.

Proof.

We begin with the following difference

[TABLE]

Replace the first integral by a trapezoidal rule and integrate the second directly to obtain

[TABLE]

as $\displaystyle\tau_{k}\rightarrow 0^{+}.$ Consider local expansions of each term. A direct comparison then leads to.

[TABLE]

and apparently, the approximation is first-order unless $\displaystyle\theta=1/2,$ in which case the approximation is second-order.

Let $\displaystyle\Phi_{x,k}\equiv E(\tau_{k}M_{x}/2)$ and $\displaystyle\Phi_{y,k}\equiv E(\tau_{k}M_{y}).$ We further denote $\displaystyle\Phi_{k}\equiv\Phi_{x,k}\Phi_{y,k}\Phi_{x,k}$ and $\displaystyle F(\varepsilon,v_{k},v_{k+1},\theta)\equiv\theta g(\varepsilon,v_{k})+(1-\theta)g(\varepsilon,v_{k+1}).$ Thus, from $\displaystyle(\ref{xxx})$ and Lemma 2.1 we acquire that

[TABLE]

where $\displaystyle\theta\in[0,1]$ is a fixed parameter, and $\displaystyle p\in\{1,2\}$ depending upon the choice of $\displaystyle\theta.$ We note that the stochastic influences are completely contained in the vector-valued function $\displaystyle F(\varepsilon,v_{k},v_{k+1},\theta).$

Before proceeding, we would like to note that, in general, $\displaystyle(\ref{xder})$, $\displaystyle(\ref{yder})$ are first order approximations [13, 20]. The following lemma demonstrates that in certain cases the accuracy can be increased to second-order.

Lemma 2.2**.**

If the nonuniform mesh $\displaystyle{\cal{D}}_{h}$ is the image of a uniform mesh on $\displaystyle{\cal{D}}$ via a smooth mapping, then $\displaystyle(\ref{xder})$ and $\displaystyle(\ref{yder})$ are in fact second-order.

Proof.

We consider the case for $\displaystyle(\ref{xder})$ and note that the result for $\displaystyle(\ref{yder})$ follows in a similar manner. Assume that $\displaystyle u\in C^{4}({\cal{D}})$ and $\displaystyle{\cal{D}}_{h}$ is generated via a smooth mapping $\displaystyle g^{(x)}$ defined as

[TABLE]

where $\displaystyle\omega^{(x)}_{i}=-1+ik_{1},\ i=0,\ldots,N+1,$ with $\displaystyle k=2/(N+1).$ Thus, $\displaystyle u(x,y,t)=u(g^{(x)}(\omega^{(x)}),y,t).$ We now consider the local truncation error, $\displaystyle\psi^{(x)}_{i,j}(t),$ of $\displaystyle(\ref{xder})$ at a point $\displaystyle(x_{i},y_{j},t)\in{\cal{D}}_{h}\times[t_{0},T_{q}):$

[TABLE]

For general nonuniform meshes, $\displaystyle(\ref{111a})$ is first-order, however, due to $\displaystyle(\ref{1111})$, we have

[TABLE]

Note that for $\displaystyle\xi_{i}=\epsilon\omega_{i}+(1-\epsilon)\omega_{i-1},\ \epsilon\in[0,1],$ we have

[TABLE]

which gives that $\displaystyle k_{1}=(g^{(x)}_{\omega}(\xi_{i}))^{-1}h_{x,i-1}={\cal{O}}(h_{x,i-1}).$ A similar process also gives $\displaystyle k_{1}={\cal{O}}(h_{x,i}).$ These facts combined with $\displaystyle(\ref{111a})$ and $\displaystyle(\ref{1122})$ gives the desired result.

3 Positivity

It has been shown that solutions to $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ remain positive within their domains of existence; i.e, for $\displaystyle t\in[t_{0},T_{q})$ [3, 7, 8, 21]. This is also a key characteristics of a healthy physical solution no matter whether an ignition, or quenching, occurs [1, 24]. The property must be properly preserved by $\displaystyle(\ref{A3})$, and likewise, by our semidiscretized problem $\displaystyle(\ref{c8})$, $\displaystyle(\ref{c8b})$.

All vector inequalities are to be interpreted componentwise in following investigations.

Lemma 3.1**.**

Let $\displaystyle K>1$ and $\displaystyle A\in{\mathbb{R}}^{K\times K}$ be nonsingular and nonnegative. Further, let $\displaystyle\beta\in{\mathbb{R}}^{K}$ be positive. Then $\displaystyle A\beta>0.$

Proof.

The proof follows directly from the definitions.

Lemma 3.2**.**

The matrix $\displaystyle E(\tau_{k}M)$ is positive for all $\displaystyle\tau_{k}\geq 0.$

Proof.

Note that $\displaystyle M=M_{x}+M_{y}$ is such that $\displaystyle M_{ij}\geq 0$ for $\displaystyle i\neq j.$ Let $\displaystyle N=M+\beta I,$ where $\displaystyle\beta=4/h_{\min}\textstyle\min\{a^{2},b^{2}\}.$ Then $\displaystyle N$ is a positive matrix and

[TABLE]

Since $\displaystyle\beta I$ is a diagonal matrix, it follows that$\displaystyle E(-\tau_{k}\beta I)$ is a diagonal matrix, with $\displaystyle\left(E(-\tau_{k}\beta I)\right)_{ij}=e^{-\tau_{k}\beta}>0$ for $\displaystyle i=j$ and $\displaystyle 0$ elsewhere. Further, $\displaystyle E(\tau_{k}N)$ is positive by the definition of the matrix exponential and the fact that $\displaystyle N$ is positive. Thus, we have the result.

Lemma 3.3**.**

The solution to $\displaystyle(\ref{c8})$, $\displaystyle(\ref{c8b})$ is positive for all $\displaystyle t_{0}<t<T_{q}$ and $\displaystyle 0\leq u<1.$

Proof.

Note that the solution to $\displaystyle(\ref{c8})$, $\displaystyle(\ref{c8b})$ is given by $\displaystyle(\ref{rs1})$, and that $\displaystyle E(tM)v_{0}$ is positive for all positive $\displaystyle v_{0}$ due to Lemmas 3.1 and 3.2. Similarly, $\displaystyle E((t-\tau)M)g(v(\tau))$ is positive since $\displaystyle g$ is a positive vector and by Lemma 3.1. Since the integral of positive vectors is positive, it follows naturally that $\displaystyle v(t)$ is positive.

Lemma 3.4**.**

The matrices $\displaystyle\Phi_{x,k}$ and $\displaystyle\Phi_{y,k}$ are positive for all $\displaystyle\tau_{k}\geq 0.$ Likewise, $\displaystyle\Phi_{k}$ is positive for all $\displaystyle\tau_{k}\geq 0.$

Proof.

The proof is similar to that of Lemma 3.2.

We study the matrix $\displaystyle M^{-1}[\Phi_{k}-I]$ in the remainder of this section. To show that $\displaystyle M^{-1}[\Phi_{k}-I]$ is positive under appropriate restrictions of $\displaystyle\tau_{k},$ we need the following definition and lemmas.

Definition 3.1**.**

Let $\displaystyle K>1$ and $\displaystyle A\in\mathbb{C}^{K\times K},$ $\displaystyle I_{K}\in\mathbb{C}^{K\times K}$ be the identity matrix. Further, let $\displaystyle\|\cdot\|$ be an induced matrix norm. Then the associated logarithmic norm $\displaystyle\mu:\mathbb{C}^{K\times K}\to{\mathbb{R}}$ of $\displaystyle A$ is defined as

[TABLE]

Remark 3.1*.*

[19] When the spectral norm is considered, we have $\displaystyle\mu_{2}(A)=\lambda_{\max}[(A+A^{*})/2].$

Remark 3.2*.*

[25] When the maximal norm is considered, we have $\displaystyle\mu_{\infty}(A)=\textstyle\sup\{\mbox{Re}(a_{ii})+\textstyle\sum_{j\neq i}|a_{ij}|\}.$

Lemma 3.5**.**

For $\displaystyle t\geq 0$ we have $\displaystyle\|E(tA)\|\leq E(t\mu(A))$ for any logarithmic norm, $\displaystyle\mu,$ and any induced matrix norm, $\displaystyle\|\cdot\|.$

Proof.

Consider the linear initial value problem,

[TABLE]

which has the solution $\displaystyle x(t)=E(tA)x_{0}.$ Recall that the upper right Dini derivative for any function $\displaystyle x(t)$ with respect to $\displaystyle t$ is defined as

[TABLE]

Using this fact and taking the norm, we obtain the differential inequality

[TABLE]

where $\displaystyle\|\cdot\|$ is any induced matrix norm. Note that

[TABLE]

is the solution to $\displaystyle(\ref{inq1})$ [26]. Dividing both sides of $\displaystyle(\ref{inq2})$ by $\displaystyle\|x_{0}\|\neq 0$ and taking the supremum over all $\displaystyle x_{0}$ such that $\displaystyle\|x_{0}\|\neq 0$ gives the desired result.

Lemma 3.6**.**

For any $\displaystyle\tau_{k}\geq 0,$ we have $\displaystyle\|\Phi_{x,k}\|_{\infty}\leq 1$ and $\displaystyle\|\Phi_{y,k}\|_{\infty}\leq 1.$

Proof.

The proof follows immediately by applying Lemma 3.5 and then utilizing Remark 3.2 to compute the logarithmic norm of each matrix involved.

Now, we wish to determine a bound for $\displaystyle\|M^{-1}\|_{\infty}$ under suitable conditions on our nonuniform spatial grids. To that end, we denote

[TABLE]

where $\displaystyle\alpha=x,y.$ Let $\displaystyle\sigma_{\max}=\textstyle\max_{1\leq i\leq N_{x};1\leq j\leq N_{y}}\sigma_{i,j}.$

Lemma 3.7**.**

Assume that $\displaystyle k(h)\geq 1.$ Then $\displaystyle\|M^{-1}\|_{\infty}\leq\sigma_{\max}/2.$

Proof.

Let $\displaystyle\tilde{M}=B^{-1}M=a^{-2}(I_{N_{y}}\otimes T_{x})+b^{-2}(T_{y}\otimes I_{N_{x}}).$ We have

[TABLE]

Thus, we only need to consider the case with $\displaystyle\tilde{M}.$ To this end, we let $\displaystyle f_{i,j}$ be any mesh function defined on $\displaystyle{\cal{D}}_{h}.$ Define $\displaystyle u_{i,j}=-\tilde{M}^{-1}f_{i,j},$ where $\displaystyle u_{i,j}$ solves the equation $\displaystyle-\tilde{M}u_{i,j}=f_{i,j}$ on $\displaystyle{\cal{D}}_{h}$ and satisfies the homogeneous Dirichlet boundary conditions $\displaystyle u_{i,j}=0,~{}(x_{i},y_{j})\in\partial{\cal{D}}_{h}.$ For $\displaystyle w(x,y)=[x^{2}+y^{2}]/4,$ we define $\displaystyle w_{i,j}$ on $\displaystyle\overline{{\cal{D}}}_{h}$ by

[TABLE]

Clearly, $\displaystyle\Delta w=w_{xx}+w_{yy}=1.$ We now consider $\displaystyle\tilde{M}w_{i,j}$ by considering $\displaystyle\overline{\Delta}_{x,h}w_{i,j}\equiv 4(I_{N_{y}}\otimes T_{x})w_{i,j}$ and $\displaystyle\overline{\Delta}_{y,h}w_{i,j}\equiv 4(T_{y}\otimes I_{N_{x}})w_{i,j}$ separately. Thus, we have

[TABLE]

A similar argument indicates that

[TABLE]

Combining these results gives that $\displaystyle\tilde{M}w_{i,j}\geq k(h)\boldsymbol{x},$ where $\displaystyle\boldsymbol{x}=(1,1,\ldots,1)^{\intercal}\in{\mathbb{R}}^{N_{x}N_{y}}.$

Note that $\displaystyle w_{i,j}$ is nonnegative on $\displaystyle\overline{{\cal{D}}}_{h}$ and $\displaystyle\|w_{i,j}\|_{\infty}=1/2.$ It follows therefore

[TABLE]

based on our earlier assumptions. Hence, by the Discrete Maximum Principle, $\displaystyle\|f_{i,j}\|_{\infty}w_{i,j}+u_{i,j}$ achieves its maximum on $\displaystyle\partial{\cal{D}}_{h}.$ Since $\displaystyle u_{i,j}=0$ on $\displaystyle\partial{\cal{D}}_{h},$ we have

[TABLE]

for $\displaystyle(x_{i},y_{j})\in{\cal{D}}_{h}.$ Similarly, we have

[TABLE]

which gives

[TABLE]

A combination of the above results gives

[TABLE]

It follows naturally that $\displaystyle\|\tilde{M}^{-1}\|_{\infty}\leq 1/2$ since $\displaystyle f_{i,j}$ are arbitrary. This completes our proof.

We note that our restriction on $\displaystyle k(h)$ is reasonable and has no negative impact on desired computations. In fact, $\displaystyle k(h)\geq 1$ can be observed if uniform or symmetric nonuniform grids are utilized. Moreover, for certain nonuniform meshes we have the following corollary.

Corollary 3.1**.**

If the nonuniform mesh $\displaystyle{\cal{D}}_{h}$ is the image of a uniform mesh on $\displaystyle{\cal{D}}$ via a smooth mapping, then we have $\displaystyle\|M^{-1}\|_{\infty}\leq\sigma_{\max}/2+{\cal{O}}(\textstyle\max_{i,j}\{h_{x,i}^{2},h_{y,j}^{2}\}),~{}\textstyle\max_{i,j}\{h_{x,i}^{2},h_{y,j}^{2}\}\rightarrow 0^{+}.$

Proof.

Assume that $\displaystyle u\in C^{4}({\cal{D}})$ and $\displaystyle{\cal{D}}_{h}$ is generated via smooth mappings $\displaystyle g^{(x)}$ and $\displaystyle g^{(y)}$ with

[TABLE]

where $\displaystyle\omega^{(x)}_{i}=-1+ik_{1},\ i=0,\ldots,N_{x}+1$ and $\displaystyle\omega^{(y)}_{j}=-1+ik_{2},\ j=0,\ldots,N_{y}+1.$ Using the Taylor expansions of the functions $\displaystyle g^{(x)}$ and $\displaystyle g^{(y)}$ in $\displaystyle(\ref{kval})$ yields $\displaystyle k(h)=1+{\cal{O}}(\max\{k_{1}^{2},k_{2}^{2}\}).$ Just as in Lemma 2.2, we may deduce that $\displaystyle k_{1}={\cal{O}}(h_{x,i})$ and $\displaystyle k_{2}={\cal{O}}(h_{y,j}),$ which gives $\displaystyle k(h)=1+{\cal{O}}(\textstyle\max_{i,j}\{h_{x,i}^{2},h_{y,j}^{2}\}).$ Using this fact in the proof of Lemma 3.7 gives the desired result.

Lemma 3.8**.**

If $\displaystyle k(h)\geq 1$ then $\displaystyle M^{-1}[\Phi_{k}-I]$ is positive for sufficiently small $\displaystyle\tau_{k}\geq 0.$

Proof.

Note that

[TABLE]

where

[TABLE]

with

[TABLE]

In a manner similar to that in [27], it follows

[TABLE]

where

[TABLE]

Note that $\displaystyle\int_{t_{k}}^{t_{k+1}}E_{1}(s)E_{2}(s)E_{1}(s)\,ds$ is positive according to Lemmas 3.1 and 3.4, so it is sufficient to show that the norm of the second term in $\displaystyle(\ref{w2})$ can be made sufficiently small. To this end,

[TABLE]

Thus, our claim is true for sufficiently small $\displaystyle\tau_{k}>0.$

Remark 3.3*.*

Lemma 3.8 holds for all $\displaystyle\tau_{k}\geq 0$ in many realistic situations where $\displaystyle\|[M_{x},M_{y}]\|_{\infty}=0.$ The equality is apparently true for cases with a symmetric degeneracy term on symmetric nonuniform grids. Numerical experiments show that $\displaystyle\|[M_{x},M_{y}]\|_{\infty}$ remains well-bounded.

4 Monotonicity

As has been discussed in [3, 7, 8, 14, 21], another fundamental property of solutions of $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ is their monotonicity. Solutions to Kawarada problems are monotonically increasing with respect to time on their domain of existence. Thus, it is important that this is rigorously observed in $\displaystyle(\ref{A3})$ as well as $\displaystyle(\ref{c8})$,$\displaystyle(\ref{c8b})$.

In the following results, there are additional restrictions place on the stochastic nonlinear reaction term in order to guarantee monotonicity. Further, our results are similar to those in the non-stochastic case. We first consider $\displaystyle(\ref{c8})$,$\displaystyle(\ref{c8b})$.

Lemma 4.1**.**

If $\displaystyle Mv_{0}+g(\varepsilon,v_{0})>0,$ then the solution of $\displaystyle(\ref{c8})$,$\displaystyle(\ref{c8b})$ is monotonically increasing with respect to $\displaystyle t\in(t_{0},T_{q})$ while maintaining $\displaystyle 0\leq u<1.$

Proof.

Recall $\displaystyle(\ref{c8})$, $\displaystyle(\ref{c8b})$. Let $\displaystyle w(t)=v^{\prime}(t)$ and consider

[TABLE]

The solution to the above is

[TABLE]

where $\displaystyle g_{v}(\varepsilon,v(t))$ is the Jacobian matrix. It is also clear that $\displaystyle(g_{v}(\varepsilon,v(t)))$ is a diagonal matrix. It follows that $\displaystyle g_{v}$ is a nonnegative matrix since $\displaystyle f_{u}\geq 0.$ Denote

[TABLE]

Hence $\displaystyle A_{ij}\geq 0$ for $\displaystyle i\neq j.$ Consequently, by Lemma 3.2, we conclude that $\displaystyle E(A(t))$ is positive for $\displaystyle t_{0}\leq t<T_{q}.$ By Lemma 3.1, we know that $\displaystyle w(t)>0$ for $\displaystyle t_{0}<t<T_{q}$ since $\displaystyle w_{0}>0.$ Thus, $\displaystyle v^{\prime}(t)>0$ for $\displaystyle t_{0}<t<T_{q}$ and our solution is monotonically increasing with respect to time $\displaystyle t.$

We now consider $\displaystyle(\ref{A3})$.

Lemma 4.2**.**

If $\displaystyle Mv_{0}+g(\varepsilon,v_{0})>0,$ then $\displaystyle Mv_{k}+g(\varepsilon,v_{k})>0$ for all $\displaystyle k\geq 0.$

Proof.

Let $\displaystyle Mv_{0}+g(\varepsilon,v_{0})>0.$ We note for $\displaystyle(\ref{A3})$ that

[TABLE]

which implies that

[TABLE]

by inductive assumptions, Lemmas 3.1 and 3.4, and the fact that $\displaystyle f(\varepsilon,u)$ is monotone.

Lemma 4.3**.**

Let $\displaystyle Mv_{0}+g(\varepsilon,v_{0})>0$ and $\displaystyle\tau_{k}\geq 0$ be sufficiently small. Then the solution sequence of $\displaystyle(\ref{A3})$, that is, $\displaystyle\{v_{k}\}_{k\geq 0},$ is monotonically increasing while maintaining $\displaystyle 0\leq v_{k}<1,~{}k\geq 0.$

Proof.

From $\displaystyle(\ref{A3})$ we have

[TABLE]

according to Lemmas 3.8 and 4.2.

The above results can be generalized readily to obtain the following.

Theorem 4.1**.**

Let $\displaystyle Mv_{k_{0}}+g(\varepsilon,v_{k_{0}})>0$ for any $\displaystyle k_{0}\geq 0$ and $\displaystyle\tau_{k}\geq 0$ be sufficiently small for $\displaystyle k\geq k_{0}.$ Then the sequence $\displaystyle\left\{v_{k}\right\}_{k\geq{k_{0}}}$ generated by the semi-discretized nonuniform scheme $\displaystyle(\ref{A3})$ increases monotonically until unity is reached or exceeded by one or more components of the solution vector, i.e., until quenching occurs.

5 Stability

Stability of different schemes for solving singular and stochastic Kawarada equations such as $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ has been a highly nontrivial and challenging concern [5, 6, 7, 10, 13, 14]. This difficulty is mostly attributed to the fact that near quenching the temporal derivative of the solution may blow up faster than an exponential rate [12, 24]. Most existing investigations simply consider a linear stability analysis, and later the results are reconfirmed through rigorous computational experiments [10, 11, 13, 14]. Since such considerations have been proven to be highly effective in localized computations [28, 12], in this paper, we will first carry out a standard linear stability analysis for $\displaystyle(\ref{A3})$. The study will then be extend to a nonlinear case by including contributions from the singular nonlinear stochastic reaction term. This extension of the analysis ensures that the singular and stochastic nature of the differential equation is incorporated near the quenching point as quenching is approached.

Lemma 5.1**.**

$\displaystyle T_{\alpha},~{}\alpha=x,y,$* utilized in $\displaystyle M$ are congruent to symmetric matrices. *

Proof.

The proof is similar to that in [24].

Lemma 5.2**.**

$\displaystyle\rho(T_{\alpha})\in(-\infty,0],$* where $\displaystyle\rho(\cdot)$ is the spectral radius and $\displaystyle\alpha=x,y.$*

Proof.

We may consider $\displaystyle T_{x}$ since the other case is similar. By Lemma 5.1 we have that $\displaystyle T_{x}$ is congruent to a symmetric matrix, hence all eigenvalues of $\displaystyle T_{x}$ are real. Further, since $\displaystyle T_{x}$ is diagonally dominant with all negative diagonal elements, the result follows immediately.

Lemma 5.3**.**

All eigenvalues of $\displaystyle M_{x}~{}\mbox{and}~{}M_{y}$ are real and negative. Further $\displaystyle\mu(B^{1/2}T_{x}B^{1/2})<0$ and $\displaystyle\mu(B^{1/2}T_{y}B^{1/2})<0.$

Proof.

The proof is similar to that in [29].

Lemma 5.4**.**

We have $\displaystyle\left\|\textstyle\prod_{j=0}^{k}\Phi_{j}\right\|_{2}\leq\sqrt{\kappa(B)}=\sqrt{\max\sigma_{i,j}/\min\sigma_{i,j}}$ for $\displaystyle\tau_{j}>0.$

Proof.

We begin with $\displaystyle\Phi_{y,j}=B^{1/2}E(\tau_{j}B^{1/2}T_{y}B^{1/2})B^{-1/2},~{}\tau_{j}>0.$ For it, we have

[TABLE]

By the same token,

[TABLE]

It then follows from $\displaystyle(\ref{xx1})$, $\displaystyle(\ref{xx2})$, and Lemmas 3.5 and 5.3 immediately that

[TABLE]

where $\displaystyle\kappa(B)$ is the condition number of the matrix $\displaystyle B.$

By combining the above results, we arrive at the following theorem.

Theorem 5.1**.**

*For $\displaystyle\tau_{k}\geq 0$ sufficiently small, the semi-adaptive nonuniform method $\displaystyle(\ref{A3})$ with the stochastic nonlinear reaction term frozen is unconditionally stable in the von Neumann sense under the spectral norm. *

Proof.

Recall $\displaystyle(\ref{A3})$. When the stochastic nonlinear reaction term is frozen, the error $\displaystyle z_{k+1}$ takes the form of

[TABLE]

Iterating $\displaystyle(\ref{a})$ repeatedly to yield

[TABLE]

The above leads to

[TABLE]

Our proof is thus accomplished by letting $\displaystyle C\equiv\sqrt{\kappa(B)}$ in the above.

In the circumstance when the stochastic nonlinear reaction term in $\displaystyle(\ref{A3})$ is not frozen, we denote $\displaystyle t_{Q}$ as the time at which numerical quenching occurs, that is, for the first time for which $\displaystyle\|v_{Q}\|_{\infty}\geq 1.$

Theorem 5.2**.**

For $\displaystyle\tau_{k}\geq 0$ sufficiently small, the semi-adaptive nonuniform method $\displaystyle(\ref{A3})$ is unconditionally stable in the von Neumann sense, that is, for every $\displaystyle t_{m}<t_{Q}$ there exists a constant $\displaystyle C(t_{m})>0$ such that

[TABLE]

where $\displaystyle z_{0}=v_{0}-\tilde{v}_{0}$ is an initial error, $\displaystyle z_{k+1}=v_{k+1}-\tilde{v}_{k+1}$ is the $\displaystyle(k+1)$th perturbed error vector, and $\displaystyle C(t_{m})>0$ is a constant independent of $\displaystyle k$ and $\displaystyle\tau_{j}$ for each $\displaystyle 0\leq j\leq k.$

Proof.

Recalling $\displaystyle(\ref{A3})$, we have

[TABLE]

It follows that

[TABLE]

where $\displaystyle\xi_{j}=\epsilon v_{j}+(1-\epsilon)\tilde{v}_{j},~{}\epsilon\in[0,1],~{}j=k,k+1,$ and $\displaystyle k\geq 0,~{}\theta\in[0,1].$ Recall $\displaystyle(\ref{w2})$. By utilizing a left-endpoint quadrature to approximate the integral in $\displaystyle(\ref{stab1})$, we observe that

[TABLE]

where $\displaystyle\Gamma_{1,k+1}=I-\tau_{k}(1-\theta)\left(I+{\cal{O}}\left(\tau_{k}\right)\right)g_{v}(\varepsilon,\xi_{k+1})$ and $\displaystyle\Gamma_{2,k}=I+\tau_{k}\theta\left(I+{\cal{O}}\left(\tau_{k}\right)\right)g_{v}(\varepsilon,\xi_{k}).$ It follows therefore that, by employing the standard [1/0] and [0/1] Padé approximations, respectively,

[TABLE]

Apply the above to $\displaystyle(\ref{stab2})$. We obtain

[TABLE]

Now, let $\displaystyle G(t_{m})=\max_{0\leq j\leq m}\left\|g_{v}(\varepsilon,\xi_{j})\right\|_{2}.$ Taking the norm of both sides of the above expression yields

[TABLE]

where $\displaystyle c_{1},~{}c_{2},~{}\mbox{and}~{}c_{3}$ are constants independent of $\displaystyle\tau_{j}$ and $\displaystyle j.$ Thus, for $\displaystyle\tau_{j}>0$ sufficiently small,

[TABLE]

where we have used the fact that $\displaystyle\textstyle\sum_{j=0}^{k}\tau_{j}\leq T_{q}.$ It follows that $\displaystyle G(t_{m})<\infty$ for all $\displaystyle t_{m}<t_{Q}.$ This completes our proof.

6 Numerical Experiments

We consider the following degenerate stochastic Kawarada model problem

[TABLE]

where $\displaystyle 0\leq t_{0}<t\leq T_{q}<\infty,$ $\displaystyle u_{0}\in C^{2}([-1,1]\times[-1,1]),$ and $\displaystyle 0\leq u_{0}\ll 1.$

Without loss of generality, we set $\displaystyle t_{0}=0,$ $\displaystyle a=b=2,$ and choose $\displaystyle u_{0}(x,y)=0.001(1-\cos(2\pi x))(1-\cos(2\pi y)),\ -1\leq x,y\leq 1.$ Temporal adaptations are initiated once

[TABLE]

since quenching is expected to occur in the situation [29, 16, 3, 10, 12]. We fix the parameter $\displaystyle\theta=1/2$ in $\displaystyle(\ref{A3})$. We primarily focus on the capability of the scheme $\displaystyle(\ref{A3})$ to recover anticipated numerical solutions and, at the same time, explore the nonlinear numerical convergence of the semi-adaptive semi-discretized scheme. Our first experiment will serve as a baseline approach, where nonlinear convergence results at different time levels in the absence of degenerate or stochastic influences will be presented. Influences of the degeneracy term on the convergence will be closely observed in the second experiment. Standard [1/1] Padé approximants will be utilized for matrix exponential evaluations [20]. Appropriate restrictions are realized for guaranteeing the required positivity and monotonicity with the particular matrices involved throughout experiments [5, 6, 29].

6.1 Example 1

Figure 1. Solution $\displaystyle u$ [LEFT] and its temporal derivative $\displaystyle u_{t}$ [RIGHT] immediately prior to quenching.

To prepare our baseline results, we set $\displaystyle\varphi(\varepsilon)=\sigma(x,y)\equiv 1$ in $\displaystyle(\ref{k1})$-$\displaystyle(\ref{k4})$. Without the concerns of stochastic and degenerate impacts, we adopt a symmetric nonuniform spatial grid generated in a similar manner to that in [29]. A typical symmetric nonuniform grid can be seen in Figure 2. Adaption in the temporal direction is implemented via $\displaystyle u_{t}$ based arc-length monitoring functions [11, 12, 29]. Quenching is observed precisely at the origin for $\displaystyle T_{q}\approx 0.519715480937553$ which agrees well with existing results [5, 6, 9, 12, 16, 10, 11]. While the peak value of $\displaystyle u$ reaches unity peacefully, the temporal derivative function $\displaystyle u_{t}$ grows at an explosive rate to reach $\displaystyle u_{t}^{*}\approx 4635.870128316449$ as $\displaystyle t\rightarrow T_{q}.$ Profiles of the numerical solution and its temporal derivative function are depicted in Figure 1.

Figure 2. [LEFT] A typical symmetric nonuniform mesh used in computations for the non-degenerate case. [RIGHT] A symmetric nonuniform mesh which can be used in degenerate computations for a possible degeneracy that is concentrated at $\displaystyle(-1,-1).$ In both plots, the red point represents the quenching point.

To examine the convergence of the numerical solution numerically, we arbitrarily choose two temporal points close to the quenching time $\displaystyle T_{q}:$

[TABLE]

Let $\displaystyle u_{h}^{\tau}$ denote the numerical solution obtained by using the set of nonuniform spatial steps represented by size $\displaystyle h$ and the particular adaptive temporal step $\displaystyle\tau.$ Suppose that $\displaystyle h/2$ represents the new set of nonuniform spatial steps obtained from dividing each nonuniform mesh step in set $\displaystyle h$ by two. Likewise, $\displaystyle h/4$ is for the new set of nonuniform spatial steps obtained from dividing each nonuniform mesh step in set $\displaystyle h/2$ by two. By the same token, we may define sets of temporal steps $\displaystyle\tau/2.~{}\tau/4$ based on set $\displaystyle\tau.$ Now, we define the following generalized Milne device for rates of point-wise convergence in space and time, respectively:

[TABLE]

where these differences are defined on mesh points corresponding to $\displaystyle\tau,h.$ The estimates can also be reformulated via any vector norm. As an example of such reformulations, we include estimates for convergence rates with respect the spectral norm, denoted by $\displaystyle p_{2}^{h}$ and $\displaystyle q_{2}^{\tau}.$

[TABLE]

Table 1. Detailed spatial and temporal convergence rates. The convergence rates in space match the expectations of Lemma 2.2. Temporal convergence rates agree well with the Courant-Friedrichs-Lewy condition.

Table 1 presents the convergence data for space and time at sampling times $\displaystyle T_{1}$ and $\displaystyle T_{2}.$ From Table 1 we may also notice that the mean, median, and maximal point-wise convergence rates and the spectral convergence rates agree well in both space and time with a slight reduction when considering spectral convergence. We note that the method used to generate the nonuniform mesh for this experiment meets the criteria of Lemma 2.2, hence, we observe second-order spatial convergence. Further, these results support the Courant-Friedrichs-Lewy condition remarkably.

Figure 3. Spatial convergence rate $\displaystyle p_{\rm PW}^{h}$ and $\displaystyle p_{2}^{h}$ [LEFT] and temporal convergence rate $\displaystyle q_{\rm PW}^{\tau}$ and $\displaystyle p_{2}^{h}$ [RIGHT] for the solution at times $\displaystyle T_{1},~{}T_{2}$ are given. The mean (blue), median (green), maximal (red), and spectral (magenta) convergence rates are included. Spatial and temporal rates are approximately two and one, respectively.

Figure 4. 3D profiles at $\displaystyle T_{2}.$ The first row of surfaces are for the solution [LEFT] and its temporal derivative function [RIGHT], whereas the second row of surfaces represent the spatial convergence rate $\displaystyle p_{\rm PW}^{h}$ [LEFT] and temporal convergence rate $\displaystyle q_{\rm PW}^{\tau}$ [RIGHT].

Figure 3 depicts the dynamic profile of convergence rates between $\displaystyle T_{1}$ and $\displaystyle T_{2}.$ The lines indicate movements of the maximal, median, and mean point-wise convergence rates as well as the spectral convergence rates in both space and time. When considering the spatial convergence, we observe that the curves are very close to two, and when considering the temporal convergence, we notice that the curves are grouped very close to one. Minimal spatial and temporal convergence rates are recorded in Table 1. The minimum spatial convergence rates most likely due to the smoothness of the solution about corners of the spatial domain [30]. The minimal temporal convergence rates are most likely due to the development of quenching singularity, which indicated by rate drops from $\displaystyle T_{1}$ to $\displaystyle T_{2}.$ The quenching singularity also affects the spatial convergence, but less violently as the gradient of the solution is bounded. These expectations follow from a generalization of the results in [31].

Further, Figure 4 depicts the solution, its temporal derivative, and convergence rates at $\displaystyle T_{2}.$ We note that the peaks of the solution and the derivative are concentrated at the origin, which is the quenching location in this particular problem [29, 12, 3, 4, 31]. In the bottom row of the figure, rates of point-wise spatial and temporal convergences are plotted as surfaces. It is easy to observe that both rates decrease near the origin. This, once again, is due to the reaction term becoming unbounded near quenching, or combustion. The singularity has an explosive impact to the overall convergence rates. It is also observed that that there are decreased spatial convergence rates near corners of the domain, which has been predicted by several authors [30, 15, 17].

Despite drops in convergence rates near the quenching point, the mean convergence rate of the scheme is well-maintained near quenching at time $\displaystyle T_{2}.$ At this time, the mean spatial and temporal convergence rates are approximately 1.998013984025239 and 0.991391747823586, respectively. The spectral convergence rates in space and time are 1.995539227366137 and 0.953687356004341, respectively. We note that the spatial convergence matches the expectations of Lemma 2.2 but the temporal convergence rate is lower than the expected rate from Lemma 2.1. However, the temporal convergence rate is appropriately proportional to the spatial rate, hence, the reduced temporal convergence is a result of the Courant-Friedrichs-Lewy condition.

Figure 5. Surface profiles representing convergence rates, $\displaystyle p_{\rm PW}^{h},~{}q_{\rm PW}^{\tau},$ of the temporal derivative of the solution, $\displaystyle u_{t},$ at $\displaystyle T_{1}$ [LEFT] and $\displaystyle T_{2}$ [RIGHT] . The minimum in the latter case is approximately 0.751476025769804.

We would like to also note that implementing uniform grids or considerations of less exotic initial data will not drastically change the convergence rates nor the overall shape of the convergence surface plots. Furthermore, similar conclusions are confirmed in our experiments when nonsymmetric grids are utilized, though these results are admitted for brevity. The effects of nonuniform grids, that are not generated in a smooth manner, on convergence rates have yet to be explored in relation to quenching-combustion problems. On the other hand, to further raise the accuracy of schemes, say, to obtain a second-order temporal accuracy, we would need to implement spatial discretizations to the fourth order. This could be fulfilled via a nine-point scheme or a compact algorithm [20]. We also observe that, while the emergence of the quenching singularity reduces the spatial convergence rate significantly near the origin, the reduction is less severe than the reduction in temporal convergence. This is probably due to the fact that $\displaystyle u$ and its spatial derivatives are well-behaved in our quenching computations [31, 23].

The temporal derivative of the solution of $\displaystyle(\ref{k1})$-$\displaystyle(\ref{k4})$ becomes unbounded as $\displaystyle t\rightarrow T_{q}^{-}$ [23]. As a reference, in Figure 5 we explore the temporal convergence rate of $\displaystyle u_{t}$ at times $\displaystyle T_{1}$ and $\displaystyle T_{2},$ respectively, in order to examine whether the explosive rate of change can significantly affect the convergence rate. The shapes of the surfaces in Figure 5 are no doubt similar to that in Figure 4. We find that the convergence rates drop significantly near the origin. However, the mean temporal convergence rate is still approximately first-order. We observe the average point-wise convergence rate is approximately 0.991674641553278 at $\displaystyle T_{1}$ and 0.986792663593675 at $\displaystyle T_{2}.$ Further, the spectral convergence rates of $\displaystyle u_{t}$ are 0.948154262871667 and 0.819817700033525 at times $\displaystyle T_{1}$ and $\displaystyle T_{2},$ respectively. This novel exploration of convergence issues provides profound insights into the full nature of singular nonlinear problems such as $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$.

6.2 Example 2

In this experiment we explore impacts of a degeneracy and a mild stochastic influence on both the solution and rates of convergence. In the following, we set $\displaystyle\varphi(\varepsilon)=\varepsilon$ for $\displaystyle\varepsilon\in[0.98,1.02],$ and $\displaystyle\sigma(x,y)=\textstyle\sqrt{(x+1)^{2}+(y+1)^{2}},\ -1\leq x,y\leq 1.$ We note that this choice for $\displaystyle\sigma(x,y)$ forces the degeneracy to occur at the spatial point $\displaystyle(x_{d},y_{d})=(-1,-1)$ [13, 24]. By choosing such a degeneracy function, we will be better able to observe its fast spreading effects on the solution. Our nonuniform grids are generated in the consideration that the quenching location will be shifted, but not too far away, from the origin [12]. See Figure 2 for a possible example of how to shift the mesh appropriately. The convergence rates are computed in a similar fashion as before.

Figure 6. Solution $\displaystyle u$ [LEFT] and its temporal derivative function $\displaystyle u_{t}$ [RIGHT] as $\displaystyle t\rightarrow T_{q}.$ We have the maximum $\displaystyle{}_{-1\leq x,y\leq 1}u_{t}(x,y,t)=u_{t}(x_{q},y_{q},t)\approx 279.5375137783287$ as $\displaystyle t\rightarrow T_{q}.$

An elevated quenching time, $\displaystyle T_{q}\approx 0.836473787737391,$ probably primarily due to the degeneracy [17], is observed. However, the vibrant stochastic feature of the reaction term may also increase quenching time [29]. Figure 6 depicts the solution and its derivative immediately prior to quenching. Just as we have anticipated, the quenching location is slightly shifted from the origin to $\displaystyle(x_{q},y_{q})\approx(-0.151533166458073,-0.151533166458073).$ The peak value of $\displaystyle u_{t},$ which is concentrated around the quenching location, reaches 279.5375137783287 prior to quenching.

Due to our slightly delayed ignition-quenching time $\displaystyle T_{q},$ we consider the following two arbitrary convergence testing times

[TABLE]

A major change in spatial convergence is expected since the degeneracy $\displaystyle\sigma(x,y)$ is a spatial function creating strong singularities on the boundary of the spatial domain. In particular, a modification of the proof of Lemma 2.2 will show that the degeneracy limits the effects of generating a nonuniform mesh via a smooth mapping. It seems that the stronger the singularity is, the more the order of spatial accuracy and convergence is reduced.

[TABLE]

Table 2. Detailed spatial and temporal convergence rates. The experimental results suggest clearly a slightly above one spatial convergence rate, and slightly below one temporal rate.

Table 2 is devoted to convergence information in space and time at sampling times $\displaystyle T_{1}$ and $\displaystyle T_{2}.$ From Table 2 we may also notice that while the mean and median convergence rates strongly suggest a slightly higher than one rate of convergence in space and slightly lower than one rate of convergence in time, minimal and maximal rate estimates differ significantly. Again, the point-wise rates of convergence and temporal spectral rate of convergence decrease as $\displaystyle t\rightarrow T_{q}^{-},$ though reductions are less significant probably also due to a stronger restriction on temporal step sizes to maintain the solution positivity when $\displaystyle\sigma$ is nontrivial [29, 6, 5, 12, 10, 11, 24]. The spatial spectral convergence rate actually increases in time, but this is due to the solution becoming more smooth near the corners of the domain as the solution is iterated. The phenomena are well documented by Figure 7 where spatial and temporal convergence rates of the numerical solution under the influences of a non-constant $\displaystyle\sigma(x,y)$ and a nonlinear stochastic reaction term are exhibited.

Figure 7. Spatial convergence rates $\displaystyle p_{\rm PW}^{h}$ and $\displaystyle p_{2}^{h}$ [LEFT] and temporal convergence rates $\displaystyle q_{\rm PW}^{h}$ and $\displaystyle q_{2}^{h}$ [RIGHT] for the solution at times $\displaystyle T_{1}$ and $\displaystyle T_{2}.$ The mean (blue), median (green), maximal (red), minimal (cyan), and spectral (magenta) rates are included. The mean and median curves almost overlap in the first picture. The negative minimal rate of convergence in space is due to spurious artifacts of the degeneracy.

Figure 8. Surface profiles at time $\displaystyle T_{2}.$ The first row is for the solution [LEFT] and its temporal derivative [RIGHT], whereas the second row is for the rates of convergences in space [LEFT] and in time [RIGHT].

Figure 8 illustrates the solution of our stochastic degenerate problem, as well as its derivative, at time $\displaystyle T_{2}.$ Surface plots of the spatial and temporal convergence are also presented. We note that there is a visible shift of the quenching location, as well as the position of the peaks of $\displaystyle u,~{}u_{t},$ to $\displaystyle(x_{q},y_{q}).$ This shift agrees well with existing predictions [12]. Slight dips at the quenching point can also be observed in Figure 8. Needless to say, the point-wise convergence rate surfaces provide us with further insights into the maximal and minimal rates noted in Table 2 and Figure 7. From the surface plot, we observe that spurious issues occur at the corner opposite of the degeneracy in the spatial domain. This phenomenon can be explained by the fact that the degeneracy simplifies the problem at $\displaystyle(-1,-1)$ by reducing the parabolic problem computationally to an elliptic problem. However, the degeneracy still creates disturbances, which are observed in the extreme convergence rates near $\displaystyle(1,1).$

Figure 9. Surface plots of the point-wise rates of convergence rates in time for the temporal derivative function of the solution, $\displaystyle u_{t},$ at $\displaystyle T_{1}$ [LEFT] and $\displaystyle T_{2}$ [RIGHT] within the spatial domain. The mean convergence rates are approximately 0.991674641553278 and 0.986792663593675, respectively.

As we have seen, surface illustrations will reflect the fact that the convergence issues can be restricted to a small number of mesh points near a corner of the domain due to the impacts from degeneracy. This is further recorded in Figure 9 which shows that the degeneracy has a noticeable effect on the temporal convergence of $\displaystyle u_{t}.$ There is an obvious ring of distinctive oscillatory waves near $\displaystyle(1,1)$ in the figure. However, these oscillations have little impact on the mean or median convergence rates which evolve from 0.991674641553278 to 0.986792663593675 in Figure 9. Further, the spectral convergence rates of $\displaystyle u_{t}$ evolve from 0.996338976463597 to 0.985326655425598. These rates of convergence for the temporal derivative function are ideal, as they remain approximately the same magnitudes as those from the previous example.

7 Conclusions

A semi-discretized method for solving degenerate stochastic Kawarada equations is derived and studied. The numerical method utilizes arbitrary grids in space and time. Various types of adaptations can thus be introduced and incorporated to more precisely capture the underlying quenching, stochastic, and degeneracy singularities. Exponential splitting strategies are employed for achieving a higher effectiveness and efficiency in computations. The semi-discretized scheme is proven to preserve key physical and mathematical properties such as the positivity and monotonicity of the solution under modest restrictions. The scheme acquired is stable in the classical von Neumann sense without additional constraints. An enriched stability analysis which incorporates the nonlinear stochastic quenching singularity is implemented. This extends the stability investigations in several recent publications [5, 6, 29], and has proven to be more effective as compared with traditional analysis (see [10, 11, 24, 12] and references therein).

Since the Lax-Richtmeyer Equivalence theorem is not applicable to Kawarada equations which are highly nonlinear, explorations of the nonlinear convergence of the aforementioned semi-discretized scheme are conducted through multiple numerical experiments. Our experiments indicate that, despite the strong quenching singularity, stochastic influence, and degeneracy, desirable rates of convergence can be expected. Although the degeneracy and quenching singularities do affect the convergence, especially within certain regions, mean and median rates of convergence remain relatively unaffected. Our experiments also depict novel presentations of global convergence rates via convergence surfaces at various times. These surfaces provide important insights into the nature of nonlinear convergence. They also help to better understand how singular components in the differential equation may affect convergence as quenching is approached.

Our continuing endeavors include a study of the nonlinear stability and convergence of degenerate stochastic Kawarada equations. While the classical stability concept has been extended in our discussions, more rigorous nonlinear stability analysis and statistical tools will be implemented in future work. We intend to extend our study to Sobolev-like spaces and norms in the near future, with the hopes that the current convergence approaches may have provided us with more valuable structural information. On the other hand, we also intend to incorporate balanced nonlocal derivatives into the model [32]. The inclusion of fractional partial differential equations may help to capture and utilize global features of the numerical combustion with solid applications to the energy industry.