An unconditionally stable implicit difference scheme for 2D porous medium equations using four-point NEGMSOR iterative method

In this paper, a numerical method has been proposed for solving several two-dimensional porous medium equations (2D PME). The method combines Newton and Explicit Group MSOR (EGMSOR) iterative method namely four-point NEGMSOR. Throughout this paper, an initialboundary value problem of 2D PME is discretized by using the implicit finite difference scheme in order to form a nonlinear approximation equation. By taking a fixed number of grid points in a solution domain, the formulated nonlinear approximation equation produces a large nonlinear system which is solved using the Newton iterative method. The solution vector of the sparse linearized system is then computed iteratively by the application of the four-point EGMSOR method. For the numerical experiments, three examples of 2D PME are used to illustrate the efficiency of the NEGMSOR. The numerical result reveals that the NEGMSOR has a better efficiency in terms of number of iterations, computation time and maximum absolute error compared to the tested NGS, NEG and NEGSOR iterative methods. The stability analysis of the implicit finite difference scheme used on 2D PME is also provided.

(3) with some real parameters: 1 m and 0 .Equation (1) can be known as the generalized form of the heat equation.By referring to equation (1), the heat equation can be obtained by taking 1 m . In a particular case like 2 m , the Boussinesq equation can be obtained.Equation ( 1) is a type of partial differential equation in the nonlinear degenerate parabolic class which is widely used to describe models for gases in porous media, heat diffusion, underground infiltration, population dynamics as well as high energy physics [1].The presence of a nonlinear term m u in equation (1) makes the exact solution difficult to be obtained.
Generally, the numerical method can be considered as an alternative option in solving the partial differential equation especially for the case where the exact solution is almost impossible to be obtained.The numerical method yields an approximate solution that is differing from the exact solution by less than a specified tolerance.Furthermore, in terms of the implementation of the numerical method, efficiency is one of the important characteristics to be observed.The selected numerical method should be the one that requires less computation time and complexity in obtaining the desired result [2].
The numerical solution of multidimensional PME has received considerable interest from a number of researchers.For instance, Borana et al. [3] studied the numerical solution of the one-dimensional (1D) case of PME that is modelled based on the instability phenomenon in the double phase flow of two immiscible fluids through the inclined homogeneous porous medium.The method used is a Crank-Nicolson finite difference scheme.Meanwhile, Pradhan et al. [4] applied the finite element method for solving the 1D PME model that representing the instability phenomenon in the recovery of oil from a reservoir.For the numerical solution of 2D PMEs, Ngo and Huang [5] studied the application of an adaptive moving mesh finite element method as an improvement to the existing moving mesh method.In addition, Chew and Sulaiman [6] started the investigation on the efficiency of a numerical method that combines Newton and MSOR (NMSOR) iterative method for solving 1D PMEs.Then, Chew and Sulaiman [7] extended their study with the implementation of the Explicit Group (EG) method [8] to improve the efficiency of the NMSOR iterative method in solving the 1D PMEs.
The aim of this paper is to study the efficiency of the Newton and Explicit Group MSOR iterative method or NEGMSOR for solving several examples of 2D PME.The outline of this paper is arranged as follows.Section 2 shows the discretization of 2D PME by using the implicit finite difference scheme into a nonlinear approximation equation and the formation of a sparse and large nonlinear system.Section 3 presents the solution of the nonlinear system by the implementation of the Newton and EGMSOR iterative methods.Section 4 illustrates the numerical experiments and Section 5 concludes the finding of this study.

Implicit finite difference approximation equation
To obtain the numerical solutions, we consider a rectangular domain that is marked off in the unit square using the grid lines as follows: and for the time step, we define The numerical solutions that are represented by the points of intersection between the two grid lines are denoted by . For the sack of simplicity, we use By using the implicit finite difference scheme, we obtain the implicit finite difference approximation equation of equation ( 1) that can be written as where 6) is the function at each interior grid point over the solution domain that links with neighbourhood grid points to both x-and y-directions.As a result, a sparse and large nonlinear system is formed and generally written as where .

Linear stability analysis
Generally, the application of the linear Fourier analysis to the nonlinear partial differential equation such as 2D PME cannot be rigorously justified.Nevertheless, it has been found to effective in practice [9][10][11][12].To show the stability of the implicit finite difference scheme used on 2D PME, we assume that the solution function t , y , x u is bounded in the spatiotemporal region of 1 0 y , x and 1 0 t .Now, the nonlinear term m u can be frozen at each grid point in the solution domain and gives where . Equation ( 8) can be rewritten into From here, the stability of the implicit finite difference scheme used in this paper can be shown by implementing the von Neumann method over equation (9) as shown in Theorem 3.

Theorem 3.
The fully implicit finite difference scheme of equation ( 9) on the finite domain 1 0 y , x with zero boundary conditions for all 0 t is unconditionally stable.
Proof.To begin, we define the corresponding implicit finite difference operators as follows. and By substituting equation ( 10) -( 12) into equation ( 9) and after some rearranging, the corresponding equation obtained is 4 Derivation of the four-point NEGMSOR iterative method In the context of reducing the computational complexity, the convergence rate of the solution process needs to be increased.There are several iterative methods that have been proposed in order to reduce the computational complexity and this paper considers the application of the Explicit Group (EG) method which was introduced by Evans [8].This method uses the concept of breaking the original matrix into several small fixed-size groups of points.By this way, the computational complexity can be reduced by solving several small groups of points instead of a large linear system per iteration.The EG method has been used extensively together with some efficient iterative methods such as SOR, MSOR and AOR [13][14][15].Before the EG method is applied, equation ( 17) can be transformed into a linear system by using the Newton method as follows: where is a penta-diagonal Jacobian matrix with a dimension of , and for the iterative index, we denote it as and 1 n is considered to be the next time level.Then, equation ( 17) is divided into several completed groups of four points and some ungroup that can be treated as a group of two points and/or single point [8].Based on any completed group of four grid points, the general form of the EG iterative method is given by where Besides that, Kincaid and Young [16] introduced Modified Successive Over-Relaxation (MSOR) method which can be considered as a point iterative method that has two relaxation parameters.The basic idea of the MSOR method is to associate the different parameter with each row and column of a sparse and large linear system.Therefore, to derive the NEGMSOR iterative method, two different relaxation parameters, , are added to equation ( 18) and forms Based on equations ( 19) and (20), the NEG and NEGSOR iterative methods can be obtained by choosing within the range of 2 1, respectively.The implementation of the NEGMSOR iterative method can be shown in Figure 1.By referring to Figure 1, the numerical solution of 2D PME by using NEGMSOR iterative method can be summarized into Algorithm 1.The optimum values of 1 and 2 are determined when the least number of iterations is obtained after a few executions.
During the executions, we consider a tolerance error of iii.Formulation of the square matrix and the column matrix From the numerical experiment results in Table 1-3, the reduction percentage of the NEGMSOR iteration family as compared the GS iterative method for three examples are summarized in Table 4. Referring to the results of numerical experiments for three tested examples in Table 4, it can be pointed out that the percentage reduction of number of iterations for the NEGMSOR iterative method have declined approximately by 60.29-95.94%,55.38-95.54%,and 75.51-96.09%respectively as compared with the NGS method.In fact, implementations of execution time for NEGMSOR iterative method are much faster about 40.96-94.24%,52.59-94.68%,and 41.84-97.47%respectively than the NGS method.

Conclusion
This paper proposed the unconditionally stable implicit difference scheme together with the NEGMSOR iterative method for solving several examples of 2D PME.The implicit finite difference scheme is successfully applied to discretize the general form of 2D PME.In addition to that, the sparse and large linear system is solved using the NEGMSOR iterative method and the numerical results reveal that the proposed method is more efficient in terms of number of iterations, and execution time as compared to the tested NGS, NEG and NEGSOR iterative methods.Overall, in term of and the magnitude of maximum absolute error, the accuracy of the proposed NEGMSOR iterative method and the three tested numerical methods are comparable.

1 Introduction
This paper considers the numerical solution of the following initial-boundary value problem of the two-dimensional porous medium equation (2D PME):

Table 2 .
Comparison in terms of execution time (seconds) based on all iterative methods for three examples.

Table 3 .
Comparison in terms of maximum absolute error based on all iterative methods for three examples.

Table 4 .
Reduction percentage of the number of iterations and the execution time for NMEGSOR iteration family compared to NGS method