Finite difference schemes for time-dependent convection q-diffusion problem

1.   Introduction

Currently, solving non-linear differential equations in many scientific phenomena has gravitated to the attention of many researchers. Today, $q$-calculus and partial differential equation are considered the center of attention by contemporary society. In the eighteenth century, Euler and Jacobi formulate the history of $q$-calculus. Later on, F.H Jackson ^[1] was considered a pioneer who provided extensions to Euler and Jacobi's work during the twentieth century. A profound surge in the study of $q$-calculus was observed during the mid-twentieth century because researchers found its extensive applications in several mathematical and analytical branches. Nowadays, $q$-calculus has found its concern in mathematical modeling in terms of quantum computing. Details are available in ^[2,3,4].

The $q$-calculus is a novel branch that deals with studying the link between physics and mathematics. Several mathematical areas, such as number theory, combinatorics, orthogonal polynomials, basic hypergeometric functions, and quantum theory, are considered $q$-calculus as a primary source of research. This branch has provided an interesting methodology in differential transform methods; hence, suitable to draw numerical approximations for the ordinary and partial differential equations.

This field is the origin of the $q$-differential equation, which describes several physical processes involved in quantum dynamics, discrete dynamical systems, and discrete stochastic processes. It is worth mentioning that q-differential equations are elaborated on a time scale set as a ${T}_{q}$ in which $q$ stands for scale index. Based on $q$-calculus theory, several other concepts have been introduced and familiarized, such as $q$-Laplace transform, $q$-Gamma and $q$-Beta functions, $q$-Mittag–Leffler functions, $q$-Taylor expansion, and $q$-integral transform theory. Details of $q$-calculus and $q$-partial differential equations can be seen in ^{[5,6,7,8,9,10,11,12]}. Studies related to classical fractional calculus are more advanced than fractional $q$-calculus.

Recent literature review revealed few advanced studies that provide a peculiar solution to the fractional $q$-differential equation. Abdeljawad et al. in ^[13] utilized a novel generalized version of discrete fractional $q$-Gronwall inequality to verify the rareness of an initial value problem involving a non-linear delay Caputo fractional $q$-difference system.

Considering Banach's contraction mapping principle, the writer of ^[14] provides a diverse solution for resolving Caputo q-fractional boundary value problem with the $p$-Laplacian operator. Ren et al. in ^[15] proved the uniqueness of nontrivial solutions by the contraction mapping principle. Additionally, they established the existence of many positive solutions under certain circumstances using traditional fixed point theorems. By using the Ascoli-Arzela Theorem and a $q$-analogue Gronwall inequality, Zhang et al. in ^[16] gave the existence and uniqueness of the solution of the Caputo fractional $q$-differential equations. Whereas in ^[17], Zhang et al. elaborated on the existence of a unique solution in the $q$-integral space.

The combination of $q$-calculus and physics and mathematics has resulted in a wide variety of applications in combinatorics, relativity theory, mechanics, number theory, and orthogonal polynomials ^[18,19,20]. $q$-Taylor's formula was firstly introduced by Jackson ^[1], which was later on modified to its q-remainder by Jing and Fan ^[21]. He utilized the $q$-differentiation approach and formulated results on the $q$-remainder in the $q$-Taylor's formula. Ernest ^[22] derived four variable $q$-Taylor's formulas with their integral $q$-remainder; however, ^[23] Prashant et al. examined the $q$-analogue of the iterative methods, particularly the $q$-analogue of generalized Newton Raphson method for the solution of algebraic transcendental equations and compared the precision of the results obtained by the traditional methods by using this formula. Several physical phenomena of the real world involving linear and non-linear models can be interpreted using $q$-differential equations. Jafari et al. ^[24] used the Daftardar decomposition methodology to solve the $q$-difference equations and determined the method's convergence.

A homogenization function has been proposed in ^[25] to address inverse source problems of the non-linear time-fractional wave equation. The new scheme was capable of solving 2D and 3D inverse source problems of non-linear time-fractional wave equations via resolving a linear matrix system. The new scheme did not involve numerical integration, regularization, mesh generation, and fundamental solutions. A new general variant ^[26], the useful and effective type of inequality, has been obtained on Chebyshev's inequality. It was concluded that the main findings had generalized various existing results and also iterated the Chebyshev inequality in special cases. With the help of Atangana-Baleanu integral operators, new results were generated for strongly convex functions ^[27].

In literature, numerous numerical methods exist to solve ordinary and partial differential equations. Some have been adopted to solve time-dependent problems, including Runge-Kutta methods. Runge-Kutta methods for partial differential equations are multi-stage methods that produce more stable regions than the multi-step method. But, stability regions of all Runge-Kutta methods are different. These methods only discretize one temporal variable, and any other scheme can be considered for spatial discretization. Some of the existing numerical methods have been constructed for solving differential equations. But there does exist much work in Quantum Calculus. The numerical methods for solving Quantum Calculus are also a future considered area where methods will be required to solve problems. This work also consists of the numerical scheme that can be considered to solve integer-order differential equations, and the method can be considered to tackle the problem having $q$-diffusion. Like integer-order derivatives, the finite difference formulas can be derived for the $q$-derivative. These difference formulas for discretizing $q$-derivative(s) terms are based on $q$-Taylor series approach. Also, in this study, the effects of convection and $q$-diffusion are added to the existing energy balance model of ODEs. The PDEs model is more general than the ODEs model because it also shows the behavior of temperatures on spatial variable(s). The summary of the work is given in the next paragraph.

The problem of climate change is constructed by considering non-linear force terms and $q$-diffusion in the existing energy balance model. After this, the construction of the proposed scheme is given with stability analysis. The leading error terms and order of convergence for non-linear convection time-dependent problem is also provided. Later on, some discussion is presented.

Since humans interact with the environment, this interaction is responsible for changing temperature. On the other side, environmental circumstances affect humans, posing a new threat in the form of climate change. The present climate change model is based on the rate at which the global mean surface temperature (GMST) increases. The global mean surface temperature has been increased by 0.07 ℃ per decade stated in ^[28]. For the study of variation in climate, some mathematical models have been introduced. These energy balance models are simple, but these models are quite effective models of the climate system. M. Budyko ^[29] and W.D. Sellers ^[30] introduced these models approximately 50 years ago. But in this study, these model is studied in the form of ordinary and partial differential equations. Since ordinary differential equations (ODEs) provide information over one variable, the partial differential equations (PDEs) model gives more information than ODEs models. Also, the ODEs model can be considered a special case of PDEs models. In this study, the mathematical model in the form of PDEs can be expressed as:

where $T$ is the temperature of GMST, ${T}_{D}$ represents global mean deep ocean temperature, ${c}_{1}$ and ${c}_{2}$ are coefficients of convection, ${d}_{1}$ and ${d}_{2}$ are coefficients of diffusion, $\lambda$ denotes the climate feedback parameter, $\gamma$ denotes deep-ocean heat uptake parameter, $C$ and ${C}_{D}$ denotes effective heat capacity of the upper box and effective heat capacity of the deep ocean, respectively.

Instead of classical derivatives in integer form, $q$-derivatives for diffusion are considered. The $q$-derivative has been defined in the literature.

Definition ^[31]: The $q$-derivative of $v\left(x\right)$ is expressed as:

At this stage of the solution procedure, stability conditions corresponding to the ODEs model are analyzed. So the equilibrium points of Eqs (1.1) and (1.2) using ${c}_{1} = {c}_{2} = {d}_{1} = {d}_{2} = 0$ and $m = 0.5$ can be found by solving the equations

Solving Eqs (1.4) and (1.5) yields two points

The Jacobin of the system (1.4) and (1.5) can be expressed as

The Jacobin evaluated at ${B}_{1}$ is given as

The characteristic polynomial for matrix (1.8) is expressed as

According to the Routh-Hurwitz criterion, the system of ODEs is stable. If the following inequality holds,

2.   Construction of numerical schemes

In this section, a numerical scheme is constructed, which can be used to solve any time-dependent partial differential equation. The scheme is explicit in both of its stages. The first stage of the scheme is just the forward Euler scheme. So for Eq (1.1), the first stage of the scheme is expressed as:

and the second stage of the scheme is given by

Substitution Eq (2.1) into Eq (2.2), it is obtained:

Expanding ${T}_{i}^{n+1}$ using the Taylor series expansion of the form

Substituting Eq (2.4) into Eq (2.3) and comparing coefficients of ${T}_{i}^{n}, \Delta t{\left(\frac{\partial T}{\partial t}\right)}_{i}^{n}$ and ${\left(\Delta t\right)}^{2}{\left(\frac{{\partial }^{2}T}{\partial {t}^{2}}\right)}_{i}^{n}$on both sides of the resulting equations gives

Solving Eq (2.5) gives

Consider the central difference formula for the second-order partial $q$-derivative term, which can be derived using $q$-Taylor series expansion,

Substituting numerical approximation (2.7) into Eq (1.1) and employing the first stage of the constructed scheme for time discretization, it is obtained:

and the second stage of the proposed scheme is expressed as

where

3.   Stability analysis

Before starting stability analysis of the constructed scheme, the system of Eqs (1.1) and (1.2) are expressed in a single vector-matrix equation as:

where ${\boldsymbol{U}} = {\left[T, {T}_{D}\right]}^{t}, A = \left[\begin{array}{cc}{c}_{1} & 0\\ 0 & {c}_{2}\end{array}\right], B = \left[\begin{array}{cc}{d}_{1} & 0\\ 0 & {d}_{2}\end{array}\right]$.

Employing the constructed numerical scheme to Eq (3.1) gives

and

By adopting Von Neumann's stability criteria, consider the transformations

Substituting transformation (3.4) into Eq (3.2) yields

Dividing both sides of Eq (3.5) by ${e}^{iI\psi }$, it is obtained

Eq (3.6) can be written as

where $I.D$ denotes the identity matrix and ${c}_{4} = \frac{\Delta t}{2\Delta x}, d = \frac{\Delta t}{{\left(\Delta x\right)}^{2}}$.

Similarly, substituting transformations from (3.4) into Eq (3.3) and dividing the resulting equation by ${e}^{iI\psi }$, it gives

Substituting Eq (3.7) into Eq (3.8) and rewrite the resulting equation in the form

The stability condition is expressed as

where ${\mu }_{A}$ & ${\mu }_{B}$ are maximum eigenvalues of $A$ and $B$, respectively.

4.   Error analysis

Since the proposed scheme is second-order accurate in time and space, the leading error terms in both space and time can be found by Taylor series analysis. For discretization of convection term, the Taylor series expansions are given as

Subtracting Eq (4.2) from Eq (4.1), and rewrite the resulting equation in the form

The leading error term in the discretization of convection term is the third-order partial derivative term and also multiple of ${\left(\Delta x\right)}^{2}$And therefore, the discretization is second-order accurate.

For finding leading errors in the spatial discretization of $q$-diffusion term, the $q$-Taylor series expansions are expressed as

Adding Eqs (4.4) and (4.5), and rewrite the resulting equation in the form of

The leading error term for the discretization of $q$-diffusion term can be seen in Eq (4.6), which contains a fourth-order spatial $q$-derivative term, which is multiple of ${\left(\Delta x\right)}^{2}$Therefore, the discretization of $q$-diffusion term is second-order accurate.

The leading error term in temporal discretization is given as

5.   Order of convergence

For finding the order of convergence of the proposed scheme, a non-linear parabolic equation is given as

subject to initial and boundary conditions

The exact solution of the problem (5.1) and (5.2) is given as ^[32]

The comparison of three numerical schemes for finding order of convergence of the problems (5.1) and (5.2) is shown in Table 1. The norm of the error is also shown in Table 1, which shows the accuracy of the considered numerical scheme for finding the solution to the problem. Theoretically, the order of accuracy of each employed numerical scheme for constructing Table 1 is the same, which is two. Still, the order of convergence of all three schemes is less than two. One of the three schemes is the classical finite difference scheme in which second-order central temporal discretization is adopted with second-order central spatial discretization. The norm of error produced by the central scheme is less than those given by the other two schemes. The following formula computes the order of convergence

where ${L}_{\infty }^{C}$ is the maximum norm of error vector at the current temporal step size ${\Delta t}^{c}$ and ${L}_{\infty }^{p}$ is the maximum norm of the error vector at the previous temporal step size ${\Delta t}^{p}$.

6.   Results and discussion

The constructed numerical schemes are explicit and second-order accurate. This accuracy of the scheme can be verified by observing the construction procedure of the scheme. Since the scheme is constructed by balancing second-order derivative terms in Taylor series expansion, it is second-order accurate. The consistency of the scheme can be verified from its order of accuracy. If the scheme is at least first-order accurate, then it is consistent. Since the constructed scheme is second-order accurate so it is also consistent. The stability of the scheme has been shown in the previous section. Using the condition of stability analysis, the scheme is conditionally stable, and it is consistent. Therefore, it is conditionally convergent for linear time-dependent partial differential having $q$-diffusion. The stability region of the constructed scheme depends on the choice of the parameter $a$. Different values of the parameter $a$ produce different stability regions.

The system of ordinary differential equations (1.1) and (1.2) is also solved with Matlab solver ode45. The ode45 solver can be used to solve linear and non-linear differential equations having initial conditions. According to Matlab code, the solver converges for small values of the heat source parameter $Q$ when it is applied for solving ordinary differential equations (1.1) and (1.2) using ${c}_{1} = {c}_{2} = {d}_{1} = {d}_{2} = 0$. Figures 1 and 2 compare the constructed numerical scheme with the forward Euler method. The error is found by finding the absolute difference between the constructed/Euler schemes and Matlab solver ode45. From this comparison, it can be concluded that the constructed scheme produces less error than that given by the existing forward Euler scheme using two different sets of numerical values of the contained parameters. Figure 3 shows the impact of $q$ on global mean surface temperature $T$ and global mean deep ocean temperature ${T}_{D}$. Both types of temperature have escalating behavior on the temporal variable $t$, but the global mean deep ocean temperature has dual behavior on the spatial variable$x$. Both global mean surface and deep ocean temperatures over temporal and spatial variables $t$ and $x$ increase by enhancing heat source parameter $Q$, this kind of behavior can be seen in Figure 4. Figures 5–8 show the contours of global mean surface and deep ocean temperatures using two different sets of numerical values of the parameters.

7.   Conclusions

The study was comprised of modifying the energy balance model of climate change. Two types of modifications were suggested. The model of ordinary differential equations models was extended to the model of partial differential equations having convections and $q$-diffusions effects. The numerical schemes were constructed, which were explicit and second-order accurate. Applying a numerical scheme for the ODEs and PDEs models and the graphs represented some results. The stability condition was given for the linear time-dependent differential equations. The impact of the parameter of $q$ on both types of temperatures was also shown, and it was concluded that global mean surface and deep oceanic temperatures had increasing behaviors over temporal variables. Also, it was observed that the deep ocean temperature had both increasing and decreasing behaviors on spatial variable $x$. Global mean surface and deep ocean temperatures escalated on spatial and temporal variables by the rising coefficient of linear and non-linear heat source terms. The presented numerical schemes can be further applied to different time-dependent problems having integer-order and $qth$ order spatial derivative terms. Also, the proposed numerical scheme can be further applied in epidemiological disease models and other problems ^[33,34,35] in fractional calculus.

Acknowledgments

The authors wish to express their gratitude to Prince Sultan University (PSU) for facilitating the publication of this article through the Theoretical and Applied Sciences Lab. This work was supported by the research grants Seed project; PSU; SEED-2022-CHS-100.

Conflict of interest

The authors declare no conflicts of interest to report regarding the present study.