Numerical solution for diffusion equations with distributed order in time using a Chebyshev collocation method

Abstract

In this work we present a new numerical method for the solution of the distributed order time-fractional diffusion equation. The method is based on the approximation of the solution by a double Chebyshev truncated series, and the subsequent collocation of the resulting discretised system of equations at suitable collocation points. An error analysis is provided and a comparison with other methods used in the solution of this type of equation is also performed.

Introduction

In 1827, the botanist Robert Brown observed the intriguing jittering movement of small particles such as pollen grains, when these were immersed in water. Nowadays it is well known that this motion is caused by the rapid movement of water molecules, and insight into this problem was provided by Albert Einstein in 1905, in his work regarding Brownian motion, entitled “On the motion, required by the molecular-kinetic theory of heat, of particles suspended in fluids at rest.” [22]. His work served as a definitive confirmation that atoms and molecules actually exist. Although molecules are too small to be seen directly, their presence can be inferred from their visible effect on larger grains (such as the pollen grains). In his doctoral dissertation, Einstein developed a statistical molecular theory of liquids and in his subsequent paper he took the view that Brownian motion could be explained in terms of a type of stochastic process called a “random walk” (it is worth mentioning that Louis Bachelier, a student of Henri Poincaré, developed a theory of Brownian motion in his 1900 thesis [1] regarding stock market fluctuation [39]).

The random walk theory was popularised by Karl Pearson in his letter to Nature (1905) [52], where he proposed the following problem: a man starts from a point O and walks l yards in a straight line; after this step he turns through a random angle and takes another step of length l. He repeats this process n times. What is the probability that after n steps he is at a distance between r and $r+\delta r$ from the original point O? This problem is basically the idea behind Brownian motion, but the connection between the diffusion law proposed by Einstein,

$$\frac{\partial u(x,t)}{\partial t}=D\frac{{\partial }^{2}u(x,t)}{\partial x^2}$$

and the random walk process, was due to the excellent work of M. von Smoluchowski [58] in 1906.

In Brownian motion theory, the mean squared displacement, $<x^2>$, of a particle is proportional to the time t, $<x^2>=2Dt$, where D is the diffusion coefficient and $<x^2>$ is the second moment of the Gaussian distribution that governs the probability of being at position x at time t. This Gaussian distribution was observed experimentally, for example in the work of Kappler [36], and, if we track the diffusion of a particle that follows the governing equation given by equation (1) and is at the origin at time $t=0$, then we obtain the solution as the normalised Gauss probability density function,

$$u(x,t)=\frac{1}{\sqrt{4\pi Dt}}\exp (-\frac{x^2}{4Dt}).$$

Although this linear growth of the mean squared displacement is observed in most diffusion problems, there is a considerable number of physical processes that do not show this behaviour [4], [48], [47], [32], [49], but present an anomalous diffusion, that can be subdiffusive or supperdifusive. Note that, by 1926, Lewis Fry Richardson [53] studied the relative diffusion of two tracer particles in a turbulent flow and found that $<x^2>\simeq t^3$, meaning that a supperdifusive behaviour was obtained.

In general, anomalous diffusion refers to the following power-law form,

$$<x^2>\simeq D_{\alpha }{t}^{\alpha }$$

where we have subdiffusion for $0<\alpha <1$ and superdiffusion for $\alpha >1$.

We may question how we can obtain a governing equation for a diffusion process that behaves like equation (3). The first idea is to take a look at the physical process under examination and find a suitable stochastic model. The problem is that we may find different models for different problems, leading to the difficult task of finding a universal theory for such non-universal and process dependent behaviour.

In order to understand the variety of possibilities, we can take, for example, the continuous time random walk process (CTRW) [50], [56], that generalises the random walk proposed by Pearson. In the CTRW we may have jumps of different length and a particle may wait before the next jump. Different distributions may be used for the waiting time and jump length and, therefore, different regimes of diffusion may be obtained, allowing in this way the development of a general theory that can be tuned according to the characteristics of the physical process.

For example, in the particular case of Lévy flights [45], [23], [34], [61], [26] there is no characteristic size for the random walk jumps, and therefore, the variance or second moment is infinite, and does not follow a law such as the one given by equation (3). This does not mean that the Lévy flights are not important for modelling real world diffusion phenomena. In fact they find a wide range of applications including chaotic phase diffusion in Josephson junctions, turbulent diffusion, micelle dynamics, vortex dynamics, anomalous diffusion in rotating flows, etc. [39]. It is all just a question of definition of the mean squared displacement. If we ask for the mean squared length of a jump then the answer is infinity, but, if instead we try to understand the properties of the distribution used in the Lévy flights by asking how far a Lévy walk has travelled from its starting point at time t, then the answer would be a time-dependent moment of the probability distribution [39].

As a second example [49], [15], [57], [55], let us consider a finite variance of the jump length and the wait time assumed to be given by a power-law distribution $\phi \left(t\right)={\tau }^{\alpha }/t^{1+\alpha }$ $(0<\alpha <1)$ that leads to a divergent mean waiting time. A process modelled in this way leads to subdiffusion with a mean squared displacement given by $<x^2>=(2K_{\alpha }/\mathrm{\Gamma }(1+\alpha ))t^{\alpha }$. The governing equation, that tells us the probability of finding a particle at position x at time t, is given by

$$\frac{{\partial }^{\alpha }u(x,t)}{\partial t^{\alpha }}=K_{\alpha }\frac{{\partial }^{2}u(x,t)}{\partial x^2}$$

where the fractional derivative on the left-hand-side of the equation is of Caputo type [19], [55]:

$$\frac{{\partial }^{\alpha }u(x,t)}{\partial t^{\alpha }}=\frac{1}{\mathrm{\Gamma }(1-\alpha )}\underset{0}{\overset{t}{\int }}{(t-s)}^{-\alpha }\frac{{\partial }^{n}u(x,s)}{\partial s^n}ds,$$

with n being the smallest integer greater than or equal to α. This equation is known as a time-fractional diffusion equation (TFDE) and has found application in a broad variety of engineering, biological, finance and physical processes where anomalous diffusion occurs (see, for example [31], [44] and [49]).

Regarding fractional differential calculus, it is wrong to think that this is a recent subject. For instance, the symbol $d^{n}y/d{x}^n$ was first proposed by Leibniz, and, in 1695 L'Hôpital asked Leibniz the meaning of $d^{n}y/d{x}^n$ if $n=1/2$. Several well known mathematicians tried to develop this general theory. In 1697, Leibniz sent letters to J. Wallis and J. Bernoulli, and, he mentioned the possible approach to fractional-order differentiation. In 1716 Leibniz died, but the urge to understand the derivatives of fractional order increased, and several other authors devoted their time to the subject. The well known Leonard Euler also contributed to the understanding and generalisation of fractional differential calculus. He generalised the notion of factorial n! to non-integer values, called the gamma function $\left(\mathrm{\Gamma }(.)\right)$, by Adrien-Marie Legendre around 1811. Later, Joseph Fourier (1768–1830) gave the first step to the generalisation of the notion of differentiation for arbitrary functions with his 1822 book, Theorie Analytique de la Chaleur. In 1832, Liouville presented two different definitions for the fractional derivative. Riemann also contributed to this subject and the ideas of Riemann and Liouville, led to the Riemann–Liouville definition of fractional integrals and derivatives. Grünwald (1867) and Post (1930) presented the idea of fractional derivative as the limit of a sum, using the classical definitions of a derivative. In 1927, Marchaud formulated an equivalent fractional derivative of arbitrary order α (with $0<\alpha <1$). In 1940, Erdélyi and Kober presented what we call now the Erdélyi–Kober fractional integral. This operator generalises the Riemann fractional integral and the Weyl integral. Later, Riesz also proposed a fractional integral of order α, and in 1967 Caputo developed the Caputo fractional derivative presented in equation (5). Several other works can be found in the literature that contributed to the development of fractional calculus and the interested reader should consult [54]. It remains impossible to say that one definition is better than the other one (although it is fair to say that the Caputo definition is most often used for the application of fractional calculus to physics).

Returning to equation (4), it is worth mentioning that in recent decades great attention has been devoted to the general form of equation (4) given by,

$$\frac{{\partial }^{\alpha }u(x,t)}{\partial t^{\alpha }}=D_{\alpha }\frac{{\partial }^{2}u(x,t)}{\partial x^2}+f(x,t),0<t\le a,0<x<b,$$

where $D_{\alpha }$ stands for a general diffusion coefficient with dimensions ${\left[length\right]}^2/{\left[time\right]}^{\alpha }$.

Equation (6) is then useful for diffusion processes that follow (3) but, in the real world, there are more complex systems that do not follow equation (3) and that may present more than one characteristic scale. For example we saw that in the Lévy flights there is an absence of a characteristic scale. This means that fractional differential equations are a good tool for modelling some of the processes, but, in order to model correctly more complex systems a more powerful mathematical tool should be used.

The answer to this problem of non-unique scalability was given by the creation of distributed order fractional differential equations [9], [10] that proved to be useful in modelling anomalous diffusion characterised by two or more scaling exponents in the mean squared displacement (equation (3)), more precisely, in modelling processes for which the diffusion exponent can change in the course of time, as is the case of complex systems where the morphology of the material changes along the process. A close look into the physics of some complex diffusive processes, suggests that an even more general theory for fractional derivatives should be devised, and, this will no doubt be developed in the near future [55], [47], [49], [39].

As Chechkin, Gorenflo and Sokolov [9] observed, “the development of numerical schemes for solving distributed-order kinetic equations and for modelling sample paths of the random processes governed by these equations is also of importance.”, meaning that numerical schemes for the solution of this type of equation are demanded. This is the main motivation for the present work, where we present a numerical method for the solution of the general distributed-order time fractional diffusion equation,

$$\underset{0}{\overset{1}{\int }}c(\alpha )\frac{{\partial }^{\alpha }u(x,t)}{\partial t^{\alpha }}d\alpha =\frac{{\partial }^{2}u(x,t)}{\partial x^2}+f(x,t),0<t\le a,0<x<b,$$

where the function $c(\alpha )$ acting as weight for the order of differentiation is such that [30], [43] $c(\alpha )\ge 0$ and ${\int }_0^{1}c(\alpha )d\alpha =C>0$. Obviously, if $c(\beta )=\frac{\delta (\alpha )}{D^{\alpha }}$, where $\delta ()$ is the delta Dirac function, then (7) reduces to (6). Note that the dimensions of $c(\alpha )$ are ${\left[time\right]}^{\alpha }/{\left[length\right]}^2$. For an alternative to these dimensions please consult the work of Chechkin et al. [10] where they use an extra constant that represents a characteristic time of the problem.

Here, we will be interested in the numerical approximation of this type of distributed-order equation with boundary conditions of Dirichlet type:

$$u(0,t)={\varphi }_0(t),u(b,t)={\varphi }_b(t),0<t<a,$$

and initial condition

$$u(x,0)=g_0(x),0<x<b.$$

The concept of distributed order was developed by Caputo [6] and further developed by Caputo, Bagley and Torvik [7], [8], [2], [3].

With the growing interest on this type of equation, numerical methods started being developed, for example in the work by Diethelm and Ford [16] where they present a basic framework for the numerical solution of distributed order differential equations (see also [17], [18]). A more recent increase in interest in the use of distributed order differential equations (particularly in the case where the derivatives are given in the Caputo sense) led Ford and Morgado [24] to discuss the existence and uniqueness of solutions for this type of equation, and also to propose a numerical method for their approximation in the case where the initial conditions are not known (with boundary conditions being given away from the origin). Two years later, Katsikadelis [37] devised a numerical method for the solution of the distributed order FDE approximating them with a multi-term FDE (that is solved by adjusting appropriately the numerical method developed for multi-term FDEs). In the same year, Liao et al. [41] investigated a class of modified Du Fort–Frankel-type schemes for fractional subdiffusion equations in the Jumarie's modified Riemanniouville form with constant, variable or distributed fractional order.

In 2015, some papers were published on the numerical solution of distributed order FDEs. Morgado and Rebelo [51] presented an implicit scheme for the numerical approximation of the distributed order time-fractional reaction–diffusion equation with a nonlinear source term (see also [25]), Ye et al. presented two papers, one [62] considering the time distributed-order and Riesz space fractional diffusion on bounded domains with Dirichlet boundary conditions (deriving an implicit difference method for the multi-term time-space fractional diffusion equation) and the other [63] presenting a numerical method based on a compact difference scheme for a distributed order time-fractional diffusion-wave equation. Hu et al. [33] considered a new time distributed-order and two-sided space-fractional advection–dispersion equation that was solved numerically using an implicit method for the solution the multi-term fractional equation. Gao et al. published a series of papers [27], [28], [29] where: [27] two difference schemes were derived for both one-dimensional and two-dimensional distributed-order differential equations (he proved that the schemes are unconditionally stable and convergent); [28] the Grünwald formula was used to solve the one-dimensional distributed-order differential equations (two difference schemes were derived and the extrapolation method was applied to improve the approximate accuracy); [29] two alternating direction implicit difference schemes were derived for two-dimensional distributed-order fractional diffusion equations (he proved that the schemes are unconditionally stable and convergent). Wang et al. [60] derived and analysed a second-order accurate implicit numerical method for the Riesz space distributed-order advection–dispersion equation. Duong et al. [21] proposed a novel numerical scheme for analysing the behaviour of a distributed order linear single input single output control system under random forcing. The method is based on the operational matrix technique to handle stochastic distributed order systems. Jin et al. [35] presented a numerical solution of an initial boundary value problem for the distributed order time fractional diffusion equation. They developed a space semidiscrete scheme based on the standard Galerkin finite element method, and established error estimates for both smooth and nonsmooth initial data.

Finally, Chen et al. [14] developed a mixed finite difference/spectral method, and, Li et al. [40] proposed a numerical method for solving distributed order diffusion equations by using a classical numerical quadrature formula, and the resulting multi-term time-fractional diffusion equation were solved by the reproducing kernel method.

In almost all the methods described so far, only finite difference approximations have been considered for the fractional time derivative, and these methods may become computationally heavy, due to the non-local property of fractional differential operators. Many of these authors assumed certain regularity assumptions on the solution in order to provide the convergence analysis of their numerical schemes, although it is widely known that the solution of fractional differential equations may be nonsmooth at $t=0$ even if the data is infinitely smooth. However, as Diethelm and Ford pointed out in [18], under certain simple hypotheses, solutions of the distributed order equations can be shown to be smooth. We have adapted their analysis (to be published in a sequel) for the distributed order diffusion equation and this provides a theoretical basis for our convergence analysis in this paper.

On the other hand, the idea of approximating the solution of fractional differential equations with truncated Chebyshev series has been widely exploited (see, for example, [38], [20] and the references cited therein), but to the best of our knowledge, a complete and detailed error analysis has not yet been provided.

Hence, this paper has two purposes: the first one is to develop a numerical scheme with a lower computational cost, based on the approximation of the solution with a double truncated Chebyshev series, and the subsequent collocation of the resulting discretised system of equations at suitable collocation points; the second one is to provide an error analysis for such methods, which, as far as we are aware, has not previously been presented for fractional differential equations.

The paper is organised as follows. After this Introduction, we present, in Section 2 some preliminaries results that will be useful in the derivation of the numerical method, which will be described in Section 3. In Section 4 we provide a convergence analysis of the numerical scheme. In Section 5 we test the numerical method considering different examples and we compare them with previous ones. The paper ends with the conclusions in Section 6.