On a Nonuniform Crank-Nicolson Scheme for Solving the Stochastic Kawarada Equation via Arbitrary Grids

TL;DR

This paper introduces a nonuniform finite difference scheme on arbitrary grids for the stochastic Kawarada equation, ensuring stability, positivity, and monotonicity, with numerical experiments validating the approach.

Contribution

It presents a novel nonuniform finite difference method on adaptive grids for the stochastic Kawarada equation, addressing singularities and uncertainty.

Findings

Scheme preserves positivity and monotonicity under mesh constraints

Numerical experiments confirm stability and accuracy

Adaptive grids effectively handle singularities and stochastic sources

Abstract

This paper studies a nonuniform finite difference method for solving the degenerate Kawarada quenching-combustion equation with a vibrant stochastic source. Arbitrary grids are introduced in both space and time via adaptive principals to accommodate the uncertainty and singularities involved. It is shown that, under proper constraints on mesh step sizes, the positivity, monotonicity of the solution, and numerical stability of the scheme developed are well preserved. Numerical experiments are given to illustrate our conclusions.

Figures

Equations

$$\sigma (x)u_t=u_{xx}+\phi (ϵ)f(u),(x,t)\in \mathcal{E},$$

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

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

$$f(0)=f_0>0,\underset{u\to 1^{-}}{\lim }f(u)=+\infty ,$$

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

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

$$\kappa =\underset{0<\xi <\infty }{\max }{\int }_0^{\xi }{e}^{t^2-{\xi }^2}dt.$$

$$u_t=\psi (x)u_{xx}+g(u,ϵ),(x,t)\in \mathcal{E},$$

$$u(-1,t)=u(1,t)=0,t>t_0,$$

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

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

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

$$v(t_0)$$

$$B=\mbox{diag}\left(1/\sigma_{1},\dots,1/\sigma_{N}\right),~{}{\color[rgb]{0,0,0}P}=\mbox{tridiag}\left(l_{i},~{}m_{i},~{}n_{i}\right)$$

$$l_i$$

$$m_i$$

$$v(t)=E(tM)v_0+{\int }_{t_0}^{t}E((t-\tau )M)g(v(\tau ))d\tau ,t_0<t<T,$$

$$p(t)={\left(I-\frac{t}{2}M\right)}^{-1}\left(I+\frac{t}{2}M\right),t_0<t<T.$$

$$v(t)=p(t)\left[v_0+\frac{t}{2}g(v_0)\right]+\frac{t}{2}g(v(t))+\mathcal{O}\left((t-t_0{)}^3\right),|t-t_0|\to 0^{+}.$$

$$v_{\ell +1}={\left(I-\frac{{\tau }_{\ell }}{2}M\right)}^{-1}\left(I+\frac{{\tau }_{\ell }}{2}M\right)\left(v_{\ell }+\frac{{\tau }_{\ell }}{2}g(v_{\ell })\right)+\frac{{\tau }_{\ell }}{2}g(v_{\ell +1}),$$

$$w_{\ell }=v_{\ell }+{\tau }_{\ell }(M{v}_{\ell }+g(v_{\ell })),0<{\tau }_{\ell }\ll 1,$$

$${\tau }_{\ell }<a^2{\beta }_{\min },$$

$${\Vert \frac{{\tau }_{\ell }}{2}M\Vert }_2\le \frac{{\tau }_{\ell }}{2a^2}\Vert B{\Vert }_2\Vert T{\Vert }_2\le \frac{2{\tau }_{\ell }}{a^2}\Vert B{\Vert }_2\left(1/\underset{i=0,1,\dots ,N}{\min }\left\{h_i^2\right\}\right)=\frac{2{\tau }_{\ell }}{a^2{h}^2}\Vert B{\Vert }_2<1.$$

$$M{v}_{k+1}+g(v_{k+1})$$

$$M{v}_1+g(v_1)>{\left(I-\frac{{\tau }_0}{2}M\right)}^{-1}\left(I+\frac{{\tau }_0}{2}M\right)\left[M{v}_0+g(v_0)\right]>0$$

$$M{v}_{\ell }+g(v_{\ell })>{\left(I-\frac{{\tau }_{\ell -1}}{2}M\right)}^{-1}\left(I+\frac{{\tau }_{\ell -1}}{2}M\right)\left[M{v}_{\ell -1}+g(v_{\ell -1})\right]>0,$$

$$v_{k+1}-v_k$$

$$\left(I-\frac{{\tau }_{\ell }}{2}M\right)x\ge x.$$

$$w=\left(I-\frac{{\tau }_{\ell }}{2}M\right)x=(w_1,w_2,\dots ,w_N{)}^{⊺}.$$

$$w_1$$

$$w_i$$

$$w_N$$

$$v_1={\left(I-\frac{{\tau }_0}{2}M\right)}^{-1}\left(I+\frac{{\tau }_0}{2}M\right)\frac{{\tau }_0}{2}g(0)+\frac{{\tau }_0}{2}g(v_1).$$

Full text

On a Nonuniform Crank-Nicolson Scheme for Solving the Stochastic Kawarada Equation via Arbitrary Grids

Joshua L. Padgett111Principal and corresponding author. Email address: Josh Padgett@baylor.edu and Qin Sheng

Department of Mathematics and Center for Astrophysics, Space Physics and Engineering Research, Baylor University, Waco, TX 76798-7328, USA

Abstract. This paper studies a nonuniform finite difference method for solving the degenerate Kawarada quenching-combustion equation with a vibrant stochastic source. Arbitrary grids are introduced in both space and time via adaptive principals to accommodate the uncertainty and singularities involved. It is shown that, under proper constraints on mesh step sizes, the positivity, monotonicity of the solution, and numerical stability of the scheme developed are well preserved. Numerical experiments are given to illustrate our conclusions.

Keywords. Stochastic Kawarada equation, quenching blow-up, nonuniform grids, numerical stability, monotonicity, positivity

1 Introduction

Let $\displaystyle{\cal{D}}=(-a,a),~{}{\cal{E}}={\cal{D}}\times(t_{0},T),~{}\partial{\cal{D}}=\bar{{\cal{D}}}\setminus{\cal{D}}$ and $\displaystyle\SS=\partial{\cal{D}}\times(t_{0},T),$ where $\displaystyle a>0,~{}0\leq t_{0}<T<\infty.$ We are interested in the monotonically increasing positive solution of the degenerate stochastic Kawarada problem,

[TABLE]

where the degeneracy function $\displaystyle\sigma(x)\geq 0$ for $\displaystyle x\in\bar{{\cal{D}}},$ and the equality occurs only on $\displaystyle\partial{\cal{D}}.$ The nonlinear source function, $\displaystyle f(u),$ is strictly increasing for $\displaystyle 0\leq u<1$ with

[TABLE]

and $\displaystyle\varphi(\epsilon):~{}0<\varphi_{\min}\leq\varphi\leq\varphi_{\max},$ is a stochastic inference function of the random variable, or white noise, $\displaystyle\epsilon(x).$ The existence and uniqueness of the solution of $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ can be viewed as a generalization of the results given by Chan and Levine [1, 2]. It is also observed that solutions of the stochastic modeling problem $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ are in general only fractional order Hölder continuous [3]. Further, the solution $\displaystyle u$ of $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ is said to quench if there exists a finite time $\displaystyle T_{a}>0$ such that

[TABLE]

Such a value $\displaystyle T_{a}$ is called the quenching time [1, 4, 5]. It has been shown that a necessary condition for quenching to occur is

[TABLE]

It is known that $\displaystyle T_{a}$ exists only when $\displaystyle a$ is greater than a certain critical value $\displaystyle a^{*}\ll\infty.$ The interval $\displaystyle{\cal{D}}$ associated with such an $\displaystyle a^{*}$ is defined as the critical domain and denoted as $\displaystyle{\cal{D}}^{*}.$ Therefore, $\displaystyle(\ref{a2})$ occurs only when $\displaystyle{\cal{D}}^{*}\subseteq{\cal{D}},$ otherwise the monotone positive solution of $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ exists globally for $\displaystyle T\rightarrow\infty$ [2, 4, 5]. In the particular circumstance when $\displaystyle\sigma(x),\varphi(\epsilon)\equiv 1$ and $\displaystyle f(u)=1/(1-u),$ it has been shown that $\displaystyle a^{*}=\kappa\sqrt{2}$ [5, 6], where

[TABLE]

Kawarada partial differential equations have been intensively used for modeling numerous important phenomena in nature. They characterize not only ignitions of liquid fuels in combustion chambers, but also turbulent macro or micro flows between channel walls [6, 7, 8]. The latter is particularly meaningful for predicting and preventing oil pipeline decays [9]. Though computational results for $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ can be found in numerous recent publications, most of numerical analysis presented relies heavily on the 1-norm or $\displaystyle\infty$-norm [10, 11, 12]. The numerical analysis in the current paper implements the more preferred spectral norm. In addition, the influence of white noise associated to the source is considered. It should be noted, however, that impacts of nonsmooth sources are in general different from those due to nonsmooth initial data. In fact, while smooth solutions are still possible if nonsmooth initial data are dealt with appropriately, nonsmooth solutions are almost certain when a nonsmooth reaction term is utilized [13, 14].

This paper proposes a temporally adaptive Crank-Nicolson scheme. Predetermined nonuniform grids are utilized in space and are chosen in order to ideally incorporate the effects of the quenching singularity and stochastic influences in space. The use of predetermined nonuniform spatial grids is practically preferable, especially in cases when quenching locations are predictable [8]. A further merit of such a semi-adaptive infrastructure is that it can be conveniently extended for solving multidimensional Kawarada problems. It also makes the subsequent numerical analysis much simpler and straightforward. Initial approaches of such an idea with uniform spatial grids can be found in [12].

It is crucial that the numerical solution acquired preserves fundamental features of the physical solution, such as the positivity, monotonicity, quenching time, and location. To this end, our discussions will be organized as follows. In the next section, the adaptive Crank-Nicolson scheme for solving $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ is implemented and evaluated. Its positivity is investigated. In Section 3, constraints under which the numerical solution is monotone are determined. Our stability analysis is conducted in Section 4. We first accomplish a standard stability analysis for a fully linearized scheme. Then an extended stability analysis is fulfilled for a fully nonlinear method. In Section 5, several numerical examples are provided. These examples provide interesting insights into the effects of the degeneracy and stochastic functions on not only quenching times but also quenching locations. Finally, our investigations are concluded through remarks and proposed future problems in Section 6.

2 Semi-adaptive Crank-Nicolson scheme and its positivity

Utilizing the transformation $\displaystyle\tilde{x}=x/a$ and reusing the original variable and other notations for simplicity, we may reformulate $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ as

[TABLE]

where $\displaystyle{\cal{D}}=(-1,1),~{}{\cal{E}}={\cal{D}}\times(t_{0},T),~{}\psi(x)=1/({a^{2}\sigma(x)}),$ and $\displaystyle g(u,\epsilon)={\varphi(\epsilon)f(u)}/{\sigma(x)}.$ For $\displaystyle N\gg 1,$ we inscribe over $\displaystyle\bar{{\cal{D}}}$ the variable grid: $\displaystyle{\cal{D}}_{h}=\left\{x_{i}\,:\,i=0,\dots,N+1;~{}x_{0}=-1,~{}x_{i}<x_{i+1},~{}x_{N+1}=1\right\}.$ Denote $\displaystyle h_{i}=x_{i+1}-x_{i}$ for $\displaystyle 0\leq i\leq N.$ Let $\displaystyle u_{i}=u_{i}(t)$ be an approximation of $\displaystyle u(x_{i},t)$ and adopt the nonuniform finite difference [12],

[TABLE]

where $\displaystyle{\cal{D}}_{h}^{\circ}={\cal{D}}_{h}\setminus\{x_{0},x_{N+1}\}.$ Further, denote $\displaystyle v(t)=(u_{1},u_{2},\dots,u_{N})^{\intercal}\in{\mathbb{R}}^{N}$ and let $\displaystyle g(v)$ be a discretization of the source term in $\displaystyle(\ref{c1})$. We obtain readily from $\displaystyle(\ref{c1})$-$\displaystyle(\ref{c2b})$ the following semi-discretized problem

[TABLE]

where \displaystyle M=\frac{1}{a^{2}}B{\color[rgb]{0,0,0}P}\in{\mathbb{R}}^{N\times N},

[TABLE]

and for the above

[TABLE]

The formal solution of $\displaystyle(\ref{c8})$, $\displaystyle(\ref{c8b})$ can thus be written as

[TABLE]

where $\displaystyle E(\cdot)=\exp(\cdot)$ is the matrix exponential [12]. We proceed by approximating $\displaystyle(\ref{rs1})$ via a trapezoidal rule and a [1/1] Padé approximation; that is, $\displaystyle E(tM)=p(t)+{\cal{O}}\left(t^{3}\right),$ where

[TABLE]

These lead to

[TABLE]

Based on the above, we obtain the following second-order in time semi-adaptive Crank-Nicolson scheme on variable spatial grids:

[TABLE]

where $\displaystyle v_{\ell}$ and $\displaystyle v_{\ell+1}$ are approximations of $\displaystyle v(t_{\ell})$ and $\displaystyle v(t_{\ell+1}),$ respectively, $\displaystyle v_{0}$ is the initial vector, $\displaystyle t_{\ell}=t_{0}+{\textstyle\sum_{k=0}^{\ell-1}\tau_{k}},~{}0<\tau_{\ell}\ll 1,~{}\ell=0,1,2,\ldots,$ and $\displaystyle\{\tau_{\ell}\}_{\ell\geq 0}$ is a set of adaptive temporal steps. In order to avoid a fully implicit scheme, $\displaystyle g(v_{\ell+1})$ may be approximated by $\displaystyle g(w_{\ell}),$ where $\displaystyle w_{\ell}$ is an approximation to $\displaystyle v_{\ell+1},$ such as

[TABLE]

in practical computations.

Recall $\displaystyle(\ref{a1})$. Due to the strong singularity of $\displaystyle u_{t}$ as $\displaystyle t$ approaches $\displaystyle T_{a},~{}a\geq a^{*},$ selecting the proper temporal steps $\displaystyle\tau_{\ell}$ is vital in computations. To this end, we may consider employing arc-length monitoring functions [8, 12, 15, 16], or allow the temporal steps to be proportional to the source term or its gradient[10, 11, 12, 17].

Positivity is one of the most important characteristics of the solution of Kawarada problems including $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ and $\displaystyle(\ref{c1})$-$\displaystyle(\ref{c2b})$ [1, 2, 4, 5]. In order for our numerical solutions to be valid, it is crucial that they preserve this property. To this end, we let $\displaystyle\wedge$ be one of the operations $\displaystyle<,~{}\leq,~{}>,~{}\geq$ and $\displaystyle\alpha,~{}\beta\in{\mathbb{R}}^{K_{1}\times K_{2}}.$ We assume the following notations in subsequent discussions:

$\displaystyle\alpha\wedge\beta$ means $\displaystyle\alpha_{i,j}\wedge\beta_{i,j},~{}i=1,2,\ldots,K_{1};\,j=1,2,\ldots,K_{2};$ 2. 2.

$\displaystyle a\wedge\alpha$ means $\displaystyle a\wedge\alpha_{i,j},~{}~{}i=1,2,\ldots,K_{1};\,j=1,2,\ldots,K_{2},$ for any given scalar $\displaystyle a.$

Lemma 2.1**.**

$\displaystyle\|T\|_{2}\leq\max_{i=0,1,\dots,N}\left\{{4}/{h_{i}^{2}}\right\}.$**

Proof.

The proof is similar to the one from our earlier investigations [17].

For the following we denote $\displaystyle\beta_{\min}=h^{2}/2\|B\|_{2}$ and $\displaystyle h=\textstyle\min_{i=0,1,\dots,N}\{h_{i}\}.$

Lemma 2.2**.**

If

[TABLE]

then $\displaystyle I-\frac{\tau_{\ell}}{2}M$ and $\displaystyle I+\frac{\tau_{\ell}}{2}M$ are nonsingular. Further, $\displaystyle I-\frac{\tau_{\ell}}{2}M$ is monotone and inverse-positive, and $\displaystyle I+\frac{\tau_{\ell}}{2}M$ is nonnegative.

Proof.

First, we note that

[TABLE]

Hence, $\displaystyle I+\frac{\tau_{\ell}}{2}M$ is nonsingular, and also nonnegative.

Next, we consider $\displaystyle A=I-\frac{\tau_{\ell}}{2}M.$ As $\displaystyle A_{ij}\leq 0$ for $\displaystyle i\neq j$ and the weak row sum criterion is satisfied; hence $\displaystyle A$ is monotone, and it follows that its inverse exists and is nonnegative. So, $\displaystyle A$ must be inverse-positive[18]. This ensures the proof.

Lemma 2.3**.**

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

Proof.

This is clear by the definition of the operations.

3 Monotonicity

Another fundamental feature which distinguishes a quenching solution from other blow-up type solutions is its monotonicity with respect to time $\displaystyle t\geq t_{0}$ [1, 2, 4, 5, 8]. It is therefore necessary to guarantee that a numerical solution preserves this important physical property when solving the Kawarada equation $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ or $\displaystyle(\ref{c1})$-$\displaystyle(\ref{c2b})$.

Lemma 3.1**.**

If $\displaystyle Mv_{0}+g(v_{0})>0,$ then it follows that $\displaystyle Mv_{\ell}+g(v_{\ell})>0$ for all $\displaystyle\ell\geq 0.$

Proof.

First, we proceed be computing the following:

[TABLE]

where the inequality follows from the fact that $\displaystyle f(\varepsilon,u)$ is strictly increasing and Lemma 2.3. Second, we proceed by induction. Letting $\displaystyle k=0$ we have

[TABLE]

by the assumption and Lemmas 2.2 and 2.3. Third, we assume that the inequality holds for $\displaystyle k=\ell-1.$ It follows that

[TABLE]

by the inductive assumption and Lemmas 2.2 and 2.3, which completes the induction.

Lemma 3.2**.**

If $\displaystyle(\ref{cfl})$ holds, $\displaystyle 0\leq\tau_{k}\leq 1$ for all $\displaystyle 0\leq k\leq\ell,$ and $\displaystyle Mv_{0}+g(v_{0})>0,$ then $\displaystyle v_{\ell+1}\geq v_{\ell}\mbox{ for all }\ell\geq 0.$ That is, the sequence $\displaystyle\left\{v_{\ell}\right\}_{\ell=0}^{\infty}$ is monotonically increasing.

Proof.

From $\displaystyle(\ref{c3})$ we observe that

[TABLE]

by the assumption and Lemmas 2.2 and 2.3. Since the result holds for all $\displaystyle k\geq 0,$ we have that the sequence $\displaystyle\{v_{\ell}\}_{\ell\geq 0}$ is monotonically increasing as desired.

Lemma 3.3**.**

Let $\displaystyle x=(1,1,\dots,1)^{\intercal}\in{\mathbb{R}}^{N}.$ Then for any $\displaystyle\tau_{\ell}>0$ we have

[TABLE]

Proof.

Consider

[TABLE]

First, we observe that

[TABLE]

Secondly, for $\displaystyle i=2,\dots,N-1,$ we have

[TABLE]

Finally, we observe that

[TABLE]

Hence, we may conclude that $\displaystyle w_{i}\geq 1,~{}i=1,\dots,N.$ Therefore $\displaystyle w\geq x.$

In the next lemma we show that numerical quenching, i.e., one or more components of $\displaystyle v_{\ell}$ reaching or exceeding unity, cannot occur immediately after the first time step under appropriate constraints.

Lemma 3.4**.**

Let $\displaystyle x$ be the vector defined in Lemma 3.3 and $\displaystyle v_{0}\equiv 0.$ If $\displaystyle(\ref{cfl})$ hold and $\displaystyle\bar{h}^{2}<{2\|B\|_{2}}/[{a^{2}f(\tau_{0}\varphi_{\max}f_{0}/\sigma_{\min})}],$ where $\displaystyle\bar{h}=\textstyle\max_{i=1,\ldots,N}\left\{h_{i}\right\},$ then $\displaystyle v_{1}<x.$

Proof.

Recall $\displaystyle(\ref{c3})$. For $\displaystyle v_{0}\equiv 0$ we have

[TABLE]

Using

[TABLE]

where $\displaystyle\gamma=\left(\varphi_{1}/\sigma_{1},\dots,\varphi_{N}/\sigma_{N}\right)^{\intercal}\in{\mathbb{R}}^{N},$ we have

[TABLE]

Based on $\displaystyle(\ref{37c})$ we observe that

[TABLE]

where

[TABLE]

It can be seen that

[TABLE]

and the above indicates that

[TABLE]

By Lemma 3.3 we conclude that $\displaystyle s^{-}\leq-x,$ and therefore,

[TABLE]

Since we again wish each component of the above vector to be negative, we need

[TABLE]

Hence $\displaystyle v_{1}-x\leq 0$ follows immediately from the assumptions.

Combining above results we obtain immediately the following.

Theorem 3.1**.**

Assume that for $\displaystyle\ell_{0}\geq 0,$

(i)

$\displaystyle\bar{h}^{2}<\frac{2\|B\|_{2}}{a^{2}f(\tau_{\ell_{0}}\varphi_{\max}f_{\ell_{0}}/\sigma_{\min})},$* where $\displaystyle\bar{h}=\max_{i=1,\dots,N}\{h_{i}\},$*

(ii)

$\displaystyle Mv_{\ell_{0}}+g(v_{\ell_{0}})>0,$**

If $\displaystyle(\ref{cfl})$ holds for all $\displaystyle\ell\geq\ell_{0},$ then the sequence $\displaystyle\left\{v_{\ell}\right\}_{\ell\geq{\ell_{0}}}$ produced by the semi-adaptive nonuniform scheme $\displaystyle(\ref{c3})$ increases monotonically until unity is reached or exceeded by one or more components of the solution vector, i.e., until quenching occurs.

4 Stability

Nonlinear stability has been an extremely challenging issue when Kawarada equations are concerned [4, 8, 10, 11, 12, 17]. On the other hand, linear stability analysis may uncover crucial information for underlying schemes locally and asymptotically [8, 21]. In the following study, we will carry out a linearized stability analysis in the von Neumann sense for $\displaystyle(\ref{c3})$ with its nonlinear source term frozen. The analysis will then be extended to circumstances where the nonlinear term is not frozen.

In the following, let $\displaystyle A\in\mathbb{C}^{N\times N},$ $\displaystyle I_{N}\in\mathbb{C}^{N\times N}$ be the identity matrix, and again denote $\displaystyle E(\cdot)=\exp(\cdot).$

Definition 4.1**.**

[17, 19]* Let $\displaystyle\|\cdot\|$ be an induced matrix norm. Then the associated logarithmic norm $\displaystyle\mu:\mathbb{C}^{N\times N}\to{\mathbb{R}}$ of $\displaystyle A$ is defined as*

[TABLE]

Remark 4.1*.*

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

Lemma 4.1**.**

[19]* For $\displaystyle t\geq 0$ we have $\displaystyle\|E(tA)\|\leq E(t\mu(A)).$*

Lemma 4.2**.**

Let $\displaystyle P$ be as in $\displaystyle(\ref{matrixT})$. Then $\displaystyle P$ is congruent to a symmetric matrix. In particular,

[TABLE]

where

[TABLE]

for which

[TABLE]

Proof.

The proof is clear by the definitions [20].

Lemma 4.3**.**

\displaystyle\rho({\color[rgb]{0,0,0}P})\in(-\infty,0],* where $\displaystyle\rho(A)=\sqrt{\lambda_{\max}(AA^{\intercal})}~{}$ for any $\displaystyle A\in{\mathbb{R}}^{N\times N}.$*

Proof.

By Lemma 4.2 we have that \displaystyle{\color[rgb]{0,0,0}P} is congruent to a symmetric matrix, hence all eigenvalues of \displaystyle{\color[rgb]{0,0,0}P} are real. Since \displaystyle{\color[rgb]{0,0,0}P} is diagonally dominant with negative diagonal elements, the result follows immediately from the Gershgorin circle theorem [18].

Lemma 4.4**.**

All eigenvalues of $\displaystyle M$ are real and negative. Further \displaystyle\mu(B^{1/2}{\color[rgb]{0,0,0}P}B^{1/2})<0.

Proof.

Since \displaystyle M=B{\color[rgb]{0,0,0}P}=B^{1/2}B^{1/2}{\color[rgb]{0,0,0}P}, we have \displaystyle B^{-1/2}M=B^{1/2}{\color[rgb]{0,0,0}P}. Hence

[TABLE]

because $\displaystyle B$ is diagonal. Further,

[TABLE]

is congruent to a symmetric matrix. Thus matrices $\displaystyle B^{-1/2}MB^{1/2}$ and \displaystyle{\color[rgb]{0,0,0}P} are congruent. Since the eigenvalues of \displaystyle{\color[rgb]{0,0,0}P} are real and negative, the eigenvalues of $\displaystyle B^{-1/2}MB^{1/2}$ and $\displaystyle M$ must be real and negative.

Since we have \displaystyle B^{-1/2}MB^{1/2}=B^{1/2}{\color[rgb]{0,0,0}P}B^{1/2}, then the eigenvalues of \displaystyle B^{1/2}{\color[rgb]{0,0,0}P}B^{1/2} are real and negative, which gives the result.

Lemma 4.5**.**

If $\displaystyle\tau_{k}>0,~{}0\leq k\leq\ell,$ we have \displaystyle\left\|\prod_{k=0}^{\ell}E(\tau_{k}M)\right\|_{2}\leq\sqrt{\kappa(B)}~{}{\color[rgb]{0,0,0}=\sqrt{\|B^{-1}\|_{2}\|B\|_{2}}}.

Proof.

We first show \displaystyle E(\tau_{\color[rgb]{0,0,0}k}M)=B^{1/2}E(\tau_{\color[rgb]{0,0,0}k}B^{1/2}{\color[rgb]{0,0,0}P}B^{1/2})B^{-1/2} for \displaystyle\tau_{\color[rgb]{0,0,0}k}>0. To this end, we have

[TABLE]

It then follows from $\displaystyle(\ref{xx})$,

[TABLE]

Lemma 4.6**.**

If $\displaystyle(\ref{cfl})$ holds and $\displaystyle 0\leq\tau_{k}\leq 1$ for all $\displaystyle 0\leq k\leq\ell,$ then

[TABLE]

Proof.

Recalling the [1/1] Padé approximation utilized in Section 2, we have

[TABLE]

By the definition of $\displaystyle T$ and $\displaystyle\{\tau_{\ell}\}_{\ell\geq 0},$ we have $\displaystyle\textstyle\sum_{k=0}^{\ell}\tau_{k}\leq T$ and $\displaystyle\tau=\textstyle\max_{0\leq k\leq\ell}\{\tau_{k}\}\leq 1.$ Now, based on Lemma 4.5,

[TABLE]

which yields the desired bound.

Combining the above results gives the following theorem.

Theorem 4.1**.**

If $\displaystyle(\ref{cfl})$ holds and $\displaystyle 0\leq\tau_{k}\leq 1$ for all $\displaystyle 0\leq k\leq\ell,$ then the semi-adaptive nonuniform method $\displaystyle(\ref{c3})$ with the source term frozen is unconditionally stable in the von Neumann sense under the spectral norm, that is,

[TABLE]

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

Proof.

When the nonlinear source term is frozen, $\displaystyle z_{\ell+1}$ takes the form of

[TABLE]

Iterating $\displaystyle(\ref{a})$ gives

[TABLE]

Taking the norm of both sides of $\displaystyle(\ref{st11})$, it follows that

[TABLE]

where $\displaystyle C$ is a positive constant independent of $\displaystyle\ell$ and $\displaystyle\tau_{k}$ for each $\displaystyle 0\leq k\leq\ell.$

We now consider the case without freezing the nonlinear source term in $\displaystyle(\ref{c3})$. In the following, let $\displaystyle t_{Q}$ be the time at which numerical quenching occurs, that is, $\displaystyle\|v_{Q}\|_{\infty}\geq 1,$ and recall that $\displaystyle t_{m}=t_{0}+\textstyle\sum_{k=0}^{m-1}\tau_{k}$ for any $\displaystyle m\geq 0.$ Also, let

[TABLE]

Theorem 4.2**.**

If $\displaystyle(\ref{cfl})$ holds and $\displaystyle\tau_{k}$ sufficiently small for all $\displaystyle 0\leq k\leq\ell,$ then the semi-adaptive nonuniform method $\displaystyle(\ref{c3})$ 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_{\ell+1}=v_{\ell+1}-\tilde{v}_{\ell+1}$ is the $\displaystyle(\ell+1)$th perturbed error vector, and $\displaystyle C(t_{m})>0$ is a constant independent of $\displaystyle\ell$ and $\displaystyle\tau_{k}$ for each $\displaystyle 0\leq k\leq\ell.$

Proof.

By definition we have

[TABLE]

It follows that

[TABLE]

where $\displaystyle\xi_{k}\in{\mathcal{L}}(v_{k},\tilde{v}_{k}),~{}k=\ell,\ell+1.$ Rearranging the above equality, we have

[TABLE]

Further, it follows that for $\displaystyle\tau_{k}\rightarrow 0^{+}$ being sufficiently small,

[TABLE]

Thus,

[TABLE]

Letting $\displaystyle G(t_{m})=\max_{0\leq k\leq m}\left\|g_{v}(\xi_{k}(t))\right\|_{2}$ and $\displaystyle\tilde{P}=B^{1/2}PB^{1/2}$ we have

[TABLE]

where $\displaystyle c_{1},~{}c_{2}$ are positive constants independent of $\displaystyle\ell$ and $\displaystyle\tau_{k},~{}k=0,1,\dots,\ell.$ Since $\displaystyle t_{m}<t_{Q},$ it follows that $\displaystyle G(t_{m})<\infty$ for all $\displaystyle t_{m}<t_{Q}.$ This yields the desired stability.

Remark 4.2*.*

The nonuniform constant $\displaystyle G(t_{m})$ is anticipated to be large within the interval $\displaystyle[t_{0},t_{m}]\subset[t_{0},t_{Q}).$ However, we observe in experiments that its values remain to be well-manageable as far as proper stopping criteria are adopted. This allows for the claimed stability.

5 Numerical Experiments

We consider the following stochastic Kawarada model problem in the first three experiments:

[TABLE]

where $\displaystyle T<\infty,~{}\theta>0,~{}u_{0}\in C[-1,1],$ and $\displaystyle 0\leq u_{0}\ll 1.$

Without loss of generality, we set $\displaystyle\theta=1,$ $\displaystyle t_{0}=0,$ $\displaystyle u_{0}(x)=0.001(1-\cos(2\pi x)),~{}-1\leq x\leq 1,$ and let our temporal adaptations start once $\displaystyle v=\textstyle\max_{-1\leq x\leq 1}u(x,t)$ reaches a certain value $\displaystyle v^{*}.$ In most computational procedures, we adopt $\displaystyle v^{*}=0.90$ for such a triggering criterion [8]. The main purposes of our computations are to examine the numerical method built, and to explore impacts of varying physical parameters used in the equation. In the first experiment we will focus on the consequences due to variations in domain sizes. Quenching behaviors and possible temporal blow-up times are recorded. The second experiment will illustrate impacts of the degeneracy on quenching profiles. The third experiment will explore the effects that a stochastic component plays on the numerical solution and quenching criteria. Our semi-adaptive algorithm $\displaystyle(\ref{c3})$ coupled with $\displaystyle(\ref{num2})$, $\displaystyle(\ref{num3})$ is shown to be satisfactorily reliable, effective, and accurate. Finally, we extend our pursuits with a two-dimensional stochastic Kawarada equation problem to demonstrate the usability and effectiveness of our semi-adaptive infrastructure in higher-dimensional applications. A typical LOD [22] strategy is used.

Experiment 1

Letting $\displaystyle\sigma(x)=\varphi(\epsilon)\equiv 1$ be fixed, we vary the value of $\displaystyle a$ to study its effects on the quenching phenomenon. Since quenching, if it exists, must occur at $\displaystyle x_{q}=0$ [8, 12], we may consider herein a set of symmetric grids for the simplicity of computations. We generate the nonuniform grids in the following way. Consider a parabola designed to have a prescribed minimum at $\displaystyle x_{\lfloor(N+2)/2\rfloor}$ and a prescribed maximum at $\displaystyle x_{0}$ and $\displaystyle x_{N+1}.$ Then a standard arc-length adaptation procedure based on the curvature [10] produces the nonuniform grids, which is scaled to fit $\displaystyle[-1,1].$

Figure 1. Initial function $\displaystyle u_{0}(x)=0.001(1-\cos(2\pi x))$ [LEFT] and its estimated temporal derivative via $\displaystyle(\ref{num1})$ [RIGHT]. Locations of multiple maximal values of the functions concerned are symmetric with respect to the origin.

To commence, we adopt $\displaystyle a=0.5$ and plot $\displaystyle u_{0}(x)$ and its reference temporal derivative in Fig 1. We note that the locations of the twin peaks of the initial functions chosen and triple peaks of the temporal derivative are symmetric with respect to the origin.

Figure 2. Two-dimensional thermal flow plots of the numerical solution $\displaystyle u$ [LEFT] and $\displaystyle u_{t}$ [RIGHT] for $\displaystyle t\in[0,T].$ It is observed that while the heat flows from the left to the right smoothly, symmetrically and increases monotonically in the first figure, the flow speed decreases rapidly but monotonically as time goes on. There is no quenching found in this situation.

Figure 3. Three-dimensional plots of the numerical solution curvatures of $\displaystyle u$ [LEFT] and $\displaystyle u_{t}$ [RIGHT] for $\displaystyle t\in[0,T].$ The magenta curve represents the maximal value trajectory of $\displaystyle u$ [LEFT], and the red curve is for that of $\displaystyle u_{t}$ [RIGHT] at different $\displaystyle t$-levels for $\displaystyle 0<t\leq T.$

Figure 4. Curves in the main frames are for profiles of the maximal values of $\displaystyle u$ [LEFT] and $\displaystyle u_{t}$ (RIGHT). That is, they are plots of $\displaystyle\textstyle\max_{-1\leq x\leq 1}u,~{}0\leq t\leq T,$ and $\displaystyle\textstyle\max_{-1\leq x\leq 1}u_{t},~{}0<\epsilon\leq t\leq T,$ respectively. The embedded graphics represent profiles of $\displaystyle u$ [LEFT] and $\displaystyle u_{t}$ [RIGHT] at time $\displaystyle t=T.$ Up to 200,000 temporal steps have been executed. Temporal adaptation is never activated in this circumstance since no quenching is generated.

Let $\displaystyle T=1.0479.$ In Fig. 2, the left and right graphics show the two-dimensional heat flow of $\displaystyle u$ and its velocity, $\displaystyle u_{t},$ respectively. It can be observed in the former that the heat flows smoothly and monotonically from the left to right. However, the multiple-peak maximal values quickly merge into one and the flow tends to be steady as time increases. Temporal locations of the slices in Fig. 3 are chosen by evenly dividing the arc-length of the function $\displaystyle L(t)=\textstyle\max_{-1\leq x\leq 1}u(x,t).$ These plots again indicate that the solution $\displaystyle u$ tends to be steady with limited changes, while the rate of change function quickly tapers and then diminishes when $\displaystyle t$ tends to $\displaystyle T.$ Both the numerical solution $\displaystyle u$ and temporal derivative $\displaystyle u_{t}$ preserve symmetry about $\displaystyle x=0$ as expected. The temporal adaptation is never activated as the maximum of the solution stays far below unity. The experimental results presented are consistent with existing results [2, 8, 12].

It is noticed in this experiment that $\displaystyle u$ exists globally. More detailed profiles of the maximal values of $\displaystyle u$ and $\displaystyle u_{t}$ at different time levels are shown in Fig. 4 for $\displaystyle 0\leq t\leq T.$ Embedded figures are for $\displaystyle u$ and $\displaystyle u_{t}$ in terminal positions at $\displaystyle t=T.$ It is again observed that while the solution $\displaystyle u$ increases monotonically, the temporal derivative function $\displaystyle u_{t}$ decreases after some initial disturbances. In the terminal position $\displaystyle T_{0}=1.052907287028235,$ we have

[TABLE]

In the next experiment, we choose $\displaystyle a=2$ in order to witness a quenching case for which the physical solution should exist only for finite time. We show three-dimensional profiles of the numerical solution $\displaystyle u$ and its temporal derivative $\displaystyle u_{t}$ in Fig. 5 for $\displaystyle 0\leq t\leq T^{*}$ and $\displaystyle T^{0}\leq t\leq T^{*},$ where $\displaystyle T^{0}=0.509286490538884,$ respectively. A quenching time $\displaystyle T^{*}\approx 0.509391490538887$ is recorded. The sole purpose of using the temporal interval $\displaystyle[T^{0},T^{*}]$ is to witness the explosive feature of $\displaystyle u_{t}$ immediately prior to quenching. This corresponds to the last 105 time steps in computations. The curvature functions are again selected via the arc-length of the maximal value function of $\displaystyle u.$ It is evident that while $\displaystyle u$ quenches peacefully as $\displaystyle t\rightarrow T^{*},$ the function $\displaystyle u_{t}$ blows up simultaneously.

Figure 5. Three-dimensional curvature views of $\displaystyle u,~{}0\leq t\leq T^{*},$ [LEFT] and $\displaystyle u_{t},~{}T^{0}\leq t\leq T^{*},$ [RIGHT] where $\displaystyle T^{0}=0.509286490538884$ and $\displaystyle T^{*}=0.509391490538887$ are used. The magenta and red curves represent functions $\displaystyle\textstyle\max_{-1\leq x\leq 1}u$ and $\displaystyle\textstyle\max_{-1\leq x\leq 1}u_{t},$ respectively. The temporal derivative values concentrates about the quenching point with $\displaystyle\textstyle\max_{-1\leq x\leq 1}u_{t}(x,T^{*})>985\gg 1.$

Figure 6. Maximal value profiles of the numerical solution $\displaystyle u$ [LEFT] and $\displaystyle u_{t}$ [RIGHT]. The mainframe curves are for $\displaystyle\textstyle\max_{-1\leq x\leq 1}u$ and $\displaystyle\textstyle\max_{-1\leq x\leq 1}u_{t}$ for $\displaystyle 0\leq t\leq 0.509391490538887.$ The embedded graphics represent profiles of $\displaystyle u$ and $\displaystyle u_{t}$ immediately prior to quenching. The quenching time is observed to be $\displaystyle T^{*}\approx 0.509391490538887.$

More details of the maximal value profiles of $\displaystyle u$ and $\displaystyle u_{t}$ can be found in Fig. 6. For additional information, we embed the terminal solution $\displaystyle u$ and corresponding derivative $\displaystyle u_{t}$ into the main frames. A logarithmic scale is used for $\displaystyle u_{t}$ in order to provide a better illustration of the explosive feature of the derivative function. The temporal adaptation is triggered automatically once $\displaystyle\textstyle\max_{-1\leq x\leq 1}u\approx 0.90,$ and remains activated throughout the remainder of computations.

Figure 7. Illustrations of the dynamic connection between quenching times and domain sizes $\displaystyle a$ [LEFT], and the rate of change of the quenching time with respect to $\displaystyle a$ [RIGHT]. As $\displaystyle a$ increases beyond $\displaystyle a_{1},$ quenching time declines rapidly until its minimal value $\displaystyle T^{*}_{a_{307}}\approx 0.499360935318447,$ where $\displaystyle a_{307}=3.8321581.$ Once $\displaystyle a$ increases beyond $\displaystyle a_{307}$ the quenching time increases in a slightly oscillatory manner to reach $\displaystyle T^{*}_{a_{1000}}.$

Fig. 7 shows the effect of domain size on quenching time. The numerical solutions based on 1000 different values of $\displaystyle a\in[0.7652281,10.7552281]$ are computed, compared, and analyzed. The reason for choosing $\displaystyle a_{1}=0.7652281$ is that it is slightly larger than the theoretical critical size of $\displaystyle a^{*}\approx 0.765228037955310$ [1, 8]. The experiments indicate that the quenching time is longer when $\displaystyle a$ is close to $\displaystyle a_{1},$ with a maximum $\displaystyle T_{a_{1}}^{*}=9.752350010587456.$ The quenching time then sharply decreases. These results firmly support theoretical expectations that the quenching time should approach infinity as the domain size decreases to $\displaystyle a^{*}.$ We also observe that the quenching time seems to have a lower threshold. Although the data acquired exhibit slight oscillations, the vibrations should not be caused by round-off errors, since various error-reduction measurements are utilized in our calculations. For example, we let the temporal step minimum $\displaystyle\tau_{\min}=c\times 10^{-6},$ where $\displaystyle c>0$ is an usual floating point number. When a temporal derivative value is evaluated, we reformulate the original formula through modifications such as

[TABLE]

which effectively reduces the risk of unfavorable round-off errors. Consequently, we obtain a final value of $\displaystyle T^{*}_{a_{1000}}\approx 0.515984311015508.$

Experiment 2

Let us consider $\displaystyle a=2,~{}\varphi(\epsilon)\equiv 1$ and utilize the same initial function $\displaystyle u_{0}(x)=0.001(1-\cos(2\pi x)).$ Set our degenerate function $\displaystyle\sigma(x)=(x+1)^{p}(1-x)^{1-p},~{}-1\leq x\leq 1,~{}0\leq p\leq 1.$ Note that $\displaystyle\sigma(x)$ creates a degeneracy near each of the spatial boundaries of the problem $\displaystyle(\ref{num1})$-$\displaystyle(\ref{num3})$. We commence by using the golden ratio $\displaystyle p=(\sqrt{5}-1)/2.$

Figure 8. Two-dimensional thermal flows of the numerical solution $\displaystyle u$ for $\displaystyle t\in[0,T]$ [LEFT], and $\displaystyle u_{t}$ for $\displaystyle t\in[0.905433681825884,0.905541681825887]$ [RIGHT]. In the former case, the heat flows smoothly and monotonically increases until quenching at $\displaystyle P=(-0.378707538403295,0.905541681825887).$ Numerical solutions in last 105 temporal steps immediately before quenching are used for estimating $\displaystyle u_{t}.$

Figure 9. Three-dimensional views of $\displaystyle u$ [LEFT] and $\displaystyle u_{t}$ [RIGHT]. The blue three-dimensional curvature plots represent profiles of $\displaystyle u$ and $\displaystyle u_{t}$ at different times, respectively. While the magenta curves indicate $\displaystyle\textstyle\max_{-1\leq x\leq 1}u,$ the red curves are for $\displaystyle\textstyle\max_{-1\leq x\leq 1}u_{t}.$

We may observe in Fig. 8 and 9 how the maximal values have shifted away from the center due to the degeneracy. The velocity map of $\displaystyle u_{t}$ matches that of $\displaystyle u$ in the thermal plots. To view more clearly the explosive profile of $\displaystyle u_{t},$ a time interval $\displaystyle[0.905433681825884,0.905541681825887]$ is used in the second frames. Note that the temporal derivative $\displaystyle u_{t}$ reaches its maximum at the quenching point with $\displaystyle\textstyle\max_{x,t}u_{t}\approx 475.5902863402550.$ Further, in Fig. 9 and 10, we can observe that there is a remarkable shift of the local maximal values to the left of the origin when approaching quenching, as compared to cases without a degeneracy. The temporal adaptation is activated once $\displaystyle\textstyle\max_{-1\leq x\leq 1}u\approx 0.90$ and remains active throughout the remainder of the computations.

Figure 10. Profiles of the maximal values of the numerical solution $\displaystyle u$ [LEFT] and its derivative $\displaystyle u_{t}$ [RIGHT]. The two embedded pictures are for $\displaystyle u$ and $\displaystyle u_{t}$ in the last position immediately before quenching. The quenching location is $\displaystyle x^{*}=-0.378707538403295$ and the quenching time is $\displaystyle T^{*}\approx 0.905541681825887$ in the experiments which are well-agreeable with known results [1, 4, 8, 12].

Figure 11. The left figure displays the relationship between quenching position and values of $\displaystyle p$ in the degeneracy function $\displaystyle\sigma(x).$ The right picture, on the other hand, shows the connection between the quenching time and values of $\displaystyle p$ in $\displaystyle\sigma(x).$ It is known that the quenching position and quenching time distribution should be symmetric about the value $\displaystyle p=0.5$ [7, 10, 12].

Fig. 11 is designed to show possible relations between the quenching location and $\displaystyle p,$ as well as the connections between the quenching time and $\displaystyle p.$ Values of $\displaystyle p$ are specified through the given degeneracy function $\displaystyle\sigma(x).$

It is found in our numerical experiments that the quenching position decays monotonically as $\displaystyle p$ increases. Extreme values, $\displaystyle x^{*}=-0.552238805970151,~{}0.552238805970147,$ are taken as $\displaystyle p=0,~{}1,$ respectively. As expected, the quenching location is at the origin when $\displaystyle p=1/2.$ However, the decay of quenching location function is apparently nonlinear and exhibits a pattern of antisymmetry. It seems that the impact of $\displaystyle p$ on quenching locations is relatively more significant for $\displaystyle p\in[0.4,0.6].$

Further, for the range of $\displaystyle p$ values used, we may observe that the minimal quenching time, $\displaystyle T^{*}\approx 0.394063444318618,$ occurs as $\displaystyle p$ approaches either $\displaystyle 0^{+}$ or $\displaystyle 1^{-}.$ On the other hand, the maximum quenching time, $\displaystyle T^{*}\approx 0.964575637131343$ can be witnessed at $\displaystyle p=0.5.$ The quenching time function is symmetric about $\displaystyle p=0.5$ but again nonlinear. The impact of $\displaystyle p$ on the quenching time is relatively more pronounced when $\displaystyle p$ is closer to the end of its defined interval, that is, when $\displaystyle p\in[0,0.12]$ or $\displaystyle p\in[0.88,1].$

Experiment 3

We proceed with $\displaystyle a=2$ and $\displaystyle\sigma(x)\equiv 1$ in this particular exploration. In order to study the effects of the stochastic influence we consider the function $\displaystyle\varphi(\epsilon)=\epsilon^{2},~{}0.01\leq\epsilon\leq 1.$ In this example we explore the numerical solutions resulting from two different white noise vectors $\displaystyle\epsilon(x).$ The purpose of these considerations is to investigate how slight changes in the vector $\displaystyle\epsilon(x)$ can result in drastically different solutions profiles.

Figure 12. Plots of two different stochastic function values corresponding to $\displaystyle\varphi(\epsilon)=\epsilon^{2}$ are shown. Herewith we have $\displaystyle 0.01\leq\epsilon\leq 1.$

Figure 13. Three-dimensional curvature views of $\displaystyle u$ [LEFT] and $\displaystyle u_{t}$ [RIGHT] corresponding to the first set of random values. The blue lines represent profiles of $\displaystyle u$ and $\displaystyle u_{t}$ at different times, while the magenta and red lines represent $\displaystyle\textstyle\max_{-1\leq x\leq 1}u$ and $\displaystyle\textstyle\max_{-1\leq x\leq 1}u_{t},$ respectively. The temporal derivative has its largest values concentrated about the quenching point with $\displaystyle\max u_{t}\approx 484.1416720746672.$

The two different random variable function outputs used are displayed in Fig. 12. Note that the vector $\displaystyle\epsilon(x)=\left(\epsilon_{1},\dots,\epsilon_{N}\right)^{\intercal}$ is generated randomly, with each component consisting of a uniformly distributed random number $\displaystyle\epsilon_{i}\in[0.01,1],~{}i=1,\dots,N.$ Three-dimensional curvature views of the numerical solution $\displaystyle u$ for $\displaystyle t\in[0,T_{1}]$ and and $\displaystyle u_{t}$ for $\displaystyle t\in[T^{0}_{1},T_{1}]$ corresponding the left figure are given in Fig. 13. Values of $\displaystyle T^{0}_{1}=1.848738391680962$ and $\displaystyle T_{1}=1.848843391680954$ are used in the second frame. Again, the heat monotonically increases from left to right until quenching at $\displaystyle P=(0.020653228859545,1.848843391680954).$ On the other hand, the velocity $\displaystyle u_{t}$ maximum trajectory matches that of $\displaystyle u$ and explodes at the aforementioned quenching position $\displaystyle P$ in an extremely short period of time. Numerical solutions in last 105 temporal steps immediately before quenching are used.

Figure 14. Profiles of the maximal values of the numerical solution $\displaystyle u$ [LEFT] and its derivative $\displaystyle u_{t}$ [RIGHT]. The two embedded pictures are for $\displaystyle u$ and $\displaystyle u_{t}$ in the last position immediately before quenching.

Note from Fig. 13 and 14 how the location of quenching has shifted, just as in the degeneracy case, but this solution is slightly to the right of the origin. The most notable feature is the lack of an initial smoothness in the solution profile of $\displaystyle u.$ The temporal adaptation is trigged once $\displaystyle\textstyle\max_{-1\leq x\leq 1}u\approx 0.90$ and remains active throughout the remainder of computations. The quenching is found at $\displaystyle x^{*}=0.066585749194076$ and the quenching time is $\displaystyle T^{*}_{1}\approx 1.848843391680954.$

Figure 15. Three-dimensional curvature views of $\displaystyle u$ [LEFT] and $\displaystyle u_{t}$ [RIGHT]. The blue lines represent profiles of $\displaystyle u$ and $\displaystyle u_{t}$ at different times, while the magenta and red lines represent $\displaystyle\textstyle\max_{-1\leq x\leq 1}u$ and $\displaystyle\textstyle\max_{-1\leq x\leq 1}u_{t},$ respectively. The temporal derivative has its largest values concentrated about the quenching point with $\displaystyle\max u_{t}\approx 885.0797753481299.$

Figure 16. Profiles of the maximal values of the numerical solution $\displaystyle u$ [LEFT] and its derivative $\displaystyle u_{t}$ [RIGHT]. The two embedded pictures are for $\displaystyle u$ and $\displaystyle u_{t}$ in the last position immediately before quenching. The quenching location is $\displaystyle x^{*}=-0.043597691787030$ and the quenching time is $\displaystyle T^{*}_{2}\approx 2.159108137916605.$

Fig. 15 and 16 are for the case when the second random output in Fig. 12 is selected. $\displaystyle T^{0}_{2}=2.159003137916590$ and $\displaystyle T_{2}=2.159108137916605$ are utilized for $\displaystyle u_{t}.$ The heat $\displaystyle u$ flows smoothly and increases monotonically until quenching at $\displaystyle P=(0.020653228859545,2.159108137916605).$ On the other hand, the trajectory of the maximal velocity of $\displaystyle u_{t}$ matches that of $\displaystyle u.$ It explodes at the aforementioned quenching position $\displaystyle P$ to the peak value. Numerical solutions in last 105 temporal steps immediately before quenching are used for $\displaystyle u_{t}.$ It can also be seen that the stochastic term causes the maximum value of the solution to shift drastically. The temporal adaptation is activated once $\displaystyle\textstyle\max_{-1\leq x\leq 1}u\approx 0.90$ and remains active throughout the remainder of the computations.

Figure 17. More detailed three-dimensional surface views of $\displaystyle u$ [LEFT] and $\displaystyle u_{t}$ [RIGHT]. Both surfaces increase monotonically with respect to the time. The function $\displaystyle u$ and its temporal derivative have their largest values concentrated about the quenching point with $\displaystyle\max u\approx 0.998883860103612$ and $\displaystyle\max u_{t}\approx 799.5162246152710.$ Quenching time $\displaystyle T^{*}\approx 1.902752501571197.$

To conclude this experiment, we rerun the second case for surface plots of the last 105 $\displaystyle t$-level numerical solutions before quenching in Fig. 17. Note that in all sets of figures the quenching time is considerably longer. This is due to the fact that the constraints on $\displaystyle\epsilon$ dampen the effects of the source term, thereby increasing the amount of time for a solution component to reach unity. Further note how small differences in the random function $\displaystyle\varphi(\epsilon)$ can shift the quenching location to either side of the origin. Further, it can be observed in Fig. 13-17 that the solution $\displaystyle u$ becomes quite non-smooth near the quenching point. We may also note how differently the maximum values of the solution $\displaystyle u$ flows with respect to time in Fig. 13, 15 and 17. While the maximal values converge to the appropriate spatial position fairly quickly, their paths to that point are extremely different.

Since the Lax equivalence theorem does not apply to nonlinear schemes such as $\displaystyle(\ref{c3})$, rigorous tests on the nonlinear convergence are in general important. But our experiments are carried out in cases when nonlinear source terms are linearized, similar to those in [8, 10, 20]. This linearization allows for the problem to be locally considered linear, and in this sense, convergence can be ensured through the numerical stability. Note that such an argument is only valid in the local sense, a global convergence analysis of the numerical method $\displaystyle(\ref{c3})$, or $\displaystyle(\ref{approx})$, must be conducted via chemical-physical energy conservations. Such an endeavor has been in our agenda and will be discussed carefully in our forthcoming papers.

Experiment 4

We now consider the following two-dimensional problem:

[TABLE]

where $\displaystyle T<\infty,~{}\theta>0,~{}u_{0}(x,y)\in C[[-1,1]\times[-1,1]],$ and $\displaystyle 0\leq u_{0}\ll 1.$ Without loss of generality, we choose $\displaystyle a=b=2,$ $\displaystyle\theta=1,$ $\displaystyle t_{0}=0,$ and $\displaystyle u_{0}(x,y)=0.001(1-\cos(2\pi x))(1-\cos(2\pi y)),$ $\displaystyle-1\leq x,y\leq 1.$ For the simplicity we once again set $\displaystyle\sigma(x,y)\equiv 1,~{}0\leq x,y\leq 1,$ and $\displaystyle\varphi(\varepsilon)=\varepsilon^{2},~{}0.01\leq\varepsilon\leq 1.$ We wish to demonstrate that the stochastic effects exhibited in the previous experiment remain to be significant in the two-dimensional case. To this end, we first approximate $\displaystyle(\ref{num11})$-$\displaystyle(\ref{num33})$ via a standard nonuniform LOD method [10, 20, 22].

Figure 18. Three-dimensional views of $\displaystyle u$ [LEFT] and $\displaystyle u_{t}$ [RIGHT] immediately prior to quenching. Quenching is observed to occur at $\displaystyle T\approx 2.564180137941836.$ The temporal derivative is observed to reach a maximal value $\displaystyle\textstyle\max u_{t}\approx 331.0481243466621.$ The blow-up of $\displaystyle u_{t}$ is concentrated around the observed spatial quenching point.

Fig. 18 depicts the final profiles of the solution $\displaystyle u$ and its derivative $\displaystyle u_{t},$ immediately prior to quenching. Quenching is observed to occur at $\displaystyle T\approx 2.564180137941836.$ We have experienced a slight increase in quenching times as compared to existing results for $\displaystyle\varphi(\varepsilon)\equiv 1$ [8, 10, 20].

Figure 19. Projections of $\displaystyle u$ [LEFT] and $\displaystyle u_{t}$ [RIGHT] onto the $\displaystyle xu$-plane and $\displaystyle xu_{t}$-plane immediately prior to quenching, respectively. A slight shift of the location of maximum values is observed. It is evident that the solution smoothness is once again affected slightly by the presence of the stochastic influence in reactions.

[TABLE]

Table 1. Experimented quenching locations $\displaystyle P_{k}^{*}=(x_{k}^{*},y_{k}^{*}),~{}k=1,\ldots,5,$ due to stochastic influences. Drifting effects can be clearly observed.

Fig. 19 depicts projections of $\displaystyle\textstyle\max_{-1\leq y\leq 1}u(x,y,t)$ and $\displaystyle\textstyle\max_{-1\leq x\leq 1}u_{t}(x,y,t),$ respectively. In the absence of stochastic influences, the theoretical spatial quenching location should be $\displaystyle P=(0,0).$ But we have observed $\displaystyle P^{*}=(0.034825870646764,-0.014925373134330)$ to be a shifted spatial quenching position under the influence of present stochastic source term. Repeated experiments with varying stochastic influence functions suggest the same phenomenon. We list results of five randomly selected stochastic functions in Table 1 as an illustration. We note similarities between the plots in Fig. 19 and those for one-dimensional cases in quenching location disturbances. This is to be expected, as the results demonstrate that the propagation of the nonsmooth feature is not limited to one-dimensional situations only.

6 Conclusions

A nontraditional Crank-Nicolson method for solving the nonlinear stochastic Kawarada differential equation is proposed and studied. Conventional uniform or symmetric mesh structures are replaced successfully by fully arbitrary grids in the space. Temporal adaptation is incorporated in order to effectively capture the strong quenching-combustion singularity and degeneracy built with the nonlinear partial differential equation. Key properties of the numerical method developed, including the solution positivity, monotonicity, and stability, are investigated and proven. Stability conditions determined are nonrestrictive. Stochastic impacts through the source term are examined and discussed carefully through simulation experiments.

Although linear stability analysis has been effective in the study of numerical solution of quenching problems while nonlinear source terms are frozen in implicit schemes [10, 11, 12, 17, 21], an improved semi-linear stability analysis is introduced and conducted. This modified analysis depends upon the Jacobian of the nonlinear reaction term of the Kawarada equation. It is found that the constraints required to guarantee the positivity and monotonicity of the underlying nonuniform numerical method are sufficient for ensuring the aforementioned semi-linear stability. This sheds further insights as to reasons why a linear stability analysis is often adequate in realistic computational applications.

In our numerical experiments we have studied effects of the size of spatial domain on quenching time, which seems to suggest a possible optimal domain size due to the minimum quenching time observed. Further, we have explored the effects of the stochastic source term on overall solutions. Computational experiments indicate that different nonlinear source terms may have impacts on not only quenching time and location, but also the smoothness of the quenching solution profile [2, 5, 12]. A two-dimensional experiment is also presented to verify the potential of the semi-adaptive infrastructure introduced in this study, as well as verify the effects of a stochastic source influence in higher-dimensional cases. Our future endeavors include studying stochastic influences which vary with respect to both time and space. Multi-dimensional Kawarada problems will also be approximated via the latest operator splitting strategies [8, 11]. On the other hand, exponential time differencing schemes, such as those explored in [13, 14], will be considered in the near future together with proper adaptations. We also plan to extend our investigations of literature by including balanced fractional derivatives in order to more precisely capture and explore global features of the numerical combustion [23].

Acknowledgements

The authors would like to thank the anonymous referees for their time and thorough comments. Integrating their suggestions has undoubtedly elevated the quality and presentation of this paper.

The second author particularly appreciates Wes Johnson and Mike Hutcheson of the Baylor University Academic and Research Computing Services for computational validations and technical support.