1. Introduction
Numerical methods typically require large mesh sizes to give better approximations to solutions of differential equations. Consequently, their memory requirements are huge and, therefore, computational efficiency is affected. In this work, we aim to improve the computational efficiency of the higher-order unconditionally positive finite difference method (HUPFD). The HUPFD is based on using higher-order finite difference schemes to solve differential equations, with its main attribute being the preservation of the positivity of solutions. The positivity of solutions is essential in several physical settings—for example, quantities like population sizes, chemical species concentration and neutron numbers require positive solutions in order to be meaningful.
One of the techniques researchers employ to reduce the computational span of computational methods is the proper orthogonal decomposition method (POD) [
1]. This has previously been used to reduce the dimensions of numerical methods such as non-standard finite difference and Crank–Nicolson to solve problems such as the Sobolev, hyperbolic, Burgers and Navier–Stokes equations [
2,
3,
4]. In this paper, we couple the POD and HUPFD to obtain a new enhanced hybrid method that we name the enhanced higher-order unconditionally positive finite difference method (EHUPFD). By utilizing the POD, the EHUPFD extracts a set of basis functions from the HUPFD solution called the snapshot matrix and then uses a small subset of leading basis functions to construct state variable approximations.
We test the applicability of the EHUPFD on the advection diffusion reaction equation below
where
D is the diffusion coefficient,
is the advection coefficient,
is the advection term,
is the diffusion term,
represent the rate of change of concentration of substances and
is the source or reaction term [
5]. The ADR equation governs the process of diffusion and advection simultaneously [
6]. The ADR equation is used to model the exponential traveling waves semiconductor modeling, biological system modeling, pollutant absorption in soil and neutron diffusion [
6,
7,
8,
9,
10]. The UPFD method has previously been used to solve the ADR equation, amongst other known numerical methods reported in the literature, such as the standard finite difference, Crank–Nicolson finite difference, and non-standard finite difference methods and acts. The UPFD method has been proven to be a stable numerical method, accurate and consistent when solving the ADR equation in rectangular domains. The method is also known to preserve the positivity of solution regardless of step size [
11]. For example, Appadu, A and Rao used UPFD to approximate the solution to the linear ADR equation, which models exponential traveling waves. The authors studied various scenarios in which advection or diffusion is dominating, as well as when both are equal. In this study, they discovered that the UPFD scheme is successful in solving the ADR equation when advection or diffusion factors are dominating [
12]. However, as always, there is a need to constantly refine the numerical methods in order to improve the method’s accuracy, computational time and convergence rate.
The HUPFD method was developed to solve the advection diffusion reaction equations to improve the solution’s accuracy [
13]. The higher-order numerical methods have been used to solve mathematical problems such as [
13,
14,
15,
16,
17,
18,
19,
20]. The standard finite difference methods, such as the backward, forward and central difference schemes, struggle to approach the exact solution of most partial differential models with large step sizes. As a result, relatively small step sizes are chosen to improve the accuracy of the solution to the problem at the cost of increased computational time due to the increased number of grid points. With a smaller grid size, HUPFD calculates a more accurate solution to the ADR equation’s linear, non-linear and system versions compared to the UPFD results.
The highly accurate finite difference method has gained much attention in recent years [
21,
22,
23,
24,
25]. Gemechis FD and Tesfaye AB studied the higher-order numerical method to solve the singularly perturbed delay reaction–diffusion equation [
26]. In order to achieve the six higher-order schemes, they used the Richardson extrapolation technique [
26]. They found that the developed method is stable and consistent, and produces a more accurate solution than some of the methods reported in the literature. Ranjana KM and Venn G investigated the fourth-order finite difference method based on spline in tension approximation for the solution of second-order quasi-linear hyperbolic equations in one space dimension [
27]. They proposed a three-level implicit nine-point compact finite difference formulation of order two in the time and four in the space directions based on split tension approximation in the t-direction for the numerical solution of space dimensional second-order quasi-linear hyperbolic equations with a first-order space derivative term. The proposed numerical method for the wave equation in polar coordinates is conditionally stable [
27]. In contrast, the method for the damped wave and telegraphic equations is unconditionally stable. Sari M et al. [
28] proposed a tenth-order finite difference scheme for solving a one-dimensional diffusion equation. Gurarslan et al. [
29] presented a sixth-order compact difference scheme and a fourth-order Runge–Kutta scheme in time to produce numerical solutions of the one-dimensional advection–diffusion equation, which was found to be accurate in solving the concurrent transport equation for the one-dimensional advection diffusion reaction equation.
The rest of the sections in the paper are organized as follows.
Section 2 and
Section 3 introduce the UPFD and HUPFD, respectively. Subsequently,
Section 4 and
Section 5 apply the schemes in
Section 2 and
Section 3 to solve the ADR equations.
Section 6 introduces the POD technique used to enhance the HUPFD, which results in the development of EHUPFD. The results obtained by the EHUPFD are illustrated in
Section 6 and a summary of the results is provided in
Section 7.
5. Formulation of the POD Technique
In this section, the POD basis function is constructed and used to reduce the degree of freedom to the UPFD and HUPFD scheme.
The POD technique is used for the HUPFD by extracting the snapshot matrix in order to produce the EHUPFD, which leads to the construction of the approximate state variable by developing a small subset of the leading basis function obtained from the SVD factorization. The POD technique seeks a low-dimensional basis function that can accurately approximate the state variable, which is a left singular vector from the SVD.
The first
L sequence of solutions is extracted from the HUPFD solutions as given by
,
and
, which gives the formulated
snapshot matrix given below
where
m and
L represent the size of the snapshot matrix, which is factorized by using the singular value decomposition (SVD) as given below
where
is an
matrix representing an orthogonal matrix consisting of the eigenvectors of
,
diag
represents the diagonal matrix from the SVD, given in decreasing order
, and
is an
orthogonal matrix which represent the orthogonal eigenvector of
and
;
O is a zero matrix [
32].
From identical eigenvalues of
and
, we obtain the eigenvalues of
for matrices
and corresponding eigenvectors
, which is given by the formula below
This leads to the construction of the matrix
given below
where the diagonal matrix given by
is constituted of the first
M positive singular values of the diagonal matrix.
5.1. The Implementation of the Algorithm of the EHUPFD Scheme for the Advection Diffusion Reaction Equations
Below are the steps to implement the algorithm of the EHUPFD to solve the linear, non-linear and system ADR equations [
32].
- Step 1:
Snapshot matrix for HUPFD is extracted.
- Step 2:
Formulation of the matrix and calculation of eigenvalues and eigenvectors of matrix and executed.
- Step 3:
Choice of POD basis for the error tolerance of is executed.
- Step 4:
The EHUPFD computes the reduced-order solution.
- Step 5:
Accuracy of the solution is checked for the renewal of the POD basis.
5.2. The EHUPFD Error Estimate
The following equation is used to estimate the error of the EHUPFD when solving the ADR equations.
Theorem 2 (Accuracy)
. The error estimate for the HUPFD and from the EHUPFD is calculated by using the following equationThe absolute error, calculated between the actual and approximate solutions from the POD reduced-order model satisfy the following equation [32] 6. Results and Discussion
In this section, we present solutions to the linear, non-linear and system ADR equations solved by using the EHUPFD, UPFD, NSFD and Crank–Nicolson methods. To check the performance of the methods, absolute errors, convergence rates and computational time are measured. The actual and numerical solutions are denoted by
and
, respectively. The error is represented by
and its magnitude is calculated by using norm to infinity denoted by
-norm given by the formula below,
The formula to calculate the convergence rate is given below,
Finding the best snapshot size is paramount because it influences the EHUPFD solution accuracy. The best snapshot is determined by finding the solution to the problem when
and then an increment of 1 to
L is evaluated until there is no significant change in the solution. The tolerance is calculated using the formula below
where
is a solution for a specific value of
L and
is the exact solution. The first value of
L that satisfies Equation (
106) is the optimal value of
L. We set
.
6.1. Results Obtained to Solve the Linear ADR Equation
The convergence rates of the UPFD, Crank–Nicolson, NSFD and HUPFD with respect to the time and spatial steps are shown in
Table 1. The results show that the EHUPFD converges to the solution faster than the UPFD.
Table 2 compares solutions obtained by the exact, UPFD, Crank–Nicolson, NSFD and EHUPFD solutions at
and
.
Table 3 compares the infinity norm results of the developed EHUPFD method against the Crank Nicolson and NSFD, as well as their computational times. The results show that the EHUPFD takes less computational time to converge to the solution as compared to the HUPFD methods.
Table 3 compares the infinity norms of the proposed EHUPFD methods with the UPFD, Crank–Nicolson and NSFD methods in solving the linear ADR equation. The results show that EHUPFD methods preserve solution accuracy.
Figure 1 shows the contribution made by each POD basis mode. The first five modes contain the predominant information where the EHUPFD meets the theoretical accuracy requirement. We observed that five modes are significant in optimizing the linear ADR equation.
Figure 2b compares the number of POD modes against the absolute error. We observe that only five POD modes are required in order to converge to the toleration of
.
Figure 3a to
Figure 4a show the surface plots of the exact, UPFD and EUPFD solutions of the linear ADR.
Figure 4b shows the increase in error when the POD modes increases. The figures show that the developed methods achieve positivity of the solutions and the results are comparable with well-known methods such as Crank–Nicolson and NSFD.
6.2. Results Obtained by Using the Non-Linear ADR Equation
The convergence rates of the UPFD, Crank–Nicolson, NSFD and the EHUPFD with respect to the time and spatial steps are shown in
Table 4.
Table 5 compares the exact, UPFD and EHUPFD solutions at
and
. The HUPFD has a slightly higher convergence rate than UPFD, while the EHUPFD has significantly reduced the computational time taken to converge to the solution of the ADR equation.
In
Table 6, we compare the infinity norms of the proposed HUPFD with the UPFD, Crank–Nicolson and NSFD methods in solving the linear ADR equation. For this example, the results show that the EHUPFD method preserves the accuracy and reduces the computational time taken to converge to the solution of the ADR equation.
Figure 5a shows the contribution made by each POD basis mode. The first five modes contain the predominant information where the EHUPFD meets the theoretical accuracy requirement. We observed five modes are significant in optimizing the linear ADR equation.
Figure 5b shows the increase in computational time when the number of modes increases.
Figure 6a compares the results obtained by the EHUPFD, Crank–Nicolson and NSFD methods. The results show that the EHUPFD solution is comparable with the Crank–Nicolson and NSFD solutions.
Figure 6b compares the error against the number of modes. The results show that the error of the solution decreases when the number of modes increases.
Figure 7a to
Figure 8a show the surface plots of the exact, UPFD and EUPFD solutions of the non-linear ADR equation.
Figure 8b show the absolute errors between the exact solution and the EHUPFD solution for the non-linear ADR against time.
6.3. Results Obtained When Solving the Schnakenberg Model
The convergence rates of the UPFD, HUPFD and EHUPFD with respect to the time and spatial steps are shown in
Table 7. The results show that the convergence rate of the EHUPFD methods is slightly higher than the UPFD when solving the Schnakenberg model.
Table 8 and
Table 9 compare the solutions of the UPFD, HUPFD and EHUPFD solutions at
and
, as well as the computational times. The results show that the EHUPFD takes less computational time to converge to the solution than the UPFD, HUPFD and EHUPFD methods.
Figure 9a to
Figure 10b show the surface plots of the UPFD, HUPFD and EHUPFD of the Schnakenberg model. The figures show that the developed methods achieve positivity of the solutions.
Figure 11 compares the solutions obtained by the exact HUPFD, EHUPFD and UPFD. The results show that the proposed methods are comparable to the known methods of the UPFD.
The convergence rates of the UPFD, HUPFD and EHUPFD with respect to the time and spatial steps are shown in
Table 7. The results show that the convergence rate of the EHUPFD methods is slightly higher than the UPFD when solving the Schnakenberg model.
Table 8 and
Table 9 compare the UPFD, HUPFD and EHUPFD solutions at
and
, as well as the computational times. The results show that the EHUPFD takes less computational time to converge to the solution than the UPFD, HUPFD and EHUPFD methods.
Figure 9a to
Figure 10b show the surface plots of the UPFD, HUPFD and EHUPFD of the Schnakenberg model. The figures show that the developed methods achieve positivity of the solutions.
Figure 11 compares the solutions obtained by the exact HUPFD, EHUPFD and UPFD. The results show that the proposed methods are comparable to the know methods of the UPFD.