A Newton-Modified Weighted Arithmetic Mean Solution of Nonlinear Porous Medium Type Equations

The mathematical theory behind the porous medium type equation is well developed and produces many applications to the real world. The research and development of the fractional nonlinear porous medium models also progressed significantly in recent years. An efficient numerical method to solve porous medium models needs to be investigated so that the symmetry of the designed method can be extended to future fractional porous medium models. This paper contributes a new numerical method called Newton-Modified Weighted Arithmetic Mean (Newton-MOWAM). The solution of the porous medium type equation is approximated by using a finite difference method. Then, the Newton method is applied as a linearization approach to solving the system of nonlinear equations. As the system to be solved is large, high computational complexity is regulated by the MOWAM iterative method. Newton-MOWAM is formulated technically based on the matrix structure of the system. Some initial-boundary value problems with a different type of nonlinear diffusion term are presented. As a result, the Newton-MOWAM showed a significant improvement in the computation efficiency compared to the developed standard Weighted Arithmetic Mean iterative method. The analysis of efficiency, measured by the reduced number of iterations and computation time, is reported along with the convergence analysis.


Introduction
The porous medium type equation is one of the classes of nonlinear parabolic evolution equation. This class of partial differential equations appears in the description and modeling of natural physical phenomena such as gas diffusion, fluid flow, and heat propagation. The mathematical theory behind the nonlinear porous medium type equation is well developed and produces many applications to the real world [1]. A porous medium type equation is used to investigate the diffusion of reactant gases across a thin layer of porous material. For instance, the authors of [2] developed a mathematical model to describe the flow of a mixture of ideal gases in a highly porous electrode for fuel cell engineering. In their developed model, a porous medium type equation is used to simulate the evolution of the gas mixture.
Porous medium type equation is also one of the important mathematical models used in underground oil production. The permeability of rock above the underground oil wells, which represents the porous medium for the oil diffusion, determines the maximum capacity of oil to be produced over a fixed period. Moreover, underground oil production is affected by a phenomenon called instability during the oil extracting process. The instability phenomenon appears in the form of irregular trembling fingers caused by the immiscibility of water and oil. One article has proposed a mathematical model using the formulated porous medium type equation to simulate the instability phenomenon [3]. In addition to that, porous medium type equation has its importance in life sciences such as bacteria biofilms growth model [4], cell population dynamics model [5], and cell to cell adhesion model [6].
Real-world phenomena modeled by the porous medium type equation contains arbitrary functions or constants. These functions and constants, which act as the model parameters, may be specified through several experiments or applying the well-established physical laws. However, some forms of these parameters are assumed in most numerical studies. Some numerical experiments also use random values for the parameters to investigate the accuracy and efficiency of the proposed numerical approach. Many examples of porous medium type equations with arbitrary parameters can be obtained in [7]. Besides that, the research and development of the fractional nonlinear porous medium type equation have progressed significantly in recent years and is recommended to read [8].
Many different numerical methodologies have been presented and applied to solve the mathematical models of porous medium types. The models are varied in the degree of nonlinearity, the number of derivatives terms, and the type of diffusion term. Thus, an efficient numerical method to solve porous medium models needs to be investigated so that the symmetry of the designed method can be extended to future fractional porous medium models. A new numerical method called Newton-Modified Weighted Arithmetic Mean (Newton-MOWAM) is introduced and studied in this paper. The solution of the porous medium type equation is approximated by using a finite difference method. Then, the Newton method is applied as a linearization approach to solving the system of nonlinear equations. The nonlinear equations generated by the finite difference approximation to a nonlinear partial differential equation usually possess a high computational complexity. In order to overcome the high computational complexity to solve the nonlinear problems, the complexity regulation by the Newton-MOWAM is formulated technically based on the matrix structure of the system. Therefore, the contribution of this paper is to propose the Newton-MOWAM method to solve the nonlinear porous medium type equations efficiently with rigorously analysis of convergence by developing a theorem and its corresponding proof.

Preliminary
The basic idea behind the numerical method to resolve the porous medium type equation can be represented in the form of a linear system as follows: where M is a coefficient matrix, F is a given vector, and ϕ signifies the unknown vector to be calculated. Basically, methods for solving Equation (1) can be categorized into two main streams, i.e., direct and iterative methods. In principle, the direct method is the solution of Equation (1) is determined through a finite number of arithmetic operations which without consideration of round-off errors. In contrast to direct method, iterative method generates a chronological order of approximations to the solution by recurrence application of the same computational prototype at each iteration. When the size of the linear system is large, iterative method is always preferable to be used compared to direct method. Direct methods have the complexity of order O(m 3 ), where m is the order of a matrix. In contrast, for iterative methods, the work per one iteration is essentially the matrix-vector multiplication, which for sparse matrices is sometimes much less than O(m 2 ). In addition, the number of iterations needed for the convergence, when using a suitable preconditioner, is usually much less than O(m). Thus, the overall work is much less than O(m 3 ), which means it has fast convergence. According to the book in [9], sometimes we do not need to solve the system of equations exactly. For example, a few methods for solving nonlinear system of equations exist in which each iteration involves a linear system to solve. Frequently, it is sufficient to solve the linear system within the nonlinear iteration only to a low degree of accuracy. Direct methods cannot accomplish this because, by definition, to obtain a solution, the process must be completed. There is no notion of early termination or an inexact solution.
A general technique to devise iterative method is based on an additive splitting of matrix M into M = P − Q where M must be a nonsingular matrix. Then, the corresponding iterative scheme for solving linear system as in Equation (1) is of the form with the new iterate ϕ ( +1) , previous iterate ϕ ( ) , vector P −1 F, and iteration matrix P −1 Q. Based on the iterative form shown by Equation (2), there are two main types of iterative methods, namely, stationary iterative and nonstationary iterative methods. Both types of methods are classified based on the nature of iteration matrix P −1 Q and vector P −1 F. For stationary methods, the iteration matrix P −1 Q and vector P −1 F remain constant throughout the iteration process, while a new iteration matrix P −1 Q and vector P −1 F are generated in every step of the nonstationary iterative methods. In this paper, several stationary iterative methods are reviewed, i.e., Alternating Group Explicit [10], Modified Alternating Group Explicit [11], Iterative Alternating Decomposition Explicit [12], Block SOR [13], and Weighted Mean (WM) [14,15]. Among these names of stationary iterative method, the concept from the WM method is adopted to develop a new efficient numerical method for porous medium type equation. The WM method is a group of algorithms that have been widely used to solve matrix problems efficiently. One of the WM methods to be investigated in this paper is the Arithmetic Mean (AM) iterative method [16]. In the previous literature, the AM method and its variants have been studied and scrutinized on linear and nonlinear systems arising from a variety of scientific problems such as [17][18][19][20][21][22]. In [14,16], the numerical simulations showed that the AM method compares favorably against the existing nonstationary iterative method, particularly when the iteration matrix is strongly asymmetric. Further, Aruchunan et al. [23] introduced a modified AM (MAM) method by initiating second weighted optimal value to improve the convergence of AM algorithm via solving fourth-order integro-differential equations. However, the authors did not discuss the convergence theorem and its corresponding proofs in the MAM method. Therefore, in this paper, the MAM method is extended together with Newton method for solving nonlinear porous medium type equation accompanied by the complete convergence theorem and proof. The complete procedures on developing Newton-MOWAM is described in the following segments.

Porous Medium Type Equation
This paper considers the numerical solution of a general one-dimensional porous medium type equation as follows [7]: where the function f (w) represents the type of diffusion term which is depending on the nonlinear physical phenomenon to be modeled. For example, when f (w) = w α with the real number α, the simplest form of porous medium equation used to model the instability phenomenon in oil recovery can be obtained [3]. Using a finite difference method, the solution domain of Equation (3) can be restricted about a unit square containing numerous number of well distributed unknown grid points. The restriction is usually assumed by imposing suitable initial condition w(x, 0) = w 0 (x), 0 ≤ x ≤ X and the Dirichlet boundary condition w(0, t) = a(t), w(X, t) = b(t), 0 ≤ t ≤ T.

Finite Difference Approximation to Porous Medium Type Equation
To formulate the finite difference approximation equation to Equation (3), it is more appropriate to differentiate the right hand side of Equation (3) into where f (w) is the derivative of the smooth function f (w).
To set up the network of grid points as the approximate solutions of Equation (4), we define the approximate solution by W (x, t) = W (i∆x, j∆t), i = 0, ..., m, j = 0, ..., n, where ∆x and ∆t are the spatial and temporal step sizes, respectively. In this work, a square-shaped domain is considered with the finite interval in space, [0, X] is divided into m − 1 intervals that yields m unknown grid points, and the finite interval in time [0, T] is divided into n − 1 intervals. By substituting a backward time operator and a second-order center space operator into the time and space derivatives in Equation (4), respectively, the implicit finite difference approximation to Equation (3) can be formulated as and simplified into As the discrete finite difference approximation equation shown by Equation (6) is nonlinear, it yields a system of nonlinear equations when a fixed number of unknown points is considered. The approximate solution of the porous medium type equation needs to be obtained through the solution of the system of nonlinear equations. Instead of solving the high computational complexity system directly, Newton method is adopted as an effective linearization approach to solve the system of nonlinear equations.

Newton Method
To apply the Newton method, Equation (6) needs to be rearranged into a function as follows: When a Jacobian matrix with respect to three unknown mesh points, i.e., W i−1,j+1 , W i,j+1 and W i+1,j+1 , is derived, a tridiagonal coefficient matrix can be formed. Therefore, the associated system of linear equations developed from the discretized Equation (3) is written in the form of with the following tridiagonal coefficients: We define the corrector to the approximate solutions of (3) as with the iteration index, and F = (F 1,j+1 , F 2,j+1 , . . . , F m−2,j+1 , F m−1,j+1 ) T .

Newton-Modified Weighted Arithmetic Mean
As stated in Section 1, the major contribution of this research is to establish new numerical method called Newton-MOWAM with its complete convergence theorem and proof. The Newton-MOWAM is developed by proposing two optimal parameters on the existing Weighted Arithmetic Mean (WAM) iterative method. Following that, the developed Newton-MOWAM method will be investigated by solving the system of linear equation in matrix form as in Equation (8). Therefore, in the next few subsections, the formulations of novel Newton-MOWAM method will be explained by demonstrating its convergence theorems and proofs. The iteration procedure for Newton-MOWAM method consists of two loops, where in Loop 1, an independent system,φ F is involved meanwhile in Loop 2 independent system ofφ B is involved. This work assume the number of unknown points m is an even number.

Weighted Optimal Value Identification
The weighted optimal value is crucial to accelerate the convergence process for the iterative methods [24]. Therefore, the Newton-MOWAM method has weighted optimal values for rapid convergence. In this research, another new weighted optimal value, Ω 2 is initiated for Loop 2. The optimal value Ω 2 was set based on trial with least number of iterations after the best Ω 1 is fixed. The execution of these two weighted optimal values improved the convergence rate of the Newton-MOWAM compared to the standard WAM methods. On condition that Ω 1 = Ω 2 , the Newton-MOWAM method is equivalent to the conventional WAM methods. The weighted optimal values Ω 1 and Ω 2 have values between 0 and 2 in the range.

Matrix Splitting
According to the authors of [25], the unique splittings can be used to represent a broad class of iterative methods, thus splittings can be used to represent both outer and inner iterations with its condition of generality. Therefore, in this paper, the generated monotone coefficient matrix splittings for both stationary and nonstationary methods, which can be converged for inner iterations were also demonstrated. In order to conduct a methodical analysis of convergence conditions, matrix splitting is required. Hence, let the matrix splitting of M be expressed as where H 1 and H 2 are given by respectively. Consequently, K r is expressed as The MOWAM iterative method is can be defined as follows: are nonsingular triangular matrices; Ω 1 and Ω 2 are weighted optimal values; ϕ ( ) is an initial vector at the th iteration; meanwhile, ϕ (F) and ϕ (B) are independent solution of forward and backward iterations, respectively. From Equation (13), the Newton-MOWAM iterative scheme for the linear system of Equation (8) is where is a square matrix. Meanwhile, is a coefficient of load vector F. For guaranteed convergence of the MOWAM iteration, the general condition for solving the Equation (8) is proven in the following Theorem.  (14) is convergent for 0 < Ω 1 < 2 and 0 < Ω 2 < 2.

Porous Medium Type Equation Selected Problems
In order to evaluate the efficacy of the Newton-MOWAM iterative method, two examples of one-dimensional porous medium type initial-boundary value problems are selected with different levels of difficulties. Below are the following problems that we consider in the numerical experiment.

Problem 1.
∂w ∂t where A and B are arbitrary constants. The exact solution used for the accuracy checking is given by [7] w(x, t) = ± with arbitrary constants C and D. We chose A = −0.5, B = 1, C = 1, and D = 5 as the specification of Problem 1.

Problem 2. ∂w ∂t
where A and B are arbitrary constants. Exact solution is given by [7] w(x, t) = Bx with an arbitrary constant C and the values that we have selected are A = 1, B = 5, and C = 4.
In order to test the applicability of the Newton-MOWAM method to solve a higher degree of difficulty problem, a two-dimensional porous medium type initial-boundary value problem is proposed. Note that the formulation of a finite difference approximation to a two-dimensional porous medium type equation can be made by extending the formulation in Section 3.2 to another spatial direction. The full formulation of a finite difference approximation to a two-dimensional porous medium type equation can be referred in one of the authors' previous work [27]. Below is the following two-dimensional porous medium type initial-boundary value problem.

Problem 3.
∂w ∂t Exact solution used to test the accuracy of the solution is based on the work in [28]: w(x, y, t) = 4 0.8x + 0.8y + 1.6t.
In executed numerical experiment, the efficiency of the Newton-MOWAM iterative method is compared to Newton-WAM iterative method, which was developed independently by using the combination of Newton method and the existing WAM iterative method. There are two criteria used to measure the difference in terms of efficiency between Newton-MOWAM and Newton-WAM methods: the number of iterations ( max ) and the program execution time recorded in seconds (s). Besides that, the accuracy of the two tested numerical methods is taken into account by recording the maximum values of the absolute errors.

C Program Implementation and Algorithm
The prototype of proposed Newton-MOWAM and the standard Newton-WAM methods are designed and developed by using C programming. The C programming language is used because of its good coding organization and comprehension in the field of numerical analysis. Furthermore, with the basic C programming concept, a correct calculation of the number of iterations and the execution time can be taken. For the comparison analysis, the implementation of the Newton-MOWAM and Newton-WAM on Problems 1-3 is made independently. The number of iterations, the execution time, and the maximum absolute error are recorded using five different orders of a matrix, that is m = 64, 128, 256, 512, and 1024 for Problems 1 and 2, while m × m = 16 × 16, 32 × 32, 64 × 64, 128 × 128, and 256 × 256 for Problem 3. The results from these different orders of a matrix are used for consistency verification. Furthermore, the numerical solution can only be obtained through the running C program after the iteration process was completed successfully. As the C program code is copyrighted, below is the Newton-MOWAM Algorithm 1 used for the program design and development.

Results and Discussion
The numerical results from the simulation of Newton-MOWAM and Newton-WAM methods to solve Problems 1-3 are recorded accordingly from the C program outputs. The collected results such as the number of iterations, the computation time based on the time of completion by the developed C programs, and the values of absolute errors are tabulated in Tables 1-3. Besides that, the percentage of reduction in the number of iterations and the computation time is calculated and showed in Table 4. Based on Tables 1-3, which correspond to the results of Problems 1-3, respectively, the study found that the number of iterations and the simulation execute time required by Newton-MOWAM are notably lower than Newton-WAM for all five different orders of a matrix. The use of two optimum weighted parameters, which have been added to both forward and backward iterations, successfully improved the rate of convergence of the solutions of tested problems. Moreover, the Newton-MOWAM method has minimized the number of iterations tremendously by 15.74% and the computation time by 6.85% for Problem 1, cut down the number of iterations by 63.40% while the computation time by 61.12% for Problem 2 and significantly reduced the number of iterations by 88.24% and the computation time by 80.62% for Problem 3, see Table 4.
In terms of numerical accuracy, the two tested iterative methods are comparable for the selected orders of a matrix. The study found that the magnitude of maximum absolute errors by Newton-MOWAM are smaller than by Newton-WAM, which can be seen at m = 512, 1024 for Problem 1, at m = 128, 256, 512, 1024 for Problem 2 and at m = 64 × 64, 128 × 128, 256 × 256 for Problem 3. Evidently, it can be observed that the magnitude of maximum absolute errors by Newton-MOWAM is decreased further when the orders of a matrix is increased. In the other words, the approximate solution by Newton-MOWAM is converged to the exact solution at a considerably large matrix size. Overall, these findings demonstrated that the efficacy of the proposed Newton-MOWAM method to obtain accurate solutions of nonlinear porous medium type equation with minimal computational complexity.

Conclusions
This paper presented the formulation of a new numerical method called Newton-MOWAM for solving several porous medium type equations. The analysis of efficiency and convergence of the method have also been discussed and supported by the results extracted from the C language-based simulation. The numerical results of all the tested examples showed that the proposed Newton-MOWAM iterative method required least number of iterations with minimum execution time compared to the typical Newton-WAM method. Further, the Newton-MOWAM also helped obtain the results, and it is much closer to the exact solution when the the order of matrix is large. Overall, it can be concluded that the Newton-MOWAM iterative method is certainly better approach in solving nonlinear porous medium type equation after the well-presented symmetry of the finite difference approximation and MOWAM iteration. In the future works, the applicability and efficiency of the Newton-MOWAM method will be further investigated to solve fractional porous medium type equations. Besides that, the Newton-MOWAM method will also be extended to solve more complex equations such as thermo-elastic-viscoplastic medium with pores and thermoelastic medium with double porosity.