Solution of convection-diffusion model in groundwater pollution

This research involves the development of the spectral collocation method based on orthogonalized Bernoulli polynomials to the solution of time-fractional convection-diffusion problems arising from groundwater pollution. The main aim is to develop the operational matrices for the fractional derivative and classical derivatives. The advantage of our approach is to orthogonalize the Bernoulli polynomials for the sake of creating sparse operational matrices in such a way that classical derivatives have one sub-diagonal non-zero entries only, and also creating an operational matrix for fractional derivative have diagonal matrix only. Due to these properties, the cost of computational our approach is very low and the convergence is fast. A discussion on the error analysis for the presented approach is given. Two test problems are considered to illustrate the effectiveness and applicability of our method. The absolute error in the computed solution compares with the existing method in the literature. The comparison shows that our method is more accurate and easily implemented.

) is the time fractional diffusion equation(TFDE).This type of equation governs the evolution of the probability density function that describes anomalously diffusing particles.1) is the time fractional advection diffusion equation(TFADE).This equation can be solved to determine the changes in tracer concentration with space and time.Also used for water mass and marine particle transport modeling and sediment diagenesis.
Many researchers trying to solve the above-mentioned equation.Using the Chebyshev wavelet collocation method proposed in 44 .The radial basis function (RBF) method combined with a modification of the finite integration method (FIM) derived in 45 .The sinc-Galerkin method is proposed in 46 .The Sinc-Legendre collocation method proposed in 47 .The Chebyshev collocation method was developed in 48 .Wherever the Eq.(1) be the constant coefficients, in 49 collocation method based on RBF is developed.In the case, when d(ξ ) = −1 and b(ξ ) = 0 , authors used finite difference and finite element 50 .
In this research, we propose a new approach for the solution of TFCDEs.Our approach is to develop the operational matrices based on the orthogonalized Bernoulli polynomials.Our approach is based on the spectral collocation method based on orthogonalized Bernoulli polynomials and developing operational matrices for derivatives which is very sparse and have one sub-diagonal non-zero and also developing an operational matrix for fractional derivative which is only diagonal matrix.Due to these properties, the method is very fast, and computational time is low.Numerical experiments demonstrate the accuracy and efficiency of the proposed method.
The layout of this work is as follows: In Section "Preliminary definitions", some preliminary definitions are given.In Section "The proposed method", we developed an approach for the solution of the Eq. ( 1)-(3).In Section "Error analysis", we study the error analysis of the proposed method.In Section "Numerical results", the presented method tested on two problems for verification, applicability and to show the nature of accuracy of the proposed method.We compare our results with the references 45,47,48,[51][52][53] .

Fractional derivative
By following 54 , we recall the essential concepts: Definition 1 The Riemann-Liouville fractional integral operator of order γ,(γ ≥ 0) , is defined as with I 0 t u(t) = u(t).

Definition 2
The fractional derivative of u(t) in the Caputo derivative is defined as These polynomials have useful properties that we do not express here; for further information see 55,56 .Despite that, these polynomials have interested properties but are not orthogonal.To overcome such disadvantages of Bernoulli, by using the Gram -Schmidt process, we try to obtain an explicit form of orthogonal Bernoulli polynomials.

Definition 5
The explicit form of orthogonal Bernoulli polynomials (OBPs) of p th degree is as follows 57 : so that where δ ij denotes the Kronecker delta function.The form of operational matrix based on the orthogonal Bernoulli polynomial is as follows: where and Since A is a lower triangular matrix with nonzero diagonal elements.Therefore A is nonsingular.Thus, we have Any function u(t) defined on (0, 1] can be expand by: (5) where the coefficients In the application, we consider only the first (p + 1)-terms OPBs, so that, we could write where B T = β 0 , β 1 , . . ., β p and � p (t) = ω 0 (t), ω 1 (t), . . ., ω p (t) T .
In same manner, any two variables function u(ξ , t) can be expand by the OPBs series as: where H = h ij (p+1)×(q+1) and and According to Eq. ( 10), we have where and according to Eq. ( 13), we get

The proposed method
We develop the spectral collocation scheme based on orthogonal Bernoulli polynomials in both the space and time direction.

Error analysis
We aim to obtain an estimation of the error bound for the approximation of the function u(ξ , t) ∈ � = [0, 1] × (0, 1] with orthogonal Bernoulli polynomials.
Proof Using Theorem 1, we have: Now, by using the useful properties of the Caputo derivative and Riemann-Liouville integral defined in 58 and Eq. ( 41), we obtain: Now, we know that I m t is the operator integral Riemann -Liouville, Then this For proof, we need to introduce an upper bound for I m t 2 .Using the definition of the left Riemann -Liouville integral operator and Schwarz's inequality, we have: Therefore, we have: Finally, by using ( 42) and ( 45), we obtain: .

Numerical results
Here we consider two problems to demonstrate the reliability and efficiency of our method.The absolute errors of the solution for different nodes (ξ i , t j ) ∈ � defined as this: also, the maximum absolute errors(MAEs) are The presented method applied to these problems for various p, q and γ .All the programming in Matlab 2016 software is run by a computer with core i5.
Consider the following TFCDE 48,51 : subjected to the initial and boundary Conds.
The computed solution and absolute errors in the solution for p = q = 4 and γ = 0.5 were obtained.We plot the graph of the computed solution, the exact solution, and absolute errors in solution in Fig. 2a-c respectively.In Fig. 3, we plot the computed solution and the exact solution for p = q = 3 with γ = 0.7 at t = 0.5 .We show that the computed solution coincides with the exact solution.In Fig. 4a-c the contour plots of the computed solution, the exact solution and absolute errors are represented for p = q = 4 with γ = 0.5 .Moreover in Table 1, we compared the MAE of the proposed method with the methods in 45,47,48,51,52 .Table 1, shows that our method is more efficient and accurate than the compared methods.Our method for p = q = 2 yields the solution of a system of 3 × 3 but in 51 for m = 5, J = 3 needs to solve a system of 110 × 110 .The reference 52 needs to solve the system of 32 × 32 .The reference 47 needs to solve the system of 31 × 8 .The reference 48 needs to solve the system of 6 × 6 .Finally, the method in 45 needs to use 10 grid points in time and space but our method uses 3 grid points in time and space.The absolute errors in the solution for this problem for p = q = 3 and p = q = 4 with γ = 0.5 for t = 0.1 and the run time of our method is tabulated to show the effectiveness of our method in Table 2.
The absolute errors are computed.The comparison of absolute errors of the proposed method with the method in 53 tabulated in the Table 3, shows that our method is more efficient and accurate.Actually our method for p = q = 3 yields the solution of system of 4 × 4 but in 53 for k = 4, M = 1 needs to solve a system of 48 × 21 .In Table 4, we report the absolute errors at t = 0.1 for various γ , p = q = 4 , and the run time of our method is tabulated to show the effectiveness of our method.
Figure 5a,b represent the computed and the exact solutions for p = q = 5 with γ = 0.5 and γ = 0.9 respec- tively.We plot the graph of the exact solution, the computed solution, and absolute errors in solution for p = q = 5 and γ = 0.9 . in Fig. 6a-c respectively.In Fig. 7a-c the contour plots of the computed solution, the exact solution and absolute errors are represented for p = q = 5 with γ = 0.9 , respectively.

Conclusion
An efficient, accurate spectral collocation method based on orthogonalized Bernoulli polynomials has been developed for time fractional convection-diffusion problems.Our approach contains operational matrices for approximate derivatives as well as fractional derivative.Operational matrices for derivatives are sparse having one sub-diagonal non-zero entries only, and for the fractional derivatives operational matrix is diagonal only.
Absolute errors for values of p = q = 2 , p = q = 3 and p = q = 4 , γ = 0.5,t = 0.1 for Problem 1.    Due to these properties and using spectral methods convergence our presented method is spectral and fast with low computation cost.The comparing numerical results justifies the effectiveness and accuracy of our proposed scheme.

Figure 1 .Definition 3
Figure 1.Landfill is one of the sources of groundwater pollution.

Figure 3 .
Figure 3.The comparison with the computed and exact solutions at t = 0.5 with γ = 0.7and p = q = 3 for Problem 1.
by solving the above system of equations, we obtain the coefficients matrix H . • Step 6.Finally, by replacing matrix H in Eq. (17), we approximate the solution of TFCDE.