Positive Fitted Finite Volume Method for Semilinear Parabolic Systems on Unbounded Domain

Abstract

This work deals with a semilinear system of parabolic partial differential equations (PDEs) on an unbounded domain, related to environmental pollution modeling. Although we study a one-dimensional sub-model of a vertical advection–diffusion, the results can be extended in each direction for any number of spatial dimensions and different boundary conditions. The transformation of the independent variable is applied to convert the nonlinear problem into a finite interval, which can be selected in advance. We investigate the positivity of the solution of the new, degenerated parabolic system with a non-standard nonlinear right-hand side. Then, we design a fitted finite volume difference discretization in space and prove the non-negativity of the solution. The full discretization is obtained by implicit–explicit time stepping, taking into account the sign of the coefficients in the nonlinear term so as to preserve the non-negativity of the numerical solution and to avoid the iteration process. The method is realized on adaptive graded spatial meshes to attain second-order of accuracy in space. Some results from computations are presented.

1. Introduction

Environmental problems are becoming increasingly important for our world, and their significance is expected to grow in the future. High levels of pollution in the air, water, and soil can cause damage to plants, animals, and humans.

Many models in environmental modeling are generally described by systems of parabolic PDEs for calculating the concentrations of a number of chemical species (pollutants and components of the air that interact with the pollutant) in a large 3D domain (part of the atmosphere above the studied geographical region). For example, the Danish Euler Model (DEM) is one of the most frequently used atmosphere pollution model and is mathematically represented by the following system PDEs [1,2,3,4]:

$\begin{array}{ccc}\hfill \frac{\partial c_s}{\partial t}& =& -\frac{\partial \left(u{c}_s\right)}{\partial x}-\frac{\partial \left(v{c}_s\right)}{\partial y}-\frac{\partial \left(w{c}_s\right)}{\partial z}+\frac{\partial }{\partial x}\left(K_s^x\frac{\partial c_s}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_s^y\frac{\partial c_s}{\partial y}\right)\hfill \\ & & +\frac{\partial }{\partial z}\left(K_s^z\frac{\partial c_s}{\partial z}\right)+f_s+r_s(c_1,c_2,\dots ,c_S)-(k_{1s}+k_{2s})c_s,s=1,2,\dots ,S,\hfill \end{array}$

where $c_s$ is the concentration of the chemical species; u, v, and w are wind velocities and $K_s^x$, $K_s^y$, and $K_s^z$ are the diffusion components; $f_s$ is the emissions; $k_{1s}$ and $k_{2s}$ are dry/wet deposition coefficients; and $r_s(c_1,c_2,\dots ,c_S)$ are nonlinear functions describing the chemical reactions between the species under consideration [1,2,3,4,5]. A typical case is as follows:

(1)$r_s(c_1,c_2,\cdots ,c_s)=\sum _{i=1}^S{\gamma }_{s,i}{c}_i+\sum _{i=1}^S\sum _{j=1}^S{\beta }_{s,i,j}{c}_i{c}_j,s=1,2,\cdots ,S,$

where ${\gamma }_{s,i}<0$ and ${\beta }_{s,i,j}<0$ are constants.

For such complex models, operator splitting [1,4,6,7,8,9] is very often applied to achieve sufficient accuracy and efficiency of the numerical solution, but it does not guarantee the non-negativity of the physical solution [10]. Although the splitting is a crucial step in the efficient numerical treatment of the model, after the discretization of the large computational domain, each sub-problem becomes a huge computational task itself.

This paper focuses on designing a numerical method for a class of reaction–convection–diffusion systems related to environmental pollution modeling. We present theoretical and computational results for a nonlinear system of parabolic PDEs, which is a vertical sub-model of air pollution transport problems.

The first part consists of the transformation of the original problem to a system of parabolic PDEs defined on a bounded domain. Subsequently, the finite volume method, initially proposed to address the degeneracy of the transformed on a finite interval Black–Scholes equation [11], is applied to the obtained parabolic system with boundary degeneration.

Here, we will concentrate on a non-stationary sub-model of a vertical advection–diffusion with chemistry, emissions, and deposition (see [4,12])

(2)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{\partial c_s}{\partial t}-\frac{\partial }{\partial z}\left(k_s\left(z\right)\frac{\partial c_s}{\partial z}\right)+w\frac{\partial c_s}{\partial z}-r_s(c_1,c_2,\dots ,c_S)=f_s(t,z),z\in (0,\infty ),t\in (0,T],\hfill \\ & \frac{\partial c_s}{\partial z}(t,0)={\nu }_s{c}_s(t,0),{\nu }_s=const\ge 0,t\in [0,T],\hfill \\ & \underset{z\to \infty }{lim}{c}_s(t,z)=0,t\in [0,T],\hfill \\ & c_s(0,z)=c_{s,0},z\in [0,\infty ),s=1,2,\dots ,S,\hfill \end{array}\end{array}$

with a more general representation of the nonlinear right-hand-side $r_s(c_1,c_2,\cdots ,c_S)$ in (1), including sign-changing coefficients. There exists many theoretical and numerical methods for parabolic PDEs (see, e.g., [13,14,15]) and weakly coupled parabolic and elliptic systems (see, e.g., [16]). However, the specific nonlinear right-hand sides of systems (1) and (2) require special considerations for providing the important positivity-preserving property of the solution [10,17,18].

Since, in general, numerical methods are applied to a finite domain, one must consider an artificial truncation of the domain—in our case, the semi-infinite intervals. A possible solution to this issue is the construction of an artificial boundary condition to confine the computational domain appropriately.

If the solution on the computational domain coincides with the exact solution on the unbounded domain (restricted to the finite domain), one refers to these boundary conditions as a transparent boundary condition (TBC). A discretization strategy of TBC to solve a linear air pollution model equation is presented in [19].

Various numerical methods for solving linear advection–diffusion single equations, which describe the processes of pollutant transport and diffusion in the atmosphere, are proposed in the literature (see, e.g., [2,9,12,19,20,21,22,23,24,25,26,27]).

Considerably fewer numerical results exist for systems of nonlinear parabolic advection–diffusion equations modeling air pollution processes. The book [4] selects efficient numerical methods and discusses their implementation for solving such problems.

In the paper [1], the authors apply sequential splitting for numerically solving an advection–diffusion system with a nonlinear chemistry module. They derive two sub-problems: linear and nonlinear systems of ordinary differential equations (ODEs). Four different numerical algorithms, including one based on a preconditioned sparse matrix approach, are developed and compared for efficiently solving these ODE systems.

A new meshless method is proposed in [28] for solving a system of linear and nonlinear advection–diffusion-reaction equations describing transport in anisotropic media. The method is based on representing the analytical solution as a series over a basis system that satisfies boundary conditions for any free parameter. Then, radial basis functions are used, and the free parameter is determined by collocation inside the solution domain.

Fourth-order compact finite difference schemes and fourth-order schemes based on Richardson extrapolation are constructed in [29] for solving a system of ten parabolic partial differential equations describing the dispersion of pollutants in the air.

The finite volume method is often used for numerically solving linear and nonlinear PDEs and PDE systems with applications in many fields, including physics, biology, finance, etc., since it is conservative and can be easily applied to unstructured meshes. A fundamental theoretical review paper is [30], where new theoretical and computational results for linear convection–diffusion problems are reported. The fitted finite volume method is first proposed by S. Wang [31,32] for handling the degeneration at $x=0$ in the Black–Scholes equation.

In [33], the problem (2) is transformed on a finite interval and then a fitted finite volume is applied. In the present work, we extend these results, providing theoretical investigation, improving the order of convergence of the numerical approach and proposing efficient realization based on implicit–explicit time stepping and decoupling of the discrete equations to avoid iteration processes and preserve the non-negativity of the solution.

Although our method works for any number of spatial dimensions, in this paper, we consider a one-dimensional problem to simplify the computations.

The rest of the paper is organized as follows. In Section 2, we formulate and investigate the transformed differential problem and discuss its well-posedness and the properties of the solution. In Section 3, we develop a fitted finite volume semidiscretization and study the positivity-preserving property. In Section 4, we construct a full discretization and prove that the numerical solution is non-negative. Results from numerical tests are presented and discussed in Section 5. Finally, we provide some conclusions.

2. The Differential Problem on Bounded Domain

In the numerical simulations, the use of boundary condition at infinity is not suitable. In [19], discrete transparent boundary conditions are developed and examined for the basic scenario of a single linear advection–diffusion equation.

In this work, we use the transformation [33]

$z=\frac{1}{2a}log\left(\frac{1+\xi }{1-\xi }\right),\xi \in \Omega =(0,1)\iff \xi =\frac{e^{2az}-1}{e^{2az}+1},z\in (0,\infty ),$

where a is a stretching factor. Applying this transformation, the system (2) becomes

(3)$$\begin{gathered} L_s{C}_s\equiv \frac{\partial C_s}{\partial t}-a^2{\left(1-{\xi }^2\right)}^2{k}_s\left(\xi \right)\frac{{\partial }^2{C}_s}{\partial {\xi }^2} \\ +a\left(1-{\xi }^2\right)\left(2a\xi k_s\left(\xi \right)+w-a\left(1-{\xi }^2\right)\frac{\partial k_s\left(\xi \right)}{\partial \xi }\right)\frac{\partial C_s}{\partial \xi } \\ =r_s(C_1,C_2,\dots ,C_S)+f_s(t,\xi ),\xi \in \Omega ,t\in (0,T], \end{gathered}$$

(4)$\begin{array}{ccc}& & a\frac{\partial C_s}{\partial \xi }(t,0)={\nu }_s{C}_s(t,0),t\in [0,T],\hfill \end{array}$

(5)$\begin{array}{ccc}& & C_s(t,1)=0,t\in [0,T],\hfill \end{array}$

(6)$\begin{array}{ccc}& & C_s(0,\xi )=C_s^0\left(\xi \right),\xi \in [0,1],s=1,2,\dots ,S,\hfill \end{array}$

where

$\begin{array}{ccc}& & C_s(t,\xi )\equiv c_s(t,z\left(\xi \right)),k_s\left(\xi \right)=K_s\left(z\left(\xi \right)\right),\hfill \\ & & r_s(C_1,C_2,\cdots ,C_s)=r_s(c_1(t,z\left(\xi \right)),c_2(t,z\left(\xi \right)),\cdots ,c_S(t,z\left(\xi \right)).\hfill \end{array}$

It can be easily seen that at $\xi =1$, the system (3) degenerates to the ODE system:

(7)$\frac{\partial C_s(t,1)}{\partial t}-r_s(C_1(t,1),\cdots ,C_S(t,1)=f_s(t,1),s=1,2,\cdots ,S,t\in (0,T].$

In the case of $C_s^0=0$ and homogeneous boundary conditions (4) and (5), the unique solution of (7) is zero.

According to the Fichera and the Oleinik–Radkevich theory [34] for degenerate parabolic equations, at the degenerate boundary $\xi =1$, the boundary condition should not be given. However, based on physical motivation, we have imposed the boundary condition (5). It is easy to check that if the functions $C_s(t,\xi )$ satisfy (5), then they also satisfy (7). Therefore, (5) is a particular case of (7).

One of the advantages of transformation is the small-sized domain (here, the interval (0, 1)), but one may also introduce an additional parameter controlling the size, which implies a lower amount of computational effort needed to obtain sufficient accuracy.

Following the physical motivation, we will discuss the non-negativity of the solution.

The idea is to reformulate the parabolic system (3)–(7) to the initial-value problem (IVP) for ODE systems and to apply the following well-known result (see, e.g., [35]).

Let us introduce the IVP

(8)$\frac{\partial W}{\partial t}=\mathcal{F}(t,W),t\ge 0,W\left(0\right)=W_0,W_0\in {\mathbb{R}}^{\mu },$

where $\mathcal{F}:R^{+}\times R^{\mu }\to R^{\mu }$ is continuous and (8) has a unique solution for all $W_0$.

In the following, we will write $v\ge 0$ for a vector $v\in R^{\mu }$ if all the components are non-negative.

We consider the ODE system (8) as positive (short for “non-negative preserving”), if $W\left(t\right)\ge 0$ holds for all $t\ge 0$ whenever $W_0\ge 0$. An often used criterion is given by the following lemma:

Lemma 1

(Theorem 7.1 in [35]). The IVP (8) is positive if and only if

(9)$W_i=0,W_j\ge 0forallj\ne i⟹{\mathcal{F}}_i(\tau ,W)\ge 0$

holds for all τ and any vector $W\in R^{\mu }$ and $i=1,2\cdots ,\mu$. Then, the system (9) is positive.

The next assertion provides new sufficient conditions for the positivity of the initial-value problem (3)–(7).

Theorem 1.

Let $\mathbf{C}\equiv (C_1,C_2,\cdots ,C_S)$ be a solution of this problem (3)–(6) (or ((3), (4), (6) and (7)). Assume that for all $s=1,\cdots ,S$, the function $r_s(C_1,\cdots ,C_S)$ is Lipshitz-continuous with respect to concentration $C_1,\cdots ,C_s$ and it satisfies

(10) $r_s(C_1,\cdots ,C_{s-1},0,C_{s+1},\cdots ,C_S)\ge 0\forall \mathbf{C}\in R_{+}^S\equiv \{C_s\ge 0\}.$

Let $f_s(t,\xi )\ge 0,s=1,\cdots ,S$. Finally, assume for the solution $\mathbf{C}$ of the problem (3)–(6) (or ((3), (4), (6) and (7)) that $C_s\in C\left({\overline{Q}}_T\right)\cap C^{1.3}\left(Q_T\right)$ and for the initial concentration

$C_s(0,\xi )=C_s^0\left(\xi \right)\ge 0,s=1,\cdots ,S.$

Then,

$C_s(t,\xi )\ge 0,s=1,\cdots ,S.$

Proof. 

For ease of exposition, we will present the proof for a model of three regime concentrations, i.e., $S=3$ in (3).

$\begin{array}{ccc}& & L_1{C}_1=r_1(C_1,C_2,C_3)+f_1(t,\xi ),\hfill \\ & & L_2{C}_2=r_2(C_1,C_2,C_3)+f_2(t,\xi ),\hfill \\ & & L_3{C}_3=r_3(C_1,C_2,C_3)+f_3(t,\xi ).\hfill \end{array}$

As usual [16], we use the linearization procedure for $r_1(C_1,C_2,C_3)$ (and similarly for $r_2(C_1,C_2,C_3)$ and $r_3(C_1,C_2,C_3)$) to obtain

$\begin{array}{ccc}\hfill L_1{C}_1& =& r_1(C_1,C_2,C_3)-r_1(0,C_2,C_3)+r_1(0,C_2,C_3)+f_1(t,\xi )\hfill \\ & \ge & r_1(C_1,C_2,C_3)-r_1(0,C_2,C_3).\hfill \end{array}$

We have used condition (10) and assumption $f_1(t,\xi )\ge 0$. In addition, the Lipschitz continuity of $r_i$ guarantees that there exists ${\widehat{C}}_1$ with $0\le {\widehat{C}}_1\le C_1$, such that

$L_1{C}_1\ge \left(\frac{\partial r_1}{\partial C_1}{\left(\mathbf{C}\right)|}_{C_1={\widehat{C}}_1}\right)C_1.$

Moreover, the Lipschitz continuity of $r_1$ guarantees that for $t\in [0,T]$, there exists a positive constant M that

$\left|\frac{\partial r_1}{\partial C_1}{\left(\mathbf{C}\right)|}_{C_1={\widehat{C}}_1}\right|\le M^{*}.$

Therefore,

(11)$L_1{C}_1\ge -M^{*}{C}_{1}fort\in [0,T].$

Let ${\mathbf{C}}^{\epsilon }\equiv (C_1^{\epsilon },C_2^{\epsilon },C_3^{\epsilon })$ be the solution of (3)–(6) (or ((3), (4), (6) and (7)) for $S=3$ with the initial condition

$C_s^{\epsilon }(0,\xi )=C_{\epsilon }^0\left(\xi \right)+\epsilon ,s=1,2,3$

for $\epsilon >0$. It is easy to see that, if $\mathbf{C}$ is the solution on $[0,T]$, then ${\mathbf{C}}^{\epsilon }$ is a solution on $[0,T_{\epsilon }]$, with $T_{\epsilon }$ as $\epsilon \to 0$. For each fixed $\epsilon >0$, since initially $C_s^{\epsilon }>0$, we have $C_s^{\epsilon }(t,\xi )>0$ for a sufficiently small t. Furthermore, we show that this is also true for all $[0,T_{\epsilon }]$.

Suppose otherwise, for some s, there exists $t_0<T_{\epsilon }$, such that

(12)$C_s^{\epsilon }(t_0,{\xi }_0)=0\mathrm{and}{\mathbf{C}}^{\epsilon }\in {\mathbb{R}}_{+}^{\epsilon }\mathrm{for}0\le t\le t_0.$

Without loss of generality, we take $s=1$, and in view of (11), we have

$\frac{\partial C_1^{\epsilon }}{\partial t}\ge L_1{C}_1^{\epsilon }-M^{*}{C}_1^{\epsilon }.$

Let $G=e^{M^{*}t}{C}_s^{\epsilon }$. Then, it follows from the last inequality that

(13)$\begin{array}{c}\hfill \begin{array}{ccc}\hfill & & \frac{\partial G}{\partial t}\ge L_{1}G,\\ \\ & & G(0,\xi )=C_1^0\left(\xi \right)+\epsilon \ge \epsilon .\end{array}\end{array}$

Now, for the linear parabolic problem (13), we can apply the classical maximum principle [17,36] to arrive at

$G(t,\xi )\ge \epsilon .$

Consequently,

$C_1^{\epsilon }\ge \epsilon e^{-M^{*}t}>0\forall t\in [0,T_{\epsilon }].$

This is a contradiction with (12). Therefore, we must have

$C_1^{\epsilon }(t,\xi )>0\forall t\in [0,T_{\epsilon }].$

Finally, letting $\epsilon \to 0$, we obtain

$C_i(t,\xi )\ge 0,\forall t\in [0,T].$

Let us now suppose that $G(t,\xi )$ attains negative minimum at a point $P(t,0),0\le t\le t_0$. Then, we have

$\frac{\partial G}{\partial \xi }(t,0)=\underset{\xi \to 0+}{lim}\frac{G(t,\xi )-G(t,0)}{\xi }\ge 0,$

which contradicts the boundary condition (4).

Finally, if we work with boundary condition (5), the proof is completed. Let us take (7). Then, again using linearization, for $s=1$, we derive

$\frac{\partial C_1(t,1)}{\partial t}\ge -M^{*}{C}_1(t,1),C_1(0,1)=C_1^0.$

Solving this differential inequality, we find

$C_1(t,1)=C_1^0\left(1\right)e^{-M^{*}t}\ge 0.$

□

Example 1.

We present, as an example of system (2), a simplified but somewhat realistic example of a three-component system that models the chemical processes occurring in the atmosphere model, based on Chapman’s cycle [3,9]. The components of the system are nitric oxide (NO), nitrogen dioxide $\left(N{O}_2\right)$, and ozone $\left(O_3\right)$ denoted by $C_1$, $C_2$, and $C_3$, respectively. A simplified model of chemical reactions in the system is

$NO+O_3\stackrel{v_1}{\to }N{O}_2,N{O}_2+O_3\stackrel{v_2}{\to }NO+O_3,$

with constant rates $v_1>0,v_2>0$. Then, reaction terms are given by

$\begin{array}{ccc}& & r_1\left(\mathbf{C}\right)=r_1(C_1,C_2,C_3)=v_2{C}_2-v_1{C}_1{C}_3,\hfill \\ & & r_2\left(\mathbf{C}\right)=-v_2{C}_2+v_1{C}_1{C}_3,\hfill \\ & & r_3\left(\mathbf{C}\right)=v_2{C}_2-v_1{C}_1{C}_3.\hfill \end{array}$

Now, it is easy to check that the condition (10) is fulfilled and Theorem 1 holds.

3. Fitted Finite Volume Difference Scheme

In the following section, we perform the spatial semidiscretization by the method of vertical lines and investigate the key property of the positivity for the discrete solution [10,17,36].

3.1. Spatial Discretization

This subsection presents the spatial discretization of (3)–(7). Further, for simplicity, we assume that $k_s\left(z\right)>0$ and $k_s\left(\xi \right)>0$. This condition is satisfied by many space-dependent functions, for example, $k_s\left(z\right)=e^{pz}+q$, where $p\le 0$ and $q>0$ are constants, including constant $k_s>0$ ($p=0$) (see, e.g., [4,12]).

We write the system (3) in divergent form:

(14)$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial C_s}{\partial t}=& \frac{\partial }{\partial \xi }\left(p_s\left(\xi \right)\frac{\partial C_s}{\partial \xi }+q_s\left(\xi \right)C_s\right)+B_s(\xi ,C_1,\dots ,C_S)\hfill \\ & +f_s(t,\xi ),(t,\xi )\in Q_T=[0,T]\times \Omega ,s=1,2,\dots ,S,\hfill \end{array}\end{array}$

where

$\begin{array}{ccc}& & p_s\left(\xi \right)=a^2{(1-{\xi }^2)}^2{k}_s\left(\xi \right),q_s\left(\xi \right)=a(1-{\xi }^2)(2a\xi k_s\left(\xi \right)-w),\hfill \\ & & B_s(\xi ,C_1,\dots ,C_S)=r_s(C_1,\cdots ,C_S)-d_s\left(\xi \right)C_s,\hfill \\ & & d_s\left(\xi \right)=2a^2(1-3{\xi }^2)k_s\left(\xi \right)+2a^2\xi (1-{\xi }^2)\frac{\partial k_s\left(\xi \right)}{\partial \xi }+2aw\xi ,\hfill \end{array}$

and $f_s(t,\xi )$, $s=1,2,\dots ,S$ is a regularization of the Dirac delta function.

We construct fitted finite volume semidiscretization [31,32] for the problem (14).

Let the interval $[0,1]$ be subdivided into N intervals $I_i=[{\xi }_i,{\xi }_{i+1}]$, $i=0,1,\dots ,N-1$ with $0={\xi }_0<{\xi }_1<\dots <{\xi }_N-1<{\xi }_N=1$ and $h_i={\xi }_{i+1}-{\xi }_i$.

We consider a dual mesh $0={\xi }_{-1/2}={\xi }_0<{\xi }_{1/2}<{\xi }_1<x_{3/2}<\cdots <{\xi }_{N-1/2}<{\xi }_N={\xi }_{M+1/2}=1$ and denote ${\hslash }_i={\xi }_{i+\frac{1}{2}}-{\xi }_{i-\frac{1}{2}}$ for $i=0,1,\dots ,N$.

A. Internal nodes $1\le i\le N-1$. We integrate Equation (14) on the cell $[{\xi }_{i-\frac{1}{2}},{\xi }_{i+\frac{1}{2}}]$:

(15)$\begin{array}{c}\underset{{\xi }_{i-\frac{1}{2}}}{\overset{{\xi }_{i+\frac{1}{2}}}{\int }}\frac{\partial C_s}{\partial t}d\xi =\underset{{\xi }_{i-\frac{1}{2}}}{\overset{{\xi }_{i+\frac{1}{2}}}{\int }}\frac{\partial }{\partial \xi }\left(p_s\left(\xi \right)\frac{\partial C_s}{\partial \xi }+q_s\left(\xi \right)C_s\right)d\xi \hfill \\ \hfill +\underset{{\xi }_{i-\frac{1}{2}}}{\overset{{\xi }_{i+\frac{1}{2}}}{\int }}\left[B_s(\xi ,C_1,\dots ,C_S)+f_s(t,\xi )\right]d\xi ,i=1,2,\dots ,N-1.\end{array}$

Applying the mid-point quadrature rule to all the integrals in (15) with the exception of the second one, we obtain

(16)$\begin{array}{c}\hfill \begin{array}{ccc}\hfill & & {\hslash }_i{\left.\frac{\partial C_s}{\partial t}\right|}_{(t,{\xi }_i)}={\left.\left(p_s\left(\xi \right)\frac{\partial C_s}{\partial \xi }+q_s\left(\xi \right)C_s\right)\right|}_{\left(t,{\xi }_{i+\frac{1}{2}}\right)}-{\left.\left(p_s\left(\xi \right)\frac{\partial C_s}{\partial \xi }+q_s\left(\xi \right)C_s\right)\right|}_{\left(t,{\xi }_{i-\frac{1}{2}}\right)}\\ & & +{\hslash }_i{\left.\left[B_s(\xi ,C_1,\dots ,C_S)+f_s(t,\xi )\right]\right|}_{(t,{\xi }_i)}.\end{array}\end{array}$

Further, to derive the semidiscrete equations, we follow the methodology of [32]. Introducing the notation $v_{s,i}=v_s({\xi }_i,t)$, we rewrite Equation (16) in the form

(17)$\frac{\partial C_{s,i}}{\partial t}{\hslash }_i=(1-{\xi }_{i+\frac{1}{2}}^2){\rho }_{s,i+\frac{1}{2}}-(1-{\xi }_{i-\frac{1}{2}}^2){\rho }_{s,i-\frac{1}{2}}+{\hslash }_i\left(B_{s,i}+f_{s,i}\right),i=1,2,\cdots ,N-1,$

where

(18)${\rho }_s\equiv {\rho }_s\left(C_s\right)=a\left[a(1-{\xi }^2)k_s\left(\xi \right)\frac{\partial C_s}{\partial \xi }+(2a\xi k_s\left(\xi \right)-w)C_s\right]$

and

$B_{s,i}=B({\xi }_i,C_{1,i},C_{2,i},\cdots ,C_{S,i}).$

We need to derive an approximation of the continuous flux ${\rho }_s$ in the point ${\xi }_{i+1/2}$, $i=0,1,\dots ,N-1$. To perform this, we consider the two-point BVP:

(19)$\begin{array}{ccc}& & {\left((1-\xi )l_{s,i+\frac{1}{2}}{V_s}^{\prime }+m_{s,i+\frac{1}{2}}{V}_s\right)}^{\prime }=0,\xi \in I_i,\hfill \end{array}$

(20)$\begin{array}{ccc}& & V_s\left({\xi }_i\right)=C_{s,i},V_s\left({\xi }_{i+1}\right)=C_{s,i+1},\hfill \end{array}$

where $l_s=a^2{k}_s\left(\xi \right)(1+\xi )$, $m_s=a(2a\xi k_s\left(\xi \right)-w)$, $l_{s,i+\frac{1}{2}}=l_s({\xi }_{i+\frac{1}{2}})$, $m_{s,i+\frac{1}{2}}=m_s({\xi }_{i+\frac{1}{2}})$. Integrating (19) yields the first-order linear equation. Solving this equation and using the boundary conditions (20), we obtain

(21)${\rho }_{s,i}=m_{s,i+\frac{1}{2}}\frac{{\left(1-{\xi }_{i+1}\right)}^{{\alpha }_{s,i}}{C}_{s,i}-{\left(1-{\xi }_i\right)}^{{\alpha }_{s,i}}{C}_{s,i+1}}{{\left(1-{\xi }_{i+1}\right)}^{{\alpha }_{s,i}}-{\left(1-{\xi }_i\right)}^{{\alpha }_{s,i}}},$

where ${\alpha }_{s,i}=\frac{m_{s,i+\frac{1}{2}}}{l_{s,i+\frac{1}{2}}}$.

B. Flux approximation close to the degenerate boundary $\xi =1$.

In order to find ${\rho }_{s,N}$, we consider its asymptotic behavior for $\xi \to 1^{-}$ [32]. If ${\alpha }_N<0$, from (21), for $i=N-1$, we derive

(22)$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{\xi \to 1^{-}}{lim}{\rho }_{s,N-1}=& m_{s,N-1/2}\underset{\xi \to 1^{-}}{lim}\frac{C_{s,N-1}-{\left(1-{\xi }_N\right)}^{-{\alpha }_{s,N-1}}{\left(1-{\xi }_{N-1}\right)}^{{\alpha }_{s,N-1}}{C}_{s,N}}{1-{\left(1-{\xi }_N\right)}^{-{\alpha }_{s,N-1}}{\left(1-{\xi }_N-1\right)}^{{\alpha }_{s,N-1}}}\hfill \\ \hfill =& m_{s,N-1/2}{C}_{s,N-1}.\hfill \end{array}\end{array}$

Analogically, if ${\alpha }_N>0$, from (21), we obtain

(23)$\underset{\xi \to 1^{-}}{lim}{\rho }_{s,N-1}=m_{s,N-1/2}{C}_{s,N}.$

From (22) and (23), we have

(24)${\rho }_{s,N-1}=m_{s,N-1/2}\left[\frac{1+\mathrm{sign}\left(m_{s,N-1/2}\right)}{2}{C}_{s,N}+\frac{1-\mathrm{sign}\left(m_{s,N-1/2}\right)}{2}{C}_{s,N-1}\right].$

C. Boundary $\xi =0$. At the vertical boundary $\xi =0$, we derive the approximation by integrating Equation (14) on the interval $[{\xi }_0,{\xi }_{\frac{1}{2}}]$. Thus, we have

$\frac{\partial C_{s,0}}{\partial t}\frac{h_1}{2}={(1-{\xi }_{\frac{1}{2}})}^2{\rho }_{s,\frac{1}{2}}-{\rho }_{s,0}+\frac{h_1}{2}\left(B_{s,0}+f_{s,0}\right),$

where ${\rho }_{s,\frac{1}{2}}$ is approximated from (21) for $i=0$. For ${\rho }_{s,0}$, from (18), taking into account the boundary condition (4), we derive

${\rho }_{s,0}=a(a{\nu }_s{k}_s\left({\xi }_0\right)-w)C_{s,0}.$

Combining all the above results, we obtain the spatial semidiscretization. It is a nonlinear ODE system of equations with unknown $C_{s,i}\left(t\right)$, $s=1,2,\dots ,S$

(25)$\begin{array}{ccc}\hfill \frac{\partial C_{s,i}}{\partial t}{\hslash }_i={\hslash }_i{\mathcal{L}}_i\left({\mathbf{C}}_s\right)& =& e_{s,i,i-1}{C}_{s,i-1}+e_{s,i,i}{C}_{s,i}+e_{s,i,i+1}{C}_{s,i+1}\hfill \\ & & +{\hslash }_i\left[B_{s,i}+f_{s,i}\left(t\right)\right],i=0,1,\dots ,N-1,\hfill \end{array}$

(26)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \left\{\begin{array}{cc}{C}_{s,N}=0,\hfill & \mathrm{for}\mathrm{b}.\mathrm{c}.\left(5\right),\hfill \\ \frac{\partial C_{s,N}}{\partial t}=r_{s,N}+f_{s,N},\hfill & \mathrm{for}\mathrm{b}.\mathrm{c}\left(7\right),\hfill \end{array}\right.\hfill \\ & C_{s,i}\left(0\right)=C_s^0\left({\xi }_i\right),i=0,1,\cdots ,N.\hfill \end{array}\end{array}$

where

$\begin{array}{ccc}& & e_{s,0,0}=-\frac{2m_{s,1/2}(1-{\xi }_{1/2}^2){(1-{\xi }_1)}^{{\alpha }_{s,0}}}{{\left(1-{\xi }_0\right)}^{{\alpha }_{s,0}}-{\left(1-{\xi }_1\right)}^{{\alpha }_{s,0}}}-2a(a{\nu }_s{k}_s\left({\xi }_0\right)-w),\hfill \\ & & e_{s,0,1}=\frac{2m_{s,1/2}(1-{\xi }_{1/2}^2){(1-{\xi }_0)}^{{\alpha }_{s,0}}}{{\left(1-{\xi }_0\right)}^{{\alpha }_{s,0}}-{\left(1-{\xi }_1\right)}^{{\alpha }_{s,0}}}\hfill \\ & & e_{s,0,-1}=0,\hfill \\ & & e_{s,i,i}=-\frac{m_{s,i+1/2}(1-{\xi }_{i+1/2}^2){(1-{\xi }_{i+1})}^{{\alpha }_{s,i}}}{{\left(1-{\xi }_i\right)}^{{\alpha }_{s,i}}-{\left(1-{\xi }_{i+1}\right)}^{{\alpha }_{s,i}}}-\frac{m_{s,i-1/2}(1-{\xi }_{i-1/2}^2){(1-{\xi }_{i-1})}^{{\alpha }_{s,i-1}}}{{\left(1-{\xi }_{i-1}\right)}^{{\alpha }_{s,i-1}}-{\left(1-{\xi }_i\right)}^{{\alpha }_{s,i-1}}}\hfill \\ & & e_{s,i,i+1}=\frac{m_{s,i+1/2}(1-{\xi }_{i+1/2}^2){(1-{\xi }_i)}^{{\alpha }_{s,i}}}{{\left(1-{\xi }_i\right)}^{{\alpha }_{s,i}}-{\left(1-{\xi }_{i+1}\right)}^{{\alpha }_{s,i}}}\hfill \\ & & e_{s,i,i-1}=\frac{m_{s,i-1/2}(1-{\xi }_{i-1/2}^2){(1-{\xi }_i)}^{{\alpha }_{s,i-1}}}{{\left(1-{\xi }_{i-1}\right)}^{{\alpha }_{s,i-1}}-{\left(1-{\xi }_i\right)}^{{\alpha }_{s,i-1}}},i=1,2,\cdots ,N-2,\hfill \\ & & e_{s,N-1,N-2}=m_{s,N-3/2}\frac{(1-{\xi }_{N-3/2}^2){(1-{\xi }_{N-1})}^{{\alpha }_{s,N-2}}}{{(1-{\xi }_{N-2})}^{{\alpha }_{s,N-2}}-{(1-{\xi }_{N-1})}^{{\alpha }_{s,N-2}}},\hfill \end{array}$

$\begin{array}{ccc}& & e_{s,N-1,N-1}=m_{s,N-1/2}(1-{\xi }_{N-1/2}^2)\frac{1-\mathrm{sign}\left(m_{s,N-1/2}\right)}{2}\hfill \\ & & -m_{s,N-3/2}\frac{(1-{\xi }_{N-3/2}^2){(1-{\xi }_{N-2})}^{{\alpha }_{s,N-2}}}{{(1-{\xi }_{N-2})}^{{\alpha }_{s,N-2}}-{(1-{\xi }_{N-1})}^{{\alpha }_{s,N-2}}},\hfill \\ & & e_{s,N-1,N}=m_{s,N-1/2}(1-{\xi }_{N-1/2}^2)\frac{1+\mathrm{sign}\left(m_{s,N-1/2}\right)}{2}.\hfill \end{array}$

3.2. Positivity-Preserving

Now, we will discuss the non-negativity of the semidiscrete solution.

Theorem 2.

Let the conditions of Theorem 1 be fulfilled, $a{\nu }_s{k}_s\left({\xi }_0\right)\ge w$ and $C_s^0\left(\xi \right)\ge 0$, $s=1,2,\cdots ,S$. Then, the solution of (25) and (26) is non-negative:

$C_{s,i}\left(t\right)\ge 0,t\in [0,T],s=1,2,\cdots ,S,i=0,1\cdots ,N.$

Proof. 

First, we establish that coefficients in the system (25) has the following property:

(27)$e_{s,i,i\pm 1}\ge 0,e_{s,i,i}\le 0,i=0,\cdots ,N-1,s=1,2,\cdots ,S.$

Indeed, since $\mathrm{sign}\left({\alpha }_{s,i}\right)=\mathrm{sign}\left(m_{s,i+1/2}\right)$ and in view of the assumptions of the theorem, we conclude that (27) is fulfilled for $i=0,1,\cdots ,N-2$. Next, for $i=N-1$, by also taking into account

$m_{s,N-1/2}\frac{1-\mathrm{sign}\left(m_{s,N-1/2}\right)}{2}=\left\{\begin{array}{cc}0,\hfill & m_{s,N-1/2}\ge 0,\hfill \\ m_{s,N-1/2}<0,\hfill & m_{s,N-1/2}<0,\hfill \end{array}\right.$

$m_{s,N-1/2}\frac{1+\mathrm{sign}\left(m_{s,N-1/2}\right)}{2}=\left\{\begin{array}{cc}{m}_{s,N-1/2}>0,\hfill & m_{s,N-1/2}\ge 0,\hfill \\ 0,\hfill & m_{s,N-1/2}<0,\hfill \end{array}\right.$

we deduce that (27) holds.

Further, we apply Lemma 1, taking

$\begin{array}{ccc}& & W=[C_{1,0},C_{1,1},\cdots ,C_{1,N},C_{2,0},C_{2,1},\cdots ,C_{2,N},\cdots ,C_{S,0},C_{S,1},\cdots ,C_{S,N}],\hfill \\ & & \mathcal{F}(t,W)={\mathcal{L}}_i\left({\mathbf{C}}_s\right),i=0,1,\cdots ,N,s=1,2,\cdots ,S.\hfill \end{array}$

Consequently, for $C_{s,i}=0$ and $C_{s,j}\ge 0$, $j=i\pm 1$, taking into account that in this case $B_{s,i}=r_{s,i}$ and in view of (27) and Theorem 1, we obtain

${\mathcal{L}}_i\left(C_s\right)=\frac{1}{{\hslash }_i}\left(e_{s,i,i-1}{C}_{s,i-1}+e_{s,i,i+1}{C}_{s,i+1}\right)+B_{s,i}+f_{s,i}\left(t\right)\ge 0,i=0,1,\cdots ,N-1.$

For $i=N$, from Theorem 1, we have $r_{s,N}+f_{s,N}\ge 0$. Therefore, from Lemma 1, we conclude that systems (25) and (26) are positive. □

4. Full Discretization

Positivity-preserving properties impose restrictions on time integration methods. Our aim is to construct an unconditionally (without restrictions on the mesh step sizes) positivity-preserving numerical scheme, which can be realized in an efficient manner, namely to achieve optimal accuracy and save computational time.

For concreteness, we consider the representation (1) for the function $r_s$. We suppose that the assumptions of Theorem 1 are fulfilled and rewrite the function $B_s$ in the form

$\begin{array}{ccc}\hfill B_s(\xi ,C_1,\dots ,C_S)& =& {\gamma }_{s,s}{C}_s+\sum _{j=1}^S{\beta }_{s,s,j}{C}_s{C}_j+\sum _{i=1}^S{\beta }_{s,i,s}{C}_i{C}_s-d_s(\xi )C_s\hfill \\ & & +r_s(C_1,C_2\cdots ,C_{s-1},0,C_{s+1},C_S)\hfill \\ & =& C_s\left[\sum _{j=1}^S{\beta }_{s,s,j}{C}_j+\sum _{\genfrac{}{}{0pt}{}{j=1}{j\ne s}}^S{\beta }_{s,j,s}{C}_j\right]+({\gamma }_{s,s}-d_s(\xi ))C_s\hfill \\ & & +r_s(C_1,C_2,\cdots C_{s-1},0,C_{s+1},\cdots ,C_S)\hfill \\ & =& C_s\sum _{j=1}^S{\eta }_{s,j}{C}_j+({\gamma }_{s,s}-d_s\left(\xi \right))C_s\hfill \\ & & +r_s(C_1,C_2,\cdots C_{s-1},0,C_{s+1},\cdots ,C_S),\hfill \end{array}$

where

${\eta }_{s,j}=\left\{\begin{array}{cc}{\beta }_{s,s,j}+{\beta }_{s,j,s},\hfill & j=1,2,\cdots ,S,j\ne s,\hfill \\ {\beta }_{s,s,s},\hfill & j=s,\hfill \end{array}\right.$

In order to discretize the problem with respect to time, we introduce uniform mesh with grid nodes $t_n=n\tau$, $n=0,1,\cdots ,M$, $\tau =T/M$. The mesh function $v_s$ at grid node $(x_i,t_n)$ is denoted by $v_{s,i}^n$.

We also introduce the notation $v^{+}=max\{0,v\}$, $v^{-}=max\{0,-v\}$ and therefore $v=v^{+}-v^{-}$. Furthermore, we apply explicit–implicit time stepping and design the numerical scheme to decouple discrete equations, using the most recent computed solution (similar to the Gauss–Seidel procedure), while simultaneously preserving positivity. Thus, the full discretization of (25)–(26) is

(28)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{C_{s,i}^{n+1}-C_{s,i}^n}{\tau }-{\mathcal{L}}_i^e\left({\mathbf{C}}_s^{n+1}\right)+\left({\mathcal{B}}_i^{-}\left({\mathbf{C}}_s\right)+d_{s,i}^{+}\right)C_{s,i}^{n+1}=\left({\mathcal{B}}_i^{+}\left({\mathbf{C}}_s\right)+d_{s,i}^{-}\right)C_{s,i}^n+f_{s,i}\left(t_{n+1}\right)\hfill \\ & +r_s(C_{1,i}^{n+1},C_{2,i}^{n+1},\cdots C_{s-1,i}^{n+1},0,C_{s+1,i}^n,\cdots ,C_{S,i}^n),i=0,1,\dots ,N-1,\hfill \\ & \left\{\begin{array}{cc}{C}_{s,N}^{n+1}=0,\hfill & \mathrm{for}\mathrm{b}.\mathrm{c}.\left(5\right),\hfill \\ \frac{C_{s,N}^{n+1}-C_{s,N}^n}{\tau }+{\mathcal{B}}_N^{-}\left({\mathbf{C}}_s\right)C_{s,N}^{n+1}={\mathcal{B}}_N^{+}\left({\mathbf{C}}_s\right)C_{s,N}^n+f_{s,N}\left(t_{n+1}\right)\hfill & \\ +r_s(C_{1,N}^{n+1},C_{2,N}^{n+1},\cdots C_{s-1,N}^{n+1},0,C_{s+1,N}^n,\cdots ,C_{S,N}^n),\hfill & \mathrm{for}\mathrm{b}.\mathrm{c}\left(7\right),\hfill \end{array}\right.\hfill \\ & C_{s,i}^0=C_s^0\left({\xi }_i\right),i=0,1,\cdots ,N,\hfill \end{array}\end{array}$

where $s=1,2,\cdots ,S$ and

$\begin{array}{ccc}& & {\mathcal{L}}_i^e\left({\mathbf{C}}_s^n\right)=\frac{1}{{\hslash }_i}\left(e_{s,i,i-1}{C}_{s,i-1}^n+e_{s,i,i}{C}_{s,i}^n+e_{s,i,i+1}{C}_{s,i+1}^n\right),\hfill \\ & & {\mathcal{B}}_i^{\pm }\left({\mathbf{C}}_s\right)={\mathcal{B}}_i^{\pm }\left({\mathbf{C}}_1^n,{\mathbf{C}}_2^n,\cdots ,{\mathbf{C}}_S^n,{\mathbf{C}}_s^{n+1}\right)=\sum _{j=1}^{s-1}{\eta }_{s,j}^{\pm }{C}_{j,i}^{n+1}+\sum _{j=s}^S{\eta }_{s,j}^{\pm }{C}_{j,i}^n+{\gamma }_{s,s}^{\pm }\left({\xi }_i\right).\hfill \end{array}$

At each time layer, we compute each discrete equation in (28) separately. Namely, let the solution ${\mathbf{C}}_s^n$, $s=1,2,\cdots ,S$ be known (already computed). To find ${\mathbf{C}}_s^{n+1}$, $s=1,2,\cdots ,S$, we consequently solve S discrete problems:

(1) Solve (28) for $s=1$ to find ${\mathbf{C}}_1^{n+1}$;

(2) Solve (28) for $s=2$ and known ${\mathbf{C}}_1^{n+1}$ to find ${\mathbf{C}}_2^{n+1}$;

⋯⋯

(S) Solve (28) for $s=S$ and known ${\mathbf{C}}_1^{n+1},{\mathbf{C}}_2^{n+1},\cdots ,{\mathbf{C}}_{S-1}^{n+1}$ to find ${\mathbf{C}}_S^{n+1}$.

Thus, instead of solving one system of $S\times N$ equations, we compute S number systems (discrete problems) of N equations.

In the next statement, we establish an unconditionally positivity-preserving property of the numerical scheme (28).

Theorem 3.

Let the assumptions of Theorem 2 hold. Then,

$C_{s,i}^n\ge 0,i=0,1,\cdots ,N,n=0,1,\cdots ,M,s=1,2,\cdots ,S.$

Proof. 

We prove the statement of the theorem by induction. By the conditions of the theorem, we have ${\mathbf{C}}_s^0\ge 0$, $s=1,2,\cdots ,S$. Suppose that ${\mathbf{C}}_s^n\ge 0$, $s=1,2,\cdots ,S$. We will show that ${\mathbf{C}}_s^{n+1}\ge 0$, $s=1,2,\cdots ,S$.

Denote by ${\mathcal{D}}_s={\left\{{\tilde{d}}_{i,j}^s\right\}}_{i,j=0,0}^{N,N}$, $s=1,2,\cdots ,S$ the coefficient matrix of the s-th discrete problem in (28) with entries

$\begin{array}{ccc}& & {\tilde{d}}_{i,i}^s=\frac{1}{\tau }-\frac{e_{s,i,i}}{{\hslash }_i}+{\mathcal{B}}_i^{\pm }\left({\mathbf{C}}_s\right)+d_{s,i}^{+},\hfill \\ & & {\tilde{d}}_{i,i+1}^s=-\frac{e_{s,i,i+1}}{{\hslash }_i},{\tilde{d}}_{i,i-1}^s=-\frac{e_{s,i,i+1}}{{\hslash }_i},i=0,1,\cdots ,N-1,\hfill \\ & & {\tilde{d}}_{N,N}^s=\left\{\begin{array}{cc}1,\hfill & \mathrm{for}\mathrm{b}.\mathrm{c}.\left(5\right),\hfill \\ \frac{1}{\tau }+{\mathcal{B}}_N^{\pm }\left({\mathbf{C}}_s\right),\hfill & \mathrm{for}\mathrm{b}.\mathrm{c}.\left(7\right).\hfill \end{array}\right.\hfill \end{array}$

Consider $s=1$. From the induction assumption and definition of ${\eta }_{s,j}^{-}\ge 0$, $s_{s,i}^{+}\ge 0$ and ${\gamma }_{s,s}^{-}\ge 0$, we obtain ${\mathcal{B}}_i^{\pm }\left({\mathbf{C}}_1\right)\ge 0$, $i=0,1,\cdots ,N$, and in view of sign conditions (27), we derive

(29)${\tilde{d}}_{i,i}^1>0,{\tilde{d}}_{i,i\pm 1}^1\le 0,i=0,1,\cdots ,N-1,{\tilde{d}}_{N,N}^1>0.$

Furthermore,

(30)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\tilde{d}}_{i,i}^1-|{\tilde{d}}_{i,i+1}^1|-|{\tilde{d}}_{i,i-1}^1|& \ge \frac{1}{\tau }-\frac{e_{1,i,i}}{{\hslash }_i}-\frac{e_{1,i,i+1}}{{\hslash }_i}-\frac{e_{1,i,i-1}}{{\hslash }_i}+d_{1,i}^{+},\hfill \\ & \ge \frac{1}{\tau }-\frac{1}{{\hslash }_i}\left(m_{1,i+1/2}{(1-{\xi }_{i+1/2})}^2-m_{1,i-1/2}{(1-{\xi }_{i-1/2})}^2\right)+d_{1,i}^{+}\hfill \\ & =\frac{1}{\tau }-(d_{1,i}^{+}-d_{1,i}^{-})+d_{1,i}^{+}+O\left({\hslash }_i^2\right)>0,i=0,1,\cdots ,N-1,\hfill \end{array}\end{array}$

taking into account that $d_s\left(\xi \right)=\frac{d{q}_s\left(\xi \right)}{d\xi }$. For $i=N$, the diagonal domination is obvious.

Again, using the induction assumption, the definition of ${\eta }_{s,j}^{+}\ge 0$, $d_s^{-}\ge 0$, ${\gamma }_{s,s}^{+}\ge 0$ and the assumptions of the theorem, we conclude that the right-hand side of the discrete problem (28) for $s=1$ is non-negative. Since the inverse of the coefficient matrix ${\mathcal{D}}_1$ is non-negative [37], we deduce that ${\mathbf{C}}_1^{n+1}\ge 0$.

We proceed with the next problem (28) for $s=2$. We apply the same considerations as above, taking into account that ${\mathcal{B}}_2^{\pm }\left({\mathbf{C}}_2\right)\ge 0$, $i=0,1,\cdots ,N$ and prove that ${\mathbf{C}}_2^{n+1}\ge 0$. We continue in the same manner with discrete problems (28) for $s=3,4,\cdots ,S$. Since ${\mathbf{C}}_{s-1}^{n+1}\ge 0$, by the induction assumption, the definition of ${\eta }_{s,j}^{\pm }\ge 0$, $d_s^{\pm }\ge 0$, ${\gamma }_{s,s}^{\pm }\ge 0$, inequalities (27), and taking into account that the same conditions as (29) and (30) are valid for each $s=1,2,\cdots ,S$

${\tilde{d}}_{i,i}^s>0,{\tilde{d}}_{i,i\pm 1}^s\le 0,{\tilde{d}}_{i,i}^s>|{\tilde{d}}_{i,i+1}^s|+|{\tilde{d}}_{i,i-1}^s|i=0,1,\cdots ,N-1,{\tilde{d}}_{N,N}^s>0,$

we obtain ${\mathbf{C}}_s^{n+1}\ge 0$, $s=1,2,\cdots ,S$. □

5. Numerical Tests

In this section, we illustrate the efficiency of the developed numerical method for solving the problem (2), using (3), (4), (6), and (7) by the numerical scheme (28). For the numerical tests, we take some of the data from [5,33]. We set $S=3$, $a=0.005$ and

$\begin{array}{cc}& r_s={\gamma }_{s,s}{c}_s+{\beta }_{s,2,2}{c}_2^2,s=1,2,3,\mathrm{for}\mathrm{Examples}2,4,\\ & r_s={\gamma }_{s,s}{c}_s+\sum _{\genfrac{}{}{0pt}{}{i=1}{i\ne s}}^3{\beta }_{s,i,i}{c}_i{c}_i,s=1,2,3,\mathrm{for}\mathrm{Example}3,\\ & K_1=5,K_2=2,K_3=3,w=1.\end{array}$

In order to achieve second-order convergence in space, we use the mesh grading technique for the primal mesh, and the nodes of the dual mesh are chosen so that the approximate flux has second-order accuracy. Following [21,32], superconvergence points can be found in each subinterval of the primal mesh and used to construct the dual mesh:

(31)${\xi }_{j+1/2}=\left\{\begin{array}{cc}1-(1-{\xi }_i)(1-{\xi }_{i+1})\left(ln\frac{1-{\xi }_i}{1-{\xi }_{i+1}}\right)/h_i,\hfill & {\alpha }_{s,i}=-1,\hfill \\ 1-h_i/ln\frac{1-{\xi }_j}{1-{\xi }_{j+1}},\hfill & {\alpha }_{s,i}=0,\hfill \\ ({\xi }_{j+1}+{\xi }_j)/2,\hfill & {\alpha }_{s,i}=1,\hfill \\ 1-\frac{{\alpha }_{s,i}}{{\alpha }_{s,i}+1}\frac{{(1-{\xi }_{i+1})}^{{\alpha }_{s,i+1}}-{(1-{\xi }_i)}^{{\alpha }_{s,i+1}}}{{(1-{\xi }_{i+1})}^{{\alpha }_{s,i}}-{(1-{\xi }_i)}^{{\alpha }_{s,i}}},\hfill & {\alpha }_{s,i}\ne \{-1,0,1\}.\hfill \end{array}\right.$

Therefore, the dual mesh is different for each discrete equation. The grading primal mesh with nodes concentrated close to $\xi =1$ is

(32)${\xi }_i=1-{\left(\frac{N-i}{N}\right)}^2,i=0,1,\cdots ,N.$

Example 2.

(Convergence) In order to estimate the accuracy and convergence rate, we deal with the exact solution of the problem (2). Function $f_s(t,z)$ is chosen such that $c_s=e^{-t/s}cos\left(\pi (e^{2az}-1)/\left[2(e^{2az}-1)\right]\right)$; $s=1,2,3$ to be the exact solution of (2); ${\nu }_s=0$, $s=1,2,3$. Let ${\tilde{C}}_s$ be the transformed numerical solution $C_s$ on the original z-domain. Then, the errors (${\mathcal{E}}_s$) and the corresponding order of convergence (${\mathcal{CR}}_s$) in maximal discrete norm are given by

${\mathcal{E}}_s={\mathcal{E}}_s\left(N\right)=\underset{0\le i\le N}{max}|{\tilde{C}}_{s,i}-c_s(t,z_i)|,{\mathcal{CR}}_s={log}_2\frac{{\mathcal{E}}_s\left(N\right)}{{\mathcal{E}}_s\left(2N\right)},s=1,2,3.$

We consider the following set of parameters:

SP 1: ${\gamma }_{1,1}=-4$, ${\gamma }_{2,2}=2$, ${\gamma }_{3,3}=-1$, ${\beta }_{1,2,2}=3$, ${\beta }_{2,2,2}=-2$, ${\beta }_{3,2,2}=2$,

SP 2: ${\gamma }_{1,1}={\gamma }_{3,3}=-1000$, ${\gamma }_{2,2}=0$, ${\beta }_{1,2,2}={\beta }_{3,2,2}=2000$, ${\beta }_{2,2,2}=-2000.$

We give the computational results for uniform and nonuniform meshes (31) and (32). For uniform primal mesh, the nodes of the dual mesh are chosen to be in the middle of the subintervals $I_i$, $i=0,1,\cdots ,N-1$. The computations are performed for $\tau =h^2$, $h=\underset{0\le i\le N-1}{max}{h}_i$, and $T=1$.

In Table 1 and Table 2, we present results for uniform and nonuniform spatial meshes, respectively, with different numbers of grid nodes and parameters from SP 1. We observe that the spatial order of convergence on uniform mesh is close to one, while on the nonuniform meshes, it is close to two. In Table 3, we give computational results on the nonuniform spatial mesh with different numbers of grid nodes and parameters from SP 2. It is evident that the accuracy is $O(\tau +h^2)$.

Example 3.

(Comparison) The main advantage of the proposed numerical method (28) is that it allows us to compute the solution of (2) on a small interval. A commonly used approach for computing the solution of PDEs, defined on an unbounded domain is to truncate the domain and impose the boundary conditions the same as on the infinity on the cut boundary. In this case, it is well known that we have to compute the solution on a large enough domain, which leads to increasing the number of grid nodes and a more time-consuming computational process. We illustrate this for our model problem (2). Specifically, we compare the performance of the transformation method (28) with the finite volume method (FVM) on the original z spatial interval, where the nonlinear term is approximated as in (28), for solving (2).

Let ${\gamma }_1=-4$, ${\gamma }_2=-2$, ${\gamma }_3=-1$, ${\beta }_{1,3,3}=3$, ${\beta }_{2,1,1}={\beta }_{3,2,2}=2$, ${\beta }_{1,2,2}={\beta }_{2,3,3}={\beta }_{3,1,1}=0$, $T=0.5$, $L=100$, $a=0.005$.

We determine the right-hand side, the left boundary condition, and initial condition in (2), such that $c_1(t,z)=e^t(1-\mathrm{erfc}\left(1/\sqrt{2+z})\right)$, $c_2(t,z)=e^{t/2}(1-\mathrm{erfc}\left(2/\sqrt{1+z})\right)$, $c_3(t,z)=e^{t/3}(1-\mathrm{erfc}\left(2/\sqrt{1+z})\right)$ are exact solutions.

The numerical method of (28) is computed on the spatial interval $[0,0.999]$ for $N=40$, while for the FVM, the spatial interval is $[0,100]$ and 400 space grid nodes are used. For all computations, we set $\tau =0.0025$.

Exact and numerical solutions computed by (28) and FVM at the final time are ploted on Figure 1.

The solution, obtained by the FVM, is relevant only in a quarter of the spatial interval, and the CPU time for the computational process is 6.82 s. In contrast, computing (2) by the numerical scheme (28), we obtain an accurate solution over the whole spatial interval, and the CPU time is 32.34 s.

Example 4.

(Point sources of pollution) We consider the original problem, with the forcing formed by the source at the point $z_s^{*}$, i.e., $f_s(t,z)=Q_s\left(t\right){\delta }_s(z-z_s^{*})$, where ${\delta }_s(·)$ is the Dirac delta function. In the problem (3)–(7), we consider the following regularization of the Dirac delta function.

${\delta }_s^h(t,\xi )=\left\{\begin{array}{cc}\frac{h_i{h}_{i-1}-|{\xi }^{*}-\xi |}{4h_i{h}_{i-1}},\hfill & \xi \in ({\xi }^{*}-h_{i-1}-h_{i-2},{\xi }^{*}+h_i-h_{i+1}),\hfill \\ 0,\hfill & \mathit{otherwise}.\hfill \end{array}\right.$

Let $Q_1\left(t\right)=5t^2(1+t)$, $Q_2\left(t\right)=4+t^2$, $Q_3\left(t\right)=10t(t+1)$, $c_{1,0}=c_{2,0}=0$, $c_{3,0}=5$, $z_1^{*}=127$, $z_2^{*}=389$, $z_3^{*}=465$, ${\nu }_1=0.2$, ${\nu }_2=0.4$, ${\nu }_3=0.5$ and consider model parameters as in Example 2, SP 1. On Figure 2, we plot numerical solutions $C_s$, $s=1,2,3$ in the original z space variable, computed on nonuniform mesh for $N=80$, $\tau =h$ and $T=1$. Next, we perform the same experiment for $Q_1=1$, $Q_2=4$, and $Q_3=8$. The numerical solution is given in Figure 3.

The numerical solution is non-negative, which confirms the statement of Theorem 3. The Figure 2 and Figure 3 show the increasing concentrations of the chemical particles over time for point source pollution.

6. Conclusions

In this paper, we considered a one-dimensional nonlinear problem of air pollution described by a parabolic system with a non-standard nonlinear right-hand side.

• Applying log-transformation, the original problem defined on a semi-infinite spatial interval is reformulated as an equivalent degenerate problem on an arbitrary small finite interval, specified in advance.

• The well-posedness and non-negativity of the solution of the transformed problem are proved.

• We develop a second-order spatially accurate, positivity-preserving iteration-free numerical method, based on a fitted finite volume method and implicit–explicit time stepping.

• The numerical scheme is implemented efficiently by decoupling the discrete equations and computing them on graded primal and superconvergence dual meshes, ensuring that both the numerical fluxes and solutions have second-order accuracy.

• Numerical experiments on uniform and nonuniform meshes confirm the theoretical results.

An extension of the present theoretical and computational results to more complex nonlinear air pollution problems, as well as to the two- and three-dimensional cases, will be addressed in our future work.

Author Contributions

Conceptualization, L.G.V.; methodology, M.N.K. and L.G.V.; investigation, M.N.K. and L.G.V.; resources, M.N.K. and L.G.V.; writing—original draft preparation, M.N.K. and L.G.V.; writing—review and editing, M.N.K. and L.G.V.; visualization, M.N.K. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data are contained within the article.

Acknowledgments

The authors are very grateful to the anonymous reviewers whose valuable comments and suggestions improved the quality of the paper.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Exact solution (solid red line), numerical solution computed by (28) (line with blue circles), and numerical solution computed by FVM (dashed green line), Example 3.

Figure 2. Numerical solutions [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.]; nonuniform meshes (31) and (32), Example 4.

Figure 3. Numerical solutions [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.]; nonuniform meshes (31) and (32), Example 4.

Table 1

Errors and spatial order of convergence, uniform mesh, SP 1, Example 2.

N	${\mathcal{E}}_1$	${\mathcal{CR}}_1$	${\mathcal{E}}_2$	${\mathcal{CR}}_2$	${\mathcal{E}}_3$	${\mathcal{CR}}_3$
20	1.6007 $\times {10}^{-3}$		1.3068 $\times {10}^{-3}$		1.0612 $\times {10}^{-3}$
40	4.7523 $\times {10}^{-4}$	1.7520	4.9387 $\times {10}^{-4}$	1.4038	4.3916 $\times {10}^{-4}$	1.2729
80	1.3223 $\times {10}^{-4}$	1.8456	2.3522 $\times {10}^{-4}$	1.0701	1.4510 $\times {10}^{-4}$	1.5977
160	3.0249 $\times {10}^{-5}$	2.1281	1.1757 $\times {10}^{-4}$	1.0005	3.0562 $\times {10}^{-5}$	2.2472
320	6.2336 $\times {10}^{-6}$	2.2787	5.9179 $\times {10}^{-5}$	0.9903	1.1084 $\times {10}^{-5}$	1.4633
640	1.1946 $\times {10}^{-6}$	2.3836	2.9742 $\times {10}^{-5}$	0.9926	5.4960 $\times {10}^{-6}$	1.0120

Table 2

Errors and spatial order of convergence, nonuniform meshes (31) and (32), SP 1, Example 2.

N	${\mathcal{E}}_1$	${\mathcal{CR}}_1$	${\mathcal{E}}_2$	${\mathcal{CR}}_2$	${\mathcal{E}}_3$	${\mathcal{CR}}_3$
20	5.7089 $\times {10}^{-5}$		4.1037 $\times {10}^{-3}$		3.1898 $\times {10}^{-3}$
40	1.6444 $\times {10}^{-3}$	1.7956	1.0919 $\times {10}^{-3}$	1.9101	1.2076 $\times {10}^{-3}$	1.4013
80	5.0567 $\times {10}^{-4}$	1.7013	2.8823 $\times {10}^{-4}$	1.9216	5.1056 $\times {10}^{-4}$	1.2420
160	1.4807 $\times {10}^{-4}$	1.7719	1.0954 $\times {10}^{-4}$	1.3957	1.8199 $\times {10}^{-4}$	1.4882
320	3.5136 $\times {10}^{-5}$	2.0753	2.9775 $\times {10}^{-5}$	1.8793	4.1742 $\times {10}^{-5}$	2.1243
640	7.4282 $\times {10}^{-6}$	2.2419	6.0162 $\times {10}^{-6}$	2.3072	7.5300 $\times {10}^{-6}$	2.4708

Table 3

Errors and spatial order of convergence, nonuniform meshes (31) and (32), SP 2, Example 2.

N	${\mathcal{E}}_1$	${\mathcal{CR}}_1$	${\mathcal{E}}_2$	${\mathcal{CR}}_2$	${\mathcal{E}}_3$	${\mathcal{CR}}_3$
20	3.5182 $\times {10}^{-3}$		1.4422 $\times {10}^{-3}$		3.4972 $\times {10}^{-3}$
40	8.9976 $\times {10}^{-4}$	1.9672	3.6957 $\times {10}^{-4}$	1.9643	8.9677 $\times {10}^{-4}$	1.9634
80	2.2778 $\times {10}^{-4}$	1.9819	9.3529 $\times {10}^{-5}$	1.9824	2.2695 $\times {10}^{-4}$	1.9823
160	5.7353 $\times {10}^{-5}$	1.9897	2.3542 $\times {10}^{-5}$	1.9906	5.7163 $\times {10}^{-5}$	1.9892
320	1.4393 $\times {10}^{-5}$	1.9945	5.9050 $\times {10}^{-6}$	1.9952	1.4341 $\times {10}^{-5}$	1.9950
640	3.5921 $\times {10}^{-6}$	2.0025	1.4798 $\times {10}^{-6}$	1.9965	3.5911 $\times {10}^{-6}$	1.9976

References

1. Georgiev, K.; Zlatev, Z. Implementation of sparse matrix algorithms in an advection–diffusion–chemistry module. J. Comput. Appl. Math.; 2011; 236, pp. 342-353. [DOI: https://dx.doi.org/10.1016/j.cam.2011.07.026]

2. Karátson, J.; Kurics, T. A preconditioned iterative solution scheme for nonlinear parabolic systems arising in air pollution modeling. Math. Model. Anal.; 2013; 18, pp. 641-653. [DOI: https://dx.doi.org/10.3846/13926292.2013.868841]

3. Zlatev, Z. Computer Treatment of Large Air Pollution Models; Kluwer Academic Publishers: Dordrecht, The Netherlands, Boston, MA, USA, London, UK, 1995.

4. Zlatev, Z.; Dimov, I. Computational and Numerical Challenges in Environmental Modeling; Elsevier: Amsterdam, The Netherlands, 2006; 373.

5. Mamonov, A.; Tsai, Y.-H.R. Point source identification in nonlinear advection–diffusion–reaction systems. Inverse Probl.; 2013; 29, 035009. [DOI: https://dx.doi.org/10.1088/0266-5611/29/3/035009]

6. Faragó, I. New Operator Splitting Methods and Their Analysis. Numerical Methods and Applications. NMA 2006; Boyanov, T.; Dimova, S.; Georgiev, K.; Nikolov, G. Lecture Notes in Computer Science Springer: Berlin/Heidelberg, Germany, 2007; Volume 4310.

7. Samarskii, A.; Matus, P.; Vabishchvich, P. Difference Schemes with Operator Factors; Kluwer Academic Publishers: Boston, MA, USA, Dordrecht, The Netherlands, London, UK, 2002.

8. Sportise, B. An analysis of operator splitting techniques in the stiff case. J. Comput. Phys.; 2000; 161, pp. 140-168. [DOI: https://dx.doi.org/10.1006/jcph.2000.6495]

9. Stokie, M. The mathematics of atmospheric dispersion modeling. SIAM Rev.; 2011; 53, pp. 349-372. [DOI: https://dx.doi.org/10.1137/10080991X]

10. Faragoó, I.; Horvath, R. Discrete maximum principle and adequate discretizations of linear parabolic problems. SIAM J. Sci. Comput.; 2006; 28, pp. 2313-2336. [DOI: https://dx.doi.org/10.1137/050627241]

11. Valkov, R. Convergence of a finite volume element method for a generalized Black-Scholes equation transformed on finite interval. Numer. Algor.; 2015; 68, pp. 61-80. [DOI: https://dx.doi.org/10.1007/s11075-014-9838-8]

12. Kim, J.; Cho, S.Y. Computation accuracy and efficiency of the time-splitiing method in solving atmospheric transport/chemistry equations. Atmos. Environ.; 1997; 31, pp. 2215-2224. [DOI: https://dx.doi.org/10.1016/S1352-2310(97)88636-0]

13. Djellab, N.; Boureghda, A. A moving boundary model for oxygen diffusion in a sick cell. Comput. Methods Biomech. Biomed. Eng.; 2022; 25, pp. 1402-1408. [DOI: https://dx.doi.org/10.1080/10255842.2021.2024168]

14. Boureghda, A.; Djellab, N. Du Fort-Frankel Finite Difference Scheme for Solving of Oxygen Diffusion Problem inside One Cell. J. Comput. Theor. Transp.; 2023; 52, pp. 363-373. [DOI: https://dx.doi.org/10.1080/23324309.2023.2271229]

15. Boureghda, A. Solution to an ice melting cylindrical problem. J. Nonlinear Sci. Appl.; 2016; 9, pp. 1440-1452. [DOI: https://dx.doi.org/10.22436/jnsa.009.04.04]

16. Pao, C.V. Nonlinear Parabolic and Elliptic Equations; Springer Science & Business Media: Petaling Jaya, Malaysia, 1992; 777p.

17. Matus, P.P. The maximum principle and some of its applications. Comput. Meth. Appl. Math.; 2002; 2, pp. 50-91. [DOI: https://dx.doi.org/10.2478/cmam-2002-0004]

18. Mitideri, E.; Sweers, G. Weakly coupled elliptic systems and positivity. Math. Nachrichten; 1995; 173, pp. 256-286. [DOI: https://dx.doi.org/10.1002/mana.19951730115]

19. Dang, Q.A.; Ehrhardt, M. Adequate numerical solution of air pollution problems by positive difference schemes on unbounded domains. Math. Comp. Model.; 2006; 44, pp. 834-856. [DOI: https://dx.doi.org/10.1016/j.mcm.2006.02.016]

20. Fatahillah, A.; Masyhudi, M.A.; Setiawan, T.B. Numerical analysis of air pollutant dispersion in steam power plant area using the finite volume method. J. Phys. Conf. Ser.; 2020; 1490, 012002. [DOI: https://dx.doi.org/10.1088/1742-6596/1490/1/012002]

21. Koleva, M.N.; Vulkov, L.G. Positivity-preserving finite volume difference schemes for atmospheric dispersion models with degenerate vertical diffusion. Comp. Appl. Math.; 2022; 41, 406. [DOI: https://dx.doi.org/10.1007/s40314-022-02083-y]

22. Lanser, D.; Verver, J.G. Analysis of operator splitting for advection-diffusion problems from air pollution modelling. J. Comput. Appl. Math.; 1999; 111, pp. 201-206. [DOI: https://dx.doi.org/10.1016/S0377-0427(99)00143-0]

23. Luoc, N.V.; Dang, G.A.; Dieu, N.C. Analytic and numerical solution of some problems of air pollution. SEA Bull. Math. Special Issue; 1993; pp. 105-117.

24. Marchuk, G.I. Mathematical modelling in environmental problems. Studies in Mathematics and Its Applications; Elsevier: Amsterdam, The Netherlands, 1986; Volume 16.

25. Morton, K.W. Numerical Solution of Convection-Diffusion Problems; Chapman and Hill: London, UK, 1995.

26. Zaia, F.N.; Gunawanb, A.Y. Effects of inversion layer on the atmospheric pollutant dispersion from a high chimney. J. Indones. Math. Soc.; 2023; 29, pp. 299-310. [DOI: https://dx.doi.org/10.22342/jims.29.3.1597.299-310]

27. Alexandrov, V.N.; Sameh, A.; Siddique, Y.; Zlatev, Z. Numerical integration of chemical ODE problems arising in air pollution models. Environ. Model. Assesment; 1997; 2, pp. 365-377. [DOI: https://dx.doi.org/10.1023/A:1019086016734]

28. Lin, J.; Reutskiy, S.Y.; Lu, J. A novel meshless method for fully nonlinear advection–diffusion-reaction problems to model transfer in anisotropic media. Appl. Math. Comput.; 2018; 339, pp. 459-476. [DOI: https://dx.doi.org/10.1016/j.amc.2018.07.045]

29. Todorov, V.; Kandilarov, J.; Dimov, I.; Vulkov, L. High-accuracy numerical methods for a parabolic system in air pollution modeling. Neural Comput. Appl.; 2020; 32, pp. 6025-6040. [DOI: https://dx.doi.org/10.1007/s00521-019-04088-x]

30. Lazarov, R.D.; Mishev, I.D.; Vassilevski, P.S. Finite volume methods for convection-diffusion problems. SIAM J. Numer. Anal.; 1996; 33, pp. 31-55. [DOI: https://dx.doi.org/10.1137/0733003]

31. Wang, S. A novel finite volume method for Black-Sholes equation governing option pricing. IMA J. Numer. Anal.; 2004; 24, pp. 699-720. [DOI: https://dx.doi.org/10.1093/imanum/24.4.699]

32. Wang, S.; Shang, S.; Fang, Z. A superconvergence fitted finite volume method for Black-Sholes equation governing European and American options. Numer. Methods Part. Differ. Equ.; 2014; 31, pp. 1190-1208. [DOI: https://dx.doi.org/10.1002/num.21941]

33. Chernogorova, T.; Dimov, I.; Vulkov, L. Positive fitted finite volume method for semi-linear parabolic systems transformed on finite interval. AIP Conf. Proc.; 2018; 2048, 030013.

34. Oleinik, O.A.; Radkevich, E.V. Second Order Equations with Nonnegative Characteristic Form; Plenum Press: New York, NY, USA, 1973; 267.

35. Hundsdorfer, W.; Verwer, J.G. Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equation; Springer: Berlin, Germany, 2003.

36. Samarskii, A. The Theory of Difference Schems; Marcel Dekker: New York, NY, USA, 2001.

37. Varga, R.S. Matrix Iterative Analysis, 2nd Revised and Expanded ed.; Springer: Berlin/Heidelberg, Germany, 2000.