On the positivity, monotonicity, and stability of a semi-adaptive LOD method for solving three-dimensional degenerate Kawarada equations

TL;DR

This paper develops a semi-adaptive LOD method for solving complex 3D degenerate Kawarada equations, ensuring positivity, monotonicity, and stability despite nonlinearities and singularities.

Contribution

It introduces a semi-adaptive LOD scheme with a novel arc-length based temporal adaptation and stability analysis for challenging degenerate PDEs.

Findings

Positivity and monotonicity criteria are established.

Stability of the splitting method is proven in the von Neumann sense.

The method effectively handles highly nonlinear source terms and singularities.

Abstract

This paper concerns the numerical solution of three-dimensional degenerate Kawarada equations. These partial differential equations possess highly nonlinear source terms, and exhibit strong quenching singularities which pose severe challenges to the design and analysis of highly reliable schemes. Arbitrary fixed nonuniform spatial grids, which are not necessarily symmetric, are considered throughout this study. The numerical solution is advanced through a semi-adaptive Local One-Dimensional (LOD) integrator. The temporal adaptation is achieved via a suitable arc-length monitoring mechanism. Criteria for preserving the positivity and monotonicity are investigated and acquired. The numerical stability of the splitting method is proven in the von Neumann sense under the spectral norm. Extended stability expectations are proposed and investigated.

Equations

$$s(x,y,z)u_t=u_{xx}+u_{yy}+u_{zz}+f(u),(x,y,z,t)\in \Omega ,$$

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

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

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

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

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

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

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

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

$${\frac{{\partial }^{2}u}{\partial x^2}|}_{i,j,k}$$

$${\frac{{\partial }^{2}u}{\partial y^2}|}_{i,j,k}$$

$${\frac{{\partial }^{2}u}{\partial z^2}|}_{i,j,k}$$

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

$$v(t_0)$$

$$M_1=\frac{1}{a^2}B(I_{N_3}\otimes I_{N_2}\otimes T_1),M_2=\frac{1}{b^2}B(I_{N_3}\otimes T_2\otimes I_{N_1}),M_3=\frac{1}{c^2}B(T_3\otimes I_{N_2}\otimes I_{N_1}),$$

$$B$$

$${\varphi }_{i,j,k}$$

$$T_{\sigma }$$

$$l_{\sigma ,j}$$

$$m_{\sigma ,j}$$

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

$$p(t)=\prod _{\sigma =1}^3{\left(I-\frac{t}{2}{M}_{\sigma }\right)}^{-1}\left(I+\frac{t}{2}{M}_{\sigma }\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{)}^2\right),t\to t_0.$$

$$v_{\ell +1}=\left[\prod _{\sigma =1}^3{\left(I-\frac{{\tau }_{\ell }}{2}{M}_{\sigma }\right)}^{-1}\left(I+\frac{{\tau }_{\ell }}{2}{M}_{\sigma }\right)\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 }(C{v}_{\ell }+g(v_{\ell })),0<{\tau }_{\ell }\ll 1,$$

$$m\left(\frac{\partial v}{\partial t},t\right)=\sqrt{1+{\left(\frac{{\partial }^{2}v}{\partial t^2}\right)}^2},t_0<t<T.$$

$${\tau }_{\ell }^2={\tau }_{\ell -1}^2+{\left(\frac{\partial v_{\ell -1}}{\partial t}-\frac{\partial v_{\ell -2}}{\partial t}\right)}^2-{\left(\frac{\partial v_{\ell }}{\partial t}-\frac{\partial v_{\ell -1}}{\partial t}\right)}^2,\ell =1,2,3,\dots ,$$

$$\displaystyle T_{1}^{\intercal}T_{1}=\left(\begin{array}[]{ccccccc}\tilde{m}_{1,1}&\tilde{n}_{1,1}&\tilde{N}_{1,1}&&&&\\ \tilde{l}_{1,1}&\tilde{m}_{1,2}&\tilde{n}_{1,2}&\tilde{N}_{1,2}&&&\\ \tilde{L}_{1,1}&\tilde{l}_{1,2}&\tilde{m}_{1,3}&\tilde{n}_{1,3}&\tilde{N}_{1,3}&&\\ &&\cdots&\cdots&\cdots&\cdots&\\ &&\tilde{L}_{1,N_{1}-4}&\tilde{l}_{1,N_{1}-3}&\tilde{m}_{1,N_{1}-2}&\tilde{n}_{1,N_{1}-2}&\tilde{N}_{1,N_{1}-2}\\ &&&\tilde{L}_{1,N_{1}-3}&\tilde{l}_{1,N_{1}-2}&\tilde{m}_{1,N_{1}-1}&\tilde{n}_{1,N_{1}-1}\\ &&&&\tilde{L}_{1,N_{1}-2}&\tilde{l}_{1,N_{1}-1}&\tilde{m}_{1,N_{1}}\end{array}\right)\in{\mathbb{R}}^{N_{1}\times N_{1}},$$

$${\tilde{N}}_{1,j}$$

$${\tilde{n}}_{1,j}$$

$${\tilde{m}}_{1,j}$$

$$|{\lambda }_{1,j}-{\tilde{m}}_{1,j}|\le |{\tilde{L}}_{1,j-2}|+|{\tilde{l}}_{1,j-1}|+|{\tilde{n}}_{1,j}|+|{\tilde{N}}_{1,j}|,j=3,\dots ,N_1-2,$$

$$-|m_{1,j-1}{n}_{1,j-1}+m_{1,j}{l}_{1,j-1}|-|m_{1,j}{n}_{1,j}+m_{1,j+1}{l}_{1,j}|$$

$$-|l_{1,j-2}{n}_{1,j-1}|-|l_{1,j}{n}_{1,j+1}|+n_{1,j-1}^2+m_{1,j}^2+l_{1,j}^2$$

$${\lambda }_{1,j}$$

$$-{\lambda }_{1,j}$$

Full text

AMS Subject Claasifications - 65K20, 65M50, 35K65, 35B40

1 Introduction

Let $\displaystyle{\cal{D}}=(0,a)\times(0,b)\times(0,c)\subset{\mathbb{R}}^{3},$ where $\displaystyle a,b,c>0,$ and $\displaystyle\partial{\cal{D}}$ be its boundary. Denote $\displaystyle\Omega={\cal{D}}\times(t_{0},~{}T),~{}\SS=\partial{\cal{D}}\times(t_{0},~{}T)$ for given $\displaystyle 0\leq t_{0}<T<\infty.$ We consider the following degenerate Kawarada problem,

[TABLE]

where $\displaystyle s(x,y,z)=\left(x^{2}+y^{2}+z^{2}\right)^{q/2},~{}q\in[0,2].$ The nonlinear source function, $\displaystyle f(u),$ is strictly increasing for $\displaystyle 0\leq u<1$ with

[TABLE]

In idealized thermal combustion applications [2, 3, 17], $\displaystyle u$ represents the temperature in the combustion channel, and the $\displaystyle x$-, $\displaystyle y$-, and $\displaystyle z$-coordinates coincide with the channel walls. The initial temperature $\displaystyle 0\leq u_{0}\ll 1$ is typically chosen to be small. The function $\displaystyle s(x,y,z)$ represents certain singularities in the temperature transportation speed within the channel, which causes the degeneracy in the differential equation $\displaystyle(\ref{b1})$ [4, 15, 18, 21]. The solution $\displaystyle u$ of $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ is said to quench if there exists a finite time $\displaystyle T>0$ such that

[TABLE]

The value $\displaystyle T$ is then defined as the quenching time [1, 2, 14]. It has been shown that a necessary condition for quenching to occur is

[TABLE]

Further, such a $\displaystyle T$ exists only when certain spatial references, such as the size and shape of $\displaystyle{\cal{D}},$ reach their critical limits. A domain $\displaystyle{\cal{D}}^{*}$ is called the critical domain if the solution of $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ exists for all time when $\displaystyle{\cal{D}}\subseteq{\cal{D}}^{*},$ and $\displaystyle(\ref{a2})$ occurs when $\displaystyle{\cal{D}}^{*}\subseteq{\cal{D}}$ for a finite $\displaystyle T$ [14].

Systematic mathematical investigations of quenching phenomena can be traced back to Karawada’s original work involving the one-dimensional model equation [12]. It was observed that for any spatial domain $\displaystyle[0,a],$ there exists a unique value $\displaystyle a^{*}>0$ such that for $\displaystyle a<a^{*},$ the solution of the equation exists globally; and for $\displaystyle a\geq a^{*},$ there exists a finite time $\displaystyle T(a),$ such that $\displaystyle\lim_{t\rightarrow T(a)}\max_{0\leq x\leq a}u(x,t)=1.$ In the latter case, $\displaystyle u$ stops existing in finite time and this phenomenon is referred to as quenching [12, 14, 20]. There have been considerable developments in the study of Karawada equations, although discussions of multidimensional problems were extremely limited until recently. In 1994, Chan and Ke proved that for a domain $\displaystyle{\cal{D}}=(0,a)\times(0,b)\subset{\mathbb{R}}^{2},$ if $\displaystyle f,~{}f_{u}$ are nonnegative, then, for any fixed ratio $\displaystyle a/b,$ there exists a unique critical domain $\displaystyle{\cal{D}}^{*}$ for $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$, and the solution of the differential equation problem is unique before quenching [6]. A numerical approximation of the relationship between $\displaystyle a/b$ and the areas of $\displaystyle{\cal{D}}^{*}$ of a nondegenerate ($\displaystyle q=0$) problem was given. These results have been well supported by realistic physical processes, in particular in solid fuel combustion [3, 4, 20].

Numerous computational procedures, including moving mesh adaptive methods, have been constructed for solving blow-up and Kawarada problems in the past decades (interested readers are referred to [1, 7, 8, 18, 21] and references therein). Though in the former case, adaptations are frequently achieved via monitoring functions on the arc-length of the function $\displaystyle u;$ in the latter situation, adaptations are more likely to be built upon the arc-length of $\displaystyle u_{t},$ since it is directly proportional to $\displaystyle f(u),$ which blows up as $\displaystyle u$ quenches [6, 14, 19].

As reported in several recent investigations, when quenching locations can be predetermined, it is preferable to use nonuniform spatial grids throughout the computations [4, 13, 21]. In this case, key quenching characteristics such as the quenching time and critical domain, are more easily observed; Also important numerical properties of underlying algorithms, including the monotonicity, stability and convergence, can be more precisely studied. To that end, this paper develops a temporally adaptive splitting scheme utilizing predetermined nonuniform spatial grids. The positivity, monotonicity, and stability of the method will be investigated. It is also observed that the impact of degeneracy is limited for our implicit scheme. Our discussions will be organized as follows. In the next section, the semi-adaptive LOD scheme for solving $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ will be constructed and discussed. Then, in Section 3, criteria to guarantee the positivity of the numerical scheme will be determined. In Section 4, appropriate criteria for guaranteeing the monotonicity will be obtained. These two sections together serve as the platform for carrying out investigations of stability. Section 5 is devoted to the stability analysis of the semi-adaptive LOD scheme. The analysis will first be carried out for a fully linearized scheme, and then a more realistic stability analysis is proposed without freezing the source term. Finally, concluding remarks and proposed future work will be given in Section 6. For now, no numerical studies of the three-dimensional degenerate Karawada problem will be given.

2 Semi-adaptive LOD scheme

Utilizing the transformations $\displaystyle\tilde{x}=x/a,~{}\tilde{y}=y/b,~{}\tilde{z}=z/c,$ and reusing the original variables for simplicity, we may reformulate $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ as

[TABLE]

where $\displaystyle g(u)=f(u)/\phi,~{}\phi=\phi(x,y,z)=\left(a^{2}x^{2}+b^{2}y^{2}+c^{2}z^{2}\right)^{q/2},$ and $\displaystyle{\cal{D}}=(0,1)\times(0,1)\times(0,1)\subset{\mathbb{R}}^{3}.$

Let $\displaystyle N_{1},N_{2},N_{3}\gg 1.$ We inscribe over $\displaystyle\bar{{\cal{D}}}$ the following variable grid: $\displaystyle{\cal{D}}_{h}=\left\{(x_{i},y_{j},z_{k})|i=0,\dots,N_{1}+1;~{}j=0,\dots,N_{2}+1;~{}k=0,\dots,N_{3}+1;~{}x_{0}=y_{0}=\right.$

$\displaystyle\left.z_{0}=0,~{}x_{N_{1}+1}=y_{N_{2}+1}\right.$ $\displaystyle\left.=z_{N_{3}+1}=1\right\}.$ Denote $\displaystyle h_{1,i}=x_{i+1}-x_{i}>0,~{}h_{2,j}=y_{j+1}-y_{j}>0,$ and $\displaystyle h_{3,k}=z_{k+1}-z_{k}>0$ for $\displaystyle 1\leq i\leq N_{1},~{}1\leq j\leq N_{2},~{}1\leq k\leq N_{3}.$ Let $\displaystyle u_{i,j,k}(t)$ be an approximation of the solution of $\displaystyle(\ref{c1})$-$\displaystyle(\ref{c2b})$ at $\displaystyle(x_{i},y_{j},z_{k},t)$ and consider the following first-order finite differences [21],

[TABLE]

Further, denote $\displaystyle v(t)=(u_{1,1,1},u_{2,1,1},\dots,u_{N_{1},1,1},u_{1,2,1},u_{2,2,1},\dots,u_{N_{1},2,1},\ldots,$

$\displaystyle u_{1,N_{2},1},u_{2,N_{2},1},\dots,u_{N_{1},N_{2},1},\ldots,u_{1,N_{2},N_{3}},u_{2,N_{2},N_{3}},\dots,u_{N_{1},N_{2},N_{3}})^{\intercal}\in{\mathbb{R}}^{N_{1}N_{2}N_{3}}$ and let $\displaystyle g(v)$ be a discretization of the nonhomogeneous term of $\displaystyle(\ref{c1})$. We obtain readily from $\displaystyle(\ref{c1})$-$\displaystyle(\ref{c2b})$ the following semi-discretized system

[TABLE]

where

[TABLE]

$\displaystyle\otimes$ stands for the Kronecker product, $\displaystyle I_{N_{\sigma}}\in{\mathbb{R}}^{N_{\sigma}\times N_{\sigma}},~{}\sigma=1,2,3,$ are identity matrices, and

[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 and $\displaystyle C=\sum_{\sigma=1}^{3}M_{\sigma}$ [17].

In principle, different approximation techniques can be used to yield different splitting methods based on $\displaystyle(\ref{rs1})$ [10, 17, 19]. Yet, we are particularly interested in approximating $\displaystyle(\ref{rs1})$ via a trapezoidal rule and a [1/1] Padé approximation, $\displaystyle E(tC)=p(t)+{\cal{O}}\left(t^{2}\right),$ where

[TABLE]

The above leads to

[TABLE]

The above LOD algorithm provides a highly efficient way to compute numerical solutions of multidimensional problems such as $\displaystyle(\ref{c1})$-$\displaystyle(\ref{c2b})$ [10, 16, 18, 21]. Based on $\displaystyle(\ref{cc3})$, we obtain the following first order in space and time semi-adaptive LOD scheme:

[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}+\sum_{k=0}^{\ell-1}\tau_{k},~{}\ell=0,1,2,\ldots,$ and $\displaystyle\{\tau_{\ell}\}_{\ell\geq 0}$ is a set of variable temporal steps determined by an adaptive procedure. 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.

Figure 1.  Numerical solution (left) and its temporal derivative (right) immediately before quenching. It is observed that as $\displaystyle\max_{x}v(x)\rightarrow 1^{-},$ we have $\displaystyle\max_{x}v_{t}\gg 600.$ The computed quenching time is $\displaystyle T\approx 0.780265747310047.$

Due to a strong quenching singularity, the selection of proper nonuniform temporal steps $\displaystyle\tau_{\ell}$ is vital. As an illustration, in Figure 1, we show the numerical solution and its temporal derivative of a typical one-dimensional Kawarada problem over the interval $\displaystyle[0,\pi].$ The initial function $\displaystyle u_{0}(x)=0.001\sin(x),~{}f(u)=1/(1-u),$ and homogeneous Dirichlet boundary condition are employed. The degeneracy function utilized is $\displaystyle s(x)=x^{p}(\pi-x)^{1-p}$ with $\displaystyle p=(\sqrt{5}-1)/2.$ It is evident that $\displaystyle v_{t}$ changes dramatically when compared with $\displaystyle v.$ Recalling $\displaystyle(\ref{a1})$ and $\displaystyle(\ref{a2})$, we consider the following arc-length monitoring function on $\displaystyle v_{t},$

[TABLE]

Setting the two maximal arc-lengths in neighboring intervals $\displaystyle[t_{\ell-2},~{}t_{\ell-1}]$ and $\displaystyle[t_{\ell-1},~{}t_{\ell}]$ equal [9, 13, 20, 21], we acquire the following quadratic equations from the above,

[TABLE]

with $\displaystyle\tau_{0}$ given.

In the above temporal adaptation procedures, we may consider a minimal temporal step size controller $\displaystyle\tilde{\tau}_{0},~{}0<\tilde{\tau}_{0}\ll\tau_{0},$ to avoid sudden changes in grid movements or unnecessarily large numbers of computations. Further, let $\displaystyle\wedge$ be one of the operations $\displaystyle<,~{}\leq,~{}>,~{}\geq$ and $\displaystyle\alpha,\beta\in{\mathbb{R}}^{N}.$ We assume the following notations in subsequent discussions:

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

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

3 Positivity

The positivity property is one of the most profound characteristics of the solution of the Kawarada problem $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ or $\displaystyle(\ref{c1})$-$\displaystyle(\ref{c2b})$ [1, 2, 6, 14]. Since positive computational solutions preserve the correct physical features of quenching phenomena, it is crucial that our numerical solution also possesses this feature.

Lemma 3.1**.**

$\displaystyle\|T_{\sigma}\|_{2}\leq\max_{j=1,\dots,N_{\sigma}}\left\{{4}/{h_{\sigma,j}^{2}}\right\},~{}~{}~{}\sigma=1,2,3$**

Proof.

Due to the similarity in structure, we only consider $\displaystyle T_{1}$ since the other two cases follow by similar arguments. Note that, in general, $\displaystyle T_{1}$ is not symmetric. However, $\displaystyle\|T_{1}\|_{2}^{2}=\rho(T_{1}^{\intercal}T_{1})$ and $\displaystyle T_{1}^{\intercal}T_{1}$ is symmetric with a bandwidth of five. Thus,

[TABLE]

where

[TABLE]

We may determine a bound on the spectral radius of $\displaystyle T_{1}^{\intercal}T_{1}$ by using Gers̆chgorin’s circle theorem. In fact, only rows containing five nontrivial elements, i.e., $\displaystyle j=3,\dots,N_{1}-2,$ need to be considered. To this end,

[TABLE]

which gives

[TABLE]

and

[TABLE]

Let $\displaystyle h_{1}\equiv\min_{j=1,\dots,N_{1}}\{h_{1,j}\}.$ From $\displaystyle(\ref{g2})$ we acquire that

[TABLE]

Now, reverse $\displaystyle(\ref{g1})$ and by the same token,

[TABLE]

where, once again, $\displaystyle h_{1}\equiv\min_{j=1,\dots,N_{1}}\{h_{1,j}\}.$ Thus, combining the bounds we have $\displaystyle\|T_{1}\|_{2}\leq\max_{i=1,\dots,N_{1}}\left\{{4}/{h_{1,j}^{2}}\right\}.$ The other bounds follow similarly.

Lemma 3.2**.**

Let

[TABLE]

If

[TABLE]

then the matrices

[TABLE]

are nonsingular. Further, $\displaystyle I-\frac{\tau_{\ell}}{2}M_{\sigma},~{}\sigma=1,2,3,$ are monotone and inverse-positive, and $\displaystyle I+\frac{\tau_{\ell}}{2}M_{\sigma},~{}\sigma=1,2,3,$ are nonnegative.

Proof.

First, note that

[TABLE]

Hence, $\displaystyle I+\frac{\tau_{\ell}}{2}M_{1}$ is nonsingular, and also nonnegative. Similar arguments give that $\displaystyle I+\frac{\tau_{\ell}}{2}M_{2}$ and $\displaystyle I+\frac{\tau_{\ell}}{2}M_{3}$ are nonsingular and nonnegative.

Now, consider $\displaystyle A=I-\frac{\tau_{\ell}}{2}M_{1}.$ As $\displaystyle A_{ij}\leq 0$ for $\displaystyle i\neq j$ and the weak row sum criterion is satisfied, $\displaystyle A$ is monotone, and hence an inverse exists and is nonnegative. So, $\displaystyle A$ must be inverse-positive [11]. Similar arguments can be given for $\displaystyle I-\frac{\tau_{\ell}}{2}M_{2}$ and $\displaystyle I-\frac{\tau_{\ell}}{2}M_{3}.$ This ensures the proof.

We also need the following lemma.

Lemma 3.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.

The proof is a straightforward application of the definitions.

4 Monotonicity

Another key characteristic which distinguishes a solution to a quenching problem from a solution to most blow-up problems is its monotonicity with respect to time $\displaystyle t\geq t_{0}$ [1, 2, 6, 14, 18]. Thus, it is necessary to guarantee that our numerical solution preserves this property strictly while solving the Kawarada equation $\displaystyle(\ref{b1})$-$\displaystyle(\ref{b3})$ or $\displaystyle(\ref{c1})$-$\displaystyle(\ref{c2b})$.

Lemma 4.1**.**

If $\displaystyle(\ref{cfl})$ holds for all $\displaystyle\ell\geq k\geq 0,$ and

• (a)

$\displaystyle Cv_{0}+\frac{1}{2}g(v_{0})>0;$**

• (b)

$\displaystyle\left(I-\frac{\tau_{0}}{2}g_{v}(\xi_{0})\right)^{-1}>0$**

hold, 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.

By $\displaystyle(\ref{cfl})$ we have $\displaystyle\left\|\frac{\tau_{k}}{2}M_{\sigma}\right\|<1,$ and thus,

[TABLE]

From $\displaystyle(\ref{c3})$ and the above, we have

[TABLE]

as $\displaystyle\tau_{k}\rightarrow 0.$ Note that $\displaystyle g(v_{k+1})=g(v_{k})+g_{v}(\xi_{k})(v_{k+1}-v_{k})$ for some $\displaystyle\xi_{k}\in{\mathcal{L}}(v_{k+1};v_{k}),$ where $\displaystyle{\mathcal{L}}(v_{k+1};v_{k})$ is the line segment connecting $\displaystyle v_{k+1}$ to $\displaystyle v_{k}$ in $\displaystyle{\mathbb{R}}^{N_{1}N_{2}N_{3}}.$ Using this fact and rearranging terms in $\displaystyle(\ref{m1})$ we have

[TABLE]

and thus,

[TABLE]

We now proceed by induction. Letting $\displaystyle k=0,$ we have

[TABLE]

Thus, if $\displaystyle\tau_{0}$ is sufficiently small, we have $\displaystyle v_{1}-v_{0}>0$ by our assumption and then Lemma 3.3. For the sake of induction, assume that the monotonicity holds for $\displaystyle k=\ell-1.$ Then we have

[TABLE]

Note that $\displaystyle g(v)$ is strictly increasing since $\displaystyle f(v)$ is strictly increasing. Utilizing Lemmas 3.2-3.3 we find that $\displaystyle v_{\ell+1}-v_{\ell}>0$ if $\displaystyle v_{\ell}-v_{\ell-1}>0,$ which completes the induction.

It is not uncommon to set $\displaystyle v_{0}\equiv 0$ in practical combustion simulations. The following corollary shows that in this case conditions in Lemma 4.1 are satisfied for $\displaystyle\ell=0.$

Corollary 4.1**.**

If $\displaystyle v_{0}\equiv 0$ and $\displaystyle\tau_{0}<\min\left\{\beta_{\min}\min\{a^{2},b^{2},c^{2}\},\min_{i,j,k}\frac{2\phi_{i,j,k}}{f_{v}(\xi_{0}(x_{i},y_{j},z_{k}))}\right\},$ then conditions (a), (b) are true.

Proof.

We first consider (a):

[TABLE]

which follows from $\displaystyle f(0)=f_{0}>0.$

We now consider (b), and under these circumstances we need to show

[TABLE]

First, we note that $\displaystyle g_{v}(\xi_{0})$ is diagonal by definition, since

[TABLE]

and

[TABLE]

Let us denote

[TABLE]

It follows readily that

[TABLE]

and (b) holds if

[TABLE]

Denote $\displaystyle d^{(0)}\equiv\max_{i,j,k}\left\{d_{i,j,k}^{(0)}\right\},$ then $\displaystyle\tau_{0}d^{(0)}<2$ which leads to (b).

Lemma 4.2**.**

For any $\displaystyle\tau_{\ell}>0$ we have

[TABLE]

where $\displaystyle x=(1,1,\dots,1)^{\intercal}.$

Proof.

We only need to show the case with

[TABLE]

First, we observe that

[TABLE]

Second, for $\displaystyle i=2,\dots,N_{1}-1$ we have

[TABLE]

Third, we have

[TABLE]

Hence, we conclude that $\displaystyle w_{i,1,1}\geq 1,~{}i=1,\dots,N_{1}.$ Similar arguments may show that all remaining elements of $\displaystyle w$ are also bounded below by 1. Therefore we have $\displaystyle w\geq x.$ Similar discussions may be utilized for the cases involving $\displaystyle M_{2}$ or $\displaystyle M_{3}.$

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. To this end, we denote $\displaystyle h=\max_{j=1,\ldots,N_{\sigma},~{}\sigma=1,2,3}\left\{h_{\sigma,j}\right\}.$

Lemma 4.3**.**

If $\displaystyle(\ref{cfl})$ holds and $\displaystyle h^{2}<\frac{1}{2\min\{a^{2},b^{2},c^{2}\}}\min\left\{\frac{1}{f_{0}},\frac{4}{f(\tau_{0}f_{0}/\phi_{\min})}\right\},$ then for given $\displaystyle v_{0}\equiv 0$, we have that all components of $\displaystyle v_{1}<1.$

Proof.

If $\displaystyle v_{0}\equiv 0,$ then from $\displaystyle(\ref{c3})$ we have

[TABLE]

Using

[TABLE]

where $\displaystyle\xi=\left(\phi_{1,1,1}^{-1},\dots,\phi_{N_{1},N_{2},N_{3}}^{-1}\right)^{\intercal}\in{\mathbb{R}}^{N_{1}N_{2}N_{3}},$ we have following decomposed connections

[TABLE]

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

[TABLE]

for which

[TABLE]

The above indicates that

[TABLE]

On the other hand, according to Lemma 4.3 we have

[TABLE]

and thus,

[TABLE]

Since we wish each component of $\displaystyle s_{1}^{+}+s_{1}^{-}$ to be negative, we require

[TABLE]

Now, recall $\displaystyle(\ref{37b})$. It follows that

[TABLE]

Note that

[TABLE]

which implies that $\displaystyle s_{2}^{+}<\frac{1}{2}x.$ Therefore we arrive at

[TABLE]

By the same token, based on $\displaystyle(\ref{37c})$ we observe that

[TABLE]

It can be seen that

[TABLE]

and the above indicates that

[TABLE]

By Lemma 4.2 we conclude that $\displaystyle s_{3}^{-}\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 $\displaystyle(\ref{v1a})$ and the above.

Remark 4.1*.*

We could generalize the above lemma to include nonzero initial vectors, if desired. Let $\displaystyle 0<v_{0}<1/8$ be given. If $\displaystyle(\ref{cfl})$ holds and $\displaystyle h^{2}<\left(1-8\max\{v_{0}\}\right)/\left(2F\min\{a^{2},b^{2},c^{2}\}\right),$ where $\displaystyle F=f\left(\|v_{0}\|_{2}+\tau_{0}(C\|v_{0}\|_{2}\right.$

$\displaystyle\left.+f(\|v_{0}\|_{2}\phi_{\min}^{-1}))\right),$ then all components of $\displaystyle v_{1}$ generated by $\displaystyle(\ref{c3})$ are bounded above by unity. This follows by modifying the proof of Lemma 4.4.

Combining above results we obtain the following theorem.

Theorem 4.1**.**

For any beginning step $\displaystyle\ell_{0}\geq 0$ if $\displaystyle\tau_{\ell}$ is sufficiently small for $\displaystyle\ell\geq\ell_{0}$ and

(i)

$\displaystyle(\ref{cfl})$* holds for all $\displaystyle\ell\geq\ell_{0},$*

(ii)

$\displaystyle h^{2}<\frac{1}{2\min\{a^{2},b^{2},c^{2}\}}\min\left\{\frac{1}{f_{0}},\frac{4}{f(\tau_{0}f_{0}/\phi_{\min})}\right\},$* where $\displaystyle h=\max_{j=1,\dots,N_{\sigma},~{}\sigma=1,2,3}\{h_{\sigma,j}\},$*

(iii)

$\displaystyle Cv_{\ell_{0}}+\frac{1}{2}g(v_{\ell_{0}})>0$* and $\displaystyle\left(I-\frac{\tau_{\ell_{0}}}{2}g_{v}(\xi_{\ell_{0}})\right)^{-1}>0,$*

then the sequence $\displaystyle\left\{v_{\ell}\right\}_{\ell\geq{\ell_{0}}}$ produced by the semi-adaptive LOD 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.

5 Stability

Nonlinear stability has been an extremely difficult issue when nonlinear Kawarada equations are concerned [2, 4, 5, 18, 19, 21]. However, when the numerical solution varies relatively slowly, that is, before reaching a certain neighborhood of quenching, instability may be detected through a linear stability analysis of the nonlinear scheme utilized [7, 13, 22]. Although the application of such an analysis to nonlinear problems cannot be rigorously justified, it has been found to be remarkably informative in practical computations. In the following study, we will first carry out a linearized stability analysis in the von Neumann sense for $\displaystyle(\ref{c3})$ with its nonlinear source term frozen. This is equivalent to assuming that the source term is effectively accurate. The analysis will then be extended to circumstances where the nonlinear term is not frozen. In the later case, the boundedness of the Jacobian of the source term, $\displaystyle\|g_{v}(v)\|_{2},$ which is equivalent to assuming that we are some neighborhood away from quenching, is assumed.

In the following, let $\displaystyle A\in\mathbb{C}^{n\times n}$ and again denote $\displaystyle E(\cdot)=\exp(\cdot)$ for $\displaystyle n>1.$

Definition 5.1**.**

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]

where $\displaystyle I_{n}\in\mathbb{C}^{n\times n}$ is the identity matrix.

Remark 5.1*.*

If the matrix norm being considered is the spectral norm, then $\displaystyle\mu(A)=\max\left\{\lambda:\lambda\ \text{is an eigenvalue of}\ (A+A^{*})/2\right\}=\frac{1}{2}\lambda_{\max}(A+A^{*}).$

Lemma 5.1**.**

For $\displaystyle\alpha\in\mathbb{C}$ we have

[TABLE]

Proof.

See [11].

For the semi-adaptive LOD method $\displaystyle(\ref{c3})$ with its nonlinear source term frozen, regularity conditions need to be imposed upon the nonuniform spatial grids for a linear stability analysis. For this purpose, let us denote $\displaystyle h_{\sigma}=\min_{j=1,\dots,N_{\sigma}}\{h_{\sigma,j}\},~{}\sigma=1,2,3.$

Lemma 5.2**.**

If

[TABLE]

where the constant $\displaystyle K>0$ is independent of $\displaystyle h_{\sigma,j},~{}j=1,\dots,N_{\sigma},~{}\sigma=1,2,3.$ then

[TABLE]

Proof.

We only need to consider the case involving $\displaystyle M_{1}$ since the other cases are similar. Note that $\displaystyle\mu(M_{1})=\frac{1}{2}\lambda_{\max}\left(M_{1}+M_{1}^{\intercal}\right)$ and

[TABLE]

where

[TABLE]

We apply Gers̆chgorin’s circle theorem to an arbitrary $\displaystyle X_{j,k}$ and note that a similar argument works for each $\displaystyle X_{j,k},~{}j=1,\dots,N_{2},~{}k=1,\dots,N_{3}.$ Further, notice that we only need to consider circumstances where the bandwidth of $\displaystyle M_{1}+M_{1}^{\intercal}$ is three. Thus,

[TABLE]

We then see that $\displaystyle(\ref{sss1})$ follows immediately from the above and the fact that

[TABLE]

Lemma 5.3**.**

If $\displaystyle(\ref{sss1})$-$\displaystyle(\ref{sss3})$ hold then

[TABLE]

for sufficiently small $\displaystyle\tau_{\ell}>0.$

Proof.

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

[TABLE]

Now, based on Lemmas 5.1 and 5.2,

[TABLE]

which is the desired bound.

Combining the above results gives the following theorem.

Theorem 5.1**.**

If $\displaystyle(\ref{cfl})$ and $\displaystyle(\ref{sss1})$-$\displaystyle(\ref{sss3})$ hold, then the semi-adaptive LOD 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_{\ell}.$

Proof.

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

[TABLE]

Recall that $\displaystyle\sum_{k=0}^{\ell}\tau_{k}\leq T,~{}\ell>0.$ It follows by taking the norm on both sides of $\displaystyle(\ref{st1})$ that

[TABLE]

where $\displaystyle c_{1},c_{2},c_{3},c_{4}$ and $\displaystyle c$ are positive constants independent of $\displaystyle\ell,~{}\tau_{k},~{}0\leq k\leq\ell.$ Therefore the theorem is clear.

We now consider the case without freezing the nonlinear source term in $\displaystyle(\ref{c3})$. In this situation, restrictions upon the Jacobian matrix $\displaystyle g_{v}(v)$ become necessary.

Theorem 5.2**.**

Let $\displaystyle\tau_{k},~{}0\leq k\leq\ell,$ be sufficiently small and $\displaystyle(\ref{cfl})$, $\displaystyle(\ref{sss1})$-$\displaystyle(\ref{sss3})$ hold. If there exists a constant $\displaystyle G<\infty$ such that

[TABLE]

then the semi-adaptive LOD method $\displaystyle(\ref{c3})$ is unconditionally stable in the von Neumann sense, 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\tilde{c}>0$ is a constant independent of $\displaystyle\ell$ and $\displaystyle\tau_{\ell}.$

Proof.

By definition we have

[TABLE]

where

[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, recall $\displaystyle(\ref{cond1})$. When $\displaystyle\tau_{k}$ is sufficiently small we may claim that

[TABLE]

Thus,

[TABLE]

It follows therefore

[TABLE]

where $\displaystyle c_{1,k},~{}k=1,2,\ldots,\ell,$ are positive constants and $\displaystyle c,c_{1},\tilde{c}$ are positive constants independent of $\displaystyle\ell$ and $\displaystyle\tau_{\ell},~{}\ell>0.$ Thus giving the desired stability.

The above theorem provides further insight as to why the standard linear analysis can be useful in estimating the nonlinear stability. The extra cost paid, however, is assuming the boundedness of $\displaystyle\|g_{v}(\xi)\|_{2}.$ Nevertheless, this is an improvement upon the traditional methodology of having the nonlinear source term frozen. In fact, the aforementioned bound is well-maintained in numerical experiments until certain neighborhoods of quenching are reached. This serves as an indication that the new analysis is valid and effective.

6 Conclusions

A semi-adaptive LOD scheme is developed for solving degenerate Kawarada equations possessing a strong quenching nonlinearity and singularity. While a temporal adaptation is performed via an arc-length monitoring mechanism of the temporal derivative of the solution, fixed nonuniform spatial grids are adopted. The novel splitting method is implicit and the impact of the degeneracy is found to be limited. Rigorous analysis is given for key computational features, including the positivity, monotonicity, and stability, of the numerical solution. Important criteria to guarantee these properties, which depend upon the variable steps and degenerate function, are obtained.

Under much weaker requirements (see the latest results in [4]), the temporal step restriction for guaranteeing monotone numerical solutions of our LOD scheme has been reduced to only one-half of those in uniform spatial mesh cases [18]. Furthermore, a realistic method of targeting the realization of nonlinear stability analysis is proposed and shown to be successful. Though this new strategy needs the boundedness of $\displaystyle\|g_{v}(\xi)\|_{2},$ the requirement is well-justified before quenching is reached. This improved methodology not only provides further insight into the stability, but also offers explanations as to why the linear stability analysis must be valid before quenching. On the other hand, simulations of real three-dimensional solutions still remain as one of the most challenging tasks. In anticipated future work we plan to utilize the latest High Performance Computing tools with large data computations for this purpose. More rigorous and generalized analysis, as well as non-exponential splitting based higher order splitting methods [17, 19] will also be be investigated, studied, and experimented with.

References