# Solution of the Porous Media Equation by  a Compact Finite Difference Method

Authors: Murat Sari
Journal: Mathematical Problems in Engineering
 (2009)

---

## Abstract

Accurate solutions of the porous media equation that usually occurs in nonlinear problems of heat and mass transfer and in biological systems are obtained using a compact finite difference method in space and a low-storage total variation diminishing third-order Runge-Kutta scheme in time. In the calculation of the numerical derivatives, only a tridiagonal band matrix algorithm is encountered. Therefore, this scheme causes to less accumulation of numerical errors and less use of storage space. The computed results obtained by this way have been compared with the exact solutions to show the accuracy of the method. The approximate solutions to the equation have been computed without transforming the equation and without using linearization. Comparisons indicate that there is a very good agreement between the numerical solutions and the exact solutions in terms of accuracy. This method is seen to be a very good alternative method to some existing techniques for such realistic problems.


---

## Body

## 1. Introduction

Most physical phenomena and processes encountered in various fields of
science are governed by partial differential equations. The nonlinear heat equation describing various physical phenomena,
for instance see references [1–4], is called the porous media equation (1.1)ut=∂∂x(umux)+bu with the
initial condition (1.2)u(x,0)=f(x) and the
boundary conditions (1.3)u(0,t)=g1(t),u(1,t)=g2(t).Here m is a rational number while b is real constant. In (1.1)–(1.3), x and t denote derivatives with respect to space and
time, respectively, when they are used as subscripts.Computing the solutions that have a physical or biological interpretation for the
nonlinear equations of the form (1.1) is of essential importance. This equation is often encountered in nonlinear
problems of heat and mass transfer,
combustion theory, and flows in porous media. For instance, it illustrates
unsteady heat transfer in a quiescent medium with the heat diffusivity being a power-law function of temperature [5].Equation (1.1)
has also applications to many physical systems including the fluid dynamics of
thin films [6]. How this model has been used
to represent population pressure in biological systems was analyzed by Murray [7].
Equation (1.1) is known as a degenerate
parabolic differential equation as the
diffusion term D(u)=um does not satisfy the condition for
classical diffusion equations, D(u)>0 [6].For the motion of thin viscous films,
(1.1) with m=3 can be obtained from the Navier-Stokes equations.
Lacking a physical law to describe the
complex behavior in a system, an appropriate value for the parameter m can be determined by comparing known solutions
with empirical data [6]. A number of analytical and numerical methods are
available in the literature for the investigation of the solution of the
corresponding equation [4, 8–13].In this work, it is aimed to effectively employ a
combination of a sixth-order compact finite difference (CFD6) scheme in space
[14–16] and a low-storage total variation diminishing third-order Runge-Kutta
(TVD-RK3) scheme in time [17, 18] to obtain accurate solutions of (1.1). Since
the TVD-RK3 is an explicit time integration scheme, there is no need for linearization of the nonlinear terms.


## 2. The Compact
Finite Difference Scheme

The compact
finite difference schemes can be dealt within two essential categories:
explicit compact and implicit compact approaches. Whilst the first category
computes the numerical derivatives directly at each grid by using large
stencils, the second one obtains all the numerical derivatives along a grid
line using smaller stencils and solving a linear system of equations. Because
of the reasons given previously, the present work uses the second approach. To
attain the solutions of the equation,
discretizations are needed in both space and time.


### 2.1. Discretization

Spatial
derivatives are evaluated by the compact finite difference scheme [14]. For any
scalar pointwise value u,
the derivatives of u can be obtained by solving a tridiagonal
or pentadiagonal system. Much work was done in deriving such formulae [14, 19].
More details on the formulae can be found in [14, 19]. A uniform one-dimensional
mesh is considered, consisting of N points: x1,x2,…,xi−1,xi,xi+1,…,xN. The mesh size is denoted by h=xi+1−xi. The first derivatives can be given at interior nodes as follows [14, 19]: (2.1)αui−1′+ui′+αui+1′=bui+2−ui−24h+aui+1−ui−12h which gives
rise to an α-family of fourth-order tridiagonal schemes
with (2.2)a=23(α+2),b=13(4α−1),where α=0 leads to the explicit fourth-order scheme for
the first derivative. A sixth-order tridiagonal scheme is obtained by α=1/3, (2.3)13ui−1′+ui′+13ui+1′=19ui+2−ui−24h+149ui+1−ui−12h.For the nodes
near the boundary, approximation formulae for the derivatives of nonperiodic
problems can be derived with the consideration of one-sided schemes. More
details on the derivations for the first- and second-order derivatives can be
found in [14, 19]. The derived  formulae at boundary points 1, 2, N−1, and N,
respectively, are as follows:(2.4)ui′+5ui+1′=1h(−19760ui−512ui+1+5ui+2−53ui+3+512ui+4−120ui+5),211ui−1′+ui′+211ui+1′=1h(−2033ui−1−35132ui+3433ui+1−733ui+2+233ui+3−1132ui+4),211ui−1′+ui′+211ui+1′=1h(2033ui+1+35132ui−3433ui−1+733ui−2−233ui−3+1132ui−4),5ui−1′+ui′=1h(19760ui+512ui−1−5ui−2+53ui−3−512ui−4+120ui−5). In order to
obtain the formulae presented above, the procedure of Gaitonde and Visbal [19]
was followed. The formulae can be re-expressed as (2.5)BU′=AU, where U=(u1,…,un)T.
The second-order derivative terms are obtained by applying the first-order
operator twice, that is, (2.6)BU′′=AU′, where
(2.7)A=1h[−19760−5125−53512−120−2033−351323433−733233−1132−136−79079136⋱⋱⋱⋱⋱−136−790791361132−233733−3433351322033120−51253−551219760]N×N,B=[15211121113113⋱⋱⋱13113211121151]N×N.   In the calculation
of the numerical derivatives, a tridiagonal band matrix algorithm is
encountered in both (2.5) and (2.6).In the current work, the equations are integrated in time with
consideration of the low-storage TVD-RK3 scheme [17]. Application
of the CFD6 technique to (1.1) leads to the following first-order
ordinary differential equation in time: (2.8)duidt=Lui,   where L stands for the differential operator. The
spatial nonlinear terms are approximated by the CFD6 scheme as follows: (2.9)Lui=(∂∂x(umux)+bu)i=(mum−1ux2+umuxx+bu)i. The choice of
time integration is influenced by several factors such as desired accuracy,
available memory, computer speed, and stability. A class of high-order TVD time
discretization technique was developed by Gottlieb and Shu [17] for solving
hyperbolic conservation laws with stable spatial discretizations. As was
pointed out by them, unlike in the case of linear operators, there is no stability
criterion for fully discrete methods where Lu is nonlinear. However, the TVD methods
guarantee the stability properties expected of the forward Euler method [18]. The
total variation (TV) of the numerical solution (2.10)TV(u)=∑i|ui+1−ui| does not
increase in time, that is, the TVD property holds: (2.11)TV(uk+1)≤TV(uk).More details
can be found in [17, 18]. The scheme integrates equation (2.8) from time t0 (step k) to t0+Δt (step k+1) through the operations [17, 18](2.12)ui(1)=uik+ΔtLuik,ui(2)=34uik+14ui(1)+14ΔtLui(1),uik+1=13uik+23ui(2)+23ΔtLui(2).The CFD6 scheme will be applied to
physical models to illustrate the strength of the present method. Each spatial
derivative on the right-hand side of (2.8) was computed with the use of
the present method, and then the semidiscrete equation (2.8) was solved with the
help of the low-storage explicit TVD-RK3 scheme. Thus, the solution is obtained
without requiring neither linearization nor transformation. For the approximate[a,b]
solution of the problem (1.1) with the taken boundary and initial conditions
using the current scheme, first the interval  is discretized such that a=x1<x2<⋯<xN=b, where N is the number of grid points.


## 2.1. Discretization

Spatial
derivatives are evaluated by the compact finite difference scheme [14]. For any
scalar pointwise value u,
the derivatives of u can be obtained by solving a tridiagonal
or pentadiagonal system. Much work was done in deriving such formulae [14, 19].
More details on the formulae can be found in [14, 19]. A uniform one-dimensional
mesh is considered, consisting of N points: x1,x2,…,xi−1,xi,xi+1,…,xN. The mesh size is denoted by h=xi+1−xi. The first derivatives can be given at interior nodes as follows [14, 19]: (2.1)αui−1′+ui′+αui+1′=bui+2−ui−24h+aui+1−ui−12h which gives
rise to an α-family of fourth-order tridiagonal schemes
with (2.2)a=23(α+2),b=13(4α−1),where α=0 leads to the explicit fourth-order scheme for
the first derivative. A sixth-order tridiagonal scheme is obtained by α=1/3, (2.3)13ui−1′+ui′+13ui+1′=19ui+2−ui−24h+149ui+1−ui−12h.For the nodes
near the boundary, approximation formulae for the derivatives of nonperiodic
problems can be derived with the consideration of one-sided schemes. More
details on the derivations for the first- and second-order derivatives can be
found in [14, 19]. The derived  formulae at boundary points 1, 2, N−1, and N,
respectively, are as follows:(2.4)ui′+5ui+1′=1h(−19760ui−512ui+1+5ui+2−53ui+3+512ui+4−120ui+5),211ui−1′+ui′+211ui+1′=1h(−2033ui−1−35132ui+3433ui+1−733ui+2+233ui+3−1132ui+4),211ui−1′+ui′+211ui+1′=1h(2033ui+1+35132ui−3433ui−1+733ui−2−233ui−3+1132ui−4),5ui−1′+ui′=1h(19760ui+512ui−1−5ui−2+53ui−3−512ui−4+120ui−5). In order to
obtain the formulae presented above, the procedure of Gaitonde and Visbal [19]
was followed. The formulae can be re-expressed as (2.5)BU′=AU, where U=(u1,…,un)T.
The second-order derivative terms are obtained by applying the first-order
operator twice, that is, (2.6)BU′′=AU′, where
(2.7)A=1h[−19760−5125−53512−120−2033−351323433−733233−1132−136−79079136⋱⋱⋱⋱⋱−136−790791361132−233733−3433351322033120−51253−551219760]N×N,B=[15211121113113⋱⋱⋱13113211121151]N×N.   In the calculation
of the numerical derivatives, a tridiagonal band matrix algorithm is
encountered in both (2.5) and (2.6).In the current work, the equations are integrated in time with
consideration of the low-storage TVD-RK3 scheme [17]. Application
of the CFD6 technique to (1.1) leads to the following first-order
ordinary differential equation in time: (2.8)duidt=Lui,   where L stands for the differential operator. The
spatial nonlinear terms are approximated by the CFD6 scheme as follows: (2.9)Lui=(∂∂x(umux)+bu)i=(mum−1ux2+umuxx+bu)i. The choice of
time integration is influenced by several factors such as desired accuracy,
available memory, computer speed, and stability. A class of high-order TVD time
discretization technique was developed by Gottlieb and Shu [17] for solving
hyperbolic conservation laws with stable spatial discretizations. As was
pointed out by them, unlike in the case of linear operators, there is no stability
criterion for fully discrete methods where Lu is nonlinear. However, the TVD methods
guarantee the stability properties expected of the forward Euler method [18]. The
total variation (TV) of the numerical solution (2.10)TV(u)=∑i|ui+1−ui| does not
increase in time, that is, the TVD property holds: (2.11)TV(uk+1)≤TV(uk).More details
can be found in [17, 18]. The scheme integrates equation (2.8) from time t0 (step k) to t0+Δt (step k+1) through the operations [17, 18](2.12)ui(1)=uik+ΔtLuik,ui(2)=34uik+14ui(1)+14ΔtLui(1),uik+1=13uik+23ui(2)+23ΔtLui(2).The CFD6 scheme will be applied to
physical models to illustrate the strength of the present method. Each spatial
derivative on the right-hand side of (2.8) was computed with the use of
the present method, and then the semidiscrete equation (2.8) was solved with the
help of the low-storage explicit TVD-RK3 scheme. Thus, the solution is obtained
without requiring neither linearization nor transformation. For the approximate[a,b]
solution of the problem (1.1) with the taken boundary and initial conditions
using the current scheme, first the interval  is discretized such that a=x1<x2<⋯<xN=b, where N is the number of grid points.


## 3. Numerical Illustrations

In order to see numerically whether the present
methodology leads to accurate solutions, the CFD6 solutions are evaluated for
some examples of the porous media equations given above.Now, numerical solutions of the
porous media equation are obtained to
validate the current numerical scheme. To verify the efficiency, measure its
accuracy and the versatility of the present scheme for our problem in
comparison with the exact solution, the absolute errors are reported which are
defined by(3.1)|unum(xi,tj)−u(xi,tj)| in the point (xi,tj), where unum(xi,tj) is the solution obtained by (2.8) solved by the combination suggested here, and u(xi,tj) is the exact solution.Consider the porous media equation
in the form (1.1) with the initial condition (1.2) and boundary conditions (1.3). The
results are compared with the analytical solutions. The numerical computations
were performed using uniform grids. All computations were carried out using
some codes produced in Visual Basic for Applications. All computations were
performed using a laptop computer, Genuine Intel(R) T2500 2.00 GHz CPU, 997 MHz 1.00 GB RAM. All computations were carried
out using double-length arithmetic. The
differences between the computed solutions and the exact solutions are shown in
Tables 1–5, 7, 8. A comparison between the CFD6 solution and the solutions
produced by a noncompact fourth-order finite difference scheme (FD4) has also
been made in Table 5. In Table 6, computational orders of the method were
calculated. In Table 7, the absolute errors were computed for various values
of x and N for t=0.05 and Δt=1.0E-06 in Example 3.2. In Table 8, computation of the
absolute errors was carried out for various values of x, Δt and N at t=0.05 in Example 3.2.Table 1
The absolute
errors for various values ofx and t for N=11, h=0.1, Δt=0.0001 in Example 3.1.
xtAbsolute error0.100.011.07E-150.101.08E-141.006.51E-140.500.011.11E-150.101.10E-141.001.32E-130.900.011.33E-150.103.11E-151.003.51E-14Table 2
The absolute
errors for various values ofx and t for N=11, h=0.1, Δt=0.0001 in Example 3.2.
xtAbsolute error0.100.011.09E-140.101.20E-141.003.70E-140.500.011.28E-150.102.32E-151.004.80E-150.900.011.35E-140.101.74E-141.004.86E-14Table 3
The absolute
errors for various values ofx and t for N=11, h=0.1, Δt=0.0001 in Example 3.3.
xtAbsolute error0.100.017.05E-130.101.01E-121.002.91E-110.500.011.07E-130.102.09E-131.003.72E-120.900.011.26E-120.101.73E-121.003.68E-11Table 4
The absolute
errors for various values ofx and t for N=11, h=0.1, Δt=0.0001 in Example 3.4.
xtAbsolute error0.100.018.10E-070.101.22E-071.001.92E-100.500.016.49E-080.104.75E-071.001.54E-110.900.011.01E-060.101.00E-061.001.85E-10Table 5
Comparison
of the CFD6 with the FD4 for various values ofx and t with Δt=0.0001, N=11,
computational domain [1,2] in Example 3.5.
xtExactNumericalAbsolute errorFD4CFD6FD4CFD61.10.101.1103661.1103651.1103671.38E-061.31E-061.006.7173246.7172896.7173313.45E-057.52E-065.0020024.05981320023.95664120024.0822341.03E-012.24E-021.50.100.8142690.8142350.8142683.33E-057.68E-071.004.9260374.9257914.9260302.47E-046.95E-065.0014684.31053014683.57475214684.2896757.36E-012.09E-021.90.100.6428440.6428300.6428421.34E-051.70E-061.003.8889773.8888813.8889649.61E-051.25E-055.0011592.87673411592.58937511592.8393342.87E-013.74E-02Table 6
Convergence rate (CR) of the present scheme:
Comparisons on various values ofh,t and Δt for  3.5.
t=0.0001,Δt=1.0E-06t=0.0003,Δt=1.0E-06t=0.0003,Δt=1.0E-05hAverage absolute
errorCRhAverage
absolute errorCRhAverage absolute
errorCR1/105.84E-091/186.82E-101/188.95E-101/156.86E-105.281/203.84E-105.451/205.06E-105.421/201.48E-105.341/222.27E-105.541/222.99E-105.511/254.36E-115.461/241.39E-105.621/241.84E-105.601/301.58E-115.571/268.80E-115.711/261.17E-105.681/356.57E-125.691/285.72E-115.801/287.60E-115.771/403.03E-125.801/303.81E-115.891/305.07E-115.861/451.51E-125.911/322.59E-115.981/323.45E-115.951/508.02E-136.021/341.79E-116.071/342.39E-116.051/554.47E-136.141/361.26E-116.161/361.68E-116.15Table 7
The absolute errors for
various values ofx and N for t=0.05 and Δt=1.0E-06 in Example 3.2.
xNAbsolute error0.1114.77E-14213.04E-14512.97E-14813.00E-141013.00E-140.5118.17E-14217.66E-14517.64E-14817.68E-141017.68E-140.9111.15E-13212.28E-14512.25E-14812.33E-141012.31E-14Table 8
The absolute errors for
various values ofx, Δt and N for t=0.05 in Example 3.2.
xΔtN=11N=21Absolute errorCPU (s)Absolute errorCPU (s)0.11.0E-041.11E-140.01.74E-040.01.0E-057.09E-142.03.04E-153.01.0E-064.77E-1416.03.04E-1430.01.0E-072.49E-13162.03.03E-13295.01.0E-083.25E-121640.03.27E-122964.00.51.0E-042.26E-150.08.19E-040.01.0E-051.19E-142.07.61E-153.01.0E-068.17E-1416.07.66E-1430.01.0E-078.26E-13162.07.84E-13295.01.0E-088.74E-121640.08.62E-122964.00.91.0E-041.65E-140.007.02E-040.01.0E-058.23E-142.002.50E-153.01.0E-061.15E-1316.002.28E-1430.01.0E-073.66E-13162.002.12E-13295.01.0E-083.07E-121640.002.76E-122964.0As various
problems of science were modeled by nonlinear partial differential equations,
and since therefore the porous media equation is of high importance, the following examples [1, 4, 5, 10, 13, 20] have been considered.
Computational domain [a,b] is taken to be [0,1] in Examples 3.1–3.4 and [1,2] in Example 3.5. In Examples 3.1–3.4 b=0,
while b=2 in Example 3.5.Example 3.1 (see [13]).
Let us takem=1 and b=0 in (1.1). Thus the equation becomes (3.2)ut=∂∂x(uux) with the initial condition (3.3)u(x,0)=x subject to boundary conditions (3.4)u(0,t)=t,u(1,t)=1+t.An exact solution
of  (1.1) is (3.5)u(x,t)=x+t.
In Table1, the absolute errors
have been shown for various values of x and t.
A comparison
has been made between the results of the present scheme and the exact results. It
is seen from the results that the current scheme are to be powerful and efficient.Example 3.2 (see [5]).
Form=−1 and b=0,  (1.1) becomes (3.6)ut=∂∂x(u−1ux). In [5] the authors give an exact solution to (3.6) as follows: (3.7)u(x,t)=1c1x−c12t+c2, where c1 and c2 are arbitrary constants. For the sake of simplicity, c1=1 and c2=10.
With these choices, the solution becomes (3.8)u(x,t)=1x−t+10.Now
(3.6) is solved with the consideration of the initial condition (3.9)u(x,0)=1x+10 and the boundary conditions (3.10)u(0,t)=110−t,u(1,t)=19−t.
The absolute
errors have been shown for various values ofx and t in Table 2. Also, a comparison
between the exact solution and the CFD6 solutions is given in Table 2. The
obtained results are seen to be very accurate.Example 3.3 (see [5]).
Whenm=−4/3 and b=0, (1.1) becomes (3.11)ut=∂∂x(u−4/3ux).An exact solution to  (3.11) is given as follows
[5]: (3.12)u(x,t)=(2c1x−3c12t+c2)−3/4,where c1 and c2 are arbitrary constants, c1=1 and c2=10.
Thus, one has (3.13)u(x,t)=(2x−3t+10)−3/4.Now (3.11) is solved with the initial condition (3.14)u(x,0)=(2x+10)−3/4 subject to the boundary conditions (3.15)u(0,t)=110−3t,u(1,t)=18−3t.
The absolute errors have been shown
for various values ofx and t in Table 3. Also, a comparison
between the exact solution and the CFD6 solution is given in Table 3. The
computed results are seen to be very accurate.Example 3.4 (see [1, 4, 10]).
whenm=−2 and b=0,
(1.1) represents different nonlinear
models and becomes  (3.16)ut=∂∂x(u−2ux). An exact
solution of  (3.16) is(3.17)u(x,t)=(x2+e2t)−1/2.Now  (3.16) is solved with the initial condition (3.18)u(x,0)=(x2+1)−1/2 subject to the boundary conditions(3.19)u(0,t)=e−t,u(1,t)=(1+e2t)−1/2
The computed results for this example have been given
for various values of x and t in Table 4. The results of the combination of the CFD6 method with the low-storage explicit TVD-RK3 have
been presented in the table. Comparisons of the current results with the exact
solution showed that the presented results are very accurate.Example 3.5 (see [20]).
Whenm=−2 and b=2,
(1.1) represents nonlinear realistic
models and becomes (3.20)ut=∂∂x(u−2ux)+2u. An exact
solution of (3.20) is [18] (3.21)u(x,t)=x−1e2t.Now (3.20) is solved with the initial condition (3.22)u(x,0)=x−1 subject to the boundary conditions(3.23)u(1,t)=e2t,u(2,t)=e2t2. In Table 5, computation of the
absolute errors was carried out for various values of x and t in Example 3.5.
To see the effect of the CFD6 solution, Example 3.5 was also solved using
the FD4. Note that the effect of the CFD6 solutions is very clear as seen in Table 5.
Numerical rate of convergence (CR)
has also been studied to know about the convergency of the scheme. To achieve this, the CR was calculated using (3.24)CR=log(E1/E2)log(h1/h2),where E1 and E2 are errors correspond to grids with mesh size h1 and h2,
respectively, (see Table 6). Also computational orders in Table 7 show the
high-order accuracy of the present method for solving such problems.In this paper,
it is shown that the
approximate solutions of the porous media equation
are very close to the exact solutions. The absolute errors have been
calculated form=1, m=−1, and m=−4/3, and m=−2,
in Tables 1, 2, 3, and 4, respectively, in the point (xi,tj).
For all four cases b=0,
but m=−2 is used twice, the second one with b=2 in Example 3.5. As seen from Tables 1–5 and 7-8,
the absolute errors in all cases are very
small. Very good approximations to
the exact solutions are achieved.The method was also applied to realistic problem sizes (see Table7), and
seen that the solution becomes stable for small number of grid points.
Increasing in the number of grid points requires small time step size for
satisfying numerical stability of the solutions. This situation makes necessary
the use of excessive CPU time and gives rise to accumulation of numerical
errors (see Tables 7 and 8). Compact high-order schemes
are closer to spectral methods [14] and at the same time maintain the freedom
to retain accuracy even though relatively small number of grid points is used. Note that the same error is
obtained for all discretizations; the error is close to the level of machine
accuracy, and this is believed to be a saturation effect (see Table 7).


### 4. Conclusions

In this paper, numerical
simulations of the porous media equation were dealt with using a combination of
the CFD6 scheme in space and a low-storage explicit TVD-RK3 scheme in time. The
method successfully worked to give very reliable and accurate solutions to the
equation. The method gives convergent
approximations and handles nonlinear problems. In this method, there is no need for linearization of
nonlinear terms. Nonlinear scientific models arise frequently in scientific problems for expressing nonlinear phenomena. For
nonlinear problems, the present method is seen to be a very good choice to
achieve a high degree of accuracy while dealing with the problems. The computed results justify the advantage of this method. The present method needs less use of storage
space.


## 4. Conclusions

In this paper, numerical
simulations of the porous media equation were dealt with using a combination of
the CFD6 scheme in space and a low-storage explicit TVD-RK3 scheme in time. The
method successfully worked to give very reliable and accurate solutions to the
equation. The method gives convergent
approximations and handles nonlinear problems. In this method, there is no need for linearization of
nonlinear terms. Nonlinear scientific models arise frequently in scientific problems for expressing nonlinear phenomena. For
nonlinear problems, the present method is seen to be a very good choice to
achieve a high degree of accuracy while dealing with the problems. The computed results justify the advantage of this method. The present method needs less use of storage
space.



---

*Source: 912541-2009-04-07.xml*