Implicit solution of 1D nonlinear porous medium equation using the fourpoint NewtonEGMSOR iterative method
J.V.L. Chew
,J. Sulaiman
Journal of Applied Mathematics and Computational Mechanics 
Download Full Text 
IMPLICIT SOLUTION OF 1D NONLINEAR POROUS MEDIUM EQUATION USING THE FOURPOINT NEWTONEGMSOR ITERATIVE METHOD
J.V.L. Chew, J. Sulaiman
Faculty of Science and Natural Resources,
University Malaysia Sabah
Sabah, Malaysia
jackelchew93@gmail.com, jumat@ums.edu.my
Abstract. The numerical method can be a good choice in solving nonlinear partial differential equations (PDEs) due to the difficulty in finding the analytical solution. Porous medium equation (PME) is one of the nonlinear PDEs which exists in many realistic problems. This paper proposes a fourpoint NewtonEGMSOR (4NewtonEGMSOR) iterative method in solving 1D nonlinear PMEs. The reliability of the 4NewtonEGMSOR iterative method in computing approximate solutions for several selected PME problems is shown with comparison to 4NewtonEGSOR, 4NewtonEG and NewtonGaussSeidel methods. Numerical results showed that the proposed method is superior in terms of the number of iterations and computational time compared to the other three tested iterative methods.
Keywords: porous medium equation, finite difference scheme, Newton method, Explicit Group, MSOR
1. Introduction
Most of the nonlinear partial differential equations (PDEs) are difficult to be solved by the analytical approach. Since the rapid advancement of computers, numerical methods have grown drastically to be the choice in solving nonlinear PDEs. This paper considers a numerical approach for the solution of one of the nonlinear PDEs which is a porous medium equation (PME). This equation has great practicality in many realistic problems such as a thin film flow through a porous medium [1], diffusion of heat beneath human skin [2] and some interesting applications such as the dispersion of miscible fluid through a porous medium [3, 4] and the instability phenomena in the oil recovery technology [57]. Currently, the last two applications have been widely studied for the development of a theoretical model and finding the approximate analytical solution to the developed model. For instance, Meher et al. [3, 4] have discussed the dispersion phenomena inside porous media that occurs in oil reservoir and applied the Adomian decompo sition method and the Backlund transformation in obtaining the analytical solutions of the problems. Meher et al. [5, 6] have applied the Adomian decomposition method and the exponential self similar solutions technique in solving PME which is formulated from the instability phenomena in double phase flow through porous media. In fact, there is much literature involving the solution of PME problems that can be found in [712]. Motivated and inspired by the ongoing research in solving PME, this paper proposes a fourpoint NewtonEGMSOR (4NewtonEGMSOR) iterative method in solving 1D nonlinear PMEs. Actually, the method is a combination of fourpoint EGMSOR iteration which is initiated by Sulaiman et al. [13] with the Newton method that is used to handle the nonlinearity of the problem. This paper dealt with an efficient numerical technique that can reduce the computational time while maintaining the accuracy of the approximate solution of PMEs.
In this paper, to secure computational stability, a PME approximation equation is developed by using the implicit finite difference scheme. A system of nonlinear equations is formed at each time level. The Newton method is then used to linearize and transform the developed system of nonlinear equations into the corresponding system of linear equations. The resultant linear system is finally solved by the fourpoint EGMSOR iterative methods. Four examples of the PMEs are chosen in order to illustrate the capability of the 4NewtonEGMSOR iterative method. The reliability of the proposed method in computing approximate solutions for the selected PME problems is shown with comparison to 4NewtonEGSOR, 4NewtonEG and NewtonGaussSeidel (GS) iterative methods. Consider a general form of the PME problem be defined as
(1) 
where K and m are real numbers.
Solution domain x can be divided uniformly into d subintervals with distance Dx. Time step Dt can also be obtained by dividing the total time T at fixed sizes of s. Both steps, Dx and Dt, can be defined as
(2) 
2. Implicit finite difference approximation equation
To formulate 4NewtonEGMSOR for solving Eq. (1), a finite grid network is built as a guide for development of the three iterative methods and facilitating the implementation of computational algorithm, refer Figure 1.
Fig. 1. Finite grid network
The implementation of the 4NewtonEGMSOR method will be applied onto the interior grid points in Figure 1, i.e. grid point 1 to n, until the convergence of approximate solutions is achieved. Before discretizing, problem (1) can be rewritten as
(3) 
Now, by using the implicit finite difference scheme, (3) is discretized to derive an approximation equation as
(4) 
for i = 1,_{ }2,_{ }3,_{ }...,_{ }n and j = 0,1,_{ }2,_{ }3,_{ }...,_{ }s, where
Eq. (4) is obviously a nonlinear finite difference approximation equation that is used to form a system of nonlinear equations at each time level j. To apply the Newton method on Eq. (4), define a nonlinear function F for each interior grid point (x_{i}_{ }, t_{j}_{+1}) at time level t_{j}_{+1} as follows:
(5) 
where .
By considering all interior grid points in Figure 1, Eq. (5) leads to nonlinear system as
(6) 
Now, the Newton method can be used to calculate the Jacobian matrix on (6) in order to construct the corresponding linear system. Indicated with kth numerical solutions at each time level t_{j}_{+1}, Eq. (6) is transformed into a linear system as
(7) 
where
The unknown vector_{ }_{ }is required to be solved so that the vector of approximate solutions_{ }_{ }can be computed iteratively by using the following expression:
(8) 
3. Formulation and implementation of the fourpoint EGMSOR method
As discussed in the second section, the coefficient matrix A in linear system (7) is sparse and largescaled, and therefore needs an efficient iterative method to be solved. A number of iterative methods can be found in [1417]. In addition to that, Evans [18] has proposed a fourpoint block iterative method which is also known as Explicit Group (EG) in solving large linear systems. Again, the Successive Over Relaxation (SOR) method was introduced by Young [15] and it is the most known and widely used iterative technique in solving a system of linear algebraic equations. Due to the advantage of the SOR method, Kincaid and Young [17] have suggested Modified Successive Over Relaxation (MSOR) method, which is classified as a point iterative method together with two weighted parameters, w_{r} and w_{b}, in order to speed up the rate of convergence in SOR. Therefore, this paper considers the application of the fourpoint EGMSOR iterative method for solving the generated linear system and this iterative method is a combination between the EG and MSOR methods. Before formulating the fourpoint EGMSOR iterative method, let the linear system in (7) be rewritten in general form as
(9) 
where
Now, the coefficient matrixin Eq. (9) needs to be decomposed as
(10) 
where D, L and V are the diagonal, lower triangular part and upper triangular part of matrix A, respectively. Then, the formulation of the SOR method is given by [15]
(11) 
When w = 1, the SOR method can be reduced to the standard GS iterative method. Apart from the concept of the SOR method, the MSOR methods can be derived from (11) with two weighted parameters, w_{r} and w_{b}. In fact, this concept is similar to the SOR method with the redblack ordering strategy. The formulation of the MSOR method is given by [17]
(12) 
And then when w_{r} = w_{b}, the MSOR method reduces to the redblack SOR method. Besides that, by setting w_{r} = w_{b} = 1, (12) will become the redblack GS method. Since this paper uses fourpoint EGMSOR to solve the linear system that is transformed by using the Newton method, the derivation of fourpoint EGMSOR method will be constructed over the linear system (9) as follows.
Referring to Figure 1, a group of four points is considered to be used to form a linear system_{ }_{ }as follows.
(13) 
where
for t = 1,_{ }2,_{ }3,_{ }4. Eq. (13) can be easily inverted and derived into a fourpoint EGMSOR formula that is similar to Eq. (12) as follows:
(14a) 
(14b) 
Hence, by using Eqs. (14a) and (14b), the fourpoint EGMSOR algorithm can be given in Algorithm 1.
Algorithm 1. Fourpoint EGMSOR
i. Initialize
ii. For_{ }, implement
a. Set_{ }
b. Calculate_{ }, A and b
c. For_{ }, compute Eq. (14a) iteratively.
d. For_{ }, compute Eq. (14b) iteratively.
e. For_{ }, compute ungroup points. For more detail, see in [19].
iii. Convergence test_{ }. If converges, go to (iv). Otherwise, back to (ii).
iv. Compute
v. Convergence test_{ }. If converges, go to next time level. Otherwise, back to (i).
vi. Display approximate solutions.
The estimate values of the two weighted parameters are determined within range ± 0.01 by running Algorithm 1 with a different combination of w_{r} and w_{b}. The combination that gives the least number of iterations will be selected as optimal values.
4. Numerical experiments
In order to verify the effectiveness of 4NewtonEGMSOR iterative method in solving (1), four selected PME problems are used together with other three tested iterative methods, i.e. 4NewtonEGSOR, 4NewtonEG and NewtonGS iterative methods, that act as comparison to the proposed iterative method. For the comparison purpose, three criteria will be considered, namely the number of iterations (Iter), execution time, which is measured in seconds (Time), and maximum absolute error (MAE). In addition to that, tolerance error for convergence is e = 10^{–10}. Below is the following for four examples of PME problems.
Example 1 [19]
Consider K = 1 and m = 1 in (1) gives
(15) 
subject to the initial condition, , and satisfies the exact solution, ,_{ }where_{ }_{ }and_{ }_{ }are arbitrary constants. For the numerical implementation, constants are set to be_{ }_{ }and_{ }.
Example 2 [19]
Consider K = 0.5 and m = –1 in (1) gives
(16) 
subject to the initial condition, , and satisfies the exact solution, ,_{ }where and are arbitrary constants. To implement the iterative methods, constants chosen are_{ }_{ }and_{ }.
Example 3 [8]
Consider K = 1 and m = 2 in (1) gives
(17) 
subject to the initial condition, ,_{ }and satisfies the exact solution, , ,_{ }where C is an arbitrary constant. For the imple mentation purpose, the constant is C = 2.
Example 4 [19]
Consider K = 0.5 and m = –2 in (1) gives
(18) 
subject to the initial condition, , and satisfies the exact solu tion, , where_{ }_{ }and_{ }_{ }are arbitrary constants. To test the numerical methods, constants are set to be_{ }_{ }and_{ }.
The numerical results obtained have been summarized in Tables 13.
Table 1
Comparison of the number of iterations (Iter), execution time in seconds (Time) and maximum absolute errors (MAE) for the iterative methods using Examples 1 and 2


Example 1 
Example 2 

n 
Method 

Iter 
Time 
MAE 

Iter 
Time 
MAE 
64 
GS 

3835 
2.38 
2.76E08 

1720 
1.13 
2.03E05 

EG 

1109 
0.34 
5.88E09 

504 
0.55 
2.03E05 

EGSOR 
w = 1.59 
263 
0.12 
3.43E11 
w = 1.47 
175 
0.25 
2.03E05 

EGMSOR 
w_{r} = 1.59 w_{b} = 1.58 
240 
0.10 
2.99E11 
w_{r} = 1.47 w_{b} = 1.48 
147 
0.16 
2.03E05 
128 
GS 

13_{ }678 
7.50 
1.22E07 

6034 
4.06 
2.02E05 

EG 

3899 
1.65 
2.64E08 

1718 
1.73 
2.03E05 

EGSOR 
w = 1.76 
513 
0.34 
4.28E11 
w = 1.68 
337 
0.41 
2.03E05 

EGMSOR 
w_{r} = 1.76 w_{b} = 1.77 
455 
0.31 
1.61E11 
w_{r} = 1.68 w_{b} = 1.69 
275 
0.38 
2.03E05 
256 
GS 

48_{ }395 
38.58 
5.33E07 

20_{ }907 
27.03 
2.00E05 

EG 

13_{ }799 
11.20 
1.10E07 

5976 
11.25 
2.02E05 

EGSOR 
w = 1.87 
1010 
1.37 
6.90E11 
w = 1.82 
656 
1.35 
2.03E05 

EGMSOR 
w_{r} = 1.87 w_{b} = 1.88 
879 
1.20 
3.62E11 
w_{r} = 1.82 w_{b} = 1.83 
532 
1.14 
2.03E05 
512 
GS 

169_{ }693 
252.94 
2.10E06 

71_{ }385 
287.34 
1.93E05 

EG 

48_{ }666 
77.31 
4.99E07 

20_{ }701 
97.75 
2.00E05 

EGSOR 
w = 1.93 
2027 
5.33 
7.69E10 
w = 1.91 
1297 
4.46 
2.03E05 

EGMSOR 
w_{r} = 1.93 w_{b} = 1.94 
1680 
4.51 
4.23E11 
w_{r} = 1.90 w_{b} = 1.91 
1050 
3.82 
2.03E05 
1024 
GS 

587_{ }031 
1712.49 
7.62E06 

239_{ }975 
1741.01 
1.72E05 

EG 

170_{ }300 
557.86 
2.08E06 

70_{ }888 
571.03 
1.94E05 

EGSOR 
w = 1.97 
4072 
22.43 
1.54E10 
w = 1.95 
2477 
24.67 
2.03E05 

EGMSOR 
w_{r} = 1.96 w_{b} = 1.97 
3417 
19.72 
1.18E10 
w_{r} = 1.95 w_{b} = 1.95 
2078 
22.87 
2.03E05 
Table 2
Comparison of the number of iterations (Iter), execution time in seconds (Time) and maximum absolute errors (MAE) for the iterative methods using Examples 3 and 4


Example 3 
Example 4 

n 
Method 

Iter 
Time 
MAE 

Iter 
Time 
MAE 
64 
GS 

1344 
1.17 
8.39E05 

2015 
1.26 
2.88E06 

EG 

402 
0.38 
8.39E05 

592 
0.53 
2.89E06 

EGSOR 
w = 1.35 
216 
0.13 
8.39E05 
w = 1.49 
191 
0.30 
2.90E06 

EGMSOR 
w_{r} = 1.35 w_{b} = 1.30 
202 
0.13 
8.39E05 
w_{r} = 1.49 w_{b} = 1.49 
165 
0.17 
2.90E06 
128 
GS 

4824 
2.84 
8.39E05 

7082 
4.90 
2.90E06 

EG 

1361 
1.00 
8.39E05 

2033 
1.85 
2.94E06 

EGSOR 
w = 1.58 
437 
0.41 
8.39E05 
w = 1.69 
380 
0.89 
2.96E06 

EGMSOR 
w_{r} = 1.59 w_{b} = 1.54 
405 
0.36 
8.39E05 
w_{r} = 1.69 w_{b} = 1.70 
316 
0.44 
2.96E06 
256 
GS 

17_{ }308 
20.03 
8.39E05 

24_{ }325 
45.42 
2.71E06 

EG 

4836 
6.77 
8.39E05 

7007 
15.11 
2.92E06 

EGSOR 
w = 1.74 
873 
1.47 
8.39E05 
w = 1.83 
734 
1.98 
2.97E06 

EGMSOR 
w_{r} = 1.74 w_{b} = 1.74 
810 
1.21 
8.39E05 
w_{r} = 1.83 w_{b} = 1.84 
609 
1.42 
2.97E06 
512 
GS 

61_{ }658 
270.11 
8.40E05 

81_{ }729 
354.79 
1.86E06 

EG 

17_{ }333 
46.85 
8.39E05 

23_{ }769 
112.83 
2.73E06 

EGSOR 
w = 1.86 
1718 
5.38 
8.38E05 
w = 1.91 
1428 
5.82 
2.98E06 

EGMSOR 
w_{r} = 1.87 w_{b} = 1.86 
1563 
4.38 
8.39E05 
w_{r} = 1.91 w_{b} = 1.91 
1202 
4.76 
2.98E06 
1024 
GS 

218_{ }147 
2008.35 
8.43E05 

265_{ }698 
2293.23 
3.33E06 

EG 

61_{ }779 
342.02 
8.40E05 

79_{ }057 
733.85 
1.89E06 

EGSOR 
w = 1.93 
3344 
20.49 
8.39E05 
w = 1.95 
2881 
26.41 
2.97E06 

EGMSOR 
w_{r} = 1.93 w_{b} = 1.92 
3066 
16.83 
8.39E05 
w_{r} = 1.95 w_{b} = 1.96 
2311 
19.75 
2.98E06 
Table 3
Reduction in percentages of the iterative methods compared with the GS method
M 
Method 
Number of iterations 
Execution time in seconds 

Example 1 
Example 2 
Example 3 
Example 4 
Example 1 
Example 2 
Example 3 
Example 4 

64 
EG 
71.08% 
70.70% 
70.09% 
70.62% 
85.71% 
51.33% 
67.52% 
57.94% 

EGSOR 
93.14% 
89.83% 
83.93% 
90.52% 
94.96% 
77.88% 
88.89% 
76.19% 

EGMSOR 
93.74% 
91.45% 
84.97% 
91.81% 
95.80% 
85.84% 
88.89% 
86.51% 
128 
EG 
71.49% 
71.53% 
71.79% 
71.29% 
78.00% 
57.39% 
64.79% 
62.24% 

EGSOR 
96.25% 
94.41% 
90.94% 
94.63% 
95.47% 
89.90% 
85.56% 
81.84% 

EGMSOR 
96.67% 
95.44% 
91.60% 
95.54% 
95.87% 
90.64% 
87.32% 
91.02% 
256 
EG 
71.49% 
71.42% 
72.06% 
71.19% 
70.97% 
58.38% 
66.20% 
66.73% 

EGSOR 
97.91% 
96.86% 
94.96% 
96.98% 
96.45% 
95.01% 
92.66% 
95.64% 

EGMSOR 
98.18% 
97.46% 
95.32% 
97.50% 
96.89% 
95.78% 
93.96% 
96.87% 
512 
EG 
71.32% 
71.00% 
71.89% 
70.92% 
69.44% 
65.98% 
82.66% 
68.20% 

EGSOR 
98.81% 
98.18% 
97.21% 
98.25% 
97.89% 
98.45% 
98.01% 
98.36% 

EGMSOR 
99.01% 
98.53% 
97.47% 
98.53% 
98.22% 
98.67% 
98.38% 
98.66% 
1024 
EG 
70.99% 
70.46% 
71.68% 
70.25% 
67.42% 
67.20% 
82.97% 
68.00% 

EGSOR 
99.31% 
98.97% 
98.47% 
98.92% 
98.69% 
98.58% 
98.98% 
98.85% 

EGMSOR 
99.42% 
99.13% 
98.59% 
99.13% 
98.85% 
98.69% 
99.16% 
99.14% 
5. Conclusion
In this paper, the effectiveness of the 4NewtonEGMSOR in solving 1D nonlinear PMEs as compared with other three tested iterative methods, i.e. 4Newton‑EGSOR, 4NewtonEG and NewtonGS iterative methods, have been demonstrated by using four PME problems. The numerical results presented in Tables 1 and 2 showed that the proposed iterative method requires a much lesser number of iterations and computational time in obtaining numerical solutions as compared to the other three tested iterative methods. By using the NewtonGS as a control method, 4NewtonEGMSOR has a reduced number of iteration of approximately 84.97 ‑99.42% and computational time approximately 85.8499.16%, see in Table 3. The performance showed by the 4NewtonEGMSOR is aided by the use of two weighted parameters, w_{r} and w_{b}. When these two weighted parameters achieved their optimal choices, the maximum speed of convergence in solving the PME problems can be reached. And in terms of the accuracy of the iterative methods, all four tested numerical methods have good agreement. Therefore, it can be concluded that the 4NewtonEGMSOR iterative method can be a promising technique for solving different types of nonlinear partial differential equation problems. In this paper, however, all four iterative methods can be classified under a family of fullsweep iterative methods. Hence, for future work, this study will be extended for a family of halfsweep iterative methods [2124].
References
[1] Mahmood T., Khan N., Thin film flow of a third grade fluid through porous medium over an inclined plane, International Journal of Nonlinear Science 2012, 14(1), 5359.
[2] Khaled A.R.A., Vafai K., The role of porous media in modelling flow and heat transfer in biological tissues, International Journal of Heat and Mass Transfer 2003, 46, 49895003.
[3] Meher R., Mehta M.N., Meher S.K., Adomian decomposition method for dispersion phenomena arising in longitudinal dispersion of miscible fluid flow through porous media, Advances in Theoretical and Applied Mechanics 2010, 3(5), 211220.
[4] Meher R.K., Meher S.K., Mehta M.N., A new approach to Backlund transformation for longitudinal dispersion of miscible fluid flow through porous media in oil reservoir during secondary recovery process, Theoretical and Applied Mechanics 2011, 38(1), 116.
[5] Meher R.K., Meher S.K., Analytical treatment and convergence of the Adomian decomposition method for instability phenomena arising during oil recovery process, International Journal of Engineering Mathematics 2013, Article ID752561.
[6] Meher R.K., Meher S.K., Mehta M.N., Exponential self similar solutions technique for instability phenomenon arising in double phase flow through porous medium with capillary pressure, Applied Mathematics Sciences 2010, 4(27), 13291335.
[7] Pradhan V.H., Mehta M.N., Patel T., A numerical solution of nonlinear equation representing onedimensional instability phenomena in porous media by finite element technique, Inter national Journal of Advanced Engineering Technology 2011, 2(1), 221227.
[8] Elsheikh A.M., Elzaki T.M., Variation iteration method for solving porous medium equation, International Journal of Development Research 2015, 5(6), 46774680.
[9] Wazwaz A.M., The variational iteration method: A powerful scheme for handling linear and nonlinear diffusion equations, Computers and Mathematics with Application 2007, 54, 933939.
[10] Bhadane P.K.G., Pradhan V.H., Elzaki transform homotopy pertubation method for solving porous medium equation, International Journal of Research in Engineering and Technology 2013, 2(12), 116119.
[11] Elzaki T.M., Hilal E.M.A., Homotopy perturbation and Elzaki transform for solving nonlinear partial differential equations, Mathematical Theory and Modeling 2012, 2(3), 3342.
[12] Pamuk S., On the solution of the porous media equation by decomposition method: A review, Physics Letters A 2005, (344), 184188.
[13] Sulaiman J., Hasan M.K., Othman M., Karim S.A.A., NewtonEGMSOR methods for solution of second order twopoint nonlinear boundary value problems, Journal of Mathematics and System Science 2012, 2, 185190.
[14] Saad Y., Iterative Methods for Sparse Linear Systems, 2nd ed., Society for Industrial and Applied Mathematics, Philadelphia 2003.
[15] Young D.M., Iterative methods for solving partial difference equations of elliptic type, Transactions of the American Mathematical Society 1954, 76, 92111.
[16] Young D.M., Iterative Solution of Large Linear Systems, Academic Press, Florida 2014.
[17] Kincaid D.R., Young D.M., The modified successive over relaxation method with fixed parameters, Mathematics of Computation 1972, 26(119), 705717.
[18] Evans D.J., Group explicit iterative methods for solving large linear systems, International Journal of Computer Mathematics 1985, 17, 81108.
[19] Evans D.J., Abdullah A.R., A new explicit method for the diffusionconvection equation, Computers and Mathematics with Applications 1985, 11, 13.
[20] Polyanin A.D., Zaitsev V.F., Handbook of Nonlinear Partial Differential Equation, Chapman and Hall, Boca Raton 2004.
[21] Akhir M.K.M., Othman M., Sulaiman J., Majid Z.A., Suleiman M., Halfsweep modified successive over relaxation for solving twodimensional Helmholtz equations, Australian Journal of Basic and Applied Sciences 2011, 5(12), 30333039.
[22] Othman M., Sulaiman J., Abdullah A.R., A parallel halfsweep multigrid algorithm on the shared memory multiprocessors, Malaysian Journal of Computer Science 2000, 13(2), 16.
[23] Muthuvalu M.S., Sulaiman J., Halfsweep geometric mean method for solution of linear Fredholm equations, Matematika 2008, 24(1), 7584.
[24] Aruchunan E., Sulaiman J., Halfsweep conjugate gradient method for solving first order linear Fredholm integrodifferential equations, Australian Journal of Basic and Applied Sciences 2011, 5, 3843.