JOURNAL OF THEORETICAL AND APPLIED MECHANICS 44, 2, pp. 393-403, Warsaw 2006 NUMERICAL SOLUTIONS TO BOUNDARY VALUE PROBLEM FOR ANOMALOUS DIFFUSION EQUATION WITH RIESZ-FELLER FRACTIONAL OPERATOR Mariusz Ciesielski Jacek Leszczynski Institute of Mathematics and Computer Science, Czestochowa University of Technology e-mail: mariusz@imi.pcz.pl; jale@imi.pcz.pl In this paper, we present a numerical solution to an ordinary differential equation of a fractional order in one-dimensional space. The solution to this equation can describe a steady state of the process of anoma- lous diffusion. The process arises from interactions within complex and non-homogeneous background. We present a numerical method which is based on the finite differencesmethod.We consider a boundary value problem (Dirichlet conditions) for an equationwith theRiesz-Feller frac- tional derivative. In the final part of this paper, some simulation results are shown.We present an example of non-linear temperature profiles in nanotubes which can be approximated by a solution to the fractional differential equation. Key words: Riesz-Feller fractional derivative, boundary value problem, Dirichlet conditions, finite difference method 1. Introduction Anomalous diffusion is a phenomenon strongly connected with interac- tions within complex and non-homogeneous background. This phenomenon is observed in transport of a fluid in porous materials, in chaotic heat baths, amorphous semiconductors, particle dynamics inside a polymer network, two- dimensional rotating flow and also in econophysics. The phenomenon of ano- malous diffusion deviates from the standard diffusion behaviour. In opposite to the standard diffusion where a linear form in the mean square displace- ment 〈x2(t)〉 ∼ k1t of the diffusing particle over time occurs, the anomalous diffusion is characterized by a non-linear one 〈x2(t)〉∼ kγt γ, for γ ∈ (0,2]. In 394 M. Ciesielski, J. Leszczynski this phenomenon, there may exist a dependence 〈x2(t)〉 → ∞ characterized by occurrence of rare but extremely large jumps of diffusing particles – well- known as the Levy motion or the Levy flights. An ordinary diffusion follows Gaussian statistics and Fick’s second law for finding the running process at time t whereas anomalous diffusion follows non-Gaussian statistics or can be interpreted as the Levy stable densities. Many authors proposedmodels which base on linear and non-linear forms of differential equations. Such models can simulate anomalous diffusion but they do not reflect its real behaviour. Several authors (Carpinteri and Ma- inardi, 1997; Hilfer, 2000; Gorenflo and Mainardi, 1998; Metzler and Klafter, 2000) apply fractional calculus to themodelling of this type of diffusion. This means that time and spatial derivatives in the classical diffusion equation are replaced by fractional ones. In comparison to derivatives of the integer order, which depend on the local behaviour of the function, the derivatives of the fractional order accumulate the whole history of this function. In our previeus works (Ciesielski and Leszczynski, 2003, 2005), we presen- ted a solution to a partial differential equation of the fractional orderwith the time fractional operator and the space ordinary operator, respectively. Those solutions were based on the Finite Difference Method (FDM) and are called the Fractional FDM (FFDM). 2. Mathematical background In this paper we consider an ordinary differential equation of fractional order in the following form dα d|x|α θ T(x)= 0 x∈ R (2.1) where T(x) is a field variable (i.e. field temperature), (dα/d|x|αθ)T(x) is the Riesz-Feller fractional operator (Metzer andKlafter, 2000; Samko et al., 1993), α is the real order of this operator and θ is the skewness parameter. Accor- ding to (Gorenflo and Mainardi, 1998), the Riesz-Feller fractional operator is defined as dα d|x|α θ T(x)= xD α θT(x)=−[cL(α,θ)−∞D α xT(x)+ cR(α,θ)xD α +∞T(x)] (2.2) Numerical solutions to boundary value problem... 395 for 0<α¬ 2, α 6=1where −∞D α xT(x)= ( d dx )m [−∞I m−α x T(x)] (2.3) xD α +∞T(x)= (−1) m ( d dx )m [xI m−α +∞ T(x)] for m ∈ N, m− 1 < α ¬ m, and the coefficients cL(α,θ), cR(α,θ) (for 0<α¬ 2, α 6=1, |θ| ¬min(α,2−α)) are defined as cL(α,θ)= sin (α−θ)π 2 sin(απ) cR(α,θ)= sin (α+θ)π 2 sin(απ) (2.4) The fractional integral operators of the order α: −∞I α xT(x) and xI α ∞ T(x) are definedas the left- and right-handofWeyl’s fractional integrals (Carpinteri andMainardi, 1997; Oldhamand Spaner, 1974; Podlubny, 1999; Samko et al., 1993) whose definitions are −∞I α xT(x)= 1 Γ(α) x∫ −∞ T(ξ) (x− ξ)1−α dξ (2.5) xI α ∞ T(x)= 1 Γ(α) ∞∫ x T(ξ) (ξ−x)1−α dξ Considering Eq. (2.1), we obtain the steady state of the classical diffusion equation for α=2, i.e. theheat transfer equation. If α=1, and theparameter of skewness θ admits extreme values in (2.4), the steady state of a transport equation is noted. Therefore we assume variations of the parameter α within the range 0<α¬ 2.Analysingbehaviour of theparameter α< 2 inEq. (2.1), we found some combination between transport and propagation processes in steady states. In this work, we will consider equation (2.1) in one dimensional domain Ω : L ¬ x ¬ R with boundary-value conditions of the first kind (Dirichlet conditions) as { x=L : T(L)= gL x=R : T(R)= gR (2.6) 396 M. Ciesielski, J. Leszczynski 3. Numerical method According to the finite difference method (Hoffman, 1992; Majchrzak and Mochnacki, 1996), we consider a discrete from of equation (2.1) in space. The problemof solving equation (2.1) lies in a properapproximation ofRiesz-Feller derivative (2.2) by a numerical scheme. 3.1. Approximation of the Riesz-Feller derivative In order to develop a discrete form of operator (2.2), we introduce a homo- genous spatial grid −∞< ... < xi−2 N (3.6) On the basis of above considerations, we modify expression (3.1) for the di- scretization of the Riesz-Feller derivative. Thus we have xiD α θT(xi)≈ 1 hα (N−i∑ k=−i Ti+kw (α,θ) k +gLsL (α,θ) i +gRsR (α,θ) N−i ) (3.7) for i=1, . . . ,N−1, where sL (α,θ) j = −j−1∑ k=−∞ w (α,θ) k =    cL(α,θ)r̃j for 0<α< 1 cL(α,θ) ˜̃rj for 1<α¬ 2 (3.8) sR (α,θ) j = ∞∑ k=j+1 w (α,θ) k =    cR(α,θ)r̃j for 0<α< 1 cR(α,θ) ˜̃rj for 1<α¬ 2 Numerical solutions to boundary value problem... 399 and r̃j = (j+2)1−αλ1+(j+1) 1−α(2−2λ1)+ j 1−α(λ1−2) 2Γ(2−α) (3.9) ˜̃rj = (j+2)2−α(2−λ2)+(j+1) 2−α(3λ2−4)+ j 2−α(2−3λ2) 2Γ(3−α) + + (j−1)2−αλ2 2Γ(3−α) Putting expression (3.7) to equation (2.1), we obtain the following finite difference scheme N−i∑ k=−i Ti+kw (α,θ) k +gLsL (α,θ) i +gRsR (α,θ) N−i =0 i=0, . . . ,N (3.10) T0 = gL TN = gR The above scheme described by expression (3.12) can be written in a matrix form as A ·T =B (3.11) where A=   1 0 0 0 . . . 0 0 0 a−1 a0 a1 a2 . . . aN−3 aN−2 aN−1 a−2 a−1 a0 a1 . . . aN−4 aN−3 aN−2 a−3 a−2 a−1 a0 . . . aN−5 aN−4 aN−2 a−4 a−3 a−2 a−1 . . . aN−6 aN−3 aN−4 ... ... ... ... ... ... ... ... a−N+2 a−N+3 a−N+4 a−N+5 . . . a0 a1 a2 a−N+1 a−N+2 a−N+3 a−N+4 . . . a−1 a0 a1 0 0 0 0 . . . 0 0 1   (3.12) B= [gL,b1,b2,b3,b4, . . . ,bN−2,bN−1,gR] > with aj =w (α,θ) j for j=−N+1, . . . ,N −1 bj = gLsL (α,θ) j +gRsR (α,θ) j for j=1, . . . ,N −1 (3.13) 400 M. Ciesielski, J. Leszczynski and T = [T0,T1,T2, . . . ,TN] > is the vector of unknown values of the func- tion T . We can observe that the boundary conditions influence all values of the function in everynode. Inopposite to the secondderivative over space,which is approximated locally, the characteristic feature of the Riesz-Feller and other fractional derivatives is the dependence on values of all domain points. For α = 2 and θ = 0, our scheme is identical as with the well known and used central difference scheme in space (Hoffman, 1992;Majchrzak andMochnacki, 1996). The skewness parameter θ has significant influence on the solution. For α→ 1+ and θ →±1+ one can obtain the classical ordinary differential equ- ation. 4. Simulation results In this section, we present results of calculation. In all presented si- mulations we assumed 0 ¬ x ¬ 1. Figure 2 shows temperature profiles over space with boundary conditions gL = 2, gR = 1 for different values α= {0.1,0.5.0.75,1.01,1.25,1.5,1.75,2} and θ=0. Figure 3 presents another example of the solution in which we assumed α=1.01 and different values of the skewness parameter θ= {0,0.1,0.3,0.5,0.7,0.9,0.99}. Fig. 2. Spatial solution to equation (2.1) for different values of α and θ=0 Numerical solutions to boundary value problem... 401 Fig. 3. Spatial solution to equation (2.1) for a steady value of the parameter α=1.01 and different values of θ=0 The last example illustrates the heat transport in nanotubes. Figure 4 presents experimental data prerformed by Zhang and Li (2005) and results of our numerical solution of equation (2.1). For such a comparison we assumed α=0.35 and θ=−0.055 to best fit the experimental data. It should be noted that the temperature profile inside the nanotube deviates from the profile obtained by solving the standard heat transfer equation. Fig. 4. Comparison of temperature profiles in nanotubes measured by Zhang and Li (2005) and obtained by numerical solution to equation (2.1) 402 M. Ciesielski, J. Leszczynski 5. Conclusions Summing up, we proposed a fractional finite difference method for the fractional diffusion equation with the Riesz-Feller fractional derivative which is an extension of the standard diffusion. We analysed a linear case of the steady state of the anomalous diffusion equation, and in the future we will work on non-linear cases. We obtained the FDM scheme which generalises classical scheme of FDM for the diffusion equation. Moreover, for α = 2, our solution is equivalent with that obtained by the classical finite difference method. Analysing plots included in this work, we can see that in the case α=2, solution of Eq. (2.1) is a linear function. In the other cases, when α < 2, solutions are non-linear functions. When we analyse the probability density function generated by fractional diffusion equation for α < 2, we observe a long tail of distribution. In this way we can simulate same rare and extreme eventswhich are characterised by arbitrary very large values of particle jumps. Analysing the changes in the skewness parameter θ, we observed intere- stingbehaviour in the solution.For α→ 1+ and for θ→±1+,we obtained the steady state of the first order wave equation. For θ∈ (0,1) (with restrictions to the order α), we generated non-symmetric solutions. It should be noted that we can good approximate the temperature profile inside the nanotubes using solution of Eq. (2.1). References 1. Carpinteri A., Mainardi F. (eds.), 1997,Fractals and Fractional Calculus in Continuum Mechanics, Springer Verlag, Vienna-New-York 2. Ciesielski M., 2005, Fractional finite difference method applied for solution of anomalous diffusion equations with initial-boundary conditions, PhDThesis, Czestochowa, (in Polish) 3. Ciesielski M., Leszczynski J., 2003, Numerical simulations of anomalous diffusion, Proc. 15th International Conference on Computer Methods in Me- chanics CMM-2003, 1-5 (proceeding on CD-ROM) 4. CiesielskiM.,Leszczynski J., 2005,Numerical solutionsofaboundaryvalue problem for the anomalous diffusion equation with the Riesz fractional deriva- tive, Proc. 16th International Conference on Computer Methods in Mechanics CMM-2005, 1-5 (proceeding on CD-ROM) Numerical solutions to boundary value problem... 403 5. GorenfloR.,Mainardi F., 1998, Fractional calculus and stable probability distributions,Archives of Mechanics, 50, 3, 377-388 6. Hilfer R., 2000,Applications of Fractional Calculus in Physics, World Scien- tific Publ. Co., Singapore 7. Hoffman J.D., 1992, Numerical Methods for Engineers and Scientists, McGraw-Hill 8. Majchrzak E., Mochnacki B., 1996,Metody numeryczne, Podstawy teore- tyczne. Aspekty praktyczne i algorytmy, WydawnictwoPolitechniki Slaskiej (in Polish), Gliwice 9. Metzler R., Klafter J., 2000, The randomwalk’s guide to anomalous dif- fusion: a fractional dynamics approach,Phys. Rep., 339, 1-70 10. OldhamK., Spanier J., 1974,The Fractional Calculus, AcademicPress,New York and London 11. Podlubny I., 1999, Fractional Differential Equations, Academic Press, San Diego 12. Samko S.G., Kilbas A.A., Marichev O.I., 1993, Integrals and Derivati- ves of Fractional Order and Same of their Applications, Gordon and Breach, London 13. Zhang G., Li B., 2005, Thermal conductivity of nanotubes revisited: Effects of chirality, isotope impurity, tube length, and temperature, J. Chem. Phys., 123, 114714 Numeryczne rozwiązanie zagadnienia brzegowego równania anomalnej dyfuzji z operatorem frakcjalnym Riesza-Fellera Streszczenie W pracy zaprezentowano numeryczne rozwiązanie jednowymiarowego równania różniczkowego zwyczajnego niecałkowitego rzędu. Rozwiązanie tego równania może opisywać stan ustalony procesu anomalnej dyfuzji. Proces ten wynika z oddziaływań zachodzącychw złożonych i niejednorodnych systemach. Zaprezentowanametoda nu- meryczna oparta jest na metodzie różnic skończonych. Rozważane było zagadnienie brzegowe z warunkami Dirichleta dla tego równania z pochodną frakcjalną Riesza- Fellera.W końcowej części przedstawiono wyniki symulacji. Jako przykład zaprezen- towano nieliniowy profil temperatury w nanorurkach, który może być przybliżony przez rozwiązanie frakcjalnego równania różniczkowego. Manuscript received December 15, 2005; accepted for print January 23, 2006