**Applications of Radial Basis Function Schemes to Fractional Partial Differential Equations**

Carlos Alberto Torres Martínez and Carlos Fuentes

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/67892

#### Abstract

In modeling using diffusion equations, the relationship between fractal geometry and fractional calculus arises by modeling the conditions of the medium as a fractal whose fractional dimension determines the order of this equation. For this reason, it is very useful to have numerical methods that solve them and discretization of the domain is not determinant for the efficiency of the algorithm. In this work, it is proposed to show that meshless methods, in particular methods with radial basis functions (RBF), are an alternative to schemes in differences or structured meshes. We show that we can obtain numerical solutions to some fractional partial differential equations using collocation and RBF, over non equally distributed data.

Keywords: fractional partial differential equations (FPDE), meshless methods, radial basis functions (RBF), Caputo derivative, Riemann-Liouville derivative, Riesz derivative, diffusion-convection

## 1. Background

#### 1.1. Radial basis function methodology

The Hardy-based radial-based functions (RBF) methodology [1] arises from the need to apply multivariate interpolation to cartography problems, with randomly dispersed data (also known as collocation nodes). Micchelli [2], Powell et al. [3] gave it a great boost by proving non-singularity theorems. Later, Kansa [4, 5] proposed to consider the analytical derivatives of the FBR to develop numerical schemes that deal with partial differential equations (PDE).

Regarding PDE over spaces of dimension greater than one, we generally opt for finite element (FE) type discretizations on meshes, structured or not; also pseudo spectral (PS) methods through base functions such as Fourier or Chebyshev. The high degree of computational

efficiency in the procedure raises the cost in the regularity constraints on the form of the computational domain.

The FE methods involve decomposition of the domain; for example, in two dimensions rectangles are constructed with curvilinear mappings that allow the refinement of the mesh in critical areas. However, this type of implementation is complex and very close nodes are needed, mainly at the domain boundaries, which impairs stability conditions in time.

For this reason, we look for numerical techniques that do not depend to a great extent on data distribution. RBF collocation methods belong to the "meshfree" methods; that is, they only require scattered collocation nodes on the domain and the boundary. They are also an alternative for dealing with problems in larger dimensions and irregular domains. Hence, in recent decades, such methods have attracted the attention of researchers in order to solve partial differential equations. See for example: Chen et al. [6].

This type of technique approximates the solution by means of a linear combination of radial basis functions, which are globally defined and of exponential convergence, but they produce dense and badly conditioned interpolation matrices. On the other hand, some RBFs contain a shape parameter, a number that influences the precision of the numerical results [7, 8].

There are algorithms that also produce better conditioned interpolants, even when the shape parameter tends to zero, with a double precision arithmetic. Some of them are: Contour-Padé [9] and RBF-QR [10]. Related to the latter is the recent RBF-GA [11].

The RBF-QR method provides a numerically more stable alternative for converting basis that are very similar to each other, or almost linearly dependent on a set of scattered nodes, to a very well conditioned base that generates exactly the same space. This method was originally implemented for nodes on the surface of a sphere in Fornberg and Piret [10]. And recently for an arbitrary set of nodes in one, two, and three dimensions [12, 13].

## 1.2. Fractional calculus and differential equations: about some applications

The growing interest in fractional calculus has been motivated by applications of fractional equations in different fields of research.

One of the applications considered in this work is that of convection-diffusion equations, which appear not only in many applications of models in Physics and Chemistry [14–17], in problems of flow or heat transfer, but also in other fields such as financial [18, 19]. The solution of a convection-diffusion problem can be interpreted as a probability distribution of one or more underlying stochastic processes [18]. The classical diffusion equation (or heat equation) and its Gaussian solution existed long before Einstein established a connection with random walks.

The anomalous diffusion equations, on the other hand, were originally developed from stochastic walks models [20]. The parameters of these equations are uniquely determined by the fractal dimension of the underlying object [21].

Recent applications of these anomalous diffusion-convection models—those related to oil extraction and hydrological models for aquifers, food production and water distribution in large cities, are quite important [22]. Determining the behavior of the fluid inside the reservoir and the loss of permeability of the medium assists in the investigation of the mechanisms of oil migration. In this context, the geometric and physical interpretation of fractional derivatives in these differential equations is a "ultra slow" or "ultra fast" diffusion.

efficiency in the procedure raises the cost in the regularity constraints on the form of the

The FE methods involve decomposition of the domain; for example, in two dimensions rectangles are constructed with curvilinear mappings that allow the refinement of the mesh in critical areas. However, this type of implementation is complex and very close nodes are

For this reason, we look for numerical techniques that do not depend to a great extent on data distribution. RBF collocation methods belong to the "meshfree" methods; that is, they only require scattered collocation nodes on the domain and the boundary. They are also an alternative for dealing with problems in larger dimensions and irregular domains. Hence, in recent decades, such methods have attracted the attention of researchers in order to solve partial

This type of technique approximates the solution by means of a linear combination of radial basis functions, which are globally defined and of exponential convergence, but they produce dense and badly conditioned interpolation matrices. On the other hand, some RBFs contain a

There are algorithms that also produce better conditioned interpolants, even when the shape parameter tends to zero, with a double precision arithmetic. Some of them are: Contour-Padé [9]

The RBF-QR method provides a numerically more stable alternative for converting basis that are very similar to each other, or almost linearly dependent on a set of scattered nodes, to a very well conditioned base that generates exactly the same space. This method was originally implemented for nodes on the surface of a sphere in Fornberg and Piret [10]. And recently for

The growing interest in fractional calculus has been motivated by applications of fractional

One of the applications considered in this work is that of convection-diffusion equations, which appear not only in many applications of models in Physics and Chemistry [14–17], in problems of flow or heat transfer, but also in other fields such as financial [18, 19]. The solution of a convection-diffusion problem can be interpreted as a probability distribution of one or more underlying stochastic processes [18]. The classical diffusion equation (or heat equation) and its Gaussian solution existed long before Einstein established a connection with random walks.

The anomalous diffusion equations, on the other hand, were originally developed from stochastic walks models [20]. The parameters of these equations are uniquely determined by the

Recent applications of these anomalous diffusion-convection models—those related to oil extraction and hydrological models for aquifers, food production and water distribution in large cities, are quite important [22]. Determining the behavior of the fluid inside the reservoir

shape parameter, a number that influences the precision of the numerical results [7, 8].

needed, mainly at the domain boundaries, which impairs stability conditions in time.

differential equations. See for example: Chen et al. [6].

4 Fractal Analysis - Applications in Physics, Engineering and Technology

and RBF-QR [10]. Related to the latter is the recent RBF-GA [11].

an arbitrary set of nodes in one, two, and three dimensions [12, 13].

equations in different fields of research.

fractal dimension of the underlying object [21].

1.2. Fractional calculus and differential equations: about some applications

computational domain.

The relationship between fractal geometry and fractional calculus is given from the consideration that a particle moves in a porous medium with fractal structure [23]. Thus, when considering the derivative over that variable, a fractional derivative is obtained.

Fuentes et al. [24, 25] proposed a model that represents anomalous diffusion of petroleum (very fast or very slow) in three types of medium with different porosities: fractures, vugs and matrix. This model has a dimensionless form, which is seen as

$$\left(1 - \omega\_{\boldsymbol{f}} - \omega\_{\boldsymbol{v}}\right) \frac{\partial^{a\_1} p\_{\boldsymbol{m}}}{\partial t^{a\_1}} = (1 - \kappa\_{\boldsymbol{f}} - \kappa\_{\boldsymbol{v}}) \frac{1}{r} \frac{\partial}{\partial r} \left(r \frac{\partial^{\boldsymbol{\beta}\_1} p\_{\boldsymbol{m}}}{\partial r^{\boldsymbol{\beta}\_1}}\right) + \lambda\_{\boldsymbol{m} \boldsymbol{f}} (p\_{\boldsymbol{f}} - p\_{\boldsymbol{m}}) + \lambda\_{\boldsymbol{m} \boldsymbol{v}} (p\_{\boldsymbol{v}} - p\_{\boldsymbol{m}}), \qquad (1)$$

$$
\lambda \omega\_{f} \frac{\partial^{a\_{2}} p\_{f}}{\partial t^{a\_{2}}} = \kappa\_{f} \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial^{\theta\_{2}} p\_{m}}{\partial r^{\theta\_{2}}} \right) - \lambda\_{mf} (p\_{f} - p\_{m}) + \lambda\_{f\upsilon} (p\_{\upsilon} - p\_{f}), \tag{2}
$$

$$
\omega\_v \frac{\partial^{\alpha\_3} p\_v}{\partial t^{\alpha\_3}} = \kappa\_v \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial^{\beta\_3} p\_m}{\partial r^{\beta\_3}} \right) - \lambda\_{mv} (p\_v - p\_m) + \lambda\_{fv} (p\_v - p\_f). \tag{3}
$$

This model must be solved for p's (pressure) at time t in the matrix of the medium pm, the fractured medium pf, and the vugular medium pv .

κ represents each permeability tensor, which is assumed to be constant; r is the dimensionless variable related to the ratios of the wellbore. Values λ are dimensionless values and they are defined in terms of transfer coefficients at each interface, the radius of the wellbore, the dynamic viscosity of the fluid and the permeability tensor.

If time orders α's are less than one we say that the process is subdiffusive and if it is greater than one it is called superdiffusive. Spatial orders β's result from generalizing the classic Darcy law, which corresponds to β = 1.

The above equations, like the flow equation that considers the medium as a whole, are modified Bessel's equations.

Different porosities imply different fractal dimensions and therefore, fractional derivatives with different orders.

In Fractional Calculus there is not a single definition of derivative and between these definitions, in general, not equivalent to each other. Definitions such as Riemann-Liouville, Riesz, Caputo, Weyl or Grünwald-Letnikov are the most used to model anomalous diffusion, from a porous medium with fractal dimension.

But, fractional equations present serious numerical and mathematical difficulties in the context of diffusion equations. From the computational point of view, the challenge to overcome is to attenuate the numerical cost of the matrices that result in discretizing the problem. The proposal is to consider radial basis interpolators [26], taking into account that the domain geometry does not determine the efficiency of the algorithm and sounds like an immediate alternative for generalizing larger dimensions.

## 2. Introduction

The objective in this work is to deal with fractional differential operators, using radial basis functions (RBF) and optimizing discretization processes of such fractional operators, through QR matrix decomposition and to attenuate the bad condition due to the shape parameter. With RBF-QR, a very high precision and convergence are obtained without the need to increase a polynomial term to the interpolator [13]. As the computational cost increases, so does the range of problems to which such a technique can be applied.

## 3. Radial basis function method for interpolation

Definition 1. Let s be a positive integer, a function <sup>ϕ</sup> : <sup>R</sup><sup>s</sup> ! <sup>R</sup> is called radial whenever there is a function in one variable φ : ½0, ∞Þ ! R such that

$$
\phi(\mathbf{x}) = \phi(r), \quad \text{where } r = ||\mathbf{x}|| \tag{4}
$$

and jj � jj is a norm in <sup>R</sup><sup>s</sup> (usually the Euclidean standard norm). Table (1) shows some examples, about the more used real functions φðrÞ.

A standard interpolator in terms of radial basis functions, given the data uk, of some real function u, in the corresponding collocation nodes xk, k ¼ 1, …, N, has the form [13]

$$s\_{\varepsilon}(\mathbf{x}) = \sum\_{k=1}^{N} \lambda\_k \phi(\varepsilon || \mathbf{x} - \mathbf{x\_k}||) \tag{5}$$

where φ is one real variable function and the constant value ε is called shape parameter.


Table 1. Common elections for φðrÞ.

Note 1. The radial basis function (FBR) we will take for the approach and examples is the Gaussian function:

$$
\phi(r) = e^{-r^2} \quad r \ge 0. \tag{6}
$$

Because it is globally smooth and for which the RBF method can be applied together with QR matrix decomposition [26].

Notation 1. Considering one dimensional case, we abbreviate φkðxÞ :¼ φðεjjx � xkjjÞ. Now the interpolator (5) is written as

$$s\_\ell(\mathbf{x}) \equiv \sum\_{k=1}^N \lambda\_k \phi\_k(\mathbf{x}). \tag{7}$$

The unknown coefficients λ<sup>k</sup> can be determined by interpolation conditions

$$s\_{\varepsilon}(\mathfrak{x}\_{i}) = \mathfrak{u}\_{i\prime} \quad i = 1, \ldots, N.$$

producing the linear equations system

$$\begin{array}{ccccc} \lambda\_1 \phi\_1(\mathbf{x}\_1) + \dots + \lambda\_N \phi\_N(\mathbf{x}\_1) & = & \mu\_1 \\ \vdots & & \vdots & \vdots \\ \lambda\_1 \phi\_1(\mathbf{x}\_N) + \dots + \lambda\_N \phi\_N(\mathbf{x}\_N) & = & \mu\_N \end{array}$$

which can be described in matrix form

$$
\underbrace{\begin{bmatrix}\phi\_1(\mathbf{x}\_1) & \cdots & \phi\_N(\mathbf{x}\_1) \\ \vdots & \vdots & \vdots \\ \phi\_1(\mathbf{x}\_N) & \cdots & \phi\_N(\mathbf{x}\_N)\end{bmatrix}}\_A \underbrace{\begin{bmatrix}\lambda\_1 \\ \vdots \\ \lambda\_N\end{bmatrix}}\_A = \underbrace{\begin{bmatrix}u\_1 \\ \vdots \\ \boldsymbol{u}\_N\end{bmatrix}}\_{u\_N}.
$$

Thus

2. Introduction

The objective in this work is to deal with fractional differential operators, using radial basis functions (RBF) and optimizing discretization processes of such fractional operators, through QR matrix decomposition and to attenuate the bad condition due to the shape parameter. With RBF-QR, a very high precision and convergence are obtained without the need to increase a polynomial term to the interpolator [13]. As the computational cost increases, so does the

Definition 1. Let s be a positive integer, a function <sup>ϕ</sup> : <sup>R</sup><sup>s</sup> ! <sup>R</sup> is called radial whenever there is a

and jj � jj is a norm in <sup>R</sup><sup>s</sup> (usually the Euclidean standard norm). Table (1) shows some examples, about

A standard interpolator in terms of radial basis functions, given the data uk, of some real

function u, in the corresponding collocation nodes xk, k ¼ 1, …, N, has the form [13]

k¼1

where φ is one real variable function and the constant value ε is called shape parameter.

<sup>s</sup>εðxÞ ¼ <sup>X</sup> N

RBF name RBF φðrÞ

Compact support ('Wendland') ð1 � εrÞ

Gaussian (GA) e�ðεr<sup>Þ</sup>

Inverse multiquadric (IMQ) <sup>1</sup><sup>=</sup>

Multiquadric (MQ) ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

Inverse quadric (IQ) <sup>1</sup>=ð<sup>1</sup> þ ðεr<sup>Þ</sup>

Bessel (BE) (<sup>d</sup> <sup>¼</sup> <sup>1</sup>; <sup>2</sup>;…) Jd=2�<sup>1</sup>ðεrÞ=ðεr<sup>Þ</sup>

Polyharmonic spline rm, <sup>m</sup> <sup>¼</sup> <sup>1</sup>; <sup>3</sup>; <sup>5</sup>;…

ϕðxÞ ¼ φðrÞ, where r ¼ jjxjj ð4Þ

λkφðεjjx � xkjjÞ ð5Þ

rmlnðrÞ, <sup>m</sup> <sup>¼</sup> <sup>2</sup>; <sup>4</sup>; <sup>6</sup>;…

2 Þ

d=2�1

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ ðεrÞ 2

<sup>þ</sup>pðεrÞ, <sup>p</sup> polynomial.

m

2

q

q

1 þ ðεrÞ 2

range of problems to which such a technique can be applied.

6 Fractal Analysis - Applications in Physics, Engineering and Technology

3. Radial basis function method for interpolation

function in one variable φ : ½0, ∞Þ ! R such that

the more used real functions φðrÞ.

Piecewiese smooth, global

Table 1. Common elections for φðrÞ.

Smooth, global

$$A\boldsymbol{\lambda} = \boldsymbol{u}\_{\boldsymbol{\lambda}^\*} \quad \text{where } \boldsymbol{A} = \begin{bmatrix} \boldsymbol{\Phi}(\mathbf{x}\_1)^T \\ \vdots \\ \boldsymbol{\Phi}(\mathbf{x}\_N)^T \end{bmatrix} \quad \text{and} \quad \boldsymbol{\Phi}(\mathbf{x}) := \begin{bmatrix} \boldsymbol{\phi}\_1(\mathbf{x}) \\ \vdots \\ \boldsymbol{\phi}\_N(\mathbf{x}) \end{bmatrix}. \tag{8}$$

A is called Gram matrix.

#### 4. RBF methodology for the discretization of differential operators

Suppose we want to apply a differential operator L to the function uðxÞ in every evaluation node from the set Z ¼ {z1, …, zNe }, given the values of the function at collocation nodes X ¼ {x1, …, xN}

$$
\mathcal{L}u\_Z \approx \mathcal{D}u\_X. \tag{9}
$$

$$\text{where } \mathcal{L}u\_{\mathcal{Z}} = \begin{bmatrix} \mathcal{L}u(z\_1) \\ \vdots \\ \mathcal{L}u(z\_{N\_\varepsilon}) \end{bmatrix} \text{ and } u\_{\mathcal{X}} = \begin{bmatrix} u(\mathbf{x}\_1) \\ \vdots \\ u(\mathbf{x}\_N) \end{bmatrix}.$$

Assuming that L is linear, we apply it to the RBF interpolator, for each of the evaluation nodes, as an approximation to the values of L in u:

$$\begin{aligned} \mathcal{L}s\_{\boldsymbol{\varepsilon}}(\boldsymbol{z}\_{i}) &= \sum\_{k=1}^{N} \lambda\_{k} \mathcal{L}\phi\_{k}(\boldsymbol{z}\_{i}), \quad \boldsymbol{i} = 1, 2, \ldots, N\_{\boldsymbol{\varepsilon}};\\ \mathcal{L}(\boldsymbol{s}\_{\boldsymbol{\varepsilon}}(\boldsymbol{z}\_{1})) &= \quad \lambda\_{1} \mathcal{L}\phi\_{1}(\boldsymbol{z}\_{1}) + \cdots + \lambda\_{N} \mathcal{L}\phi\_{N}(\boldsymbol{z}\_{1})\\ \vdots & \quad \vdots\\ \mathcal{L}(\boldsymbol{s}\_{\boldsymbol{\varepsilon}}(\boldsymbol{z}\_{N\_{\boldsymbol{\varepsilon}}})) &= \quad \lambda\_{1} \mathcal{L}\phi\_{1}(\boldsymbol{z}\_{N\_{\boldsymbol{\varepsilon}}}) + \cdots + \lambda\_{N} \mathcal{L}\phi\_{N}(\boldsymbol{z}\_{N\_{\boldsymbol{\varepsilon}}}) \end{aligned} \tag{10}$$

which in the matrix form is

$$\mathcal{L}\mathbf{s}\_{\ell} = B\lambda \tag{11}$$

where

$$
\mathcal{L}s\_{\varepsilon} = \begin{bmatrix}
\mathcal{L}(\mathbf{s}\_{\varepsilon}(\mathbf{z}\_{1})) \\
\vdots \\
\mathcal{L}(\mathbf{s}\_{\varepsilon}(\mathbf{z}\_{N\_{\varepsilon}}))
\end{bmatrix},
\mathcal{B} = \begin{bmatrix}
\mathcal{L}\phi\_{1}(\mathbf{z}\_{1}) & \cdots & \mathcal{L}\phi\_{N}(\mathbf{z}\_{1}) \\
\vdots & \vdots & \vdots \\
\mathcal{L}\phi\_{1}(\mathbf{z}\_{N\_{\varepsilon}}) & \cdots & \mathcal{L}\phi\_{N}(\mathbf{z}\_{N\_{\varepsilon}})
\end{bmatrix}.
$$

The coefficient matrix λ ¼ λ1 ⋮ λ<sup>N</sup> 2 4 3 5 is calculated by solving for λ in Eq. (8). In this way, from

Eq. (11) we obtain

$$\mathcal{L}\mathbf{s}\_{\varepsilon} = (BA^{-1})\mu\_X. \tag{12}$$

Thus, the matrix BA�<sup>1</sup> is that matrix D which is mentioned in Eq. (9) and that approximates the values of the differential operator applied to the function.

#### 5. Fractional calculus definitions

The definitions shown in this section are based on the texts of Podlubny [27] and Samko [28]. Generalizing the formula of Cauchy, we obtain

$$\mathcal{J}^n\left(f(t)\right) := \int \cdots \int\_0^t f(\tau)d\tau = \frac{1}{(n-1)!} \int\_0^t (t-\tau)^{n-1} f(\tau)d\tau$$

Definition 2. The left-sided Riemann-Liouville fractional integral of order α of function f(x) is defined as

Applications of Radial Basis Function Schemes to Fractional Partial Differential Equations http://dx.doi.org/10.5772/67892 9

$$\mathbf{r}\_a^{RL} I\_x^{\alpha} f(\mathbf{x}) = \frac{1}{\Gamma(\alpha)} \int\_a^{\infty} (\mathbf{x} - \tau)^{\alpha - 1} f(\tau) d\tau,\ \mathbf{x} > a. \tag{13}$$

Definition 3. The right-sided Riemann-Liouville fractional integral of order α of function f(x) is defined as

$$\, \_x^{RL}I\_b^\alpha f(\mathbf{x}) = \frac{1}{\Gamma(\alpha)} \int\_{\mathbf{x}}^b (\pi - \mathbf{x})^{\alpha - 1} f(\pi) d\tau,\ \mathbf{x} < b. \tag{14}$$

Based on the definition of fractional integral, the following fractional derivatives are constructed, which we will use in the numerical examples

Definition 4. The left-sided Riemann-Liouville fractional derivative of order α of function f(x) is defined as

$$\mathbf{P}\_a^{RL} D\_x^a f(\mathbf{x}) = \frac{1}{\Gamma(m-a)} \frac{d^m}{d\mathbf{x}^m} \int\_a^\mathbf{x} (\mathbf{x} - \tau)^{m-a-1} f(\tau) d\tau,\ \mathbf{x} > a,\tag{15}$$

where m ¼ ⌈α⌉.

LuZ ≈ DuX: ð9Þ

Ls<sup>ε</sup> ¼ Bλ ð11Þ

3 5:

ÞuX: ð12Þ

Lφ1ðz1Þ ⋯ LφNðz1Þ ⋮⋮⋮ Lφ1ðzNe Þ ⋯ LφNðzNe Þ

5 is calculated by solving for λ in Eq. (8). In this way, from

ð10Þ

where LuZ ¼

where

Luðz1Þ ⋮ LuðzNe Þ

as an approximation to the values of L in u:

Ls<sup>ε</sup> ¼

2 4

3

8 Fractal Analysis - Applications in Physics, Engineering and Technology

5 and uX ¼

<sup>L</sup>sεðziÞ ¼ <sup>X</sup>

N

k¼1

Lðsεðz1ÞÞ ⋮ LðsεðzNe ÞÞ

3

λ1 ⋮ λ<sup>N</sup>

2 4

values of the differential operator applied to the function.

5. Fractional calculus definitions

J <sup>n</sup> � fðtÞ � :¼ Z ⋯ Z <sup>t</sup> 0

Generalizing the formula of Cauchy, we obtain

3 5, B ¼

uðx1Þ ⋮ uðxNÞ 3 5.

Assuming that L is linear, we apply it to the RBF interpolator, for each of the evaluation nodes,

Lðsεðz1ÞÞ ¼ λ1Lφ1ðz1Þ þ ⋯ þ λNLφNðz1Þ

LðsεðzNe ÞÞ ¼ λ1Lφ1ðzNe Þ þ ⋯ þ λNLφNðzNe Þ

2 4

<sup>L</sup>s<sup>ε</sup> ¼ ðBA�<sup>1</sup>

Thus, the matrix BA�<sup>1</sup> is that matrix D which is mentioned in Eq. (9) and that approximates the

The definitions shown in this section are based on the texts of Podlubny [27] and Samko [28].

<sup>f</sup>ðτÞd<sup>τ</sup> <sup>¼</sup> <sup>1</sup>

Definition 2. The left-sided Riemann-Liouville fractional integral of order α of function f(x) is

ðn � 1Þ!

Z <sup>t</sup> 0

ðt � τÞ n�1 fðτÞdτ

⋮⋮ ⋮

λkLφkðziÞ, i ¼ 1; 2, …, Ne;

2 4

2 4

which in the matrix form is

The coefficient matrix λ ¼

Eq. (11) we obtain

defined as

Definition 5. The right-sided Riemann-Liouville fractional derivative of order α of function f(x) is defined as

$$\mathbf{u}^{RL}\_{\mathbf{x}} D\_{b}^{\alpha} f(\mathbf{x}) = \frac{(-1)^{m}}{\Gamma(m-a)} \frac{d^{m}}{d\mathbf{x}^{m}} \int\_{\mathbf{x}}^{b} (\boldsymbol{\tau} - \mathbf{x})^{m-a-1} f(\boldsymbol{\tau}) d\boldsymbol{\tau}, \ \mathbf{x} < b,\tag{16}$$

where m ¼ ⌈α⌉.

Definition 6. [28–30] The Riesz fractional operator for α on a finite interval 0 ≤ x ≤ L is defined as

$$\frac{\partial^{\alpha}}{\partial |\mathbf{x}|^{\alpha}} f(\mathbf{x}, t) = -c\_{\alpha} \left( {}^{RL}\_{a}D\_{x}^{\alpha} + {}^{RL}\_{x}D\_{b}^{\alpha} \right) f(\mathbf{x}, t), \tag{17}$$

where

$$\mathcal{L}\_{\alpha} = \frac{1}{2\cos\left(\frac{\alpha\alpha}{2}\right)}, \ \alpha \neq 1,\tag{18}$$

$${}^{RL}\_{a}D\_{x}^{a}f(\mathbf{x},t) = \frac{1}{\Gamma(m-\alpha)}\frac{d^{m}}{d\mathbf{x}^{m}}\int\_{a}^{\mathbf{x}}(\mathbf{x}-\tau)^{m-\alpha-1}f(\tau,t)d\tau,\tag{19}$$

$${}^{RL}\_{x}D\_{b}^{\alpha}f(\mathbf{x},t) = \frac{(-1)^{m}}{\Gamma(m-\alpha)}\frac{d^{m}}{d\mathbf{x}^{m}}\int\_{\mathbf{x}}^{b} (\tau-\mathbf{x})^{m-\alpha-1}f(\tau,t)d\tau,\tag{20}$$

where m ¼ ⌈α⌉.

Definition 7. The left-sided Caputo fractional derivative of order α of function f(x) is defined as

$$\, \, \_a^c D\_x^a f(\mathbf{x}) = \frac{1}{\Gamma(m-a)} \int\_a^\mathbf{x} (\mathbf{x} - \tau)^{m-a-1} f^{(m)}(\tau) d\tau,\ \mathbf{x} > a,\tag{21}$$

and m ¼ ⌈α⌉.

The fact that the derivative of integer order appears within the integral in Definition 7, makes Caputo derivative the most suitable for dealing with initial conditions of the FPDEs.

Definition 8. The right-sided Caputo fractional derivative of order α of function f(x) is defined as

$$\, \, \_x^{\mathbb{C}}D\_b^a f(\mathbf{x}) = \frac{(-1)^m}{\Gamma(m-a)} \int\_{\mathbf{x}}^b (\tau - \mathbf{x})^{m-a-1} f^{(m)}(\tau) d\tau, \; \mathbf{x} < b,\tag{22}$$

and m ¼ ⌈α⌉.

#### 6. Numerical examples

#### 6.1. Fractional partial differential equation with Riesz space fractional derivatives

The equation that we take in this part is the advection-diffusion equation (see [30, 6]).

$$\begin{cases} \frac{\partial u(\mathbf{x},t)}{\partial t} = -K\_a \frac{\partial^a u(\mathbf{x},t)}{\partial |\mathbf{x}|^a}, & \mathbf{x} \in [0, \pi], \ t \in (0, T],\\ u(\mathbf{x},0) = u\_0(\mathbf{x}), \\ u(0,t) = u(\pi, t) = 0, \end{cases} \tag{23}$$

where u can be, for example, the concentration of a dissolute substance, K<sup>α</sup> the dispersion coefficient and the fractional derivative of Riesz is given with fractional order 1 < α ≤ 2.

Taking advantage of the fact that for RBFs we can consider non-equispaced collocation nodes, we divide the interval ½0, π� into N nodes, using Chebyshev distribution (see Figure 1):

$$\text{If } \mathbf{x}\_i = \frac{\pi}{2} \cos \left(\theta\_i\right) + \frac{\pi}{2}, \text{ } \theta\_i = \pi - i \frac{\pi}{N - 1}, \text{ } i = 0, 1, \dots, N - 1.$$

The FPDE is solved by the method of lines based on the spatial trial spaces spanned by the Lagrange basis associated to RBFs. The Lagrange basis L1ðxÞ, …, LNðxÞ is generated by radial functions φ<sup>j</sup> ðxÞ ¼ φðjx � xjjÞ, j ¼ 1; 2,…, N, taking into account the collocation nodes. This is done by solving the system

Figure 1: *Chebyshev nodes distribution over* [0, π] *interval.*

Figure 1. Chebyshev nodes distribution over ½0; π� interval.

Applications of Radial Basis Function Schemes to Fractional Partial Differential Equations http://dx.doi.org/10.5772/67892 11

$$L(\mathbf{x})^T = \Phi(\mathbf{x})^T \mathbf{A}^{-1} \tag{24}$$

where

C aD<sup>α</sup>

10 Fractal Analysis - Applications in Physics, Engineering and Technology

C xD<sup>α</sup>

and m ¼ ⌈α⌉.

and m ¼ ⌈α⌉.

functions φ<sup>j</sup>

done by solving the system

6. Numerical examples

<sup>x</sup> <sup>f</sup>ðxÞ ¼ <sup>1</sup>

<sup>b</sup> <sup>f</sup>ðxÞ ¼ ð�1<sup>Þ</sup>

∂uðx, tÞ

xi <sup>¼</sup> <sup>π</sup>

Figure 1. Chebyshev nodes distribution over ½0; π� interval.

<sup>∂</sup><sup>t</sup> ¼ �K<sup>α</sup>

uðx; 0Þ ¼ u0ðxÞ, uð0, tÞ ¼ uðπ, tÞ ¼ 0,

<sup>2</sup> cos <sup>ð</sup>θiÞ þ <sup>π</sup>

Γðm � αÞ

Z <sup>x</sup> a

Caputo derivative the most suitable for dealing with initial conditions of the FPDEs.

Z <sup>b</sup> x

6.1. Fractional partial differential equation with Riesz space fractional derivatives The equation that we take in this part is the advection-diffusion equation (see [30, 6]).

> <sup>∂</sup><sup>α</sup>uðx, t<sup>Þ</sup> ∂jxj

where u can be, for example, the concentration of a dissolute substance, K<sup>α</sup> the dispersion coefficient and the fractional derivative of Riesz is given with fractional order 1 < α ≤ 2.

Taking advantage of the fact that for RBFs we can consider non-equispaced collocation nodes,

The FPDE is solved by the method of lines based on the spatial trial spaces spanned by the Lagrange basis associated to RBFs. The Lagrange basis L1ðxÞ, …, LNðxÞ is generated by radial

Figure 1: *Chebyshev nodes distribution over* [0, π] *interval.*

N � 1

ðxÞ ¼ φðjx � xjjÞ, j ¼ 1; 2,…, N, taking into account the collocation nodes. This is

, i ¼ 0; 1, …, N � 1:

we divide the interval ½0, π� into N nodes, using Chebyshev distribution (see Figure 1):

<sup>2</sup> , <sup>θ</sup><sup>i</sup> <sup>¼</sup> <sup>π</sup> � <sup>i</sup> <sup>π</sup>

m

Γðm � αÞ

ðx � τÞ

The fact that the derivative of integer order appears within the integral in Definition 7, makes

Definition 8. The right-sided Caputo fractional derivative of order α of function f(x) is defined as

ðτ � xÞ

m�α�1 f ðmÞ

m�α�1 f ðmÞ

<sup>α</sup> , x ∈½0, π�, t∈ ð0, T�,

ðτÞdτ, x > a, ð21Þ

ðτÞdτ, x < b, ð22Þ

ð23Þ

$$L(\mathbf{x})^T = \begin{bmatrix} L\_1(\mathbf{x}) & \cdots & L\_N(\mathbf{x}) \end{bmatrix} \ \Phi(\mathbf{x})^T = \begin{bmatrix} \phi\_1(\mathbf{x}) & \cdots & \phi\_N(\mathbf{x}) \end{bmatrix} \tag{25}$$

and the Gram matrix

$$A = \begin{bmatrix} \phi\_1(\mathbf{x}\_1) & \cdots & \phi\_N(\mathbf{x}\_1) \\ \vdots & \vdots & \vdots \\ \phi\_1(\mathbf{x}\_N) & \cdots & \phi\_N(\mathbf{x}\_N) \end{bmatrix}.$$

If L is a differential operator and RBF φ is sufficiently smooth, then the application of such operator to the Lagrange base is calculated through the relation

$$(\mathcal{L}L)(\mathfrak{x}) = (\mathcal{L}\phi)A^{-1}.$$

Because of Lagrange's standard conditions, the zero boundary conditions in x<sup>1</sup> ¼ 0 and xN ¼ π are de facto fulfilled if we use an approximation generated by the functions L2, …, LN�1. This approximation is then represented as

$$u(\mathbf{x}, t) = \sum\_{j=2}^{N-1} \beta\_j(t) L\_j(\mathbf{x})\_j$$

with unknown vector

$$\beta(t) = \begin{bmatrix} \beta\_2(t) \\ \vdots \\ \beta\_{N-1}(t) \end{bmatrix}.$$

Evaluating the interpolator at the PDE in each node xi, we obtain

$$\sum\_{j=2}^{N-1} \beta\_j'(t) L\_j(\mathbf{x}\_i) = -K\_\alpha \sum\_{j=2}^{N-1} \beta\_j(t) \frac{\partial^\alpha}{\partial |\mathbf{x}|^\alpha} L\_j(\mathbf{x}\_i)$$

and initial conditions

$$\beta\_j(0) = \mu\_0(\mathfrak{x}\_j), \ 2 \le j \le N - 1.$$

From the latter two equations we obtain the following system of ordinary differential equations

$$\beta'(t) = -K\_{\alpha} \left( \frac{\partial^{\alpha}}{\partial |\mathbf{x}|^{\alpha}} \mathbf{L} \right) \cdot \beta(t), \ \beta(0) = \mathcal{U}\_{0\prime}$$

where

$$\frac{\partial^{\alpha}}{\partial |\mathbf{x}|^{\alpha^{\prime}}} \mathbf{L} = \begin{bmatrix} \frac{\partial^{\alpha}}{\partial |\mathbf{x}|^{\alpha^{\prime}}} L\_{2}(\mathbf{x}\_{2}) & \cdots & \frac{\partial^{\alpha}}{\partial |\mathbf{x}|^{\alpha^{\prime}}} L\_{N-1}(\mathbf{x}\_{2})\\ \vdots & \cdots & \vdots\\ \frac{\partial^{\alpha}}{\partial |\mathbf{x}|^{\alpha^{\prime}}} L\_{2}(\mathbf{x}\_{N}) & \cdots & \frac{\partial^{\alpha}}{\partial |\mathbf{x}|^{\alpha^{\prime}}} L\_{N-1}(\mathbf{x}\_{N}) \end{bmatrix} \text{and } \mathcal{U}\_{0} = \begin{bmatrix} u\_{0}(\mathbf{x}\_{2})\\ \vdots\\ u\_{0}(\mathbf{x}\_{N-1}) \end{bmatrix}. \tag{26}$$

According to the equation data in articles [30, 6], the problem (23) is taken with parameters <sup>α</sup> <sup>¼</sup> <sup>1</sup>:8, <sup>K</sup><sup>α</sup> <sup>¼</sup> <sup>0</sup>:25 and <sup>u</sup>0ðxÞ ¼ <sup>x</sup><sup>2</sup>ð<sup>π</sup> � <sup>x</sup>Þ. For the numerical solution, Gaussian RBF is used with a shape parameter ε ¼ 0:8, applying the scheme shown in Section 6.1. Graphs are shown in Figure 2.

Figure 2. RBF approximation to solution, for Eq. (23), with α ¼ 1:8, ε ¼ 0:8 and K<sup>α</sup> ¼ 0:25.

#### 6.2. Riemann-Liouville space-fractional diffusion equation

In this part we consider the problem introduced by Sousa [31] and taken up in Ref. [26] by Cécile and Hanert, the fractional diffusion equation in one dimension

$$\frac{\partial f(\mathbf{x},t)}{\partial t} = d(\mathbf{x})^{RL}\_{\ 0} D\_{x}^{\alpha} f(\mathbf{x},t) + q(\mathbf{x},t), \ \mathbf{x} \in [0,1] \text{ and } t > 0 \tag{27}$$

with RL 0D<sup>α</sup> <sup>x</sup> Riemannn-Liouville derivative, 1 < α ≤ 2,

$$d(\mathbf{x}) = \frac{\Gamma(5-\alpha)}{24} \mathbf{x}^{\alpha} \ y \ q(\mathbf{x}, t) = -2e^{-t} \mathbf{x}^{4}.$$

Initial conditions are considered

$$f(\mathbf{x}, \mathbf{0}) = \mathbf{x}^4, \ \mathbf{x} \in (a, b) \tag{28}$$

and Dirichlet boundary conditions

$$f(0,t) = 0, \; f(1,t) = e^{-t}.\tag{29}$$

The exact solution to Eq. (27) is

$$f(\mathbf{x}, t) = e^{-t} \mathbf{x}^4.$$

We apply RBF and the method of lines with Lagrange polynomials, in a similar way as it was applied in Section 6.1. With N ¼ 21 Chebyshev collocation nodes and Ne ¼ 100 evaluation nodes, various values for the order of the derivative α and for shape parameters.

The results obtained by Sousa [31] were obtained by applying implicit Crank-Nicholson schemes for time. The discretization of the fractional derivative was done using splines, with second-order precision. The results of both Sousa and RBF are shown in Tables 2 and 3, taking into account the error

$$||\boldsymbol{\mu}\_{\text{exact}} - \boldsymbol{\mu}\_{\text{approx}}\boldsymbol{\mu}\_{\text{ $\boldsymbol{\sigma}$ }}||\boldsymbol{\epsilon}\_{\text{\textquotedblleft}}\tag{30}$$

where jj � jj<sup>∞</sup> is the <sup>ℓ</sup><sup>∞</sup> norm.

∂α ∂jxj <sup>α</sup> L ¼

in Figure 2.

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

**u**

with RL 0D<sup>α</sup>

∂α ∂jxj

12 Fractal Analysis - Applications in Physics, Engineering and Technology

∂α ∂jxj

Approximation with alpha = 1.8 Ka = 0.25 and t = 0.5

6.2. Riemann-Liouville space-fractional diffusion equation

<sup>0</sup> 0.5 <sup>1</sup> 1.5 <sup>2</sup> 2.5 <sup>3</sup> 3.5 <sup>0</sup>

**x**

<sup>∂</sup><sup>t</sup> <sup>¼</sup> <sup>d</sup>ðx<sup>Þ</sup>

<sup>x</sup> Riemannn-Liouville derivative, 1 < α ≤ 2,

∂fðx, tÞ

Initial conditions are considered

and Dirichlet boundary conditions

The exact solution to Eq. (27) is

Cécile and Hanert, the fractional diffusion equation in one dimension

Figure 2. RBF approximation to solution, for Eq. (23), with α ¼ 1:8, ε ¼ 0:8 and K<sup>α</sup> ¼ 0:25.

RL <sup>0</sup> D<sup>α</sup>

<sup>d</sup>ðxÞ ¼ <sup>Γ</sup>ð<sup>5</sup> � <sup>α</sup><sup>Þ</sup>

<sup>f</sup>ðx; <sup>0</sup>Þ ¼ <sup>x</sup><sup>4</sup>

fð0, tÞ ¼ 0, fð1, tÞ ¼ e

<sup>α</sup> <sup>L</sup>2ðx2<sup>Þ</sup> <sup>⋯</sup> <sup>∂</sup><sup>α</sup>

<sup>α</sup> <sup>L</sup>2ðxN<sup>Þ</sup> <sup>⋯</sup> <sup>∂</sup><sup>α</sup>

⋮⋯ ⋮

∂jxj

∂jxj

According to the equation data in articles [30, 6], the problem (23) is taken with parameters <sup>α</sup> <sup>¼</sup> <sup>1</sup>:8, <sup>K</sup><sup>α</sup> <sup>¼</sup> <sup>0</sup>:25 and <sup>u</sup>0ðxÞ ¼ <sup>x</sup><sup>2</sup>ð<sup>π</sup> � <sup>x</sup>Þ. For the numerical solution, Gaussian RBF is used with a shape parameter ε ¼ 0:8, applying the scheme shown in Section 6.1. Graphs are shown

In this part we consider the problem introduced by Sousa [31] and taken up in Ref. [26] by

<sup>24</sup> <sup>x</sup><sup>α</sup> y qðx, t޼�2<sup>e</sup>

<sup>α</sup> LN�<sup>1</sup>ðx2Þ

and U<sup>0</sup> ¼

<sup>x</sup> fðx, tÞ þ qðx, tÞ, x∈ ½0; 1� and t > 0 ð27Þ

**t**

, x∈ða, bÞ ð28Þ

: ð29Þ

�t x4 :

�t

u0ðx2Þ ⋮ u0ðxN�<sup>1</sup>Þ 3

<sup>0</sup> 0.5 <sup>1</sup> 1.5 <sup>2</sup> 2.5 <sup>3</sup> 3.5 <sup>4</sup>

**x**

5: ð26Þ

2 4

RBF, Lagrange polynomial, 16 collocation nodes and 100 evaluation nodes

<sup>α</sup> LN�<sup>1</sup>ðxNÞ

**u**



1/30 0:2067 · 10�<sup>3</sup> 1.9 0:1150 · 10�<sup>3</sup> 2.0

Table 2. Comparison of results for FPDE (27): Sousa results [31].


Table 3. Comparison of results for FPDE (27): RBF results.

#### 6.3. Caputo time fractional partial differential equations

In this part we consider examples of fractional partial differential equations (see [32]) of the type

$$\frac{\partial^{\alpha}u(\mathbf{x},t)}{\partial t^{\alpha}} + \delta \frac{\partial u(\mathbf{x},t)}{\partial \mathbf{x}} + \gamma \frac{\partial^{2}u(\mathbf{x},t)}{\partial \mathbf{x}^{2}} = f(\mathbf{x},t), \tag{31}$$

where t > 0, x∈ ½a, b�, 0 < α ≤ 1, δ and γ are real parameters, bounded initial condition u(x,0) ¼ u0ðxÞ and Dirichlet boundary conditions uða, tÞ ¼ g1ðtÞ and uðb, tÞ ¼ g2ðtÞ, t ≥ 0. Caputo Derivative will be considered and rbf-qr routines from www.it.uu.se/research/scientific\_computing/ software/rbf\_qr

#### 6.3.1. Example 1

Putting <sup>δ</sup> <sup>¼</sup> <sup>1</sup>, <sup>γ</sup> ¼ �1 and <sup>f</sup>ðx, tÞ ¼ <sup>2</sup><sup>t</sup> 2�α <sup>Γ</sup>ð3�α<sup>Þ</sup> <sup>þ</sup> <sup>2</sup><sup>x</sup> � 2 in Eq. (31), we obtain a non-homogeneous, fractional and linear Burger equation

$$\frac{\partial^{\alpha}u(\mathbf{x},t)}{\partial t^{\alpha}} + \frac{\partial u(\mathbf{x},t)}{\partial \mathbf{x}} - \frac{\partial^{2}u(\mathbf{x},t)}{\partial \mathbf{x}^{2}} = \frac{2t^{2-\alpha}}{\Gamma(3-\alpha)} + 2\mathbf{x} - \mathbf{2}.\tag{32}$$

Using initial condition

$$
\mu(\mathbf{x},0) = \mathbf{x}^2,\tag{33}
$$

Dirichlet boundary conditions

$$
\mu(0, t) = t^2, \ u(1, t) = 1 + t^2. \tag{34}
$$

The exact solution (see [33]) is

$$
\mu(\mathbf{x}, t) = \mathbf{x}^2 + t^2,\tag{35}
$$

This problem is solved using the method described in section 6.1 along with QR decomposition for spatial part and an implicit scheme for time, on domain ½0; 1�, fractional derivative of order α ¼ 0:5; several shape parameters are considered and the one that gives the best approximation is taken. Figure 3 shows a 3D image of the solution and maximal error for the solution at time t ¼ 1.

In Figure 4 we compare errors that result when choosing uniform collocation nodes (with a fixed step size) and Chebyshev, on definition interval ½0; 1�. The election of Chebyshev nodes is due to attenuating instability, which manifests as oscillations in the graph, called the Gibbs phenomenon [34, 35].

[41 uniform collocation nodes and 100 evaluation nodes, maximal error for t ¼ 0:5, α ¼ 0:5, <sup>ε</sup> <sup>¼</sup> <sup>3</sup>:6 is 0:<sup>83595</sup> · <sup>10</sup>�<sup>4</sup> .] [41 Chebyshev collocation nodes and 100 evaluation nodes, maximal error for <sup>t</sup> <sup>¼</sup> <sup>0</sup>:5, <sup>α</sup> <sup>¼</sup> <sup>0</sup>:5, <sup>ε</sup> <sup>¼</sup> <sup>3</sup>:6 is 0:<sup>14385</sup> · <sup>10</sup>�<sup>4</sup> .]

Applications of Radial Basis Function Schemes to Fractional Partial Differential Equations http://dx.doi.org/10.5772/67892 15

Figure 3. Numerical solution for Eq. (32) with α ¼ 0:5.

Figure 4. Comparison between results for Eq. (32) with α ¼ 0:5 and t ¼ 0:5, due to election of nodes. (a) 41 uniform collocation nodes and 100 evaluation nodes, maximal error for <sup>t</sup> <sup>¼</sup> <sup>0</sup>:5, <sup>α</sup> <sup>¼</sup> <sup>0</sup>:5, <sup>ε</sup> <sup>¼</sup> <sup>3</sup>:6 is 0:<sup>83595</sup> · <sup>10</sup>�4. (b) 41 Chebyshev collocation nodes and 100 evaluation nodes, maximal error for <sup>t</sup> <sup>¼</sup> <sup>0</sup>:5, <sup>α</sup> <sup>¼</sup> <sup>0</sup>:5, <sup>ε</sup> <sup>¼</sup> <sup>3</sup>:6 is 0:<sup>14385</sup> · <sup>10</sup>�4.

#### 6.3.2. Example 2

6.3. Caputo time fractional partial differential equations

14 Fractal Analysis - Applications in Physics, Engineering and Technology

<sup>∂</sup><sup>α</sup>uðx, t<sup>Þ</sup> ∂t

<sup>α</sup> þ δ

type

software/rbf\_qr

6.3.1. Example 1

Using initial condition

Dirichlet boundary conditions

The exact solution (see [33]) is

at time t ¼ 1.

phenomenon [34, 35].

<sup>ε</sup> <sup>¼</sup> <sup>3</sup>:6 is 0:<sup>83595</sup> · <sup>10</sup>�<sup>4</sup>

error for <sup>t</sup> <sup>¼</sup> <sup>0</sup>:5, <sup>α</sup> <sup>¼</sup> <sup>0</sup>:5, <sup>ε</sup> <sup>¼</sup> <sup>3</sup>:6 is 0:<sup>14385</sup> · <sup>10</sup>�<sup>4</sup>

Putting <sup>δ</sup> <sup>¼</sup> <sup>1</sup>, <sup>γ</sup> ¼ �1 and <sup>f</sup>ðx, tÞ ¼ <sup>2</sup><sup>t</sup>

fractional and linear Burger equation

<sup>∂</sup><sup>α</sup>uðx, t<sup>Þ</sup> ∂t <sup>α</sup> þ

In this part we consider examples of fractional partial differential equations (see [32]) of the

þ γ

where t > 0, x∈ ½a, b�, 0 < α ≤ 1, δ and γ are real parameters, bounded initial condition u(x,0) ¼ u0ðxÞ and Dirichlet boundary conditions uða, tÞ ¼ g1ðtÞ and uðb, tÞ ¼ g2ðtÞ, t ≥ 0. Caputo Derivative will be considered and rbf-qr routines from www.it.uu.se/research/scientific\_computing/

<sup>∂</sup><sup>2</sup>uðx, t<sup>Þ</sup>

<sup>∂</sup>x<sup>2</sup> <sup>¼</sup> <sup>f</sup>ðx, tÞ, <sup>ð</sup>31<sup>Þ</sup>

þ 2x � 2: ð32Þ

, ð33Þ

: ð34Þ

, ð35Þ

<sup>Γ</sup>ð3�α<sup>Þ</sup> <sup>þ</sup> <sup>2</sup><sup>x</sup> � 2 in Eq. (31), we obtain a non-homogeneous,

2�α Γð3 � αÞ

2

.] [41 Chebyshev collocation nodes and 100 evaluation nodes, maximal

∂uðx, tÞ ∂x

2�α

<sup>∂</sup><sup>x</sup> � <sup>∂</sup><sup>2</sup>uðx, t<sup>Þ</sup>

<sup>u</sup>ðx; <sup>0</sup>Þ ¼ <sup>x</sup><sup>2</sup>

<sup>u</sup>ðx, tÞ ¼ <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>t</sup>

This problem is solved using the method described in section 6.1 along with QR decomposition for spatial part and an implicit scheme for time, on domain ½0; 1�, fractional derivative of order α ¼ 0:5; several shape parameters are considered and the one that gives the best approximation is taken. Figure 3 shows a 3D image of the solution and maximal error for the solution

In Figure 4 we compare errors that result when choosing uniform collocation nodes (with a fixed step size) and Chebyshev, on definition interval ½0; 1�. The election of Chebyshev nodes is due to attenuating instability, which manifests as oscillations in the graph, called the Gibbs

[41 uniform collocation nodes and 100 evaluation nodes, maximal error for t ¼ 0:5, α ¼ 0:5,

.]

2

<sup>∂</sup>x<sup>2</sup> <sup>¼</sup> <sup>2</sup><sup>t</sup>

, uð1, tÞ ¼ 1 þ t

2

∂uðx, tÞ

uð0, tÞ ¼ t

Putting <sup>δ</sup> <sup>¼</sup> 1, <sup>γ</sup> <sup>¼</sup> 0 and <sup>f</sup>ðx, tÞ ¼ <sup>t</sup> 1�α <sup>Γ</sup>ð2�α<sup>Þ</sup> sin <sup>ð</sup>xÞ þ <sup>t</sup> cos <sup>ð</sup>x<sup>Þ</sup> in Eq. (31), we obtain

$$\frac{\partial^a u(\mathbf{x},t)}{\partial t^a} + \frac{\partial u(\mathbf{x},t)}{\partial \mathbf{x}} = \frac{t^{1-a}}{\Gamma(2-a)} \sin(\mathbf{x}) + t \cos(\mathbf{x}).\tag{36}$$

Initial condition

$$
\mu(\mathbf{x}, \mathbf{0}) = \mathbf{0},
\tag{37}
$$

Next function is the exact solution (see [33]), which is used to set boundary conditions.

$$u(\mathbf{x},t) = t\sin\mathbf{x}.\tag{38}$$

The problem is solved by Lagrange (like Section 6.1) and RBF-QR, for α ¼ 0:6, N ¼ 121 collocation nodes on interval ½�3; 3�. In Figure 5 we notice that when the number of uniform

Figure 5. Unstability due to node election, for Eq. (36). (a) 121 uniform collocation nodes and 101 evaluation nodes, maximal error for <sup>t</sup> <sup>¼</sup> 2, <sup>α</sup> <sup>¼</sup> <sup>0</sup>:6, <sup>ε</sup> <sup>¼</sup> <sup>1</sup>:5 is 0:<sup>10558</sup> · 1019. (b) 121 Chebyshev collocation nodes and 101 evaluation nodes, maximal error for <sup>t</sup> <sup>¼</sup> 2, <sup>α</sup> <sup>¼</sup> <sup>0</sup>:6, <sup>ε</sup> <sup>¼</sup> <sup>0</sup>:2 is 0:<sup>71054</sup> · <sup>10</sup>�15.

Figure 6. Numerical solution for equation (39), order α ¼ 0:5 (superdiffusive phenomenon).

Figure 7. Numerical solution for Eq. (39), order α ¼ 0:7 (subdiffusive phenomenon).

Applications of Radial Basis Function Schemes to Fractional Partial Differential Equations http://dx.doi.org/10.5772/67892 17

Figure 8. Numerical solution for Eq. (39), order α ¼ 1:2 (superdiffusive phenomenon).

collocation nodes grows, there are large oscillations which cause unstability. For this reason, Chebyshev nodes are also preferred, as they attenuate this misbehavior.

#### 6.3.3. Example 3

Putting δ ¼ 0, γ ¼ �1 and fðx, tÞ ¼ 0 in Eq. (31), we obtain

$$\frac{\partial^{\alpha}u(\mathbf{x},t)}{\partial t^{\alpha}} = \frac{\partial^{2}u(\mathbf{x},t)}{\partial \mathbf{x}^{2}}.\tag{39}$$

using initial condition

$$
\mu(\mathbf{x}, \mathbf{0}) = 4\mathbf{x}(1-\mathbf{x}),
\tag{40}
$$

boundary conditions

$$
\mu(0,t) = \mu(1,t) = 0.\tag{41}
$$

The exact solution of this problem is not known, but is shown as an example of subdiffusive phenomenon which was discussed in the introduction. The problem is solved using the method described in section 6.1 along with QR decomposition for spatial part and an implicit scheme for time, on domain ½0; 1�, fractional derivative of order α ¼ 0:5 and α ¼ 0:7, according to the results shown in Refs. [36, 32]. Figures 6 and 7 show these results. Figure 8 shows a what if situation when α ¼ 1:2 (superdiffusive phenomenon), mentioned in Section 1.2.

#### 7. Discussion

0 0.2 0.4 0.6 0.8 1

(a)

RBF QR, Lagrange, alpha = 0.7, 51 collocation nodes and 100 evaluation nodes

(a)

0 0.2 0.4 0.6 0.8 1

**x**

Figure 6. Numerical solution for equation (39), order α ¼ 0:5 (superdiffusive phenomenon).

0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.2 0.4 0.6 0.8 1

**t**

**x**

Figure 7. Numerical solution for Eq. (39), order α ¼ 0:7 (subdiffusive phenomenon).

0 0.2 0.4 0.6 0.8 1

**x**

<sup>4</sup> <sup>2</sup> <sup>0</sup> <sup>2</sup> <sup>4</sup> <sup>2</sup>

Caputo approximation Exact solution

x

1.5 1 0.5 0 0.5 1 1.5 2

Figure 5. Unstability due to node election, for Eq. (36). (a) 121 uniform collocation nodes and 101 evaluation nodes, maximal error for <sup>t</sup> <sup>¼</sup> 2, <sup>α</sup> <sup>¼</sup> <sup>0</sup>:6, <sup>ε</sup> <sup>¼</sup> <sup>1</sup>:5 is 0:<sup>10558</sup> · 1019. (b) 121 Chebyshev collocation nodes and 101 evaluation nodes,

u

0 0.2 0.4 0.6 0.8 1

0 0.2 0.4 0.6 0.8 1

**u**

**u**

0.05 0.06

(b)

RBF QR, Lagrange, alpha = 0.7, 51 collocation nodes and 100 evaluation nodes

(b)

0 0.2 0.4 0.6 0.8 1

**x**

0.05 0.06

<sup>0</sup> 0.01

0

0.01 0.02

0.03 0.04

**t**

0.02 0.03

<sup>0</sup> 0.5 <sup>1</sup> 1.5 <sup>2</sup> <sup>0</sup>

epsilon

0.5

(b)

RBF QR, Lagrange, alpha = 0.5, 51 collocation nodes and 100 evaluation nodes

1

1.5


2

2.5 x 1087 MAXIMAL ERROR

0.04

**t**

0 0.01 0.02 0.03 0.04 0.05 0.06 0 0.2 0.4 0.6 0.8 1

**t**

<sup>4</sup> <sup>2</sup> <sup>0</sup> <sup>2</sup> <sup>4</sup> <sup>12</sup>

x

<sup>0</sup> 0.5 <sup>1</sup> 1.5 <sup>2</sup> <sup>0</sup>

epsilon

<sup>18</sup> x 10148 MAXIMAL ERROR

(a)

RBF QR, Lagrange, alpha = 0.5, 51 collocation nodes and 100 evaluation nodes

maximal error for <sup>t</sup> <sup>¼</sup> 2, <sup>α</sup> <sup>¼</sup> <sup>0</sup>:6, <sup>ε</sup> <sup>¼</sup> <sup>0</sup>:2 is 0:<sup>71054</sup> · <sup>10</sup>�15.


16 Fractal Analysis - Applications in Physics, Engineering and Technology

Caputo approximation Exact solution

u

**u**

**u**

The idea was to show that radial basis schemes are efficient and are on par with schemes like Finite Differences. They are an option to deal with multidimensional and irregular domain problems. The challenge is to adapt them to deal with diffusive problems, particularly with multidimensional systems of equations that consider that the medium does not have a single characteristic.

Fuentes et al. show in the work "The fractal models of saturated and unsaturated flows (capillary, diffusion and matrix) at micro, macro and mega scales", for the PEMEX oil company, how to model the pressure which the oil must leave in a wellbore, from a system of triple porosity and permeability of the fractured medium.

By converting the model into a dimensionless system, in terms of radial distances, non-mesh schemes such as RBFs sound like a viable option and it is a work that is still being addressed.

## Author details

Carlos Alberto Torres Martínez<sup>1</sup> \* and Carlos Fuentes<sup>2</sup>

\*Address all correspondence to: inocencio3@gmail.com

1 Universidad Nacional Autónoma de México, Universidad Autónoma de la Ciudad de México (UACM), Mexico

2 Instituto Mexicano de Tecnología del Agua, Coordinación de Tecnología de Riego y Drenaje, Cuernavaca, Mexico

## References


[7] Marjan Uddin. On the selection of a good value of shape parameter in solving timedependent partial differential equations using rbf approximation method. Applied Mathematical Modelling, 38(1):135–144, 2014.

multidimensional systems of equations that consider that the medium does not have a single

Fuentes et al. show in the work "The fractal models of saturated and unsaturated flows (capillary, diffusion and matrix) at micro, macro and mega scales", for the PEMEX oil company, how to model the pressure which the oil must leave in a wellbore, from a system of triple

By converting the model into a dimensionless system, in terms of radial distances, non-mesh schemes such as RBFs sound like a viable option and it is a work that is still being addressed.

\* and Carlos Fuentes<sup>2</sup>

1 Universidad Nacional Autónoma de México, Universidad Autónoma de la Ciudad de

2 Instituto Mexicano de Tecnología del Agua, Coordinación de Tecnología de Riego y Drenaje,

[1] Rolland L Hardy. Multiquadric equations of topography and other irregular surfaces.

[2] Charles A Micchelli. Interpolation of scattered data: distance matrices and conditionally positive definite functions. In Approximation theory and spline functions pp. 143–145.

[3] Michael JD Powell. The theory of radial basis function approximation in 1990. Department of Applied Mathematics and Theoretical Physics, University of Cambridge. pp.

[4] Edward J Kansa. Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics—i surface approximations and partial derivative

[5] Edward J Kansa. Multiquadrics—a scattered data approximation scheme with applications to computational fluid-dynamics—ii solutions to parabolic, hyperbolic and elliptic partial differential equations. Computers & mathematics with applications, 19(8):147–161, 1990. [6] Wen Chen, Zhuo-Jia Fu, and Ching-Shyang Chen. Recent advances in radial basis func-

estimates. Computers & Mathematics with applications, 19(8):127–145, 1990.

characteristic.

Author details

Carlos Alberto Torres Martínez<sup>1</sup>

Springer Netherlands. 1984.

105–209. 1990.

México (UACM), Mexico

Cuernavaca, Mexico

References

porosity and permeability of the fractured medium.

18 Fractal Analysis - Applications in Physics, Engineering and Technology

\*Address all correspondence to: inocencio3@gmail.com

Journal of geophysical research, 76(8):1905–1915, 1971.

tion collocation methods. Heidelberg: Springer, 2014.

