Jie Yao^1*, Cameron L. Williams², Fazle Hussain¹, Donald J. Kouri^2,3*
1Department of Mechanical Engineering, Texas Tech University, Lubbock, TX, USA.
²Department of Mathematics, University of Houston, Houston, TX, USA.
³Departments of Mechanical Engineering and Physics, University of Houston, Houston, TX, USA.
DOI: 10.4236/am.2019.1012072   PDF    HTML   XML   660 Downloads   2,101 Views   Citations

Abstract

The solution of a nonlinear diffusion equation is numerically investigated using the generalized Fourier transform method. This equation includes fractal dimensions and power-law dependence on the radial variable and on the diffusion function. The generalized Fourier transform approach is the extension of the Fourier transform method used for the normal diffusion equation. The feasibility of the approach is validated by comparing the numerical result with the exact solution for a point-source. The merit of the numerical method is that it provides a way to calculate anomalous diffusion with an arbitrary initial condition.

Keywords

Generalized Fourier Transform, Anomalous Diffusion, Nonlinear

Share and Cite:

1. Introduction

In the last few decades, anomalous diffusion has been extensively studied in a variety of physical applications, such as turbulent diffusion [1], surface growth [2], transport of fluid in porous media [2], hydraulics problems [3], etc. The diffusion is usually characterized by the time dependence of mean-square displacement (MSD) viz., $\langle r^2\rangle \propto t^{\sigma }$. The MSD grows linearly with time ( $\sigma =1$ ) for the normal diffusion case, and nonlinearly with time for the anomalous diffusion. The process is called sub-diffusion for $0<\sigma <1$ and super-diffusion for $\sigma >1$. The standard normal diffusion described by the Gaussian distribution can be obtained from the usual Fokker-Planck equation with a constant diffusion coefficient and zero drift [4]. Extensions of the conventional Fokker-Planck equation have been used to study anomalous diffusion. For example, anomalous diffusion can be obtained by the usual Fokker-Planck equation, but with a variable diffusion coefficient [5] [6]. It can also be achieved by incorporating nonlinear terms in the diffusion term, or external forces [7] [8] [9] [10]. In some approaches, fractional equations have been employed to analyze anomalous diffusion and related phenomena [11] [12] [13] [14] [15].

In this paper, we study the generalized nonlinear diffusion equation including a fractal dimension d and a diffusion coefficient which depends on the radial variable and the diffusion function $\rho$ [16] [17] [18]

$\frac{\partial \rho }{\partial t}=K_0\frac{1}{r^{d-1}}\frac{\partial }{\partial r}\left(r^{d-1-\theta }\frac{\partial }{\partial r}{\rho }^{\nu }\right)\mathrm{,}$ (1)

with the initial and boundary conditions

$\rho \left(r\mathrm{,}{t}_0\right)={\rho }_0\left(r\right)\mathrm{,}$ (2)

$\rho \left(\infty \mathrm{,}t\right)=\mathrm{0,}$ (3)

where r is the radial coordinate, and $\theta$ and $\nu$ are real parameters. When the diffusion coefficient is a function of r only, it is a generalization of the diffusion equation for fractal geometry [19]. It is the traditional nonlinear diffusion equation when the diffusion coefficient depends on $\rho$ only [20] [21]. Analytical solutions of Equation (1) with a point source have been reported in [16] [17], where an ansatz for $\rho$ is proposed as a general stretched Gaussian function. In [18], the same analytic solutions were also obtained by using Lie group symmetry analysis.

Motivated by the research on generalized nonlinear diffusion, we propose here a numerical method for solving Equation (1) using a generalized Fourier transform. The generalized Fourier transform (also called the ${\Phi }_n$ transform) is a new family of integral transforms developed by Willams et al. [22] [23]. These transforms share all the properties of the Fourier transform; hence can be employed to perform more general frequency and time-frequency analysis [24] [25].

In Section 2, a brief introduction of the generalized Fourier transform is provided. The procedure of using the generalized Fourier transform for solving the generalized nonlinear diffusion equation is discussed in 3.1 and 3.2. The method is validated by comparison between analytical and numerical results. Then, some numerical results for a non-Delta function initial condition are given in 3.3. Conclusions are drawn in 4.

2. Generalized Fourier Transform

The generalized Fourier transform ${\Phi }_n$ is defined as

${\Phi }_{n}f\left(k\right)={\displaystyle \int {\phi }_n}\left(kx\right)f\left(x\right)\text{d}x\mathrm{,}$ (4)

where the integral kernel ${\phi }_n\left(\omega x\right)=c_n\left(kx\right)+\text{i}{s}_n\left(kx\right)$ is,

$c_n\left(\eta \right)=\frac{1}{2}{\left|\eta \right|}^{n-1/2}{J}_{-1+\frac{1}{2n}}\left(\frac{{\left|\eta \right|}^n}{n}\right)\mathrm{,}$ (5)

and

$s_n\left(\eta \right)=\frac{1}{2}\text{sgn}\left(\eta \right){\left|\eta \right|}^{n-1/2}{J}_{1-\frac{1}{2n}}\left(\frac{{\left|\eta \right|}^n}{n}\right)\mathrm{,}$ (6)

where $J_{\nu }\left(\eta \right)$ is the cylindrical Bessel function, and n is the transform order, i.e., $n=1,2,\cdots$.

The Fourier transform $\mathcal{F}$ is the special case with $n=1$. The ${\Phi }_n$ transform shares many properties of the Fourier transform. Here we focus on two properties which will be used later. It is well-known that the Fourier transform preserves the functional form of a Gaussian; particularly, $\mathcal{F}\left[g_1\right]=g_1$ if $g_1\left(x\right)=\exp \left(-x^2/2\right)$. For the generalized Fourier transform, we have ${\Phi }_n{g}_n=g_n$, if $g_n\left(x\right)={\text{e}}^{-\frac{x^{2n}}{2n}}$.

In addition, the generalized Fourier transform also has the following derivative property:

${\Phi }_n\left[\frac{\partial }{\partial x}{x}^{2-2n}\frac{\partial }{\partial x}f\right]=k^{-2n}{\Phi }_{n}f\mathrm{.}$ (7)

In [22], the ${\Phi }_n$ transform is developed for integer order case. However, it can be easily extended to the non-integer case; see [22] for additional discussion of the properties of the ${\Phi }_n$ transform.

3. Solving the Generalized Diffusion Equation with Generalized Fourier Transform

It is well known that the Fourier transform can be used to find the solution for the standard diffusion equation [26]. Motivated by this idea, here we explore employing the generalized Fourier transform for solving the generalized nonlinear diffusion equation.

3.1. The O’Shaugnessy-Procaccia Anomalous Diffusion Equation on Fractals

Let us first consider the generalization of the diffusion equation for fractal geometry, where the diffusion coefficient is a function of r only (i.e. $\nu =1$ ) [19]. Equation (1) can be reduced to

$\frac{\partial \rho \left(r\mathrm{,}t\right)}{\partial t}=\frac{K_0}{r^{d-1}}\left(\frac{\partial }{\partial r}{r}^{d-1-\theta }\frac{\partial }{\partial r}\rho \left(r\mathrm{,}t\right)\right)\mathrm{.}$ (8)

In order to perform the ${\Phi }_n$ transform, we apply the following scaling relationship

$\frac{\partial }{\partial r}(\cdot )=d{r}^{d-1}\frac{\partial }{\partial r^d}(\cdot ),$ (9)

to Equation (8); and with some simplification, we obtain

$\frac{\partial \rho \left(\tilde{r}\mathrm{,}t\right)}{\partial t}={\tilde{K}}_0\frac{\partial }{\partial \tilde{r}}{\tilde{r}}^{2-\lambda /d}\frac{\partial }{\partial \tilde{r}}\rho \left(\tilde{r}\mathrm{,}t\right)\mathrm{,}$ (10)

where ${\tilde{K}}_0=K_0{d}^2$, $\tilde{r}=r^d$, and $\lambda =2+\theta$.

By applying the ${\Phi }_n$ transform to both sides and employing the derivative identity (Equation (7)), we obtain the diffusion equation in the wavenumber domain

$\frac{\partial \tilde{\rho }}{\partial t}=-{\tilde{K}}_0{k}^{\lambda /d}\tilde{\rho },$ (11)

with $\tilde{\rho }={\Phi }_n\rho$.

Equation (11) can be exactly solved as

$\tilde{\rho }\left(k\mathrm{,}t\right)={\text{e}}^{-{\tilde{K}}_0{k}^{\lambda /d}t}{\tilde{\rho }}_0\mathrm{.}$ (12)

The solution to Equation (8) is then obtained by applying the inverse ${\Phi }_n$ transform to $\tilde{\rho }\left(k\mathrm{,}t\right)$.

We validate the ${\Phi }_n$ transform method by comparing the numerical results with the analytical solution or a point source at the origin (i.e. $\rho \left(r\mathrm{,}{t}_0\right)=\delta \left(r\right)$ ), which is given as [19]

${\rho }_a\left(r\mathrm{,}t\right)=\frac{\lambda }{d\Gamma \left(d/\lambda \right)}{\left(\frac{1}{K_0{\lambda }^{2}t}\right)}^{d/\lambda }\exp \left(-\frac{r^{\lambda }}{K_0{\lambda }^{2}t}\right)\mathrm{.}$ (13)

Figure 1 shows the analytical and the numerical solution for $K_0=1$, $D=1$, and $\theta =2.5$ at different times. According to the classification discussed in [16], this example is a subdiffusion case with $\theta >D\left(1-\nu \right)$. From Figure 1, it can be seen that the numerical solution is in good agreement with the analytical solution. In addition, we observe the short tail behaviours of the solution $\rho \left(r\mathrm{,}t\right)$.

3.2. Generalized Nonlinear Equation

Now we consider the generalized nonlinear diffusion equation with $\nu \ne 1$. For the point source (or Dirac delta initial condition), Equation (1) was analytically solved using a generalized stretched Gaussian function approach in [16] :

$\rho \left(r\mathrm{,}t\right)=(\begin{array}{l}{\left[1-\left(1-q\right)\beta \left(t\right)r^{\lambda }\right]}^{1/\left(1-q\right)}/Z\left(t\right)\mathrm{,}\text{if}1-\left(1-q\right)\beta \left(t\right)r^{\lambda }\ge \mathrm{0\; <mo>\; ;\; </mo>}\\ \mathrm{0,}\text{otherwise}\mathrm{.}\end{array}$

Here $q=2-\nu$, and $\beta \left(t\right)$ and $Z\left(t\right)$ are functions given in Equation (12) in [16]. The same solution is derived in [17] using Lie group symmetry method.

In order to solve the generalized nonlinear diffusion equation numerically, we follow the procedure in 3.1, transforming the spatial domain equation to the wavenumber domain using the ${\Phi }_n$ transform. Instead of Equation (11), the wavenumber domain diffusion equation becomes

$\frac{\partial \tilde{\rho }}{\partial t}=-{\tilde{K}}_0{k}^{\lambda /d}\widetilde{{\rho }^{\nu }},$ (14)

with $\widetilde{{\rho }^{\nu }}={\Phi }_n\left({\rho }^{\nu }\right)$.

Due to the presence of nonlinearity term in the right hand side of Equation (13) (i.e. ${\rho }^{\nu }$ ), an analytical solution in the form of Equation (12) is difficult to be obtained. However, Equation (13) can be numerically solved by employing certain types of time-stepping discretization methods for the time derivative. Here, the simple forward Euler finite difference scheme is employed for time discretization with $\Delta t=0.01\text{s}$. An equally spaced mesh with $N_r=1001$ is used over the domain $\tilde{r}=\left[\mathrm{0,30}\right]$. The comparisons between the exact [16] [17] and numerical solution for the point source (Dirac delta function $\delta \left(r\right)$ ) initial condition in scaled ( $\tilde{r}$ ) and original (r) coordinates are shown in Figure 2(a) and Figure 2(b), respectively. The parameters used here are $K_0=1$, $D=3$, $\theta =2.5$, and $\nu =0.8$. Note that to avoid performing a ${\Phi }_n$ transform for the fractional order of Dirac delta function (as the definition of ${\delta }^{\nu }$ is also an ongoing research topic [27] [28], we use ${\rho }_a\left(r,t_0=0.1\right)$ as the initial condition for our numerical simulation. Again, good agreement between the numerical and analytical solutions can be observed.

3.3. Generalized Diffusion for Arbitrary Initial Condition

The merit of the numerical approach using the generalized Fourier transform is that it provides a way for solving the generalized diffusion equation with arbitrary initial condition. In Figure 3, we present the numerical solution of the generalized diffusion equation for $K_0=1$, $D=1$ and $\theta =2.5$ with the Gaussian initial condition

${\rho }_0\left(r\mathrm{,}{t}_0\right)=\frac{1}{\sqrt{4\pi t_0}}\exp \left(-\frac{r^2}{4K_0{t}_0}\right)\mathrm{,}$ (14)

where $t_0=0.1$.

As we can see, the diffusion process finally approaches the same generalized

Gaussian shape as in the point source case (Figure 1). In [29], it was analytically shown that the normal diffusion equation, when initialized with a generalized Gaussian distribution will asymptotically approach its final solution, i.e. a Gaussian distribution. Here, we present a numerical example of what amounts to the “generalized central limit’’ behaviour in which the diffusion process will finally transform the arbitrary initial distribution to the corresponding generalized Gaussian distribution [30] [31]. A rigorous proof of the existence of the attractor of the generalized Gaussian diffusion has been done [24] ] for linear diffusion. This is a consequence of the negative semi-definite spectrum of the usual Laplacian. However, as mentioned in [31], the diffusion procedure, initialized with different distribution, may take very long time to reach its asymptotic behaviour. In addition, by comparing with the solution for normal diffusion with the same initial condition, the sub-diffusion process clearly exhibits the short tail behaviour.

4. Conclusion

In this paper, a numerical method for solving the generalized nonlinear diffusion equation has been presented and validated. The method is based on the generalized Fourier transform and has been validated by comparing the numerical solution with analytical solution for the point source. The presented method may serve as a useful tool to study a variety of systems involving the anomalous diffusion. Currently, no fast transform algorithm has yet been developed for the ${\Phi }_n$ transform. This issue will be investigated in future study.

Acknowledgements

Discussions with Bernhard G. Bodmann are appreciated. Partial support for this work was provided by resources of the uHPC cluster managed by the University of Houston under NSF Award Number 1531814.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Gavrilov, V., Klepikova, N. and Rodean, H. (1995) Trial of a Nonlinear Diffusion Equation as a Model of Turbulent Diffusion. Atmospheric Environment, 29, 2317-2322. https://doi.org/10.1016/1352-2310(95)00148-R
[2]	Spohn, H. (1993) Surface Dynamics below the Roughening Transition. Journal de Physique I, 3, 69-81. https://doi.org/10.1051/jp1:1993117
[3]	Daly, E. and Porporato, A. (2004) Similarity Solutions of Nonlinear Diffusion Problems Related to Mathematical Hydraulics and the Fokker-Planck Equation. Physical Review E, 70, Article ID: 056303. https://doi.org/10.1103/PhysRevE.70.056303
[4]	Risken, H. (1984) Fokker-Planck Equation. In: The Fokker-Planck Equation, Springer, Berlin, 63-95. https://doi.org/10.1007/978-3-642-96807-5_4
[5]	Fa, K.S. (2005) Exact Solution of the Fokker-Planck Equation for a Broad Class of Diffusion Coefficients. Physical Review E, 72, Article ID: 020101.
[6]	Fa, K.S. (2011) Solution of Fokker-Planck Equation for a Broad Class of Drift and Diffusion Coefficients. Physical Review E, 84, Article ID: 012102.
[7]	Bologna, M., Tsallis, C. and Grigolini, P. (2000) Anomalous Diffusion Associated with Nonlinear Fractional Derivative Fokker-Planck-Like Equation: Exact Time-Dependent Solutions. Physical Review E, 62, 2213. https://doi.org/10.1103/PhysRevE.62.2213
[8]	Assis Jr., P., da Silva, P., da Silva, L., Lenzi, E. and Lenzi, M. (2006) Nonlinear Diffusion Equation and Nonlinear External Force: Exact Solution. Journal of Mathematical Physics, 47, Article ID: 103302. https://doi.org/10.1063/1.2354334
[9]	Lenzi, E., Lenzi, M., Gimenez, T. and da Silva, L. (2010) Some Results for an N-Dimensional Nonlinear Diffusion Equation with Radial Symmetry. Journal of Engineering Mathematics, 67, 233-240. https://doi.org/10.1007/s10665-009-9351-6
[10]	Zola, R., Lenzi, M., Evangelista, L., Lenzi, E., Lucena, L. and da Silva, L. (2008) Exact Solutions for a Diffusion Equation with a Nonlinear External Force. Physics Letters A, 372, 2359-2363. https://doi.org/10.1016/j.physleta.2007.12.007
[11]	Metzler, R. and Klafter, J. (2000) The Random Walk’s Guide to Anomalous Diffusion: A Fractional Dynamics Approach. Physics Reports, 339, 1-77. https://doi.org/10.1016/S0370-1573(00)00070-3
[12]	Sokolov, I.M., Klafter, J. and Blumen, A. (2002) Fractional Kinetics. Physics Today, 55, 48. https://doi.org/10.1063/1.1535007
[13]	Tsallis, C. and Lenzi, E. (2002) Anomalous Diffusion: Nonlinear Fractional Fokker-Planck Equation. Chemical Physics, 284, 341-347. https://doi.org/10.1016/S0301-0104(02)00557-8
[14]	Liu, F., Anh, V. and Turner, I. (2004) Numerical Solution of the Space Fractional Fokker-Planck Equation. Journal of Computational and Applied Mathematics, 166, 209-219. https://doi.org/10.1016/j.cam.2003.09.028
[15]	Tateishi, A.A., Ribeiro, H.V. and Lenzi, E.K. (2017) The Role of Fractional Time-Derivative Operators on Anomalous Diffusion. Frontiers in Physics, 5, 52. https://doi.org/10.3389/fphy.2017.00052
[16]	Malacarne, L., Mendes, R., Pedron, I. and Lenzi, E. (2001) Nonlinear Equation for Anomalous Diffusion: Unified Power-Law and Stretched Exponential Exact Solution. Physical Review E, 63, Article ID: 030101. https://doi.org/10.1103/PhysRevE.63.030101
[17]	Pedron, I.T., Mendes, R., Malacarne, L.C. and Lenzi, E.K. (2002) Nonlinear Anomalous Diffusion Equation and Fractal Dimension: Exact Generalized Gaussian Solution. Physical Review E, 65, Article ID: 041108. https://doi.org/10.1103/PhysRevE.65.041108
[18]	Abraham-Shrauner, B. (2005) Lie Symmetry Solutions for Anomalous Diffusion. Journal of Physics A: Mathematical and General, 38, 2547. https://doi.org/10.1088/0305-4470/38/12/001
[19]	O’Shaughnessy, B. and Procaccia, I. (1985) Analytical Solutions for Diffusion on Fractal Objects. Physical Review Letters, 54, 455-458. https://doi.org/10.1103/PhysRevLett.54.455
[20]	Crank, J. (1979) The Mathematics of Diffusion. Oxford University Press, Oxford.
[21]	Bluman, G. and Kumei, S. (1980) On the Remarkable Nonlinear Diffusion Equation (∂/∂x)[a (u + b) - 2(∂u/∂x)] - (∂u/∂t) = 0. Journal of Mathematical Physics, 21, 1019. https://doi.org/10.1063/1.524550
[22]	Williams, C.L., Bodmann, B.G. and Kouri, D.J. (2017) Fourier and Beyond: Invariance Properties of a Family of Integral Transforms. Journal of Fourier Analysis and Applications, 23, 660-678. https://doi.org/10.1007/s00041-016-9482-x
[23]	Williams, C.L. (2017) From Generalized Fourier Transforms to Coupled Supersymmetry. PhD Thesis, University of Houston, Houston.
[24]	Kouri, D., Pandya, N., Williams, C.L., Bodmann, B.G. and Yao, J. (2018) Point Transformations and Relationships among Linear Anomalous Diffusion, Normal Diffusion and the Central Limit Theorem. Applied Mathematics, 9, 178-197. https://doi.org/10.4236/am.2018.92013
[25]	Kouri, D.J., Zhang, M.M. and Zhang, D.S. (2017) On a New Nonlinear Image Filtering Technique. 11th International Conference on Signal Processing and Communication Systems, Gold Coast, Australia, 13-15 December 2017, 1-5.
[26]	Haberman, R. (1983) Elementary Applied Partial Differential Equations. Vol. 987, Prentice Hall, Englewood Cliffs.
[27]	Li, C. (2016) The Powers of the Dirac Delta Function by Caputo Fraction Derivatives. Journal of Fractional Calculus and Applications, 7, 12.
[28]	Jarad, F., Adjabi, Y., Baleanu, D. and Abdeljawad, T. (2018) On Defining the Distributions δr and (δ’)r by Conformable Derivatives. Advances in Difference Equations, 2018, Article No. 407. https://doi.org/10.1186/s13662-018-1865-7
[29]	Anteneodo, C., Dias, J. and Mendes, R. (2006) Long-Time Behavior of Spreading Solutions of Schrödinger and Diffusion Equations. Physical Review E, 73, Article ID: 051105. https://doi.org/10.1103/PhysRevE.73.051105
[30]	Toscani, G. (2005) A Central Limit Theorem for Solutions of the Porous Medium Equation. Journal of Evolution Equations, 5, 185-203. https://doi.org/10.1007/s00028-005-0183-1
[31]	Schwämmle, V., Nobre, F.D. and Tsallis, C. (2008) q-Gaussians in the Porous-Medium Equation: Stability and Time Evolution. The European Physical Journal B, 66, 537-546. https://doi.org/10.1140/epjb/e2008-00451-y