Table Of ContentGeneralized Fourier transform method for nonlinear anomalous diffusion equation
Jie Yao,1,∗ Cameron L. Williams,2 Fazle Hussain,1 and Donald J. Kouri2,3,†
1Texas Tech University, Department of Mechanical Engineering, Lubbock, Texas, USA, 79409
2University of Houston, Department of Mathematics, Houston, Texas, USA,77204
3University of Houston, Departments of Mechanical Engineering and Physics, Houston, Texas, USA,77204
(Dated: January 16, 2017)
The solution of a nonlinear diffusion equation is numerically investigated using the generalized
7 Fouriertransformmethod. Thisequationincludesfractaldimensionsandpower-lawdependenceon
1 theradialvariableandonthediffusionfunction. ThegeneralizedFouriertransformapproachisthe
0 extension of the Fourier transform method used for the normal diffusion equation. The feasibility
2 of theapproach is validated bycomparing thenumerical result with theexact solution for a point-
source. Themeritofthenumericalmethodisthatitprovidesawaytocalculateanomalousdiffusion
n
with an arbitrary initial condition.
a
J
2
I. INTRODUCTION of r only, it is a generalization of the diffusion equation
1
for fractal geometry [18]. It is the traditional nonlinear
] In the last few decades, anomalous diffusion has been diffusion equation when the diffusion coefficient depends
h extensively studied in a variety of physical applications, on ρ only [19, 20]. Analytical solutions of eq. (1) with
p
- such as turbulent diffusion [1], surface growth [2], trans- a point source have been reported in [15, 16], where an
p port of fluid in porous media [2], hydraulics problems ansatz for ρ is proposed as a general stretched Gaussian
m [3],etc. Theanomalousdiffusionisusuallycharacterized function. In [17], the same analytic solutions were also
o by the time dependence of mean-square displacement obtained by using Lie group symmetry analysis.
c (MSD) viz., r2 tσ. The MSD grows linearly with Motivatedbytheresearchongeneralizednonlineardif-
s. time (σ = 1)hforit∝he normal diffusion case. The process fusion, we propose here a numerical method for solving
c is called sub-diffusion for 0 < σ < 1 and super-diffusion eq. (1)usingageneralizedFouriertransform. Thegener-
si for σ > 1. The standard normal diffusion described by alizedFouriertransform(alsocalledtheΦntransform),is
y theGaussiandistributioncanbeobtainedfromtheusual anewfamilyofintegraltransformsdevelopedbyWillams
h Fokker-Planck equation with a constant diffusion coeffi- et. al. [21, 22]. These transforms share many properties
p cient and zero drift [4]. Extensions of the conventional of the Fourier transform. It can be employed to perform
[
Fokker-Planckequationhavebeen usedto study anoma- a more general frequency and time-frequency analysis.
1 lous diffusion. For example, anomalous diffusion can be In section II, a brief introduction of the generalized
v obtainedbytheusualFokker-Planckequation,butwitha Fourier transform is provided. The procedure of using
8 variablediffusioncoefficient[5,6]. Itcanalsobeachieved the generalized Fourier transform for solving the gener-
8
by incorporating nonlinear terms in the diffusion term, alized nonlinear diffusion equation is discussed in IIIA
4
3 or external forces [7–10]. In some approaches, fractional and IIIB. The method is validated by comparison be-
0 equations havebeenemployedto analyzeanomalousdif- tween analytical and numerical results. Finally, some
. fusion and related phenomena [11–14]. numerical results for a non-Delta function initial condi-
1
0 In this paper, we study the generalized nonlinear dif- tion are given in IIIC.
7 fusion equation including a fractal dimension d and a
1 diffusion coefficientwhichdepends onthe radialvariable
: and the diffusion function ρ [15–17] II. GENERALIZED FOURIER TRANSFORM
v
i
rX ∂∂ρt =K0rd1−1∂∂r(cid:16)rd−1−θ∂∂rρν(cid:17), (1) The generalized Fourier transform Φn is defined as
a
with the initial and boundary conditions Φnf(k)=Z ϕn(kx)f(x)dx, (4)
ρ(r,t )=ρ (r), (2)
0 0 where the integral kernel ϕ (ωx)=c (kx)+is (kx) is,
n n n
ρ( ,t)=0, (3)
∞ 1 η n
where r is the radial coordinate, and θ and ν are real cn(η)= 2|η|n−1/2J−1+21n(cid:16)|n| (cid:17), (5)
parameters. When the diffusion coefficient is a function
and
1 η n
sn(η)= 2sgn(η)|η|n−1/2J1−21n(cid:16)|n| (cid:17), (6)
∗ PreviouslyatUniversityofHouston, DepartmentofMechanical
Engineering;Email:[email protected],[email protected] where J (η) is the cylindrical Bessel function, and n is
ν
† Email:[email protected] the transform order, i.e., n=1,2, .
···
2
TheFouriertransform isthespecialcasewithn=1.
F
The Φ transformshares many properties of the Fourier
n t=0.2 Ex. Sol.
transform. Herewefocusontwopropertieswhichwillbe 0.9 t=0.2 Num. Sol.
used later. It is well-known that the Fourier transform 0.8 t=0.4 Ex. Sol.
t=0.4 Num. Sol.
preservesthefunctionalformofaGaussian;particularly, 0.7 t=0.6 Ex. Sol.
[g1] = g1 if g1(x) = exp( x2/2). For the generalized 0.6 t=0.6 Num. Sol.
FFouriertransform,wehaveΦ−ngn =gn,ifgn(x)=e−x22nn. ρ 0.5
Inaddition,the generalizedFouriertransformalsohas 0.4
the following derivative property: 0.3
∂ ∂ 0.2
Φ [ x2−2n f]=k−2nΦ f. (7)
n n 0.1
∂x ∂x
0
In [21], the Φ transform is developed for integer or-
n −0.1
der case. However,it can be easily extended to the non- 0 1 2 3 4 5
r
integercase;see[21]foradditionaldiscussionoftheprop-
erties of the Φ transform.
n FIG. 1: Comparison between exact (line) and numerical
(symbols) solutions with K =1, D =1 and θ =2.5.
0
III. SOLVING THE GENERALIZED DIFFUSION
EQUATION WITH GENERALIZED FOURIER
The solution to eq. (8) is then obtained by applying the
TRANSFORM
inverse Φ transform to ρ˜(k,t).
n
The exact solution to eq. (8) for a point source at the
It is well known that the Fourier transform can be
origin (i.e., ρ(r,t )=δ(r)) is given by [18]
usedto find the solutionfor the standarddiffusionequa- 0
tion [23]. Motivated by this idea, we explore employing λ 1 d/λ rλ
tnhoenlgienneearradliizffeudsFioonureiqeurattriaonns.formtosolvethegeneralized ρa(r,t)= dΓ(d/λ))(cid:16)K0λ2t(cid:17) exp(cid:16)− K0λ2t(cid:17).(13)
We validate the Φ transform method by comparing
n
the numerical results with the analytical solution. Fig.
A. The O’Shaugnessy-Procaccia anomalous (1) shows the analytical and the numerical solution for
diffusion equation on fractals
K =1,D =1,andθ =2.5atdifferenttimes. According
0
to the classification discussed in [15], this example is a
Let us first consider the generalizationof the diffusion subdiffusioncasewithθ >D(1 ν). Fromfig. (1),itcan
−
equation for fractal geometry, where the diffusion coeffi- be seenthatthe numericalsolutionis ingoodagreement
cient is a function of r only ( or ν = 1)[18]. Eq. (1) can with the analytical solution. In addition, we observe the
be reduced to short tail behaviours of the solution ρ(r,t).
∂ρ(r,t) K ∂ ∂
= 0 rd−1−θ ρ(r,t) . (8)
∂t rd−1(cid:16)∂r ∂r (cid:17)
B. Generalized nonlinear equation
InordertosatisfytheformoftheΦ transform,weapply
n
the following scaling relationship
Now we consider the generalized nonlinear diffusion
∂ ∂ equation with ν = 1. In [15], eq. (1) is analytically
()=drd−1 (). (9) 6
∂r · ∂rd · solved using a generalized stretched Gaussian function
approach:
Then eq. (8) becomes
ρ(r,t)=[1 (1 q)β(t)rλ]1/(1−q)/Z(t), (14)
∂ρ(r˜,t) ∂ ∂ − −
=K˜ r˜2−λ/d ρ(r˜,t), (10)
∂t 0∂r˜ ∂r˜ if1 (1 q)β(t)rλ 0,andρ(r,t)=0if1 (1 q)β(t)rλ <
− − ≥ − −
0. Here q =2 ν, and β(t) and Z(t) are functions given
where K˜0 =K0d2, r˜=rd, and λ=2+θ. in eq. (12) in−[15]. The same solution is derived in [16]
By applying the Φ transform to both sides and em-
n using Lie group symmetry method.
ploying the derivative identity (eq. (7)), we obtain the
Inordertosolvethegeneralizeddiffusionequationnu-
diffusion equation in the wavenumber domain
merically, we follow the procedure in IIIA, transform-
∂ρ˜ ing the spatial domain equation to the wavenumber do-
= K˜ kλ/dρ˜. (11)
∂t − 0 main using the Φn transform. Instead of eq. (11), the
wavenumber domain diffusion equation becomes
Eq. (11) can be exactly solved as
∂ρ
ρ˜(k,t)=e−K˜0kλ/dtρ˜0. (12) ∂t =−K0kλ/dρν. (15)
e
e f
3
1
0.14 t0=0.1 Initial t=t0 Initial
t=0.2 Ex. Sol. 0.9
t=t0+0.05
0.12 t=0.2 Num. Sol. 0.8 t=t0+0.15
t=0.4 Ex. Sol.
t=t0+0.25
0.1 tt==00..46 NExu.m S. oSl.ol. 0.7 Normal t=t0+0.25
t=0.6 Num. Sol. 0.6
0.08
ρ0.5
ρ 0.06
0.4
0.04 0.3
0.02 0.2
0.1
0
0
−0.02 0 0.5 1 1.5 2 2.5 3
0 5 10 15 20 25 30 r
r
FIG. 3: Numerical solution with K =1, D =1 and
(a) 0
θ =2.5 for the initial condition eq. (16).
0.14 t0=0.1 Initial
t=0.2 Ex. Sol.
C. Generalized diffusion for arbitrary initial
0.12 t=0.2 Num. Sol.
t=0.4 Ex. Sol. condition
t=0.4 Num. Sol.
0.1
t=0.6 Ex. Sol.
t=0.6 Num. Sol.
0.08 The numerical method provides us a way to solve the
generalized diffusion equation for the non-point-source
ρ 0.06
initial condition. In fig. (3), we present the numerical
0.04 solutionofthe generalizeddiffusionequationforK =1,
0
D =1 and θ =2.5 with the Gaussian initial condition
0.02
0 1 r2
ρ (r,t )= exp( ), (16)
0 0
−0.020 0.5 1 1.5 2 2.5 3 √4πt0 −4K0t0
r
where t =0.1.
(b) 0
As we can see, the diffusion process finally approaches
FIG. 2: Comparison between exact (line) and numerical the same generalized gaussian shape as in the point
(symbols) solutions in (a) scaled coordinate and (b) source case (fig. (1)). This is a numerical example of
original coordinate withK0 =1, D =3, θ =2.5 and whatamounts to a “generalizedcentrallimit” behaviour
ν =0.8. in which the diffusion process will finally transform the
arbitrary initial distribution to the corresponding gener-
alizedGaussiandistribution[24,25]. Bycomparingwith
the solution for normal diffusion with the same initial
condition, the subdiffusion process clearly exhibits the
There is no exact solution for eq. (15) except for ν =
short tail behaviour.
1. However, it can be numerically solved by employing
certaintypes oftime-stepping discretizationmethods for
the time derivative.
IV. CONCLUSION
The comparisons between the exact [15, 16] and nu-
merical solution for the Delta function initial condition
in scaled and original coordinates are shown in figs. In this paper,a numericalmethod for solving the gen-
(2a) and(2b) , respectively. To avoid performing a Φ eralized nonlinear diffusion equation has been presented
n
transform of the fractional order Delta function, we use and validated. The method is based on the generalized
ρ (r,t = 0.1) as the initial condition for our numerical Fourier transform and has been validated by comparing
a 0
algorithm. TheparametersusedhereareK =1,D =3, the numerical solution with analytical solution for the
0
θ =2.5, and ν =0.8. The forwardEuler finite difference point source. The presented method may be useful to
scheme with ∆t = 0.01s is employed for time discretiza- study a variety of systems involving the anomalous dif-
tioninthenumericalsimulation. Again,goodagreement fusion. Currently, no fast transform algorithm has yet
between the numericaland analyticalsolutions can been been developed for the Φ transform. This issue will be
n
observed. investigated in a future study.
4
ACKNOWLEDGMENTS Welch Foundation Grant E-0608 is gratefully acknowl-
edged. Partial support for this work was provided by
Discussion with Bernhard G. Bodmann are appreci- resources of the uHPC cluster managed by the Univer-
ated. Financial support of this research under R.A. sity of Houston under NFS Award Number 1531814.
[1] V.Gavrilov,N.Klepikova, andH.Rodean,Atmospheric [15] L.Malacarne,R.Mendes,I.Pedron, andE.Lenzi,Phys-
Environment 29, 2317 (1995). ical Review E 63, 030101 (2001).
[2] H.Spohn,Journal dePhysiqueI 3, 69 (1993). [16] I. T. Pedron, R. Mendes, L. C. Malacarne, and E. K.
[3] E.DalyandA.Porporato,PhysicalReviewE70,056303 Lenzi, Physical Review E 65, 041108 (2002).
(2004). [17] B.Abraham-Shrauner,JournalofPhysicsA:Mathemat-
[4] H. Risken, in The Fokker-Planck Equation (Springer, ical and General 38, 2547 (2005).
1984) pp.63–95. [18] B. O’Shaughnessy and I. Procaccia, Physical review let-
[5] K.S. Fa, Physical Review E 72, 020101 (2005). ters 54, 455 (1985).
[6] K.S. Fa, Physical Review E 84, 012102 (2011). [19] J. Crank, The mathematics of diffusion (Oxford univer-
[7] M.Bologna,C.Tsallis, andP.Grigolini,PhysicalReview sity press, 1979).
E 62, 2213 (2000). [20] G. Bluman and S. Kumei, Journal of Mathematical
[8] P. Assis Jr, P. da Silva, L. da Silva, E. Lenzi, and Physics 21, 1019 (1980).
M. Lenzi, Journal of mathematical physics 47, 3302 [21] C.L.Williams, B.G.Bodmann, andD.J.Kouri,arXiv
(2006). preprint arXiv:1403.4168 (2014).
[9] E.Lenzi,M.Lenzi,T.Gimenez, andL.daSilva,Journal [22] D. J. Kouri, C. L. Williams, N. Pandya, and H. Eisen-
of Engineering Mathematics 67, 233 (2010). berg, arXiv preprint arXiv:1610.04925 (2016).
[10] R. Zola, M. Lenzi, L. Evangelista, E. Lenzi, L. Lucena, [23] R. Haberman, Elementary applied partial differential
and L. da Silva,Physics Letters A 372, 2359 (2008). equations, Vol. 987 (Prentice Hall Englewood Cliffs, NJ,
[11] R.MetzlerandJ.Klafter,Physicsreports339,1(2000). 1983).
[12] I.M.Sokolov,J.Klafter, andA.Blumen,PhysicsToday [24] G.Toscani,JournalofEvolutionEquations5,185(2005).
55, 48 (2002). [25] V. Schw¨ammle, F. D. Nobre, and C. Tsallis, The Euro-
[13] C. Tsallis and E. Lenzi, Chemical Physics 284, 341 pean Physical Journal B 66, 537 (2008).
(2002).
[14] F.Liu,V.Anh, andI.Turner,JournalofComputational
and Applied Mathematics 166, 209 (2004).
