Solving One-Dimensional Porous Medium Equation Using Unconditionally Stable Half-Sweep Finite Difference and SOR Method

A porous medium equation is a nonlinear parabolic partial differential equation that presents many physical occurrences. The solutions of the porous medium equation are important to facilitate the investigation on nonlinear processes involving fluid flow, heat transfer, diffusion of gas-particles or population dynamics. As part of the development of a family of efficient iterative methods to solve the porous medium equation, the Half-Sweep technique has been adopted. Prior works in the existing literature on the application of Half-Sweep to successfully approximate the solutions of several types of mathematical problems are the underlying motivation of this research. This work aims to solve the one-dimensional porous medium equation efficiently by incorporating the Half-Sweep technique in the formulation of an unconditionally-stable implicit finite difference scheme. The noticeable unique property of Half-Sweep is its ability to secure a low computational complexity in computing numerical solutions. This work involves the application of the Half-Sweep finite difference scheme on the general porous medium equation, until the formulation of a nonlinear approximation function. The Newton method is used to linearize the formulated Half-Sweep finite difference approximation, so that the linear system in the form of a matrix can be constructed. Next, the Successive Over Relaxation method with a single parameter was applied to efficiently solve the generated linear system per time step. Next, to evaluate the efficiency of the developed method, deemed as the Half-Sweep Newton Successive Over Relaxation (HSNSOR) method, the criteria such as the number of iterations, the program execution time and the magnitude of absolute errors were investigated. According to the numerical results, the numerical solutions obtained by the HSNSOR are as accurate as those of the Half-Sweep Newton Gauss-Seidel (HSNGS), which is under the same family of Half-Sweep iterations, and the benchmark, Newton-Gauss-Seidel (NGS) method. The improvement in the numerical results produced by the HSNSOR is significant, and requires a lesser number of iterations and a shorter program execution time, as compared to the HSNGS and NGS methods.


Introduction
Porous medium equation (PME) is a nonlinear parabolic partial differential equation that exists in many nonlinear physical occurrences. For instance, PME is a general equation that brings up the Boussinesq equation that is used to model groundwater flow. PME is also used to describe the flow of ideal gas in a homogeneous porous medium, which is formulated by the laws such as mass balance, Darcy's law and state equation. In addition, PME is an important equation to be solved for a better understanding of the theory of heat propagation, particularly involving temperature-dependent thermal conductivity [1].
From the application of the PME side, [2] analyzed the heat transfer through human tissue, and found that the transport theory of porous media can be applied onto the biological heat transfer, as the theory reduces the number of assumptions when compared to other existing biological heat models. Then, [3] studied the qualitative properties of the PME in order to describe the dispersal processes in the dynamics of living things. The author found that the PME can be used to improve the qualitative as well as the quantitative agreement of population dynamics models. PME, without doubt, has great importance in many scientific fields, and more details about the theory and application of PME are available in [1].
The solutions of several one-dimensional PME problems via the finite difference method have been studied by many researchers [4][5][6][7][8][9]. As part of the development of a family of efficient iterative methods to solve the PME, this research adopted the Half-Sweep technique in the formulation of the finite difference method. Several researchers have discussed the success of the Half-Sweep technique in approximating the solutions of several types of mathematical problems [10][11][12][13][14][15][16]. Motivated by the unique property of Half-Sweep in securing a low computational complexity while computing the numerical solutions, this work aims to solve the one-dimensional PME using the unconditionally stable Half-Sweep finite difference approximation.
For this particular nonlinear type of partial difference equation, the finite difference discretization through the implementation of Half-Sweep yields a nonlinear type of approximation equation. Before the solution of PME is computed, the formulated nonlinear approximation equation is linearized using the Newton method to form a sparse and large linear system. A Successive Over Relaxation (SOR) iterative method with optimum parameters was applied for an efficient solution to a generated linear system.

Half-Sweep Finite Difference Method
Let us consider the general form of the one-dimensional PME [17]: where and are assumed to be any rational numbers.
It is worthy to mention that Eq. (1) can exist for all ∈ ℝ and 0 < < ∞. For our numerical study, we attempt to investigate the numerical solution of Eq. (1) in a rectangular domain, subject to the boundary and initial conditions, as follows: where 0 ( ), 1 ( )and 0 ( )are the prescribed functions based on the provided exact solutions. Before we show the formulation of Half-Sweep finite difference approximation to Eq. (1), it is best to discuss the formulation of the standard implicit finite difference approximation to Eq. (1), since our proposed method is based on the implicit finite difference method. Now, by defining the approximate solutions to Eq. (1), , = ( , ) , = 0,1,2, . . . , − 1, = 0,1,2, . . . and both spatial and temporal steps are = 1/ and = 1/ respectively, the standard implicit finite difference approximation equation to Eq. (1) can be written as follows [4,7,8,18]: where The approximation equation shown in Eq. (3) can also be known as the Full-Sweep finite difference approximation equation, because it approximates all mesh points in a bounded domain. Hence, Eq. (3) can be extended to develop our Half-Sweep finite difference approximation equation by lengthening the distance between two consecutive mesh points from to 2 , as follows [8]: where The approximation equation (4) is proven to be unconditionally stable, and the proof is at the appendix.
Using Eq. (4), we may obtain a nonlinear system for time level + 1, in the form of: where +1 = � 2, +1 , 4, +1 , … , −2, +1 � and for each function, Since solving the nonlinear system (5) deals with great computational cost, we use the Newton method to linearize the nonlinear system (5), and then apply the SOR iterative method to obtain the solution. Using the Newton method, the linear system can be written as follows [7,8,19]: where +1 The approximate solutions to Eq. (1) are computed by:

HSNSOR Iterative Method
Based on the linear system (7), we find out that the coefficient matrix +1 ( ) has the form of a tridiagonal. Thus, to apply the SOR iterative method for solving the linear system (7) [20,21], we consider the three components' decomposition of +1 ( ) as follows, where +1 ( ) is the diagonal of the matrix, +1 ( ) is the strictly lower triangular matrix, and +1 ( ) is the strictly upper triangular matrix, at the time level + 1 and -th iteration.
Hence, using the linear system (7) and the decomposition (9), the proposed method (HSNSOR) can be derived into Based on the formula shown in (10), the relaxation parameter lies within 1 < < 2 . When = 1 , the formula can be known as the HSNGS [8]. According to Figure 1, the implementation of the HSNSOR method for solving Eq. (1) can be explained as follows. After the boundary and initial conditions were imposed on the solution domain, the HSNSOR approximates the solutions on all the interior mesh points that are labelled by black dots, i.e. 2, 4, …, 14. After the iteration process is completed and the values of the black dots are obtained, the remaining mesh points that are labelled by white dots, i.e. 1, 3, …, 15 are computed directly using the approximation equation. The full algorithm for the computation using the HSNSOR method is described in Algorithm 1.
If the correctors converge, compute (8) and then the remaining mesh points, vi.
In practice, the optimum value of is determined (±0.01) by running Algorithm 1 several times, and the one that gives the least number of iterations is selected as the optimum.

Stability Analysis of the Half-Sweep Finite Difference Method on the One-Dimensional Porous Medium Equation
The application of Fourier analysis to prove the stability of the applied finite-difference on nonlinear partial differential equation (like PME) cannot be rigorously justified. Nevertheless, it is practically effective [22].

Numerical Experiment
For the numerical experiment, several criteria are observed such as the number of iterations ( ), the program execution time (seconds) and the magnitude of absolute errors ( ) . These criteria are used to evaluate the efficiency of the HSNSOR method to solve Eq. (1) subjects to both initial and boundary conditions as in Eq. (2). The efficiency of the HSNSOR method is then compared to the HSNGS and NGS [23] methods using four selected examples. Four examples used for the numerical experiment are presented hereafter.

Example 4 [24]
Given a one-dimensional PME wherein equals to -2 and the parameter is 0.5: The exact solution is ) with five different sizes of mesh points, . Additionally, Table 5 shows the percentages of reduction in the number of iterations and the program execution time by the HSNSOR and HSNGS against the control method, NGS.