Comparative Study of Finite Difference and Finite Element Methods for Solving 2D Parabolic Convection–Diffusion–Reaction Equations With Variable Coefficients

Abstract

This study investigates the numerical solution of two‐dimensional parabolic convection–diffusion–reaction (CDR) equations with variable coefficients using the finite difference method (FDM) and the finite element method (FEM). The FDM employs central differences for spatial discretization and the implicit Euler method for time integration, while the FEM uses the Galerkin approach with rectangular elements and three‐point Gauss–Legendre quadrature for spatial integrals, followed by implicit Euler discretization in time. Three test problems are considered to compare accuracy and efficiency. For small diffusion coefficients, the FEM provides higher accuracy, whereas for larger diffusion coefficients, both methods deliver nearly identical accuracy. Despite its improved accuracy in certain cases, the FEM typically involves a higher computational cost than the FDM. Based on the results, the study recommends the use of FEM for problems with boundary or interior layers.

1. Introduction

Convection–diffusion–reaction (CDR) equations are extensively used to model a wide range of physical and engineering phenomena, including elastic deformation in solid mechanics, fluid flow in porous media, heat and mass transfer, and chemical reaction processes. For example, in elasticity theory, CDR equations arise in the modeling of thermoplastic materials [1] and in analyzing stress distribution under diffusive influences [2, 3]. They are also fundamental in fluid flow applications, such as pollutant transport and groundwater contamination [4, 5], and in chemical reaction modeling, including catalytic reactors and combustion processes [6, 7]. In biological systems, CDR equations are employed to describe population dynamics [8], drug delivery and transport mechanisms [9–11], and tumor growth modeling [12].

Considering the vital role that CDR equations play in modeling diverse physical, chemical, and biological processes, obtaining accurate solutions is essential for properly describing the complex phenomena they represent. Due to the presence of variable coefficients, accurate numerical methods are essential to effectively capture the intricate behaviors associated with CDR processes. Without accurate solutions, predictions can be misleading, potentially compromising the design, analysis, and control of systems in engineering and science. Therefore, ensuring high-quality numerical approximations is fundamental to advancing our understanding and practical application of the processes governed by CDR equations.

The finite difference method (FDM) has been widely employed to solve CDR equations due to its simplicity and ease of implementation. Ndou et al. [13] developed an enhanced higher-order unconditionally positive FDM for solving advection–diffusion–reaction equations, demonstrating improved accuracy and stability while preserving the positivity of the numerical solution. Their method proved effective in both linear and nonlinear cases, particularly for problems with sharp gradients. Zhang et al. [14] developed adaptive FDMs specifically designed for CDR equations with spatially varying coefficients. Their approach dynamically adjusted the mesh in response to the solution behavior, thereby capturing boundary layers and steep fronts more accurately than uniform grids. Kumar and Singh [15] applied an implicit FDM to nonlinear CDR equations, emphasizing its ability to handle the stiff nature of the equations and suppress numerical oscillations in convection-dominated regimes. Furthermore, Lee and Park [16] introduced a robust finite difference scheme incorporating a flux-limiter strategy for high-gradient CDR systems, showing reduced numerical diffusion and improved accuracy near reaction fronts. Similarly, Li et al. [17] proposed a fourth-order compact finite difference scheme for solving the two-dimensional convection–diffusion equation, aimed at modeling groundwater pollution. The scheme achieves fourth-order spatial accuracy and second-order temporal accuracy, effectively capturing the transport behavior of pollutants in porous media.

Further developments have expanded the scope of FDM for more complex and realistic CDR models. Aniley and Duressa [18] and Tassaddiq et al. [19] developed advanced schemes specifically for fractional and classical CDR problems, demonstrating greater stability and accuracy in capturing boundary layer behavior. Li et al. [20] and Ferreira and Pena [21] proposed high-order and bound-preserving FDM schemes that effectively addressed nonlinearity while ensuring physically meaningful solutions. Tsega [22] solved two-dimensional nonlinear unsteady advection–diffusion–reaction equations with variable coefficients using a finite difference scheme designed to handle the challenges of nonlinearity, unsteadiness, and spatial variability. Gao and Liu [23] as well as Li et al. [24] introduced compact and efficient finite difference schemes applicable to multidimensional and time-fractional CDR systems, with particular relevance to environmental pollution and reactive transport modeling. Kaya [25] tackled the numerical difficulties associated with small diffusion and high Peclet numbers in multidimensional domains by implementing a specialized grid and operator-splitting approach. Collectively, these contributions highlight the adaptability and effectiveness of FDM in capturing the diverse behaviors exhibited by CDR systems across various scientific and engineering contexts.

Several authors have also used the finite element method (FEM) to solve CDR equations, which offers greater flexibility in handling complex geometries and variable coefficients inherent in these problems. Sun [26] developed a second-order characteristic mixed FEM that combines standard mixed FEM for diffusion with a characteristic approach for convection, improving accuracy in transient CDR simulations. Mekuria and Rao [27] employed an adaptive FEM with streamline diffusion stabilization and a posteriori error estimation to effectively capture boundary layer phenomena in steady-state 2D convection–diffusion problems. Codina [28] compared various FEM stabilization techniques, demonstrating enhanced stability and accuracy for convection-dominated flows. Deuring and Eymard [29] proposed a combined FEM and finite volume method using Crouzeix–Raviart elements with upwind barycentric finite volumes, ensuring L2-stability for CDR equations under mixed boundary conditions. Chowdhury and Kumar [30] introduced a subgrid multiscale stabilized FEM with rigorously derived stabilization parameters to improve reliability in solving CDR equations with variable coefficients. Wang and Hu [31] analyzed two-grid characteristic FEMs, establishing unconditional optimal error estimates and demonstrating computational efficiency in semilinear convection–diffusion problems. These contributions highlight the robustness and adaptability of FEM approaches in accurately resolving complex CDR phenomena across a range of applications. Sing et al. [32] applied a weak Galerkin finite element approach to the one-dimensional unsteady convection–diffusion equation with a nonlinear reaction term, employing a nonuniform mesh. They verified the convergence of the method specifically on layer-adapted meshes.

Several researchers have conducted comparative studies between the FDM and the FEM to evaluate their performance in solving CDR equations. For instance, Sun [26] examined the efficiency and accuracy of both methods for two-dimensional transient CDR problems and found that FEM offers superior flexibility in handling complex geometries, while FDM remains computationally efficient for regular domains. Codina [28] highlighted that FEM, especially when stabilized using streamline upwind Petrov–Galerkin (SUPG) techniques, performs better in convection-dominated problems by suppressing spurious oscillations, which FDM may struggle to eliminate without fine meshing. Similarly, Deuring and Eymard [29] discussed the integration of the FEM and finite volume method to exploit the structural advantages of both techniques, demonstrating that FEM provides better stability properties for problems with mixed boundary conditions. In contrast, Kaya [25] showed that properly designed FDM schemes can efficiently address multidimensional problems with small diffusion parameters when paired with tailored grid refinement strategies. Collectively, these studies suggest that selecting between FDM and FEM should be guided by factors such as the geometry of the domain, the nature of the boundary conditions, and the desired level of accuracy. FEM is often preferred for problems involving complex geometries, while FDM remains advantageous for its simplicity and computational efficiency in structured domains.

The objective of this study is to provide a systematic comparison between the FDM and the FEM for solving two-dimensional CDR equations in which the convection and reaction coefficients, as well as the source terms, vary with both spatial coordinates and time. Unlike most existing studies that primarily focus on one-dimensional cases or assume constant coefficients, this work addresses variable-coefficient, time-dependent problems that better reflect realistic applications. To ensure reliability, the analysis includes validation against benchmark problems with known analytical or reference solutions. The study further presents detailed convergence and efficiency results by reporting error norms, observed convergence rates, and computational cost, thereby providing a balanced assessment of the accuracy and performance of both methods. In addition, the strengths and limitations of each method are explicitly discussed, highlighting conditions under which one approach may be preferable to the other. The findings contribute a practical and scientifically grounded resource for researchers and practitioners seeking guidance on the selection of numerical methods for complex CDR problems.

2. Description of the Problem

This work addresses the numerical solution of the two-dimensional CDR equation defined as follows: 1 $\begin{array}{}{u}_t-D\left(u_{xx}+u_{yy}\right)+p\left(x,y,t\right)u_x+q\left(x,y,t\right)u_y+r\left(x,y,t\right)u=f\left(x,y,t\right),\left(x,y\right)\in \mathrm{\Omega },t\in \left(0,T\right],\end{array}$ where u(x, y, t) represents the unknown variable of interest. The coefficient D is a positive constant that characterizes the diffusion process. The functions p(x, y, t) and q(x, y, t) are the convection coefficients in the x- and y-directions, respectively, while r(x, y, t) denotes the reaction coefficient. The function f(x, y, t) is the source term, which indicates the amount of quantity being added to or removed from the system at the spatial location (x, y) and time t. In the computations, the spatial domain is $\mathrm{\Omega }=\left[a,b\right]\times \left[c,d\right]\subset {\mathbb{R}}^2$, which is bounded with a smooth boundary ∂Ω, and the time interval considered is [0, T]. To fully specify the problem, Dirichlet boundary conditions are imposed as 2 $\begin{array}{}u\left(x,y,t\right)=\varphi \left(x,y,t\right),\left(x,y\right)\in \mathrm{\partial }\mathrm{\Omega },t\in \left(0,T\right),\end{array}$ along with the initial condition 3 $\begin{array}{}u\left(x,y,0\right)=g\left(x,y\right),\left(x,y\right)\in \mathrm{\Omega }.\end{array}$

The p(x, y, t), q(x, y, t), r(x, y, t) and the source term f(x, y, t) are assumed to be positive, smooth, and bounded functions. Under these assumptions, the existence and uniqueness of the solution to the CDR equation are guaranteed [33].

This problem formulation captures a wide range of transport phenomena in which convection, diffusion, and reaction mechanisms interact under spatially and temporally varying conditions. As such, it is essential to develop and assess accurate and robust numerical methods capable of simulating these processes effectively in various scientific and engineering applications.

3. Numerical Methods

In this section, we present the numerical techniques used to approximate the solution of the two-dimensional CDR equation introduced in Section 2. For each method, we detail the spatial and temporal discretization procedures, derive the fully discrete iterative schemes, and discuss the stability and convergence properties.

3.1. FDM

In this study, FDM is applied to discretize the two-dimensional CDR equation on a uniform grid. The method involves approximating derivatives by difference quotients, converting the continuous problem into a system of algebraic equations that can be solved iteratively.

The first step in approximating the derivatives is to discretize the spatial domain. This involves creating a uniform Cartesian grid by dividing the interval [a, b] into N_x equal subintervals of length Δx, and the interval [c, d] into N_y equal subintervals of length Δy. The grid points are x_i = a + iΔx for i = 0, 1, …, N_x, and y_j = c + jΔy for j = 0, 1, …, N_y. This discretization allows us to approximate spatial derivatives using finite difference formulas at each grid point. For the temporal derivative, we discretize the time interval [0 ,  T] into steps of size Δt = T/N_t, with time levels t_n = nΔt, n = 0, 1, …, N_t.

Using the central difference approximation for the spatial derivatives and the implicit Euler method for the time derivative, equation (1) can be discretized at the grid point (x_i, y_j) and time level t_n+1 as follows: 4 $\begin{array}{c}\frac{u_{i,j}^{n+1}-u_{i,j}^n}{\mathrm{\Delta }t}-D\left(\frac{u_{i+1,j}^{n+1}-2u_{i,j}^{n+1}+u_{i-1,j}^{n+1}}{\mathrm{\Delta }{x}^2}+\frac{u_{i,j+1}^{n+1}-2u_{i,j}^{n+1}+u_{i,j-1}^{n+1}}{\mathrm{\Delta }{y}^2}\right)+p_{i,j}^{n+1}\mathrm{·}\frac{u_{i+1,j}^{n+1}-u_{i-1,j}^{n+1}}{2\mathrm{\Delta }x}+q_{i,j}^{n+1}\mathrm{·}\frac{u_{i,j+1}^{n+1}-u_{i,j-1}^{n+1}}{2\mathrm{\Delta }y}+r_{i,j}^{n+1}{u}_{i,j}^{n+1}=f_{i,j}^{n+1}.\end{array}$

This scheme is classified as fully implicit since all spatial terms and associated coefficients are evaluated at the upcoming time level. Rearranging equation (4) and assuming a uniform spatial grid in both directions, that is, Δx = Δy = h, with N_x = N_x = N, we obtain the following five-point stencil from the fully discrete scheme: 5 $\begin{array}{}{C}_1{u}_{i-1,j}^{n+1}+C_2{u}_{i,j}^{n+1}+C_3{u}_{i+1,j}^{n+1}+C_4{u}_{i,j-1}^{n+1}+C_5{u}_{i,j+1}^{n+1}=u_{i,j}^n+\mathrm{\Delta }t{f}_{i,j}^{n+1},\end{array}$ where $C_1=\mathrm{\Delta }t\left(-\left(D/h^2\right)-p_{i,j}^{n+1}/\left(2h\right)\right)$, $C_2=14+\mathrm{\Delta }t\left(D/h^2\right)+\mathrm{\Delta }t{r}_{i,j}^{n+1}$, $C_3=\mathrm{\Delta }t\left(-\left(D/h^2\right)+p_{i,j}^{n+1}/\left(2h\right)\right)$, $C_4=\mathrm{\Delta }t\left(-\left(D/h^2\right)-q_{i,j}^{n+1}/\left(2h\right)\right)$, and $C_5=\mathrm{\Delta }t\left(-\left(D/h^2\right)+q_{i,j}^{n+1}/\left(2h\right)\right)$.

This discrete formulation leads to a system of linear equations at each time level, which can be efficiently solved using iterative or direct methods.

3.1.1. Stability of the Iterative Scheme

The fully discrete finite difference scheme presented in (5) is based on the implicit Euler method for time discretization and central difference approximations for the spatial derivatives. This combination yields a fully implicit scheme, which is generally known for its unconditional stability in the context of linear parabolic partial differential equations, including CDR problems [34, 35].

Assuming the coefficients p(x, y, t), q(x, y, t), and r(x, y, t) and the source term f(x, y, t) are smooth and bounded, and under a sufficiently fine spatial grid, the stability of the scheme is ensured by the implicit treatment of the time derivative [36]. In particular, the diffusion terms benefit from the implicit time integration, which damps numerical oscillations and allows for larger time steps without compromising stability. The convection terms, discretized using central differences, may introduce spurious oscillations if the convection is dominant (i.e., high Peclet number). However, the implicit time integration provides some stabilization even in such cases. The reaction term contributes positively to the diagonal of the matrix, which also helps stabilize the solution. Mathematically, the stability of the scheme can be analyzed using matrix spectral analysis, where the resulting linear system at each time step forms an M-matrix under appropriate conditions, guaranteeing that the numerical solution does not grow unbounded in time [37].

3.1.2. Convergence of the Iterative Scheme

As demonstrated in standard references [36], the convergence of the scheme results from its consistency and stability, in accordance with the Lax–Richtmyer equivalence theorem. Since the finite difference scheme is consistent with the original PDE (i.e., the truncation error tends to zero as h ⟶ 0 and Δt ⟶ 0 and it is stable), we conclude that it is convergent. For the specific case of CDR problems, additional theoretical insights can be found in [38]. Numerical experiments involving smooth solutions demonstrate that the scheme attains first-order temporal accuracy and second-order spatial accuracy, in agreement with the accuracy expected from the employed discretization techniques.

3.2. FEM

In this study, rectangular elements are employed for the spatial discretization of the domain using the FEM. This choice facilitates a direct comparison between the FEM and FDM, as both utilize structured meshes aligned with Cartesian coordinates. The Galerkin method of the weight residual method [39–42] is used in the formulation of element equations in this paper.

3.2.1. Formulation of Element Equations

Consider a bilinear rectangular element e having local node numbers 1, 2, 3, and 4 and vertices (x₁,  y₁), (x₂,  y₂), (x₃,  y₃), and (x₄,  y₄) as shown in Figure 1. Suppose the values of u at Nodes 1, 2, 3, and 4 are denoted as $u_1^e$, $u_2^e$, $u_3^e$, and $u_4^e$, respectively. Local coordinates (ξ, η) are used to simplify interpolation and integration over the element. The element domain is converted to local coordinates by mapping it onto a standard reference square with coordinates ranging from −1 to 1 in both directions, forming the domain [−1, 1] × [−1,1]. These local coordinates are given by 6 $\begin{array}{c}\xi =\frac{2\left(x-x_c\right)}{\mathcal{l}x},\eta =\frac{2\left(y-y_c\right)}{\mathcal{l}y},\end{array}$ where (x_c,  y_c) are the coordinates of the centroid of the element, and $\mathcal{l}x$ and $\mathcal{l}y$ represent its length in x- and y-directions. The following shape functions N_i(ξ, η) are used to interpolate the solution over the element [40]: 7 $\begin{array}{c}\begin{array}{c}{N}_1\left(\xi ,\eta \right)=\frac{1}{4}\left(1-\xi \right)\left(1-\eta \right),\\ N_2\left(\xi ,\eta \right)=\frac{1}{4}\left(1+\xi \right)\left(1-\eta \right),\end{array}\begin{array}{c}{N}_3\left(\xi ,\eta \right)=\frac{1}{4}\left(1+\xi \right)\left(1+\eta \right),\\ N_4\left(\xi ,\eta \right)=\frac{1}{4}\left(1-\xi \right)\left(1+\eta \right).\end{array}\end{array}$

[IMAGE OMITTED. SEE PDF]

The spatial coordinates (x, y) in the global coordinate system and (ξ, η) ∈ [−1, 1] × [−1, 1] in the local coordinate system are related by 8 $\begin{array}{}x=x_1{N}_1\left(\xi ,\eta \right)+x_2{N}_2\left(\xi ,\eta \right)+x_3{N}_3\left(\xi ,\eta \right)+x_4{N}_4\left(\xi ,\eta \right),\end{array}$ 9 $\begin{array}{}\eta =y_1{N}_1\left(\xi ,\eta \right)+y_2{N}_2\left(\xi ,\eta \right)+y_3{N}_3\left(\xi ,\eta \right)+y_4{N}_4\left(\xi ,\eta \right).\end{array}$

Assume that within element e, u(x, y, t) is approximated using a combination of shape functions and the values at the element’s nodes as 10 $\begin{array}{}U=N_1^e{u}_1^e+N_2^e{u}_2^e+N_3^e{u}_3^e+N_4^e{u}_4^e.\end{array}$

Implementing the Galerkin FEM on equation (1) using the local coordinates with $\mathcal{l}x=\mathcal{l}y=h$ and applying the weak form of the diffusion equation yield 11 $\begin{array}{c}{\int }_{-1}^1{\int }_{-1}^1\left[N_i\frac{\mathrm{\partial }U}{\mathrm{\partial }t}+D\left(\frac{4}{h^2}\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }\xi }\frac{\mathrm{\partial }U}{\mathrm{\partial }\xi }+\frac{4}{h^2}\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }\eta }\frac{\mathrm{\partial }U}{\mathrm{\partial }\eta }\right)+\frac{2}{h}p{N}_i\frac{\mathrm{\partial }U}{\mathrm{\partial }\xi }+\frac{2}{h}q{N}_i\frac{\mathrm{\partial }U}{\mathrm{\partial }\eta }+r{N}_{i}U\right]\left|J\right|d\xi d\eta ={\int }_{-1}^1{\int }_{-1}^1{N}_{i}f\left|J\right|d\xi d\eta ,\end{array}$ where N_i is the shape function associated with node i, U is the approximation of u as indicated in (10), and |J| is the determinant of the Jacobian matrix of the transformation from the global coordinates (x, y) to the local coordinates (ξ, η). From (8) and (9), the matrix J can be obtained as 12 $\begin{array}{}J=\left[\begin{array}{cc}\frac{\mathrm{\partial }x}{\mathrm{\partial }\xi }& \frac{\mathrm{\partial }y}{\mathrm{\partial }\xi }\\ \frac{\mathrm{\partial }x}{\mathrm{\partial }\eta }& \frac{\mathrm{\partial }y}{\mathrm{\partial }\eta }\end{array}\right]=\left[\begin{array}{cccc}\frac{\mathrm{\partial }{N}_1}{\mathrm{\partial }\xi }& \frac{\mathrm{\partial }{N}_2}{\mathrm{\partial }\xi }& \frac{\mathrm{\partial }{N}_3}{\mathrm{\partial }\xi }& \frac{\mathrm{\partial }{N}_4}{\mathrm{\partial }\xi }\\ \frac{\mathrm{\partial }{N}_1}{\mathrm{\partial }\eta }& \frac{\mathrm{\partial }{N}_2}{\mathrm{\partial }\eta }& \frac{\mathrm{\partial }{N}_3}{\mathrm{\partial }\eta }& \frac{\mathrm{\partial }{N}_4}{\mathrm{\partial }\eta }\end{array}\right]\left[\begin{array}{cc}{x}_1& y_1\\ x_2& y_2\\ x_3& y_3\\ x_4& y_4\end{array}\right].\end{array}$

Equation (11) leads to a semidiscrete equation of the form 13 $\begin{array}{}{M}_e\frac{du}{dt}+A_{e}u=B_e,\end{array}$ where $u={\left[u_1^e{u}_2^e{u}_3^e{u}_4^e\right]}^T$, and 14 $\begin{array}{}{M}_{i,j}^e={\int }_{-1}^1{\int }_{-1}^1{N}_i{N}_j\left|J\right|d\xi d\eta ,\end{array}$ 15 $\begin{array}{}{A}_{i,j}^e={\int }_{-1}^1{\int }_{-1}^{-1}\left[D\left(\frac{4}{h^2}\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }\xi }\frac{\mathrm{\partial }{N}_j}{\mathrm{\partial }\xi }+\frac{4}{h^2}\frac{\mathrm{\partial }{N}_i}{\mathrm{\partial }\eta }\frac{\mathrm{\partial }{N}_j}{\mathrm{\partial }\eta }\right)+\frac{2}{h}p{N}_i\frac{\mathrm{\partial }{N}_j}{\mathrm{\partial }\xi }+\frac{2}{h}q{N}_i\frac{\mathrm{\partial }{N}_j}{\mathrm{\partial }\eta }+r{N}_i{N}_j\right]\left|J\right|d\xi d\eta ,\end{array}$ 16 $\begin{array}{}{B}_i^j={\int }_{-1}^1{\int }_{-1}^1{N}_{i}f\left|J\right|d\xi d\eta ,\end{array}$ for i, j = 1, 2, 3, 4.

The integrals (14)–(16) are approximated using the three-point Gauss–Legendre integration formula [41] as 17 $\begin{array}{}{\int }_{-1}^1{\int }_{-1}^1\mathrm{\Psi }\left(\xi ,\eta \right)d\xi d\eta \approx {\displaystyle \sum }_{i=1}^3{\displaystyle \sum }_{j=1}^3{w}_i{w}_j\mathrm{\Psi }\left({\xi }_i,{\eta }_j\right),\end{array}$ where ${\xi }_1=-\sqrt{35/}$, ξ₂ = 0, ${\xi }_3=\sqrt{35/}$, w₁ = (5/9), w₂ = (8/9), w₃ = 5/9.

3.2.2. Assembly of Element Equations and Temporal Discretization

After computing the element-level matrices M_e and A_e, and vector B_e using the shape functions and numerical integration, the next step is to assemble them into the global system. This process involves summing the contributions of each element to the global matrices and vector according to the connectivity of the mesh. After assembling the element matrices using standard finite element procedures as outlined in [41], we obtain the global system of equations as 18 $\begin{array}{}M\frac{du}{dt}+Au=B,\end{array}$ where M and A are (N + 1) × (N + 1) matrices and B is the (N + 1) × 1 vector. M is a symmetric positive definite matrix.

Using the implicit Euler method and the time discretization approach described in Section 3.1, the assembled equation (18) can be discretized in time as follows: 19 $\begin{array}{}M\frac{u^{n+1}-u^n}{\mathrm{\Delta }t}+A{u}^{n+1}=B^{n+1}.\end{array}$

Rearranging equation (19) yields 20 $\begin{array}{}\left(M+\mathrm{\Delta }tA\right)u^{n+1}=A{u}^n+\mathrm{\Delta }t{B}^{n+1}.\end{array}$

Using the boundary and initial conditions given in equations (2) and (3), equation (20) is solved using a suitable method for solving systems of equations. The matrix (M + Δt A) is generally sparse under standard assumptions on the coefficients and finite element basis functions, which makes the system well-suited for solution via direct solvers or iterative methods.

3.2.3. Stability of the Scheme

Let us use the energy method for stability analysis of scheme (20) [43]. Consider the positive symmetric definite matrix M and the discrete vector B^{(n + 1)} at time t^{(n + 1)} in the scheme. The M-inner product and the induced energy norm as B^{(n + 1)} are mentioned in equation (22). 21 $\begin{array}{}\left(X,Y\right)=X^{T}MY\text{\hspace{0.17em}and\hspace{0.17em}}\parallel X{\parallel }_M=\sqrt{X^{T}MX,}\end{array}$ for vectors X and Y. Here, the subscript T denotes the transpose. Multiplying both sides of (20) by u^n+1T, we get 22 $\begin{array}{}{u}^{n+1T}M{u}^{n+1}+\mathrm{\Delta }t{u}^{n+1T}A{u}^{n+1}=u^{n+1T}M{u}^n+\mathrm{\Delta }t{u}^{n+1T}{B}^{n+1}.\end{array}$

Using the polarization identity for the vectors uⁿ and uⁿ⁺¹ in the M-inner product space, we get 23 $\begin{array}{}{\left(u^{n+1}-u^n\right)}^{T}M\left(u^{n+1}-u^n\right)=u^{n+1T}M{u}^{n+1}+u^{nT}M{u}^n-2u^{n+1T}M{u}^n.\end{array}$

Solving for u^n+1TMuⁿ gives 24 $\begin{array}{}{u}^{n+1T}M{u}^n=\frac{1}{2}\left(u^{n+1T}M{u}^{n+1}+u^{nT}M{u}^n-{\left(u^{n+1}-u^n\right)}^{T}M\left(u^{n+1}-u^n\right)\right),\end{array}$ or 25 $\begin{array}{}{u}^{n+1T}M{u}^n=\frac{1}{2}\left(\parallel u^{n+1}{\parallel }_M^2+\parallel u^n{\parallel }_M^2-\parallel u^{n+1}-u^n{\parallel }_M^2\right).\end{array}$

Substituting (25) into (22) gives 26 $\begin{array}{}{u}^{n+1T}M{u}^{n+1}+\mathrm{\Delta }t{u}^{n+1T}A{u}^{n+1}=\frac{1}{2}\left(\parallel u^{n+1}{\parallel }_M^2+\parallel u^n{\parallel }_M^2-\parallel u^{n+1}-u^n{\parallel }_M^2\right)+\mathrm{\Delta }t{u}^{n+1T}{B}^{n+1}.\end{array}$

Rearranging (26), we get 27 $\begin{array}{}\frac{1}{2}\left(\parallel u^{n+1}{\parallel }_M^2-\parallel u^n{\parallel }_M^2+\parallel u^{n+1}-u^n{\parallel }_M^2\right)+\mathrm{\Delta }t{u}^{n+1T}A{u}^{n+1}=\frac{1}{2}\mathrm{\Delta }t{u}^{n+1T}{B}^{n+1}.\end{array}$

From (27), we have 28 $\begin{array}{}\frac{1}{2}\left(\parallel u^{n+1}{\parallel }_M^2-\parallel u^n{\parallel }_M^2\right)\le \frac{1}{2}\mathrm{\Delta }t{u}^{n+1T}{B}^{n+1}.\end{array}$

Using the Cauchy–Schwarz inequality for the SPD matrix gives 29 $\begin{array}{}{u}^{n+1T}{B}^{n+1}\le {‖u^{n+1}‖}_M{‖B^{n+1}‖}_{M^{-1}}.\end{array}$

From Young’s inequality, we obtain 30 $\begin{array}{}{‖u^{n+1}‖}_M{‖B^{n+1}‖}_{M^{-1}}\le \frac{1}{2}{‖u^{n+1}‖}_M^2+\frac{1}{2}{‖B^{n+1}‖}_{M^{-1}}^2.\end{array}$

From (29) and (30), we get 31 $\begin{array}{}{u}^{n+1T}{B}^{n+1}\le \frac{1}{2}{‖u^{n+1}‖}_M^2+\frac{1}{2}{‖B^{n+1}‖}_{M^{-1}}^2.\end{array}$

Multiplying both sides of (31) by Δt yields 32 $\begin{array}{}{u}^{n+1T}{B}^{n+1}\le \frac{\mathrm{\Delta }t}{2}{‖u^{n+1}‖}_M^2+\frac{\mathrm{\Delta }t}{2}{‖B^{n+1}‖}_{M^{-1}}^2.\end{array}$

From (28) and (32), we get 33 $\begin{array}{}\frac{1}{2}\left({‖u^{n+1}‖}_M^2-{‖u^n‖}_M^2\right)\le \frac{\mathrm{\Delta }t}{2}{‖u^{n+1}‖}_M^2+\frac{\mathrm{\Delta }t}{2}{‖B^{n+1}‖}_{M^{-1}}^2,\end{array}$ or 34 $\begin{array}{}{‖u^{n+1}‖}_M^2\le \frac{1}{1-\mathrm{\Delta }t}\left({‖u^n‖}_M^2+\mathrm{\Delta }t{‖B^{n+1}‖}_{M^{-1}}^2\right),\text{\hspace{0.17em}for\hspace{0.17em}}\mathrm{\Delta }t<1.\end{array}$

The global bound over all time steps (iterating over n = 0, 1, 2, …, N − 1) is 35 $\begin{array}{}{‖u^n‖}_M^2\le {‖u^0‖}_M^2+\mathrm{\Delta }t{\displaystyle \sum }_{k=1}^n{‖B^k‖}_{M^{-1}}^2.\end{array}$

Hence, the backward Euler scheme is stable in the M-norm.

3.2.4. Convergence of the Scheme

Let u(x, y, t) denote the exact solution of the CDR equation on the spatial domain Ω at time t such that u(x, y, t) ∈ H²(Ω) in space u(x, y, t) ∈ C¹([0, T]) in time. Let u_h(x, y, t) denote the finite element approximation in space at time t, obtained by projecting u(x, y, t) onto the finite-dimensional space V_h, spanned by bilinear rectangular basis functions 36 $\begin{array}{}{u}_h\left(x,y,t\right)\in V_h\subset H^1\left(\mathrm{\Omega }\right),u_h\left(x,y,t\right)\approx u\left(x,y,t\right),\end{array}$ where h is the mesh size of the spatial discretization.

We denoted the fully discrete solution at the time level tⁿ⁺¹ = (n + 1)Δt by Uⁿ⁺¹(x, y), obtained after applying the backward Euler method in time. 37 $\begin{array}{}{{U^n}^{+}}^1\left(x,y\right)\approx u_h\left(x,y,{{t^n}^{+}}^1\right)\approx u\left(x,y,{{t^n}^{+}}^1\right).\end{array}$

Under the regularity assumptions above, the standard finite element theory for parabolic problems [37] ensures the following spatial accuracy: 38 $\begin{array}{}\parallel u\left(x,y,t\right)-u_h\left(x,y,t\right){\parallel }_{L^2\left(\mathrm{\Omega }\right)}\le C{h}^2,\parallel u\left(x,y,t\right)-u_h\left(x,y,t\right){\parallel }_{H^1\left(\mathrm{\Omega }\right)}\le Ch,\end{array}$ where C is a constant independent of h.

For time discretization, the backward Euler method provides the first-order temporal accuracy as 39 $\begin{array}{}\parallel u^h\left(x,y,t^{n+1}\right)-U^{n+1}\left(x,y\right){\parallel }_{L^2\left(\mathrm{\Omega }\right)}\le C\mathrm{\Delta }t.\end{array}$

Combining spatial and temporal discretizations, the global error at time tⁿ⁺¹ satisfies 40 $\begin{array}{}\parallel u\left(x,y,t^{n+1}\right)-U^{n+1}\left(x,y\right){\parallel }_{L^2\left(\mathrm{\Omega }\right)}\le C\left(h^2+\mathrm{\Delta }t\right),\end{array}$ where the constant C depends on the solution regularity and the final time but is independent of h and Δt. Thus, under the smoothness assumptions on the coefficient and source functions, and with the use of a structured mesh, the fully discrete scheme is convergent, achieving second-order accuracy in space and first-order accuracy in time.

4. Numerical Results and Discussion

In this section, numerical examples with known exact solutions are presented to evaluate the performance of the FDM and FEM. All computations are carried out using programs developed in MATLAB 2015a on a PC (Intel Core i7-7500U CPU, 8 GB RAM). Computational results are displayed using surface visualizations, plotted curves, and tables. To assess the accuracy of the methods, the errors in the computed solutions are analyzed at various levels of mesh refinement within the computational domain. The same number of subintervals, N, is used in the x- and y-directions for the spatial meshes.

The maximum absolute error at all mesh points for time t = t_n is calculated by 41 $\begin{array}{}{E}^N=\underset{0\le i,\text{j}\le N}{\max }\left|u_{i,j}^n-u\left(x_i,y_j,t_n\right)\right|,\end{array}$ where $u_{i,j}^n$ and u(x_i,  y_j, t_n) represent the numerical and exact solutions of u, respectively.

In all examples, a square domain Ω = [0, 1] × [0, 1] is considered, and the initial condition function g(x, y), boundary condition function ϕ(x, y, t) , and source term f(x, y, t) are derived from the exact solution.

$$

Statement

Example 1.

Consider the CDR equation in [44] with D = 1, p(x, y, t) = 2, q(x, y, t) = 3, and r(x, y, t) = 1 with the exact solution 42 $\begin{array}{}u\left(x,y,t\right)=16\sin \left(\pi t\right)x\left(1-x\right)y\left(1-y\right)\left(\frac{1}{2}+\frac{\arctan \left(2D^{-0.5}\left({0.25}^2-{\left(x-0.5\right)}^2-{\left(y-0.5\right)}^2\right)\right)}{\pi }\right).\end{array}$

Figure 2 illustrates the exact and numerical solutions at y = 0.5, based on the mesh presented in [44], serving as a benchmark for assessing the accuracy of the employed numerical methods. The maximum absolute errors for the FDM and the FEM are found to be 0.0091 and 0.00070897, respectively. In comparison, the dual reciprocity boundary element method (DRBEM) reported in [44] yielded a maximum absolute error of 0.0013. These results indicate that the FEM provides a more accurate solution than the method used in [44], while the FDM also demonstrates a reasonable level of accuracy. Figure 3 presents the surface plot of the solutions for N = 50, Δt = 0.001, at t = 0.5. The figures show that the numerical solutions closely match the exact solution.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

$$

Statement

Example 2.

Let us consider a CDR equation as described in (1) with p(x, y, t) = 1 + x² + y² + t², q(x, y, t) = 2 + xy + t², r(x, y, t) = 1 + x²y²t² and exact solution u(x, y, t) = (1 + sin(πx)cos(πy))e^−t. In the equation, the convection and reaction coefficients, as well as the source term, are functions of space and time.

Figure 4 presents the numerical (FDM and FEM) and exact solutions of Example 2 at time t = 1, with parameters D = 1, N = 32, and Δt = 0.001. The surface plots reveal a strong agreement between the numerical and exact solutions, and the absence of steep gradients in the solution profile can be attributed to the relatively large diffusion coefficient. Both numerical methods successfully capture the smooth behavior of the solution across the domain. For a more detailed view, Figure 5 displays the solution profiles at the mid-line, x = 0.5. The FDM and FEM solution curves follow the exact solution closely, and the differences among them are minimal. This further confirms the stability and accuracy of both schemes in the presence of dominant diffusion.

A quantitative comparison is provided in Table 1, which reports the maximum absolute errors for different spatial discretizations. As the mesh is refined (i.e., as N increases), the error decreases for both methods, demonstrating clear spatial convergence. Temporal convergence is illustrated in Table 2, where the maximum absolute errors are shown for various time step sizes with a fixed spatial resolution, N = 32. As the time step Δt decreases from 0.1 to 0.001, the error for both FDM and FEM decreases significantly, confirming the schemes’ accuracy in time.

Overall, both numerical methods demonstrate high accuracy and consistency with the exact solution. Due to the smoothness of the solution caused by the relatively high diffusion parameter D = 1, the differences between FDM and FEM are negligible. These results validate the effectiveness of both schemes for solving time-dependent CDR problems under diffusion-dominated conditions.

A quantitative comparison is presented in Table 3, which summarizes the maximum absolute errors of the FDM and FEM at various mesh sizes (N = 16, 32, and 64) and diffusion coefficients ranging from 1E1 to 1E − 4. The results clearly indicate that both methods exhibit improved accuracy as the mesh is refined. At relatively high diffusion levels, both FDM and FEM yield similar levels of accuracy across all mesh sizes. However, as the diffusion coefficient decreases, noticeable differences begin to emerge. For moderately low diffusion, the methods still perform comparably, though the FEM begins to show a slight advantage. The difference is more significant as the diffusion coefficient decreases further. In these cases, the FDM exhibits a marked decline in accuracy, whereas the FEM maintains a relatively small error, demonstrating greater robustness to decreasing diffusion.

For better visualization of these effects, Figure 6 (curve plot at x = 0.5) and Figure 7 (surface plot) display the exact solution alongside the numerical solutions obtained from FDM and FEM using N = 64 with a diffusion coefficient of 1E − 4. At very low diffusion levels, the FDM’s error increases substantially, resulting in oscillatory behaviors near the layer, while the FEM continues to produce reliable and accurate results. This indicates that the FEM is more suitable for problems dominated by convection or reaction, or those involving sharp gradients, where achieving numerical stability and accuracy is more challenging.

Table 4 compares the computational time and corresponding accuracy of the FDM and FEM for Example 2 under fixed mesh size and time step. The results reveal a notable contrast between the two methods in terms of efficiency and cost. The FDM generally demands considerably less computational time than the FEM. Across all evaluated time levels, the FDM completes the simulation much faster, highlighting its computational efficiency and simplicity of implementation. In contrast, the FEM incurs a substantially higher computational cost due to the complexity of assembling and solving systems involving global basis functions.

Figure 8 presents the surface plot of the absolute maximum errors of the FDM and FEM solutions for Example 2 at t = 1 with D = 1, N = 32, and Δt = 0.001. The results show that the error distributions of both methods are nearly identical, indicating that the finite difference and finite element approaches provide comparable accuracy under the given parameters. This close agreement demonstrates the consistency and reliability of both numerical schemes in approximating the exact solution.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Table 1 Comparison of maximum absolute errors for FDM and FEM in Example 2 at different mesh sizes, with D = 1 and Δt = 0.001, evaluated at t = 1.

N	Maximum absolute error
FDM	FEM
8	4.9413E − 3	4.9935E − 3
16	1.1218E − 3	1.2191E − 3
32	2.9660E − 4	2.9795E − 4
64	6.6573E − 5	6.6907E − 5

Table 2 Comparison of maximum absolute errors for FDM and FEM in Example 2 at different time step sizes at t = 0.5, with D = 1 and N = 32.

Δt	Maximum absolute error
FDM	FEM
0.1	2.0543E − 3	2.2835E − 3
0.05	1.0218E − 3	1.0541E − 3
0.01	2.8164E − 4	2.79475E − 4
0.005	2.0912E − 4	2.0596E − 4
0.001	1.4513E − 4	1.4134E − 4

Table 3 Comparison of maximum absolute errors of FDM and FEM at different mesh sizes and diffusion coefficients for Example 2 with Δt = 0.001 at t = 1.

D	N = 16	N = 32	N = 64
FDM	FEM	FDM	FEM	FDM	FEM
1E1	4.6722E − 4	4.6920E − 4	1.1689E − 4	1.1717E − 4	4.5646E − 5	2.8786E − 5
1E0	1.2118E − 3	1.2191E − 3	2.9660E − 4	2.9795E − 4	6.6573E − 5	6.6907E − 5
1E − 1	2.6172E − 3	2.6506E − 3	6.0260E − 4	6.0608E − 4	1.2272E − 4	1.2373E − 4
1E − 2	5.9335E − 3	5.7769E − 3	9.6554E − 4	9.6309E − 4	1.7834E − 4	1.7999E − 4
1E − 3	3.2358E − 1	1.6423E − 2	2.0693E − 2	1.8180E − 3	2.4744E − 4	2.5094E − 4
1E − 4	7.4354E − 1	2.8331E − 2	5.1247E − 1	5.5075E − 3	1.5518E − 1	5.9856E − 4

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Table 4 Computational time for the FDM and FEM for Example 2 with D = 1, N = 16, and Δt = 0.001.

Time (t)	FDM	FEM
CPU time (s)	Max.abs.error	CPU time (s)	Max.abs.error
0.1	0.4550	1.5289E − 2	41.3045	2.2439E − 2
0.5	2.3887	1.7115E − 3	230.1860	1.7153E − 3
1	4.7300	1.2118E − 3	443.6938	1.2191E − 3

[IMAGE OMITTED. SEE PDF]

$$

Statement

Example 3.

Consider the CDR equation in (1) with p(x, y, t) = sin(πx)cos(πy) + t, q(x, y, t) = cos(πx)sin(πy) + e^−t, and r(x, y, t) = e^x+y−t, which has an exact solution 43 $\begin{array}{}u\left(x,y,t\right)=100{\left(1-x\right)}^2{x}^{2}y\left(1-y\right){\left(12-y\right)}^2\tan h^2\left(2x+y-t\right).\end{array}$

Figure 9 displays the numerical and exact solutions of Example 3 at time t = 0.5, with diffusion coefficient D = 1, N = 64 spatial divisions, and time step Δt = 0.001. Specifically, Figure 10 shows the solution along the line x = 0.5 obtained using the FDM, the FEM, and the exact solution. The close agreement between the numerical and exact solutions confirms the accuracy of the numerical methods for a relatively large diffusion coefficient, as also observed in the solution of Example 2.

Similar to Table 3, for Example 2, Table 5 presents a comparison of the FDM and FEM solutions for Example 3 using different diffusion coefficients and mesh sizes. As the mesh is refined, the accuracy of both numerical methods improves consistently. However, as the diffusion coefficient decreases, the error increases, indicating the growing challenge in capturing sharp gradients. The solution obtained using the FDM exhibits slight oscillations, particularly at lower diffusion levels, whereas the FEM provides results that closely match the exact solution. This behavior is consistent with the observations reported in Example 2.

Figure 11 displays the surface plots of the numerical and exact solutions of Example 3 at time t = 0.5, with a diffusion coefficient D = 1E − 4, N = 64 spatial divisions, and time step Δt = 0.001. A cross-sectional view along x = 0.5 is shown in Figure 12. These figures clearly illustrate that the FDM solution is more sensitive to low diffusion coefficients, exhibiting noticeable deviations from the exact solution. In contrast, the FEM solution remains remarkably close to the exact solution, further demonstrating its robustness and superior accuracy under challenging conditions, consistent with the trends observed in Example 2.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Table 5 Comparison of maximum absolute errors of FDM and FEM at different mesh sizes and diffusion coefficients for Example 3 with Δt = 0.001 at t = 1.

D	N = 16	N = 32	N = 64
FDM	FEM	FDM	FEM	FDM	FEM
1E1	5.4322E − 3	2.2890E − 3	1.3572E − 3	5.6480E − 4	3.3917E − 4	1.4149E − 4
1E0	4.8316E − 3	2.4572E − 3	1.2057E − 3	6.0870E − 4	3.0130E − 4	1.5175E − 4
1E − 1	3.2418E − 3	4.3830E − 3	7.9084E − 4	1.0715E − 3	3.0329E − 4	2.6141E − 4
1E − 2	2.0024E − 2	1.3745E − 2	8.5025E − 3	2.7771E − 3	7.6867E − 3	6.5036E − 4
1E − 3	4.5462E − 1	2.4804E − 2	2.8174E − 2	5.9789E − 3	2.8636E − 2	1.3165E − 3
1E − 4	7.0429E − 1	3.3705E − 2	6.1147E − 1	7.7437E − 3	4.1855E − 2	2.1061E − 3

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

5. Conclusions

This study has presented a comparative analysis of the FDM and the FEM for solving CDR equations with spatially and temporally varying convection and reaction coefficients, and a constant diffusion coefficient. Refining the mesh leads to improved accuracy and convergence in both the numerical methods. For relatively large diffusion coefficients, both FDM and FEM produce stable and accurate solutions.

However, as the diffusion coefficient decreases, the FDM tends to generate more oscillatory solutions with higher errors, while the FEM maintains greater stability and accuracy. The results indicate that FEM consistently achieves lower maximum absolute errors than FDM under low-diffusivity conditions, yielding solutions that closely match the exact solution. Although FEM involves a higher computational cost, its superior accuracy makes it a more reliable choice for high-precision simulations.

Based on these findings, this study recommends using the FEM rather than the FDM, particularly when high accuracy is required or when addressing problems with low diffusivity. Adaptive techniques, such as mesh refinement, boundary layer–fitted meshes, and tailored operators, are also likely to be more effective when implemented within the FEM framework than within the FDM.

Data Availability Statement

No data were used to support this study.

Conflicts of Interest

The author declares no conflicts of interest.

Funding

No funding was received for this study.

References

1 Wang Y., A Multifield Coupled Thermo‐Chemo‐Mechanical Theory for the Degradation of Materials, International Journal for Numerical Methods in Engineering. (2024) 125, no. 3, 567–589.

2 Zhang L. and Li M., Stochastic Process Model for Interfacial Gap of Purely Normal Elastic Contact, Journal of the Mechanics and Physics of Solids. (2024) 165.

3 Doe J. and Smith A., Analysis of Stress Distribution in Composite Materials Using CDR Models, Materials Science and Engineering A. (2023) 812.

4 Kourakos G., Bastani M., and Harter T., Simulating Nonpoint Source Pollution Impacts in Groundwater: Three-Dimensional Advection–Dispersion Versus Quasi-3D Streamline Transport Approach, Hydrology. (2025) 12, no. 3, https://doi.org/10.3390/hydrology12030042.

5 Wu Y., Yu J., Huang Z. et al., Migration of Total Petroleum Hydrocarbon and Heavy Metal Contaminants in the Soil–Groundwater Interface of a Petrochemical Site Using Machine Learning: Impacts of Convection and Diffusion, RSC Advances. (2024) 14, no. 44, 32304–32313, https://doi.org/10.1039/d4ra06060a.

6 Chen H., Zhang Y., and Wang Z., Design and Performance Investigation of a Compact Catalytic Reactor Integrated With Heat Recuperator, Energies. (2023) 15, no. 2.

7 Luo Y. and Liu Y., Development of a Reaction–Diffusion–Convection Model and Its Application to NO Oxidation in Reactors, Industrial & Engineering Chemistry Research. (2024) 63, no. 32, 14381–14390, https://doi.org/10.1021/acs.iecr.4c01864.

8 Alessio B. M. and Gupta A., A Reaction-Diffusion-Chemotaxis Model for Human Population Dynamics Over Evolving Terrains, Mathematics. (2023) 11, no. 3.

9 Mendoza M., Drug Delivery Enhanced by Ultrasound: Mathematical Modeling and Simulation, Computers in Biology and Medicine. (2021) 134.

10 Zhou Y., Theoretical Model for Diffusion-Reaction Based Drug Delivery From a Multilayer Sphere, International Journal of Heat and Mass Transfer. (2021) 174.

11 Smith J., Drug Delivery in Biological Tissues: A Two-Layer Reaction-Diffusion Model, Journal of Biomedical Engineering. (2020) 42, no. 5, 123–135.

12 Zhang Y., Tumor Growth Prediction With Reaction-Diffusion and Hyperplastic Biomechanical Models, IEEE Transactions on Medical Imaging. (2015) 34, no. 3, 586–599.

13 Ndou N., Dlamini P., and Jacobs B. A., Solving the Advection Diffusion Reaction Equations by Using the Enhanced Higher-Order Unconditionally Positive Finite Difference Method, Mathematics. (2024) 12, no. 7, https://doi.org/10.3390/math12071009.

14 Zhang Y., Li X., and Wang Z., Adaptive Finite Difference Schemes for Convection–Diffusion–Reaction Equations With Variable Coefficients, Journal of Computational Physics. (2023) 475.

15 Kumar S. and Singh R., Implicit Finite Difference Approach for Solving Nonlinear Convection–Diffusion–Reaction Equations, Applied Mathematics and Computation. (2021) 398.

16 Lee J. and Park H., Hybrid High-Order Finite Difference Methods for Convection–Diffusion–Reaction Equations With Variable Coefficients, Computational Methods in Applied Mathematics. (2025) 25, no. 2, 345–362.

17 Li L., Jiang Z., and Yin Z., Fourth-Order Compact Finite Difference Method for Solving Two-Dimensional Convection–Diffusion Equation, Advances in Difference Equations. (2018) 2018, 234–24, https://doi.org/10.1186/s13662-018-1652-5, 2-s2.0-85049783143.

18 Aniley W. T. and Duressa G. F., Nonstandard Finite Difference Method for Time-Fractional Singularly Perturbed Convection–Diffusion Problems With a Delay in Time, Results in Applied Mathematics. (2024) 21, https://doi.org/10.1016/j.rinam.2024.100432.

19 Tassaddiq A., Alzahrani E., and Khan M., Numerical Study of Convection–Diffusion–Reaction Equations Via Finite Difference Method With a New Fractional Operator, Alexandria Engineering Journal. (2021) 60, no. 5, 4677–4688.

20 Li Y., Wang Y., and Zhao H., High-Order Finite Difference Methods for 2D Nonlinear Convection–Diffusion–Reaction Equations, Applied Mathematics Letters. (2020) 105.

21 Ferreira J. A. and Pena G., A Bound-Preserving Numerical Method for Convection–Diffusion–Reaction Equations, Mathematics and Computers in Simulation. (2022) 194, 146–161.

22 Tsega E. G., Numerical Solution of Two‐Dimensional Nonlinear Unsteady Advection‐Diffusion‐Reaction Equations With Variable Coefficients, International Journal of Mathematics and Mathematical Sciences. (2024) 2024, no. 1, 1–12, https://doi.org/10.1155/2024/5541066.

23 Gao S. and Liu D., A Compact Finite Difference Scheme for Solving Time-Fractional Convection–Diffusion–Reaction Equations in Two Dimensions, Mathematics. (2021) 9, no. 15.

24 Li Z., Tan Y., and Lin C., Finite Difference Analysis of a Multidimensional Convection–Diffusion–Reaction Equation With Nonhomogeneous Boundary Conditions, Computers & Mathematics With Applications. (2020) 79, no. 6, 1726–1739.

25 Kaya A., A Finite Difference Method for Convection–Diffusion–Reaction Equations With Small Diffusion Parameters in Multidimensional Domains, Journal of Computational and Applied Mathematics. (2021) 387.

26 Sun H., A Second-Order Characteristic Mixed Finite Element Method for Transient Convection–Diffusion–Reaction Problems, Applied Numerical Mathematics. (2023) 185, 1–18.

27 Mekuria G. T. and Rao K. S., Adaptive Finite Element Methods With Streamline Diffusion Stabilization for Convection–Diffusion Problems, Computers & Mathematics With Applications. (2022) 104, 320–338, https://doi.org/10.1016/j.camwa.2021.11.021.

28 Codina R., Comparison of Finite Element Stabilization Methods for Convection–Diffusion–Reaction Problems, Computer Methods in Applied Mechanics and Engineering. (2020) 344, 529–547.

29 Deuring P. and Eymard R., Finite Element–Finite Volume Discretizations for Convection–Diffusion–Reaction Equations: L2-Stability and Mixed Boundary Conditions, Numerische Mathematik. (2021) 147, no. 4, 987–1014.

30 Chowdhury R. and Kumar A., A Stabilized Subgrid Multiscale Finite Element Method for Convection–Diffusion–Reaction Equations With Variable Coefficients, Journal of Scientific Computing. (2024) 98, no. 2, 1–25.

31 Wang J. and Hu Y., Two-Grid Characteristic Finite Element Methods for Semilinear Convection–Diffusion Problems: Unconditional Optimal Error Estimates, Numerical Methods for Partial Differential Equations. (2023) 39, no. 3, 1234–1260.

32 Singh J., Kumar N., and Jiwari R., A High-Order Numerical Method for Semilinear Convection–Diffusion Equation on Bakhvalov-Type Mesh, Mathematical Methods in the Applied Sciences. (2025) 48, no. 10, 10586–10601, https://doi.org/10.1002/mma.10904.

33 John V., Finite Element Methods for Incompressible Flow Problems, Springer Series in Computational Mathematics. (2016) 51, Springer.

34 Smith G. D., Numerical Solution of Partial Differential Equations: Finite Difference Methods, 1985, 3rd edition, Oxford University Press.

35 Hundsdorfer W. and Verwer J. G., Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations, 2003, Springer.

36 Thomée V., Galerkin Finite Element Methods for Parabolic Problems, 2006, 2nd edition, Springer.

37 Quarteroni A., Sacco R., and Saleri F., Numerical Mathematics, 2007, 2nd edition, Springer.

38 Morton K. W. and Mayers D. F., Numerical Solution of Partial Differential Equations: An Introduction, 2005, 2nd edition, Cambridge University Press.

39 Reddy J. N., An Introduction to the Finite Element Method, 2006, 3rd edition, McGraw-Hill.

40 Taylor C. and Hughes T. G., Finite Elements Programming of the Navier–Stokes Equation, 1981, Pineridge Press.

41 Smith I. M., Griffiths D. V., and Margetts L., Programming the Finite Element Method, 2014, 5th edition, John Wiley & Sons.

42 Tsega E. G., A Numerical Scheme Based on Finite Element Method for Solving Two-Dimensional Time-Dependent Convection-Diffusion-Reaction Equation, International Journal of Algorithms, Computing and Mathematics. (2025) 11, no. 5, 187–220, https://doi.org/10.1007/s40819-025-02021-6.

43 John V. and Novo J., Error Analysis of the SUPG Finite Element Discretization of Evolutionary Convection-Diffusion-Reaction Equations, SIAM Journal on Numerical Analysis. (2011) 49, no. 3, 1149–1176, https://doi.org/10.1137/100789002, 2-s2.0-79960435756.

44 Fendoğlu H., Bozkaya C., and Tezer-Sezgin M., DBEM and DRBEM Solutions to 2D Transient Convection-Diffusion-Reaction Type Equations, Engineering Analysis With Boundary Elements. (2018) 93, 124–134, https://doi.org/10.1016/j.enganabound.2018.04.011, 2-s2.0-85046878475.