Mathematical Problems in Engineering

Mathematical Problems in Engineering / 2009 / Article

Research Article | Open Access

Volume 2009 |Article ID 912541 | 13 pages | https://doi.org/10.1155/2009/912541

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

Academic Editor: Alois Steindl
Received30 Apr 2008
Revised10 Oct 2008
Accepted16 Jan 2009
Published07 Apr 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.

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 with the initial condition and the boundary conditions Here π‘š is a rational number while 𝑏 is real constant. In (1.1)–(1.3), π‘₯ and 𝑑 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 𝐷(𝑒)=π‘’π‘š does not satisfy the condition for classical diffusion equations, 𝐷(𝑒)>0 [6].

For the motion of thin viscous films, (1.1) with π‘š=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 π‘š 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 𝑒, the derivatives of 𝑒 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 𝑁 points: π‘₯1,π‘₯2,…,π‘₯π‘–βˆ’1,π‘₯𝑖,π‘₯𝑖+1,…,π‘₯𝑁. The mesh size is denoted by β„Ž=π‘₯𝑖+1βˆ’π‘₯𝑖. The first derivatives can be given at interior nodes as follows [14, 19]: which gives rise to an 𝛼-family of fourth-order tridiagonal schemes with where 𝛼=0 leads to the explicit fourth-order scheme for the first derivative. A sixth-order tridiagonal scheme is obtained by 𝛼=1/3, 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, π‘βˆ’1, and 𝑁, respectively, are as follows: In order to obtain the formulae presented above, the procedure of Gaitonde and Visbal [19] was followed. The formulae can be re-expressed as where π‘ˆ=(𝑒1,…,𝑒𝑛)𝑇. The second-order derivative terms are obtained by applying the first-order operator twice, that is, where 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: where 𝐿 stands for the differential operator. The spatial nonlinear terms are approximated by the CFD6 scheme as follows: 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 𝐿𝑒 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 does not increase in time, that is, the TVD property holds: More details can be found in [17, 18]. The scheme integrates equation (2.8) from time 𝑑0 (step π‘˜) to 𝑑0+Δ𝑑 (step π‘˜+1) through the operations [17, 18]

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[π‘Ž,𝑏] solution of the problem (1.1) with the taken boundary and initial conditions using the current scheme, first the interval is discretized such that π‘Ž=π‘₯1<π‘₯2<β‹―<π‘₯𝑁=𝑏, where 𝑁 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 in the point (π‘₯𝑖,𝑑𝑗), where 𝑒num(π‘₯𝑖,𝑑𝑗) is the solution obtained by (2.8) solved by the combination suggested here, and 𝑒(π‘₯𝑖,𝑑𝑗) 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 π‘₯ and 𝑁 for 𝑑=0.05 and Δ𝑑=1.0E-06 in Example 3.2. In Table 8, computation of the absolute errors was carried out for various values of π‘₯, Δ𝑑 and 𝑁 at 𝑑=0.05 in Example 3.2.


π‘₯ 𝑑 Absolute error

0.100.011.07E-15
0.101.08E-14
1.006.51E-14

0.500.011.11E-15
0.101.10E-14
1.001.32E-13

0.900.011.33E-15
0.103.11E-15
1.003.51E-14


π‘₯ 𝑑 Absolute error

0.100.011.09E-14
0.101.20E-14
1.003.70E-14

0.500.011.28E-15
0.102.32E-15
1.004.80E-15

0.900.011.35E-14
0.101.74E-14
1.004.86E-14


π‘₯ 𝑑 Absolute error

0.100.017.05E-13
0.101.01E-12
1.002.91E-11

0.500.011.07E-13
0.102.09E-13
1.003.72E-12

0.900.011.26E-12
0.101.73E-12
1.003.68E-11


π‘₯ 𝑑 Absolute error

0.100.018.10E-07
0.101.22E-07
1.001.92E-10

0.500.016.49E-08
0.104.75E-07
1.001.54E-11

0.900.011.01E-06
0.101.00E-06
1.001.85E-10


π‘₯ 𝑑 ExactNumericalAbsolute error
FD4CFD6FD4CFD6

1.10.101.1103661.1103651.1103671.38E-061.31E-06
1.006.7173246.7172896.7173313.45E-057.52E-06
5.0020024.05981320023.95664120024.0822341.03E-012.24E-02

1.50.100.8142690.8142350.8142683.33E-057.68E-07
1.004.9260374.9257914.9260302.47E-046.95E-06
5.0014684.31053014683.57475214684.2896757.36E-012.09E-02

1.90.100.6428440.6428300.6428421.34E-051.70E-06
1.003.8889773.8888813.8889649.61E-051.25E-05
5.0011592.87673411592.58937511592.8393342.87E-013.74E-02


𝑑 = 0 . 0 0 0 1 , Ξ” 𝑑 = 1 . 0 E - 0 6 𝑑 = 0 . 0 0 0 3 , Ξ” 𝑑 = 1 . 0 E - 0 6 𝑑 = 0 . 0 0 0 3 , Ξ” 𝑑 = 1 . 0 E - 0 5

β„Ž Average absolute errorCR β„Ž Average absolute errorCR β„Ž Average absolute errorCR

1 / 1 0 5.84E-09 1 / 1 8 6.82E-10 1 / 1 8 8.95E-10
1 / 1 5 6.86E-105.28 1 / 2 0 3.84E-105.45 1 / 2 0 5.06E-105.42
1 / 2 0 1.48E-105.34 1 / 2 2 2.27E-105.54 1 / 2 2 2.99E-105.51
1 / 2 5 4.36E-115.46 1 / 2 4 1.39E-105.62 1 / 2 4 1.84E-105.60
1 / 3 0 1.58E-115.57 1 / 2 6 8.80E-115.71 1 / 2 6 1.17E-105.68
1 / 3 5 6.57E-125.69 1 / 2 8 5.72E-115.80 1 / 2 8 7.60E-115.77
1 / 4 0 3.03E-125.80 1 / 3 0 3.81E-115.89 1 / 3 0 5.07E-115.86
1 / 4 5 1.51E-125.91 1 / 3 2 2.59E-115.98 1 / 3 2 3.45E-115.95
1 / 5 0 8.02E-136.02 1 / 3 4 1.79E-116.07 1 / 3 4 2.39E-116.05
1 / 5 5 4.47E-136.14 1 / 3 6 1.26E-116.16 1 / 3 6 1.68E-116.15


π‘₯ 𝑁 Absolute error

0.1114.77E-14
213.04E-14
512.97E-14
813.00E-14
1013.00E-14

0.5118.17E-14
217.66E-14
517.64E-14
817.68E-14
1017.68E-14

0.9111.15E-13
212.28E-14
512.25E-14
812.33E-14
1012.31E-14


π‘₯ Ξ” 𝑑 𝑁 = 1 1 𝑁 = 2 1
Absolute errorCPU (s)Absolute errorCPU (s)

0.11.0E-041.11E-140.01.74E-040.0
1.0E-057.09E-142.03.04E-153.0
1.0E-064.77E-1416.03.04E-1430.0
1.0E-072.49E-13162.03.03E-13295.0
1.0E-083.25E-121640.03.27E-122964.0

0.51.0E-042.26E-150.08.19E-040.0
1.0E-051.19E-142.07.61E-153.0
1.0E-068.17E-1416.07.66E-1430.0
1.0E-078.26E-13162.07.84E-13295.0
1.0E-088.74E-121640.08.62E-122964.0

0.91.0E-041.65E-140.007.02E-040.0
1.0E-058.23E-142.002.50E-153.0
1.0E-061.15E-1316.002.28E-1430.0
1.0E-073.66E-13162.002.12E-13295.0
1.0E-083.07E-121640.002.76E-122964.0

As 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 [π‘Ž,𝑏] 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 𝑏=0, while 𝑏=2 in Example 3.5.

Example 3.1 (see [13]). Let us take π‘š=1 and 𝑏=0 in (1.1). Thus the equation becomes with the initial condition subject to boundary conditions An exact solution of (1.1) is
In Table 1, the absolute errors have been shown for various values of π‘₯ and 𝑑. 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]). For π‘š=βˆ’1 and 𝑏=0, (1.1) becomes In [5] the authors give an exact solution to (3.6) as follows: where 𝑐1 and 𝑐2 are arbitrary constants. For the sake of simplicity, 𝑐1=1 and 𝑐2=10. With these choices, the solution becomes Now (3.6) is solved with the consideration of the initial condition and the boundary conditions
The absolute errors have been shown for various values of π‘₯ and 𝑑 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]). When π‘š=βˆ’4/3 and 𝑏=0, (1.1) becomes An exact solution to (3.11) is given as follows [5]: where 𝑐1 and 𝑐2 are arbitrary constants, 𝑐1=1 and 𝑐2=10. Thus, one has Now (3.11) is solved with the initial condition subject to the boundary conditions
The absolute errors have been shown for various values of π‘₯ and 𝑑 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]). when π‘š=βˆ’2 and 𝑏=0, (1.1) represents different nonlinear models and becomes An exact solution of (3.16) isNow (3.16) is solved with the initial condition subject to the boundary conditions The computed results for this example have been given for various values of π‘₯ and 𝑑 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]). When π‘š=βˆ’2 and 𝑏=2, (1.1) represents nonlinear realistic models and becomes An exact solution of (3.20) is [18] Now (3.20) is solved with the initial condition subject to the boundary conditions In Table 5, computation of the absolute errors was carried out for various values of π‘₯ and 𝑑 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 where 𝐸1 and 𝐸2 are errors correspond to grids with mesh size β„Ž1 and β„Ž2, 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 for π‘š=1, π‘š=βˆ’1, and π‘š=βˆ’4/3, and π‘š=βˆ’2, in Tables 1, 2, 3, and 4, respectively, in the point (π‘₯𝑖,𝑑𝑗). For all four cases 𝑏=0, but π‘š=βˆ’2 is used twice, the second one with 𝑏=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 Table 7), 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.

Acknowledgments

The author would like to thank anonymous referees of the journal of MPE for their valuable comments and suggestions to improve most of this paper. The author is very grateful to G. GΓΌrarslan (Department of Civil Engineering, Faculty of Engineering, Pamukkale University, Denizli, Turkey) for his comments and careful reading the paper.

References

  1. E. A. Saied, β€œThe non-classical solution of the inhomogeneous non-linear diffusion equation,” Applied Mathematics and Computation, vol. 98, no. 2-3, pp. 103–108, 1999. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  2. E. A. Saied and M. M. Hussein, β€œNew classes of similarity solutions of the inhomogeneous nonlinear diffusion equations,” Journal of Physics A, vol. 27, no. 14, pp. 4867–4874, 1994. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  3. L. Dresner, Similarity Solutions of Nonlinear Partial Differential Equations, vol. 88 of Research Notes in Mathematics, Pitman, Boston, Mass, USA, 1983. View at: Zentralblatt MATH | MathSciNet
  4. Q. Changzheng, β€œExact solutions to nonlinear diffusion equations obtained by a generalized conditional symmetry method,” IMA Journal of Applied Mathematics, vol. 62, no. 3, pp. 283–302, 1999. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  5. A. D. Polyanin and V. F. Zaitsev, Handbook of Nonlinear Partial Differential Equations, Chapman & Hall/CRC Press, Boca Raton, Fla, USA, 2004. View at: Zentralblatt MATH | MathSciNet
  6. T. P. Witelski, β€œSegregation and mixing in degenerate diffusion in population dynamics,” Journal of Mathematical Biology, vol. 35, no. 6, pp. 695–712, 1997. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  7. J. D. Murray, Mathematical Biology, vol. 19 of Biomathematics, Springer, Berlin, Germany, 2nd edition, 1993. View at: Zentralblatt MATH | MathSciNet
  8. J. R. King, β€œExact results for the nonlinear diffusion equations u/t=(/x)(u4/3u/x) and u/t=(/x)(u2/3u/x),” Journal of Physics A, vol. 24, no. 24, pp. 5721–5745, 1991. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  9. W. F. Ames, Nonlinear Partial Differential Equations in Engineering. Vol. II, Mathematics in Science and Engineering, Vol. 18, Academic Press, New York, NY, USA, 1972. View at: Zentralblatt MATH | MathSciNet
  10. A.-M. Wazwaz, β€œExact solutions to nonlinear diffusion equations obtained by the decomposition method,” Applied Mathematics and Computation, vol. 123, no. 1, pp. 109–122, 2001. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  11. A. Sadighi and D. D. Ganji, β€œExact solutions of nonlinear diffusion equations by variational iteration method,” Computers & Mathematics with Applications, vol. 54, no. 7-8, pp. 1112–1121, 2007. View at: Google Scholar | Zentralblatt MATH | MathSciNet
  12. A.-M. Wazwaz, β€œThe variational iteration method: a powerful scheme for handling linear and nonlinear diffusion equations,” Computers & Mathematics with Applications, vol. 54, no. 7-8, pp. 933–939, 2007. View at: Google Scholar | Zentralblatt MATH | MathSciNet
  13. S. Pamuk, β€œSolution of the porous media equation by Adomian's decomposition method,” Physics Letters A, vol. 344, no. 2-4, pp. 184–188, 2005. View at: Publisher Site | Google Scholar | MathSciNet
  14. S. K. Lele, β€œCompact finite difference schemes with spectral-like resolution,” Journal of Computational Physics, vol. 103, no. 1, pp. 16–42, 1992. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  15. M. Sari and G. Gürarslan, β€œA sixth-order compact finite difference scheme to the numerical solutions of Burgers' equation,” Applied Mathematics and Computation, vol. 208, no. 2, pp. 475–483, 2009. View at: Publisher Site | Google Scholar
  16. M. Sari, G. Gürarslan, and İ. Dağ, β€œA compact finite difference method for the solution of the generalized Burgers-Fisher equation,” Numerical Methods for Partial Differential Equations. In press. View at: Publisher Site | Google Scholar
  17. S. Gottlieb and C.-W. Shu, β€œTotal variation diminishing Runge-Kutta schemes,” Mathematics of Computation, vol. 67, no. 221, pp. 73–85, 1998. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  18. S. Gottlieb, C.-W. Shu, and E. Tadmor, β€œStrong stability-preserving high-order time discretization methods,” SIAM Review, vol. 43, no. 1, pp. 89–112, 2001. View at: Publisher Site | Google Scholar | Zentralblatt MATH | MathSciNet
  19. D. V. Gaitonde and M. R. Visbal, β€œHigh-order schemes for Navier–Stokes equations: algorithms and Implementation into FDL3DI,” Air Force Research Laboratory, Wright-Patterson AFB, Dayton, Ohio, USA, 1998. View at: Google Scholar
  20. N. Pamuk, β€œSeries solution for porous medium equation with a source term by Adomian's decomposition method,” Applied Mathematics and Computation, vol. 178, no. 2, pp. 480–485, 2006. View at: Google Scholar | Zentralblatt MATH | MathSciNet

Copyright © 2009 Murat Sari. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


More related articles

1529Β Views | 885Β Downloads | 11Β Citations
 PDF  Download Citation  Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

We are committed to sharing findings related to COVID-19 as quickly and safely as possible. Any author submitting a COVID-19 paper should notify us at help@hindawi.com to ensure their research is fast-tracked and made available on a preprint server as soon as possible. We will be providing unlimited waivers of publication charges for accepted articles related to COVID-19. Sign up here as a reviewer to help fast-track new submissions.