Article Search
eISSN 0454-8124
pISSN 1225-6951

### Article

Kyungpook Mathematical Journal 2022; 62(4): 699-713

Published online December 31, 2022

### Numerical Solution of Nonlinear Diffusion in One Dimensional Porous Medium Using Hybrid SOR Method

Jackel Vui Lung Chew∗, Elayaraja Aruchunan, Andang Sunarto, Jumat Sulaiman

Faculty of Computing and Informatics, Universiti Malaysia Sabah Labuan International Campus, Labuan 87000, Malaysia
e-mail : jackelchew93@ums.edu.my

Institute of Mathematical Sciences, Faculty of Science, University of Malaya, Kuala Lumpur 50603, Malaysia
e-mail : elayarajah@um.edu.my

Tadris Matematika, Universitas Islam Negeri (UIN) Fatmawati Sukarno, Bengkulu 38211, Indonesia
e-mail : andangs@iainbengkulu.ac.id

Faculty of Science and Natural Resource, Universiti Malaysia Sabah, Kota Kinabalu 88400, Malaysia
e-mail : jumat@ums.edu.my

Received: May 19, 2021; Revised: May 23, 2022; Accepted: September 8, 2022

### Abstract

This paper proposes a hybrid successive over-relaxation iterative method for the numerical solution of a nonlinear diffusion in a one-dimensional porous medium. The considered mathematical model is discretized using a computational complexity reduction scheme called half-sweep finite differences. The local truncation error and the analysis of the stability of the scheme are discussed. The proposed iterative method, which uses explicit group technique and modified successive over-relaxation, is formulated systematically. This method improves the efficiency of obtaining the solution in terms of total iterations and program elapsed time. The accuracy of the proposed method, which is measured using the magnitude of absolute errors, is promising. Numerical convergence tests of the proposed method are also provided. Some numerical experiments are delivered using initial-boundary value problems to show the superiority of the proposed method against some existing numerical methods.

Keywords: porous medium equation, half-sweep finite difference, Newton method, explicit group, modified successive over-relaxation method

### 1. Introduction

Nonlinear partial differential equations (NPDE) exist in various pure and applied fields. Different types of NPDE have been introduced to the literature to support research on complex physical phenomena. One such equation is called the porous medium equation (PME). PMEs can describe the nonlinear diffusion process in a one-dimensional porous mediums. For instance, it describes the water-oil instability phenomenon where water is injected into an injection well in order to push oil towards the production well [4]. Besides that, PMEs play an important role in the mathematical modelling of the water coning phenomenon in oil reservoirs [16]. The exact solutions of PMEs provide a genuine understanding of the qualitative features of these physical phenomena. Although many exact solutions of PMEs are proposed, few exact solutions give clear meaning to the physical phenomena, which can be used as test examples to verify the accuracy and efficiency of various numerical methods. In some cases, when the modelled PME is difficult to solve by using analytical methods, numerical methods become the alternative way to obtain the solution.

Many different schemes for solving PME have been developed to achieve either a high degree of accuracy or of efficiency. In [11] the authors designed two high degree accuracy finite difference schemes: the sixth and eighth-order weighted essentially non-oscillatory schemes. A linear multi-step approach in their mixed weighted with non-oscillatory finite difference scheme was then introduced in [2]. In [8], an adaptive time mesh technique to formulate a finite difference approximation to solve the PME blow-up problem was proposed. In another article, [13], an explicit finite difference scheme to solve the PME model of interface tracking and hole filling was given. Furthermore, in [7] the authors initiated the application of the half-sweep finite difference (HSFD) scheme of [1] for constructing a family of efficient iterative methods for PMEs. In the present paper, we aim to extend the application of HSFD by developing a new efficient iterative method. Due to the capability of HSFD to reduce the computational complexity of solving a large-sized system of algebraic equations, the HSFD has been applied to solve linear PDEs [12, 17], fuzzy PDEs [3] and fractional PDEs [14, 18]. All reported results were promising for linear problems and motivated the present research to investigate an efficient iterative method based on HSFD to solve nonlinear problems.

This paper extends the previous work [7] by proposing a new hybrid type of successive over-relaxation (SOR) iterative method to solve PME. The local truncation error of the HSFD scheme is investigated to show the consistency of the approximation. The stability of the HSFD scheme is analyzed to determine the stability condition. The proposed iterative method, which uses the explicit group technique [9] and modified SOR iteration [10], is formulated using the HSFD approximation equation. The numerical experiment is delivered using two initial-boundary value problems to show the efficiency and accuracy of the proposed method. The next section discusses the concept and derivation of HSFD approximation to a PME.

### 2.1 HSFD approximation to PME

Let a one-dimensional PME be mathematically generalized in the form of

δuδt=Cδδx(umδuδx),

where C,u(x,t) and um represent the arbitrary constant, the solution, and the nonlinear diffusion term, respectively. The theory and properties of PME are extended from the theory of the linear diffusion or heat equation [20]. Furthermore, the fast and slow diffusions in a one-dimensional porous medium can be characterized based on the range of the value of m. Fast diffusion is characterized by um for m<0, while the slow diffusion is characterized by um for m>0 [21]. The linear diffusion equation can be obtained by simply setting the value of m to equal zero. The solution of Eq. 2.1 can be obtained by imposing the following initial-boundary condition:

I(x)=u(x,0),0x1,

and

B0(t)=u(0,t),B1(t)=u(1,t),0t1,

where I(x),B0(t) and B1(t) are the predetermined functions.

Equation 2.1 can be approximated by the mean of finite differences in its solution domain. By partitioning the space at ph for p=0,1,2,,M and h=1/M, and ime at nk seconds for n=0,1,2,,T and k=1/T seconds, the solution domain has MT mesh points denoted by Up,n=(ph,nk). To discretize Eq. 2.1 using the implicit finite difference scheme, we firstly apply calculus to derive

δuδt=Cumδ2uδx2+Cmum1 δuδx 2.

The partial derivatives in Eq. 2.4 can be approximated using the following finite differences

δuδtUp,nUp,n1k, δuδxUp+1,nUp1,n2h,

and

δ2uδx2Up+1,n2Up,n+Up1,nh2.

Hence, the implicit approximation equation to PME is written as

Up,n1=Up,nαUp,nmUp+1,n+2αUp,nm+1αUp,nmUp1,nβmUp,nm1Up+1,n2+2βmUp,nm1Up+1,nUp1,nβmUp,nm1Up1,n2,

where α=Ck/h2,β=Ck/4h2, and p=1,2,...,M2,M1.

Since the distance between two consecutive points is at a fixed length of h, Eq. 2.8 can be called the full-sweep finite difference (FSFD) approximation. The performance of FSFD to solve PMEs was studied in [5, 6], and it was found that the use of FSFD to solve a large system of equations met with high computational complexity. Using FSFD causes a high number of iterations which prolongs the solution program time. The impact becomes more serious when finer mesh sizes or many mesh points are considered in the computation.

The formulation of the FSFD approximation can be extended to the HSFD by skipping alternate the mesh points. The distance between two consecutive points to be approximated increases from h to 2h when doing so. As a result, the mesh points that need to be computed are grouped into two different indexes: odd and even index. The even index points are approximated using the HSFD iterative method, and the odd index points can be computed using the governing equation derived from implicit finite differences. The use of the HSFD can reduce the number of iterations and solution program times by 50% compared to the FSFD for linear problems [1, 12, 17]. A comparison of the FSFD and HSFD frameworks can be seen in Figure 1.

Figure 1. The FSFD and HSFD for M mesh points

Referring to Figure 1, HSFD approximations to the partial derivatives in Eq. 2.4 can be derived into

δuδtUp,nUp,n1k, δuδxUp+2,nUp2,n4h,

and

δ2uδx2Up+2,n2Up,n+Up2,n4h2.

Thus, the HSFD approximation equation to PME is written as

Up,n1=Up,nαUp,nmUp+2,n+2αUp,nm+1αUp,nmUp2,nβmUp,nm1Up+2,n2+2βmUp,nm1Up+2,nUp2,nβmUp,nm1Up2,n2,

where α=Ck/4h2,β=Ck/16h2, and p=2,4,...,M4,M2.

This paper uses Eq. 2.12 to compute even index points and Eq. 2.8 to compute odd index points at each time level. In addition to HSFD approximation, a hybrid SOR (HSOR) iterative method is formulated to aid the computation process for both index points. A system of algebraic equations corresponding to Eq. 2.1 can be constructed as follows. Let a nonlinear function based on Eq. 2.12 be expressed as

Fp,n=Up,nαUp,nmUp+2,n+2αUp,nm+1αUp,nmUp2,nβmUp,nm1Up+2,n2+2βmUp,nm1Up+2,nUp2,nβmUp,nm1Up2,n2Up,n1,

The function shown by Eq. 2.13 approximates the mesh points p=2,4,,M2, while the remaining points, p=1,3,,M1 are calculated by using the governing equation based on the implicit scheme. By considering M-2 mesh points, a large-sized nonlinear system is formed into

F˜(U˜)=F2,n(U˜ ),F4,n(U˜ ),...,FM2,n(U˜ )T=0,

where U˜=(U2,n,U4,n,...,UM2,n).

The nonlinear system shown by Eq. 2.14 can be solved using the Newton method. First, find the Jacobian matrix from all first-order partial derivatives of Eq. 2.13 as follows:

F˜(U˜)=δ F2,nδU2,nδ F2,nδU4,nδ F2,nδUM2,nδ F4,nδU2,nδ F4,nδU4,nδ F4,nδUM2,nδ FM2,nδU2,nδ FM2,nδU4,nδ FM2,nδUM2,n,

where

δFj,nδUj2,n=αUj,nm+2βmUj,nm1Uj+2,n2βmUj,nm1Uj2,n,j=4,...,M2, δFj,nδUj,n=1αmUj,nm1Uj+2,n+2α(m+1)Uj,nmαmUj,nm1Uj2,nβm(m1)Uj,nm2Uj+2,n2+2βm(m1)Uj,nm2Uj+2,nUj2,nβm(m1)Uj,nm2Uj2,n2,j=2,4,...,M2,

and

δFj,nδUj+2,n=αUj,nm2βmUj,nm1Uj+2,n+2βmUj,nm1Uj2,n,j=2,4,...,M4.

Then, a sparse and large-sized linear system corresponds to Eq. 2.14 can be expressed as

F˜(U˜)W˜=F˜(U˜),

where

W˜= U ˜ (l+1) U ˜ (l),l=0,1,2,....

### 2.2 Local truncation error

In this section, the local truncation error of Eq. 2.12 is investigated by substituting the coefficient Up,n1,Up+2,n and Up2,n with the following Taylor’s expansions:

Up,n1=Up,nkδδtUp,n+k22δ2δt2Up,n+..., Up+2,n=Up,n+2hδδxUp,n+2h2δ2δx2Up,n+...,

and

Up2,n=Up,n2hδδxUp,n+2h2δ2δx2Up,n+....

Substituting Eq. 2.21-2.23 into Eq. 2.12 yields

Up,nαUp,nmUp,n+2hδδxUp,n+2h2δ2δx2Up,n+...+2αUp,nm+1  αUp,nmUp,n2hδδxUp,n+2h2δ2δx2Up,n+...βmUp,nm1 U p,n +2h δ δx U p,n +2 h 2 δ 2 δ x 2 U p,n +... 2+2βmUp,nm1Up,n+2hδδxUp,n+2h2δ2δx2Up,n+...  Up,n2hδδxUp,n+2h2δ2δx2Up,n+...βmUp,nm1 U p,n 2h δ δx U p,n +2 h 2 δ 2 δ x 2 U p,n +... 2  =Up,nkδδtUp,n+k22δ2δt2Up,n+...,

and is then simplified into

δδtUp,n+O(k)=CUp,nmδ2δx2Up,n+CmUp,nm1 δδx Up,n 2+O(h2).

Hence, the local truncation error of Eq. 2.12 is

limk,h0O(k)+O(h2)=0,

which shows that the HSFD scheme is consistent with Eq. 2.1, and the accuracy of the HSFD approximation is first order in time and second order in space.

### 2.3 Analysis of stability

In this section, the analysis of the stability of HSFD approximation is made by using Gerschgorin’s theorem [19] and the following stability theorem is established.

Theorem 2.1. Let Up,n be the numerical approximation to the solution of Eq. 2.1. Then, the HSFD approximation shown by Eq. 2.12 is unconditionally stable and has an eigenvalue

11+λ<1.

Proof. Let the matrix coefficient of the HSFD approximation be expressed in the form of

(I+A)U˜n=U˜n1F˜(U˜n1),

where U˜n=(U2,n,U4,n,...,UM2,n) and U˜n1=(U2,n1,U4,n1,...,UM2,n1).

Define the errors e˜n=(e2,n,e4,n,...,eM2,n) and e˜n1=(e2,n1,e4,n1,...,eM2,n1). Then, we have

e˜n=(I+A)1Ie˜n1.

We must show that every eigenvalue of the matrix A has a positive real part. According to Gerschgorin’s theorem, the eigenvalues of the matrix A are located inside the disks centred at each diagonal entry. The diagonal entry that we have is

Aj,j=1αmUj,nm1Uj+2,n+2α(m+1)Uj,nmαmUj,nm1Uj2,nβm(m1)Uj,nm2Uj+2,n2+2βm(m1)Uj,nm2Uj+2,nUj2,nβm(m1)Uj,nm2Uj2,n2>0,

where α,β>0 and m.

The eigenvalue of the matrix A has a positive real part implies that the matrix A has an eigenvalue λ if, and only if, (I+A)1I has an eigenvalue of 1/(1+λ). We observe that every eigenvalue of the matrix (I+A) has a radius larger than unity which implies that the matrix (I+A) has an inverse. Therefore, we can conclude that the eigenvalue 1/(1+λ) is always smaller than 1 which implies that the HSFD approximation shown by Eq. 2.12 is unconditionally stable.

### 2.4 Derivation of hybrid SOR method

In this section, the proposed HSOR method is derived, and the derivation begins with the explicit group technique [9]. Let's consider a general group of four mesh points as follows: At time level n ≥ 0,

bpcp00ap+2bp+2cp+200ap+4bp+4cp+400ap+6bp+6WpWp+2Wp+4Wp+6=RpRp+2Rp+4Rp+6,

where Rp=Fp(U˜)apWp2,Rp+2=Fp+2(U˜),Rp+4=Fp+4(U˜), and Rp+6=Fp+6(U˜)cp+6Wp+8.

The solution of Eq. 2.31 can be obtained by using

WpWp+2Wp+4Wp+6(l+1)=bpcp00ap+2bp+2cp+200ap+4bp+4cp+400ap+6bp+61RpRp+2Rp+4Rp+6,

wherep=2,10,....

The values of the coefficients aj, bj, and cj, j=p,p+2,p+4, and p+6, are varied at each iteration index. The values are calculated iteratively using the following equations:

aj=αUj,nm+2βmUj,nm1Uj+2,n2βmUj,nm1Uj2,n,j=4,...,M2, bj=1αmUj,nm1Uj+2,n+2α(m+1)Uj,nmαmUj,nm1Uj2,nβm(m1)Uj,nm2Uj+2,n2+2βm(m1)Uj,nm2Uj+2,nUj2,nβm(m1)Uj,nm2Uj2,n2,j=2,4,...,M2,

and

cj=αUj,nm2βmUj,nm1Uj+2,n+2βmUj,nm1Uj2,n,j=2,4,...,M4.

Then, a relaxation factor ω1 is added to the first four mesh points, and another relaxation factor ω2 is added to the next four mesh points and continues alternately for all mesh points in the solution domain. Thus, the HSOR method can be written as follows:

Wp Wp+2 Wp+4 Wp+6(l+1)=(1ω1) Wp Wp+2 Wp+4 Wp+6(l)+ω1 bpcp00 ap+2bp+2cp+20 0ap+4bp+4cp+4 00ap+6bp+61 Rp R p+2 R p+4 R p+6 ,p=2,18,...,

and

Wp Wp+2 Wp+4 Wp+6(l+1)=(1ω2) Wp Wp+2 Wp+4 Wp+6(l)+ω2 bpcp00 ap+2bp+2cp+20 0ap+4bp+4cp+4 00ap+6bp+61 Rp R p+2 R p+4 R p+6 ,p=10,26,...,

where 1<ω1,ω2<2.

The algorithm of the HSOR method to solve PME is shown as follows.

### 3. Numerical Experiment

In this section, the discussion on the efficiency and accuracy of the proposed HSOR method is presented using two different initial-boundary value problems: fast and slow diffusion problems.

Problem 1 [15] A fast diffusion equation problem,

δuδt=cδδx(1u2δuδx).

The exact general solution is given by

u(x,t)=12A1xA12t+A2.

Here we have C=0.5,A1=0.35 and A2=1.35.

Problem 2 [21] A slow diffusion equation problem,

δuδt=δδx(u2δuδx).

The exact general solution is

u(x,t)=x+A12A2t,

and we choose C=1,A1=1 and A2=4.

Table 1a shows the comparison of total iterations and program elapsed time between the proposed method HSOR, the method from [5] and the method from [6]. The total iterations are counted based on the sum of both inner and outer loops (ltotal). The elapsed time (s) is the total time required by the developed simulation program to finish the iteration process to get the final solution. In Table 1a we see that HSOR required the least total iterations to obtain the solution. Using the combination of the HSFD scheme, the Newton method, explicit group techniques, and modified successive over-relaxation, the total iterations can be reduced by 74.55% compared to the method of [5] and 49.84% compated to the method of [6]. Since the total iterations are positively correlated to the elapsed time, the HSOR method is faster than the method of [5] by 75.39% and method of [6] by 60.48%.

Problems 1 and 2 with k=0.01.

(a) Total iterations and elapsed time.

Problem 1Problem 2
M[5][6]HSOR[5][6]HSOR
ltotal25611686103171479796404
5122326119461029261563796
1024457823101194578430551563
204810401482623101151460803055
40962731210673482622449119966080
s2562.261.140.551.641.060.57
5123.842.731.733.682.821.37
102415.518.533.5415.2110.723.69
204856.3234.318.6477.3550.6411.82
4096338.54137.4238.03282.64152.9144.49

(b) Absolute errors.

M
MethodProb.256512102420484096
[5]2.9743e-62.9754e-62.9746e-62.9643e-62.8339e-6
[6]12.9725e-62.9769e-62.9772e-62.9771e-62.9761e-6
HSOR2.9576e-62.9726e-62.9764e-62.9770e-62.9751e-6
[5]8.3876e-58.3876e-58.3876e-58.3876e-58.3876e-5
[6]28.3875e-58.3875e-58.3876e-58.3876e-58.3876e-5
HSOR8.3874e-58.3875e-58.3876e-58.3876e-58.3876e-5

Table 1b compares solutions accuracy between the HSOR method, and the methods of [5] and [6]. The accuracy of the approximate solutions is measured using the magnitude of absolute errors. The paper uses the following formula to calculate the magnitude of absolute errors:

Emax=max1pM1U(xp,tn)u(xp,tn),

where U(xp,tn) and u(xp,tn) are the exact and approximate solutions at the point (xp,tn), respectively. Table 1.b shows an agreement in terms of the accuracy of the approximate solutions between the three tested numerical methods. Since the degree of accuracy between HSFD and FSFD approximations is similar, the proposed HSOR method’s numerical convergence is tested for the various step sizes of h and k. Table 2 shows that the HSOR method is numerically convergent for spatial steps h=1/M=1/256,1/512,1/1024,1/2048, and 1/4096, at time steps k=1/T=1/10,1/100,1/1000, and 1/10000, respectively.

Numerical convergence of HSOR with various h and k.

(a) Problem 1.

k
h1/101/1001/10001/10000
1/2562.9535e-52.9576e-62.7805e-79.9592e-9
1/5122.9545e-52.9726e-62.9297e-72.4825e-8
1/10242.9538e-52.9764e-62.9670e-72.8554e-8
1/20482.9512e-52.9770e-62.9765e-72.9488e-8
1/40962.9452e-52.9751e-62.9787e-72.9729e-8

(b) Problem 2.

k
h1/101/1001/10001/10000
1/2568.1407e-48.3874e-58.4129e-68.4155e-7
1/5128.1408e-48.3875e-58.4130e-68.4155e-7
1/10248.1409e-48.3876e-58.4130e-68.4155e-7
1/20488.1410e-48.3876e-58.4130e-68.4156e-7
1/40968.1423e-48.3876e-58.4130e-68.4156e-7

Table 3 compares computational complexity between HSOR, and the methods of [5] and [6]. The computational complexity of the iterative methods are measured by the number of arithmetic operations per iteration. Table 3 shows that the HSOR method is much more efficient than the other two methods, because it runs fewer arithmetic operations per iteration. By using a total of M-1 mesh points for the numerical simulation, the HSOR uses six operations of PLUS/MINUS and ten operations of MULTIPLY/DIVIDE for the entire groups of four points. Subjected to the boundary conditions, the HSOR uses five operations of PLUS/MINUS, and eight operations of MULTIPLY/DIVIDE for the remaining ungroup points per iteration. Moreover, the implementation of HSOR considers only (M/8-1) points per iteration, which is much smaller than the other two methods with (M-1) and (M/4-1) respectively.

Number of arithmetic operations per iteration.

Operations[5][6]HSOR
PLUS/MINUS4(M-1)6(M/4-1)+56(M/8-1)+5
MULTIPLY/DIVIDE5(M-1)10(M/4-1)+810(M/8-1)+8

The reliability of the numerical test using Problems 1 and 2 can be determined by comparing the exact solutions given by the literature, numerical solutions obtained by the FSFD approximation and numerical solutions obtained by the HSFD approximation. The computed solutions at time level t=1.0s for varied points x are tabulated in Table 4. Figure 2 shows the plotted solutions of Problems 1 and 2.

Computed solutions of Problems 1 and 2 at t=1.0s.

(a) Problem 1.

xExactFSFDHSFD
00.2886751350.2886751350.288675135
0.06250.3067173310.3067497210.306749714
0.1250.3247595260.3248136730.324813662
0.18750.3428017220.3428701800.342870166
0.250.3608439180.3609212850.360921270
0.31250.3788861140.3789683200.378968307
0.3750.3969283100.3970121680.397012156
0.43750.4149705060.4150534230.415053413
0.50.4330127020.4330924960.433092487
0.56250.4510548980.4511296780.451129671
0.6250.4690970940.4691651830.469165178
0.68750.4871392900.4871991740.487199170
0.750.5051814860.5052317800.505231778
0.81250.5232236810.5232631080.523263107
0.8750.5412658770.5412932480.541293247
0.93750.5593080730.5593222780.559322277
10.5773502690.5773502690.577350269

(b) Problem 2.

xExactFSFDHSFD
00.9025873650.9025873650.902587365
0.06250.8869201400.8869211080.886921025
0.1250.8720414400.8720431360.872042992
0.18750.8578872660.8578894870.857889302
0.250.8444006620.8444032380.844403027
0.31250.8315307500.8315335380.831533312
0.3750.8192319210.8192347980.819234569
0.43750.8074631510.8074660140.807465789
0.50.7961874290.7961901890.796189976
0.56250.7853712620.7853738440.785373649
0.6250.7749842580.7749866000.774986426
0.68750.7649987690.7650008140.765000665
0.750.7553895750.7553912780.755391156
0.81250.7461336200.7461349430.746134850
0.8750.7372097810.7372106900.737210627
0.93750.7285986590.7285991250.728599093
10.7202824060.7202824060.720282406

Figure 2. Plotted solutions to Problems 1 and 2

Figure 2 shows a strong agreement between the plotted exact solutions and numerical solutions by FSFD and HSFD approximations. Thus, the comparison among the plotted solutions suggests that the numerical test using Problems 1 and 2 is reliable.

### 7. Conclusion

This paper proposed a hybrid SOR method based on the HSFD scheme to solve one-dimensional nonlinear PME. The formulated HSFD approximation is consistent and accurate in first-order in time and second order in space. The HSFD approximation is also unconditionally stable according to Gerschgorin’s theorem. The numerical results illustrated that the proposed hybrid SOR method is superior to the existing methods that used FSFD approximation in total iterations and elapsed time. The numerical convergence of the proposed hybrid SOR method for various spatial and temporal step sizes is also presented.

### References

1. A. R. Abdullah, The four point explicit decoupled group (EDG) Method: a fast Poisson solver, Int. J. Comput. Math., 38(1-2)(1991), 61-70.
2. R. Abedian, H. Adibi and M. Dehghan, A high-order weighted essentially non-oscillatory (WENO) finite difference scheme for nonlinear degenerate parabolic equations, Comput. Phys. Commun., 184(8)(2013), 1874-1888.
3. L. H. Ali, J. Sulaiman and S. R. M. Hashim, Numerical solution of fuzzy Fredholm integral equations of second kind using half-sweep Gauss-Seidel iteration, J. Eng. Sci. Technol., 15(5)(2020), 3303-3313.
4. R. Borana, V. Pradhan and M. Mehta, Numerical solution of instability phenomenon arising in double phase flow through inclined homogeneous porous media, Perspect Sci., 8(2016), 225-227.
5. J. V. L. Chew and J. Sulaiman. Application of MSOR iteration with Newton scheme for solutions of 1D nonlinear porous medium equations. AIP Conference Proceedings. article no. 20017. doi: 10.1063/1.4952497.
6. J. V. L. Chew and J. Sulaiman, Implicit solution of 1D nonlinear porous medium equation using the four-point Newton-EGMSOR iterative method, J. Appl. Comput. Mech., 15(2)(2016), 11-21.
7. J. V. L. Chew and J. Sulaiman, Solution of one-dimensional porous medium equation using half-sweep Newton-MSOR iteration, Adv. Sci. Lett., 24(3)(2018), 1906-1911.
8. C. H. Cho, On the finite difference approximation for blow-up solutions of the porous medium equation with a source, Appl. Numer. Math., 65(2013), 1-26.
9. D. J. Evans, Group explicit iterative methods for solving large linear systems, Int. J. Comput. Math., 17(1)(1985), 81-108.
10. D. R. Kincaid and D. M. Young, The modified successive overrelaxation method with fixed parameters, Math. Comp., 26(119)(1972), 705-717.
11. Y. Liu, C. W. Shu and M. Zhang, High order finite difference WENO schemes for nonlinear degenerate parabolic equations, SIAM J. Sci. Comput., 33(2)(2011), 939-965.
12. N. A. Mat Ali, R. Rahman, J. Sulaiman and K. Ghazali, Solutions of reaction-diffusion equations using similarity reduction and HSSOR iteration, Indones. J. Electr. Eng. Comput. Sci., 16(3)(2019).
13. L. Monsaingeon, An explicit finite-difference scheme for one-dimensional generalized porous medium equations: interface tracking and the hole filling problem, ESAIM Math. Model. Numer. Anal., 50(4)(2016), 1011-1033.
14. F. A. Muhiddin, J. Sulaiman and A. Sunarto, Grunwald implicit solution of one-dimensional time-fractional parabolic equations using HSKSOR iteration, J. Phys. Conf. Ser., 1489(2020).
15. A. D. Polyanin and V. F. Zaitsev. Handbook of nonlinear partial differential equations. Boca Raton: Chapman & Hall; 2004.
16. M. Safari, M. J. Ameri and A. Naderifar, An efficient boundary control for porous media equation: motivated by water coning problem, Can. J. Chem. Eng., 97(2019), 888-902.
17. A. Saudi and J. Sulaiman, An efficient two-stage half-sweep modified arithmetic mean (HSMAM) method for the solution of 2D elliptic equation, Adv. Sci. Lett., 24(3)(2018), 1917-1921.
18. A. Sunarto and J. Sulaiman, Performance numerical method half-sweep preconditioned Gauss-Seidel for solving fractional diffusion equation, Math. Model. Eng. Probl., 7(2)(2020), 201-204.
19. R. S. Varga. Geršgorin and His Circles. Berlin: Springer; 2004.
20. J. L. Vazquez. The Porous Medium Equation. Oxford University Press; 2006.
21. A. M. Wazwaz, Exact solutions to nonlinear diffusion equations obtained by the decomposition method, Appl. Math. Comput., 123(1)(2001), 109-122.