- Email: [email protected]
Discretization of fractional order differentiator over Paley–Wiener space
Applied Mathematics and Computation 247 (2014) 162–168
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Discretization of fractional order differentiator over Paley–Wiener space Yan Wu Department of Mathematical Sciences, Georgia Southern University, Statesboro, GA 30460, USA
a r t i c l e
i n f o
Keywords: Bandlimited functions Extrapolation Fractional order differentiator Grünwald–Letnikov approximation Interpolation Paley–Wiener space
a b s t r a c t The Paley–Wiener space consists of functions whose Fourier transform is compactly supported in the frequency domain. In the context of signal processing, such functions are also known as bandlimited signals, which represent a large class of signals in signal processing. An analog fractional order differentiator is representable by way of the Cauchy integral formula, special functions, as well as the Fourier/Laplace transformer, while a digital differentiator of fractional order can be obtained through direct or indirect discretization techniques. In this paper, we present the design of a finite impulse response (FIR) filter that discretizes the fractional order differentiator over functions in Paley–Wiener space. The proposed FIR model has some meritorious properties that are preferred in applications: the filter coefficients are independent of the signal samples; it is capable of interpolating or extrapolating at an arbitrary point in the sampling domain; it is adaptive to uniform or non-uniform sampling scenarios. We present explicit formulas on the matrices that lead to computing the filter coefficients. A closed form formula on the error of approximation is derived to demonstrate the accuracy of the proposed discretization model of fractional derivatives. Numerical results also show that the proposed method is more computationally efficient than the well-known methods such as the Grünwald–Letnikov approximation. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction The concept of fractional order differentiation and integration came into being at the same time as the conventional integral order calculus. The early development of fractional calculus is associated with the names such as Euler, Riemann, Leibniz, Lagrange, Fourier, Laplace, and so on. A thorough historical review of fractional calculus can be found in [1]. Fractional order differentiation or a fractional derivative is considered a generalization of an integral order derivative, however, with its own characteristics. Only in recent years have fractional derivatives become important mathematical tools in modeling various physical phenomena. Some also believe that fractional order differential equations are more effective in modeling nonlinearity. Podlubny gave celebrated geometrical and physical interpretation of fractional derivative in [2]. Meanwhile, the applications of the fractional derivative are readily seen in areas such as thermodynamics, electromagnetic field, quantum mechanics [3–5], signal processing, and control [6,7], and cardiac electrophysiology [8], to name a few. In this paper, we are concerned with constructing a discrete analog of a fractional order differentiator over Paley–Wiener space, which consists of functions with compact support in the frequency domain. These types of functions are also known as bandlimited signals, representing an important class of signals in signal processing. In practice, all signals can be filtered into
E-mail address: [email protected] http://dx.doi.org/10.1016/j.amc.2014.08.089 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
Y. Wu / Applied Mathematics and Computation 247 (2014) 162–168
163
bandlimited signals, either low-passed or band-passed, based on the appropriate filters. Tseng, Pei and Hsia introduced their design of a digital FIR differentiator in [6] based on Cauchy integral formula and Fourier transform to discretize the fractional derivative kernel in the Fourier frequency domain. Their digital fractional differentiator (DFD) is capable of estimating the fractional derivative at the current sampling point (uniform sampling). The performance of their DFD is demonstrated with numerical experiments. The DFD proposed in this paper has the following novel and meritorious features: (i) it is capable of interpolating and extrapolating fractional derivatives at an arbitrary point in the sampling domain; (ii) it is adaptive to uniform or non-uniform sampling; (iii) it comprises the discretization of an integral order differentiator as a special case. We present a closed form bound on the error of approximation to prove the accuracy of the proposed model. The rest of this paper is organized as follows. In Section 2, we give necessary background on fractional derivatives along with a list of formulas that motivate the proposed design, the fractional order differentiator and underlying function space. The main results are reported in Section 3, including derivations, analysis, and numerical tests, followed by conclusive remarks and potential applications of the proposed DFD. 2. Definitions n
n
Unlike general integral order derivatives, d f =dx , where one can obtain higher order derivatives by repeating the process of first-order differentiation, the fractional order derivative, denoted by, Da , where a 2 Rþ , is defined through a fractional integral,
Da f ðtÞ ¼
1 CðaÞ
Z
t
ðt sÞa1 f ðsÞds
a
known as a Riemann–Liouville fractional integral. As such a fractional derivative is obtained as n
Da f ðtÞ ¼
d ðnaÞ f ðtÞ nD dt
ð1Þ
where n ¼ dae. The more popular Caputo fractional derivative [9] is particularly useful for solving fractional differential equations,
Da f ðtÞ ¼
1 Cðn aÞ
Z
t
a
f
ðnÞ
ðsÞ
ðt sÞaþ1n
ds;
ðn 1 6 a < nÞ:
ð2Þ
There are a few other alternative formulas for fractional order derivatives. The most relevant formula to this paper is the Grünwald–Letnikov derivative a
Da f ðtÞ ¼ lim h h!0
bðt aÞ=hc X
ð1Þj
j¼0
a f ðt jhÞ: j
Suppose we only approximate the derivative at the mesh points, t j ¼ t0 þ jh; j ¼ 0; 1; . . . ; N, the Grünwald–Letnikov approximation formula can be written as [10] a
Da f ðtN Þ¼h _
N X
xðj aÞ f ðtNj Þ
ð3Þ
j¼0 ðaÞ
ðaÞ
where x0 ¼ 1; xj
ðaÞ
¼ ð1 ða þ 1Þ=kÞxj1 ; j ¼ 1; 2; 3; . . . N. Eq. (3) is clearly a finite difference approximation of the frac-
tional derivative of order a. There are higher-order finite difference methods for numerical differentiation of integral-order derivatives, such as the ones derived from Taylor series and those based on Richardson extrapolation. Consequently, those difference formulas are generalized to approximate fractional derivatives [10–12]. Lubich also developed higher order method by way of the numerical quadrature of fractional order in [13]. Moreover, in the same paper, he also laid groundwork in the convergence and stability analysis of numerical methods for discretized fractional calculus. The proposed algorithm in this paper for a discrete fractional differentiator is motivated by (3). We observe from (3) that a fractional derivative is a process with memory. The derivative depends on the function’s historical values. Unlike the integral-order derivatives, where only local function values are needed, in the case of a fractional derivative, it depends on function values over the whole domain and so is a nonlocal operator. This is also observed from the integral definitions of fractional derivatives in (1) and (2). There are a number of known formulas, see [1], in addition to the above mentioned ones, for computing fractional order derivatives, none of which would readily produce closed form, explicit expressions for the derivative, except for some rather basic functions. This is due to the complexity of the definitions of fractional derivatives, and most of the differentiation rules for ordinary derivative cannot be carried over to the fractional derivative. Therefore, there is a need for developing computational methods that are capable of approximating fractional derivatives numerically with high accuracy, whilst the method itself has a simple structure for implementation. The model proposed in this paper meets this criterion since it is formulated in the form of a convolution between a set of optimal coefficients and a set of functional samples, and it is capable of
164
Y. Wu / Applied Mathematics and Computation 247 (2014) 162–168
approximating fractional derivatives of any non-integral order in the form of extrapolation or interpolation. We present details of the design and performance analysis in Section 3. 3. Discrete fractional differentiator In this paper, we are concerned with the fractional derivative of a class of functions from Paley–Wiener space. In signal processing, a bandlimited signal, or a function whose Fourier transform has a compact support in the frequency domain, usually does not have a closed form representation (in the time domain), hence, making it impossible to compute its fractional derivative analytically. Furthermore, in the context of signal processing, only signal samples are required for all the computational needs. It is unnecessary to acquire an algebraic expression of the signal. It is indeed a mainstream endeavor in signal processing to develop innovative and efficient algorithms for the purpose of signal transformation. In our case, we design a simple digital filter with finite impulse response for computing the fractional derivative of a bandlimited signal off its samples. As mentioned earlier, a bandlimited function belongs to the Paley–Wiener space defined in [14], Definition 1. Let r > 0 and 1 6 p < 1, the Paley–Wiener space, denoted by PW pr , consists of functions f with a Rr Rr representation f ðzÞ ¼ r gðuÞeizu du; z 2 C, for some g 2 Lp ðr; rÞ, where r jg jp du < 1. A thorough characterization of the Paley–Wiener space can be found in [14]. In signal processing, specifically, the functions in the Paley–Wiener space are also known as bandlimited functions with p ¼ 2 and r relating to the bandwidth of f, R1 whereas g is the Fourier transform of f, i.e. gðwÞ ¼ 1 f ðtÞeiwt dt. It can be deduced from Definition 1 that f satisfies f 2 L2 ðRÞC 1 ðRÞ. Obviously, the fractional derivative of a bandlimited function is still a bandlimited function. Therefore, let Da : PW r ! PW r , be the fractional differentiator of order a. Since we are only interested in approximating the discrete values ~ a under the case of uniform sampling with T the sampling of fractional derivative, we denote the discretization of Da by D T
~ a corresponds to the case of non-uniform sampling. The proposed digital fractional differentiator in the form interval, while D of a FIR filter is given as follows, n X ~ a f ðt lTÞ ¼ 1 D ak f ðt kTÞ; a T T k¼1
ð4Þ
where the ak ’s are the filter coefficients, and l 2 R. The factor 1=T a is consistent with the Grünwald–Letnikov derivative (3). It is also due to the derivation of the filter coefficients ak to be shown later. Formula (4) is an interpolator for fractional derivatives if 1 6 l 6 n; otherwise, it is an extrapolator. Let kk be an induced norm of the differentiator over the Paley–Wiener space, which is defined as Da ¼ supt2R;kf k¼1 Da f ðtÞ. The fractional derivative of a function is represented in terms of the Fourier transform,
Da f ðtÞ ¼
1 2p
Z pr
a
ðixÞ FðxÞeixt dx
ð5Þ
pr
where FðwÞ is the Fourier transform of f ðtÞ. It is understood that the bandwidth of f ðtÞ is r=2. Now, consider the error of ~ a Þf ðt lTÞ, and applying the Cauchy–Schwarz inequality as well as Parseval’s identity to get approximation ðD D T
Z " # 2 pr n 1X a ilxT ixt ikxT dx FðxÞe ðixÞ e a ak e 2 T k¼1 ð2pÞ pr Z pr 2 n 1 1 1X 2 a 6 kF k2 ak eikxT dx ¼ kf k E2l ðixÞ eilxT a 2p 2p pr T k¼1
2 ~ a Þf ðt lTÞ ¼ ðD D T
where
1
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ffi u n u 1 Z pr 1X a ilxT t El ¼ a ak eikxT dx: ðixÞ e 2p pr T k¼1
ð6Þ
~ a Hence, we have D D T ¼ El since the equality is achievable for certain f in the Cauchy–Schwarz inequality. Consequently, the discrepancy between the fractional differentiator and its discretization is measured by El . The filter coefficients ak are then obtained through minimizing El , or equivalently E2l , realizing the fact that E2l is in a quadratic form. This observation is useful for quantifying the error as in the performance analysis discussed later. The integral in (6) can be rewritten with a substitution w ¼ xT, such that 2
El ¼ ð2pT
2aþ1 1
Þ
2 Z ps n X a ilw ikw ak e dw ðiwÞ e ps k¼1
ð7Þ
where s ¼ rT known as time-bandwidth product. Since E2l is a quadratic form, there exists a set of optimal filter coefficients
165
Y. Wu / Applied Mathematics and Computation 247 (2014) 162–168
that minimize E2l satisfying the system of equations, @E2l =@aj ¼ 0; j ¼ 1; 2; . . . ; n. The result is stated in Theorem 1. ~ a that discretize the fractional order differentiator Da satisfy the system of Theorem 1. The optimal filter coefficients fak g of D T equations Ax ¼ b as a result of minimizing the performance index E2l , where
2 6 6 A¼6 6 4
sð0Þ
sðsÞ
sðsÞ
sð0Þ
. . . sððn 1ÞsÞ
7 . . . sððn 2ÞsÞ 7 7; 7 ... 5
...
sððn 1ÞsÞ sððn 2ÞsÞ . . . with sðtÞ ¼ sinðptÞ=ðptÞ , bk ¼ 2ps 1
3
R ps
ps
sð0Þ
3 b1 6 7 6 b2 7 6 7 b ¼ 6. 7 6. 7 4. 5 2
ð8Þ
bn
wa cosððk lÞw þ ap=2Þdw;
k ¼ 1; 2; . . . ; n, and x ¼ ½ a1
a2
...
a n T .
Proof. We first rewrite E2l into
E2l ¼ ð2pT 2aþ1 Þ
1
#" # Z ps " n n X X a a ðiwÞ eilw ak eikw ðiwÞ eilw ak eikw dw: ps
k¼1
k¼1
It is known that the optimal coefficients fak g that minimize (7) satisfy the system of equations @E2l =@aj ¼ 0; j ¼ 1; 2; . . . ; n, where
@E2l @aj
¼ ð2pT 2aþ1 Þ
1
" # " #) Z ps ( n n X X a a eijw ðiwÞ eilw ak eikw eijw ðiwÞ eilw ak eikw dw ps
k¼1
k¼1
Z ps X Z ps n 1 1 a a ¼ ð2pT 2aþ1 Þ ak ðeiðkjÞw þ eiðkjÞw Þdw ð2pT 2aþ1 Þ ½ðiwÞ eiðjlÞw þ ðiwÞ eiðjlÞw dw: ps k¼1
ps
The first integral in @E2l =@aj can be evaluated quickly as
Z ps n n X X ak ðeiðkjÞw þ eiðkjÞw Þdw ¼ 4 ak sinððk jÞpsÞ=ðk jÞ; k¼1
ps
k¼1
while the second integral is rewritten into a real form
Z ps
wa ei½ðjlÞwþpa=2 þ ei½ðjlÞwþpa=2 dw ¼ 2
ps
Z ps
wa cos½ðj lÞw þ pa=2dw:
ps
By setting @E2l =@aj ¼ 0, one has
Z ps n X 1 ak sinððk jÞpsÞ=ðpsðk jÞÞ ¼ wa cosððj lÞw þ pa=2Þdw; 2ps ps k¼1
j ¼ 1; 2; . . . ; n;
The integral for bk has to be evaluated numerically if a is non-integer. Numerical quadrature, such as the Romberg method and Gaussian quadrature, produces accurate result. Since the matrix A in (8) is a symmetric Toeplitz matrix, the system can be solved iteratively by using the well-known Levinson’s method [15]. This method is particularly effective when the matrix A is ill-conditioned. It is observed from (7) that El ! 0 as s ! 0 (or T ! 0Þ. This implies that the discrete fractional differentiator (4) yields more accurate approximation of a fractional derivative for smaller sampling interval T or a combined effect of the time-bandwidth product. This in turn also makes the matrix A ill-conditioned. Going through the derivations, it is readily seen that the filter coefficients in (4) are independent of function values at the sampling points. They only depend on the time-bandwidth product and the point of approximation, i,e, l, which are usually known a priori. The coefficients are fixed throughout the convolution process as the fractional derivative is computed, a standard feature of a digital filter such as the one in (4). It is clear that the discrete fractional differentiator (4) is based on uniform sampling. For the case of non-uniform sampling, the derivation is similar but with some modifications. First off, the FIR filter based on non-uniform sampling is given by
~ a f ðt t e Þ ¼ D
n X ak f ðt tk Þ; k¼1
~ a where t i < tiþ1 ; i ¼ 1; 2; . . . ; n 1. It can be shown that the error of approximation is D D ¼ Ee , where
ð9Þ
166
Y. Wu / Applied Mathematics and Computation 247 (2014) 162–168
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ffi u n u 1 Z pr X a Ee ¼ t ak eixtk dx ðixÞ eixte 2p pr k¼1
ð10Þ
The filter coefficients ak are obtained from minimizing (10), and they satisfy the system of equations Ax ¼ b where
2
sððt 2 t1 ÞrÞ . . . sððtn t1 ÞrÞ
sð0Þ
6 sððt t ÞrÞ 1 2 6 A¼6 4 ...
2
3
. . . sððtn t2 ÞrÞ 7 7 7; 5 ...
sð0Þ
sððt n t 1 ÞrÞ sððt n t 2 ÞrÞ . . .
sð0Þ
b1
3
6b 7 6 27 7 b¼6 6 .. 7 4. 5
ð11Þ
bn
R pr
1 where bk ¼ 2pr xa cosððtk te Þx þ ap=2Þdx; k ¼ 1; 2; . . . ; n. It is noted that the matrix A in (11) is symmetric but no pr longer Toeplitz. An interesting case of non-uniform sampling is the so-called interlaced sampling. The non-uniform sampling itself can be cumbersome as one has to keep track with the sampling points; however, if such sampling data are available, the proposed model (9) is ready to generate the fractional derivatives at the designated points, whereas most methods in the literature only work with uniform sampling. It is worth noting that the above derivations still hold for the case where a is a positive integer, corresponding to the conventional differentiator. Therefore, the proposed model for discretization of the fractional differentiator is also applicable to the integral order differentiator. Before presenting some numerical results, we analyze the performance of the proposed models (4) and (9). We derive an alternative expression on the error of approximation E2l , followed by an explicit formula on the minimized E2l . As mentioned earlier, the integral on E2l in (7) is in a quadratic form. To this end, we rewrite (7) into
2
El ¼ ð2pT ¼ ð2pT
2aþ1 1
Þ
Z ps " ps
2aþ1 1
Þ
Z ps
"
a ilw
ðiwÞ e
n X ak eikw
#"
a il w
ðiwÞ e
k¼1
# n X ikw dw ak e k¼1
# X XX a a a iðklÞw a iðklÞw iðjkÞw ðiwÞ ðiwÞ ak ððiwÞ e þ ðiwÞ e Þþ ak aj e dw
ps
k
k
j
After evaluating the integral over each term, we have
X XX E2l ¼ 2ðpsÞ2aþ1 =ð2a þ 1Þ 4ps ak bk þ 2ps ak aj Akj k
k
j
where bk and Akj , entries of the matrix A, are defined in Theorem 1. Consequently, the error E2l is written as
E2l ¼
r T 2a
where u ¼ ½ 1
" H¼
uT H u
ð12Þ
T
xT and x is the vector of filter coefficients ak , and the symmetric matrix H
h b b
T
#
A
ð13Þ
where h ¼ ðpsÞ2a =ð2a þ 1Þ . We observe that the matrix H is positive definite because the quadratic form (7) is positive definite. It is well known that, given a positive definite matrix H and a vector x, the following inequality holds [16], xT Hx P k0 kxk22 , and the equality is achieved if x ¼ v , where k0 is the smallest eigenvalue of H, andv is the corresponding eigenvector. Equipped with this result, T let v ¼ ½ v 0 v 1 . . . v n and kv k2 ¼ 1, the error of approximation (12) is minimized by choosing u ¼ v =v 0 and the minimal value is
min E2l ¼
k0 r T 2a v 20
ð14Þ
The implication of (14) is twofold: it first gives a quantitative measure of the error of approximation on the fractional derivative via (4); secondly, the optimal set of filter coefficients can be obtained from the eigenvector v above with ak ¼ v k =v 0 ; k ¼ 1; 2; . . . ; n. Hence, it provides an alternative approach for computing the filter coefficients. It is known that the matrix H is ill-conditioned since the Toeplitz matrix A is ill-conditioned with small time-bandwidth product. This implies that the smallest eigenvalue of H is close to zero, hence indicating accurate approximation of the fractional derivative due to (14). Numerical experiments are performed to demonstrate the accuracy of the proposed discrete fractional differentiator design and the results are consistent. We construct the Paley–Wiener functions and their fractional derivatives via Definition 1 and (5), respectively. For instance, the bandlimited function used in Fig. 1 is generated from the frequency domain with R pr FðwÞ ¼ w2 , i.e. f ðtÞ ¼ 1=ð2pÞ pr FðwÞeiwt dw, where r ¼ 6. We compute the precise fractional derivative, in this case the order is a ¼ 5=3, through (5). The sampling interval is T ¼ 0:05 and the length of the filter is n ¼ 6. The functional samples are con-
167
Y. Wu / Applied Mathematics and Computation 247 (2014) 162–168 6 Exact Approximation
4
2
0
-2
-4
-6
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
Fig. 1. Numerical approximation of fractional derivatives of a bandlimited function.
3 Exact Approximation
2.5
2
1.5
1
0.5
0
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5 1=2
Fig. 2. Numerical approximation of half-derivative of the function x with D
5
pffiffiffi pffiffiffiffi x ¼ 2 x= p.
voluted with the filter coefficients according to (4) to extrapolate the fractional derivative, i.e. l ¼ 0 in (4). The result is shown in Fig. 1. The point-wise error of approximation of the fractional derivative shown in Fig. 1 is between 104 and 103 . We also tested an interesting case with polynomial functions. Since polynomial functions have zero content in the frequency domain, we let the bandwidth parameter r approach zero when computing the filter coefficients of (4). Meanwhile, it is known that Da xn ¼ Cðn þ 1Þxna =Cðn a þ 1Þ as the exact fractional derivative of xn . When approximating the half-derivative of x, see Fig. 2, we choose T ¼ 0:1 and n ¼ 6 in (4). The approximation matches the exact half-derivative function well as shown in Fig. 2. Meanwhile, the upper bound on the error of approximation in this example is 0.082. The Grünwald–Letnikov methods for approximating fractional derivatives are frequently cited in the literature, see [10,11]. We performed comparison tests between the proposed approximation method and the Grünwald–Letnikov methods. The results show that our method yielded more accurate approximations than the Grünwald–Letnikov methods, and our method is more efficient as well, when we apply those methods to Paley–Wiener functions. A possible explanation is that the Grünwald–Letnikov methods are obtained from a generalization of the conventional numerical differentiation formulas, such as backward difference and central difference, which are conveniently derived from the Taylor series. Meanwhile, the truncation error of the Taylor series, which depends on the magnitude of the derivative, dictates the error of approximation. On the contrary, the proposed method in this paper is derived from minimizing the error of approximation directly. Hence, it is understood as the most accurate method in the sense of least squares optimization. The following table demonstrates the efficiency of our method versus Grünwald–Letnikov methods. The testing function is again obtained from the Fourier transform with FðwÞ ¼ w3 and bandwidth r ¼ 10, and its fractional derivatives are computed from (5) with
Table 1 Comparison between the proposed method and Grünwald–Letnikov methods. Error/methods
G–L backward difference
G–L central difference
Proposed
102, T = 0.1 103, T = 0.05 104, T = 0.02
16 21 28
13 16 22
7 8 10
168
Y. Wu / Applied Mathematics and Computation 247 (2014) 162–168
a ¼ 1=3. We fix the error of approximation at different levels and compare the number of samples needed to meet the error (see Table 1). In this paper, we present numerical methods for computing the fractional derivatives of functions from the Paley–Wiener space. These functions are also known as bandlimited signals in signal processing. Therefore, the proposed digital filters are applicable for generating fractional order sequences from signal samples, which can be used in signal encryption in secure communication. The filter is efficient as its coefficients are invariant with respect to functional samples and the computational process is a simple discrete convolution. The accuracy of the discrete model is demonstrated with a closed form formula on the minimized error of approximation as well as numerical experiments. All these make it possible to reconstruct a smooth fractional derivative of a rather complicated function, whereas its fractional derivative is impossible to obtain analytically. Acknowledgment The author would like to thank the anonymous reviewers for their valuable suggestions that helped improve the presentation of the original manuscript. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.amc.2014.08.089. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
K.B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, New York and London, 1974. I. Podlubny, Geometric and physical interpretation of fractional integration and fractional differentiation, Fract. Calc. Appl. Anal. 5 (2002) 367–386. L. Debnath, Recent applications of fractional calculus to science and engineering, Int. J. Math. Math. Sci. 54 (2003) 3413–3442. L. Wang, Y. Ma, Z. Meng, Haar wavelet method for solving fractional partial differential equations numerically, Appl. Math. Comput. 227 (2014) 66–76. S. Wheatcraft, M. Meerschaert, Fractional conservation of mass, Adv. Water Resour. 31 (2008) 1377–1381. C. Tseng, S. Pei, S. Hsia, Computation of fractional derivatives using Fourier transform and digital FIR differentiator, Signal Process. 80 (2000) 151–159. A. Razminia, D.F.M. Torres, Control of a novel chaotic fractional order system using a state feedback technique, Mechatronics 23 (2013) 755–763. A. Bueno-Orovio, D. Kay, V. Grau, B. Rodriguez, K. Burrage, Fraction diffusion models cardiac electrical propagation reveal structural heterogeneity effects on dispersion of repolarization, J. R. Soc. Interface 11 (2014). M. Caputo, Linear model of dissipation whose Q is almost frequency independent-II, Geophys. J. R. Astr. Soc. 13 (1967) 529–539. E. Sousa, How to approximate the fractional derivative of order 1 < a < 2, Int. J. Bifur. Chaos 22 (2012). D. Baleanu, O. Defterli, O.P. Agrawal, A central difference numerical scheme for fractional optimal control problems, J. Vibr. Control 15 (2009) 583–597. K. Diethelm, N.J. Ford, A.D. Freed, Y. Luchko, Algorithms for the fractional calculus: a selection of numerical methods, Comput. Methods Appl. Mech. Eng. 194 (2005) 743–773. C. Lubich, Discretized fractional calculus, SIAM J. Math. Anal. 17 (1986) 704–719. J.R. Higgins, Sampling Theory in Fourier and Signal Analysis Foundations, Cearendon Press, Oxford, 1996. W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes, The Art of Scientific Computing, Cambridge University Press, 1986. G.H. Golub, C.F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore and London, 1989.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Discretization of fractional order differentiator over Paley–Wiener space Yan Wu Department of Mathematical Sciences, Georgia Southern University, Statesboro, GA 30460, USA
a r t i c l e
i n f o
Keywords: Bandlimited functions Extrapolation Fractional order differentiator Grünwald–Letnikov approximation Interpolation Paley–Wiener space
a b s t r a c t The Paley–Wiener space consists of functions whose Fourier transform is compactly supported in the frequency domain. In the context of signal processing, such functions are also known as bandlimited signals, which represent a large class of signals in signal processing. An analog fractional order differentiator is representable by way of the Cauchy integral formula, special functions, as well as the Fourier/Laplace transformer, while a digital differentiator of fractional order can be obtained through direct or indirect discretization techniques. In this paper, we present the design of a finite impulse response (FIR) filter that discretizes the fractional order differentiator over functions in Paley–Wiener space. The proposed FIR model has some meritorious properties that are preferred in applications: the filter coefficients are independent of the signal samples; it is capable of interpolating or extrapolating at an arbitrary point in the sampling domain; it is adaptive to uniform or non-uniform sampling scenarios. We present explicit formulas on the matrices that lead to computing the filter coefficients. A closed form formula on the error of approximation is derived to demonstrate the accuracy of the proposed discretization model of fractional derivatives. Numerical results also show that the proposed method is more computationally efficient than the well-known methods such as the Grünwald–Letnikov approximation. Ó 2014 Elsevier Inc. All rights reserved.
1. Introduction The concept of fractional order differentiation and integration came into being at the same time as the conventional integral order calculus. The early development of fractional calculus is associated with the names such as Euler, Riemann, Leibniz, Lagrange, Fourier, Laplace, and so on. A thorough historical review of fractional calculus can be found in [1]. Fractional order differentiation or a fractional derivative is considered a generalization of an integral order derivative, however, with its own characteristics. Only in recent years have fractional derivatives become important mathematical tools in modeling various physical phenomena. Some also believe that fractional order differential equations are more effective in modeling nonlinearity. Podlubny gave celebrated geometrical and physical interpretation of fractional derivative in [2]. Meanwhile, the applications of the fractional derivative are readily seen in areas such as thermodynamics, electromagnetic field, quantum mechanics [3–5], signal processing, and control [6,7], and cardiac electrophysiology [8], to name a few. In this paper, we are concerned with constructing a discrete analog of a fractional order differentiator over Paley–Wiener space, which consists of functions with compact support in the frequency domain. These types of functions are also known as bandlimited signals, representing an important class of signals in signal processing. In practice, all signals can be filtered into
E-mail address: [email protected] http://dx.doi.org/10.1016/j.amc.2014.08.089 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.
Y. Wu / Applied Mathematics and Computation 247 (2014) 162–168
163
bandlimited signals, either low-passed or band-passed, based on the appropriate filters. Tseng, Pei and Hsia introduced their design of a digital FIR differentiator in [6] based on Cauchy integral formula and Fourier transform to discretize the fractional derivative kernel in the Fourier frequency domain. Their digital fractional differentiator (DFD) is capable of estimating the fractional derivative at the current sampling point (uniform sampling). The performance of their DFD is demonstrated with numerical experiments. The DFD proposed in this paper has the following novel and meritorious features: (i) it is capable of interpolating and extrapolating fractional derivatives at an arbitrary point in the sampling domain; (ii) it is adaptive to uniform or non-uniform sampling; (iii) it comprises the discretization of an integral order differentiator as a special case. We present a closed form bound on the error of approximation to prove the accuracy of the proposed model. The rest of this paper is organized as follows. In Section 2, we give necessary background on fractional derivatives along with a list of formulas that motivate the proposed design, the fractional order differentiator and underlying function space. The main results are reported in Section 3, including derivations, analysis, and numerical tests, followed by conclusive remarks and potential applications of the proposed DFD. 2. Definitions n
n
Unlike general integral order derivatives, d f =dx , where one can obtain higher order derivatives by repeating the process of first-order differentiation, the fractional order derivative, denoted by, Da , where a 2 Rþ , is defined through a fractional integral,
Da f ðtÞ ¼
1 CðaÞ
Z
t
ðt sÞa1 f ðsÞds
a
known as a Riemann–Liouville fractional integral. As such a fractional derivative is obtained as n
Da f ðtÞ ¼
d ðnaÞ f ðtÞ nD dt
ð1Þ
where n ¼ dae. The more popular Caputo fractional derivative [9] is particularly useful for solving fractional differential equations,
Da f ðtÞ ¼
1 Cðn aÞ
Z
t
a
f
ðnÞ
ðsÞ
ðt sÞaþ1n
ds;
ðn 1 6 a < nÞ:
ð2Þ
There are a few other alternative formulas for fractional order derivatives. The most relevant formula to this paper is the Grünwald–Letnikov derivative a
Da f ðtÞ ¼ lim h h!0
bðt aÞ=hc X
ð1Þj
j¼0
a f ðt jhÞ: j
Suppose we only approximate the derivative at the mesh points, t j ¼ t0 þ jh; j ¼ 0; 1; . . . ; N, the Grünwald–Letnikov approximation formula can be written as [10] a
Da f ðtN Þ¼h _
N X
xðj aÞ f ðtNj Þ
ð3Þ
j¼0 ðaÞ
ðaÞ
where x0 ¼ 1; xj
ðaÞ
¼ ð1 ða þ 1Þ=kÞxj1 ; j ¼ 1; 2; 3; . . . N. Eq. (3) is clearly a finite difference approximation of the frac-
tional derivative of order a. There are higher-order finite difference methods for numerical differentiation of integral-order derivatives, such as the ones derived from Taylor series and those based on Richardson extrapolation. Consequently, those difference formulas are generalized to approximate fractional derivatives [10–12]. Lubich also developed higher order method by way of the numerical quadrature of fractional order in [13]. Moreover, in the same paper, he also laid groundwork in the convergence and stability analysis of numerical methods for discretized fractional calculus. The proposed algorithm in this paper for a discrete fractional differentiator is motivated by (3). We observe from (3) that a fractional derivative is a process with memory. The derivative depends on the function’s historical values. Unlike the integral-order derivatives, where only local function values are needed, in the case of a fractional derivative, it depends on function values over the whole domain and so is a nonlocal operator. This is also observed from the integral definitions of fractional derivatives in (1) and (2). There are a number of known formulas, see [1], in addition to the above mentioned ones, for computing fractional order derivatives, none of which would readily produce closed form, explicit expressions for the derivative, except for some rather basic functions. This is due to the complexity of the definitions of fractional derivatives, and most of the differentiation rules for ordinary derivative cannot be carried over to the fractional derivative. Therefore, there is a need for developing computational methods that are capable of approximating fractional derivatives numerically with high accuracy, whilst the method itself has a simple structure for implementation. The model proposed in this paper meets this criterion since it is formulated in the form of a convolution between a set of optimal coefficients and a set of functional samples, and it is capable of
164
Y. Wu / Applied Mathematics and Computation 247 (2014) 162–168
approximating fractional derivatives of any non-integral order in the form of extrapolation or interpolation. We present details of the design and performance analysis in Section 3. 3. Discrete fractional differentiator In this paper, we are concerned with the fractional derivative of a class of functions from Paley–Wiener space. In signal processing, a bandlimited signal, or a function whose Fourier transform has a compact support in the frequency domain, usually does not have a closed form representation (in the time domain), hence, making it impossible to compute its fractional derivative analytically. Furthermore, in the context of signal processing, only signal samples are required for all the computational needs. It is unnecessary to acquire an algebraic expression of the signal. It is indeed a mainstream endeavor in signal processing to develop innovative and efficient algorithms for the purpose of signal transformation. In our case, we design a simple digital filter with finite impulse response for computing the fractional derivative of a bandlimited signal off its samples. As mentioned earlier, a bandlimited function belongs to the Paley–Wiener space defined in [14], Definition 1. Let r > 0 and 1 6 p < 1, the Paley–Wiener space, denoted by PW pr , consists of functions f with a Rr Rr representation f ðzÞ ¼ r gðuÞeizu du; z 2 C, for some g 2 Lp ðr; rÞ, where r jg jp du < 1. A thorough characterization of the Paley–Wiener space can be found in [14]. In signal processing, specifically, the functions in the Paley–Wiener space are also known as bandlimited functions with p ¼ 2 and r relating to the bandwidth of f, R1 whereas g is the Fourier transform of f, i.e. gðwÞ ¼ 1 f ðtÞeiwt dt. It can be deduced from Definition 1 that f satisfies f 2 L2 ðRÞC 1 ðRÞ. Obviously, the fractional derivative of a bandlimited function is still a bandlimited function. Therefore, let Da : PW r ! PW r , be the fractional differentiator of order a. Since we are only interested in approximating the discrete values ~ a under the case of uniform sampling with T the sampling of fractional derivative, we denote the discretization of Da by D T
~ a corresponds to the case of non-uniform sampling. The proposed digital fractional differentiator in the form interval, while D of a FIR filter is given as follows, n X ~ a f ðt lTÞ ¼ 1 D ak f ðt kTÞ; a T T k¼1
ð4Þ
where the ak ’s are the filter coefficients, and l 2 R. The factor 1=T a is consistent with the Grünwald–Letnikov derivative (3). It is also due to the derivation of the filter coefficients ak to be shown later. Formula (4) is an interpolator for fractional derivatives if 1 6 l 6 n; otherwise, it is an extrapolator. Let kk be an induced norm of the differentiator over the Paley–Wiener space, which is defined as Da ¼ supt2R;kf k¼1 Da f ðtÞ. The fractional derivative of a function is represented in terms of the Fourier transform,
Da f ðtÞ ¼
1 2p
Z pr
a
ðixÞ FðxÞeixt dx
ð5Þ
pr
where FðwÞ is the Fourier transform of f ðtÞ. It is understood that the bandwidth of f ðtÞ is r=2. Now, consider the error of ~ a Þf ðt lTÞ, and applying the Cauchy–Schwarz inequality as well as Parseval’s identity to get approximation ðD D T
Z " # 2 pr n 1X a ilxT ixt ikxT dx FðxÞe ðixÞ e a ak e 2 T k¼1 ð2pÞ pr Z pr 2 n 1 1 1X 2 a 6 kF k2 ak eikxT dx ¼ kf k E2l ðixÞ eilxT a 2p 2p pr T k¼1
2 ~ a Þf ðt lTÞ ¼ ðD D T
where
1
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ffi u n u 1 Z pr 1X a ilxT t El ¼ a ak eikxT dx: ðixÞ e 2p pr T k¼1
ð6Þ
~ a Hence, we have D D T ¼ El since the equality is achievable for certain f in the Cauchy–Schwarz inequality. Consequently, the discrepancy between the fractional differentiator and its discretization is measured by El . The filter coefficients ak are then obtained through minimizing El , or equivalently E2l , realizing the fact that E2l is in a quadratic form. This observation is useful for quantifying the error as in the performance analysis discussed later. The integral in (6) can be rewritten with a substitution w ¼ xT, such that 2
El ¼ ð2pT
2aþ1 1
Þ
2 Z ps n X a ilw ikw ak e dw ðiwÞ e ps k¼1
ð7Þ
where s ¼ rT known as time-bandwidth product. Since E2l is a quadratic form, there exists a set of optimal filter coefficients
165
Y. Wu / Applied Mathematics and Computation 247 (2014) 162–168
that minimize E2l satisfying the system of equations, @E2l =@aj ¼ 0; j ¼ 1; 2; . . . ; n. The result is stated in Theorem 1. ~ a that discretize the fractional order differentiator Da satisfy the system of Theorem 1. The optimal filter coefficients fak g of D T equations Ax ¼ b as a result of minimizing the performance index E2l , where
2 6 6 A¼6 6 4
sð0Þ
sðsÞ
sðsÞ
sð0Þ
. . . sððn 1ÞsÞ
7 . . . sððn 2ÞsÞ 7 7; 7 ... 5
...
sððn 1ÞsÞ sððn 2ÞsÞ . . . with sðtÞ ¼ sinðptÞ=ðptÞ , bk ¼ 2ps 1
3
R ps
ps
sð0Þ
3 b1 6 7 6 b2 7 6 7 b ¼ 6. 7 6. 7 4. 5 2
ð8Þ
bn
wa cosððk lÞw þ ap=2Þdw;
k ¼ 1; 2; . . . ; n, and x ¼ ½ a1
a2
...
a n T .
Proof. We first rewrite E2l into
E2l ¼ ð2pT 2aþ1 Þ
1
#" # Z ps " n n X X a a ðiwÞ eilw ak eikw ðiwÞ eilw ak eikw dw: ps
k¼1
k¼1
It is known that the optimal coefficients fak g that minimize (7) satisfy the system of equations @E2l =@aj ¼ 0; j ¼ 1; 2; . . . ; n, where
@E2l @aj
¼ ð2pT 2aþ1 Þ
1
" # " #) Z ps ( n n X X a a eijw ðiwÞ eilw ak eikw eijw ðiwÞ eilw ak eikw dw ps
k¼1
k¼1
Z ps X Z ps n 1 1 a a ¼ ð2pT 2aþ1 Þ ak ðeiðkjÞw þ eiðkjÞw Þdw ð2pT 2aþ1 Þ ½ðiwÞ eiðjlÞw þ ðiwÞ eiðjlÞw dw: ps k¼1
ps
The first integral in @E2l =@aj can be evaluated quickly as
Z ps n n X X ak ðeiðkjÞw þ eiðkjÞw Þdw ¼ 4 ak sinððk jÞpsÞ=ðk jÞ; k¼1
ps
k¼1
while the second integral is rewritten into a real form
Z ps
wa ei½ðjlÞwþpa=2 þ ei½ðjlÞwþpa=2 dw ¼ 2
ps
Z ps
wa cos½ðj lÞw þ pa=2dw:
ps
By setting @E2l =@aj ¼ 0, one has
Z ps n X 1 ak sinððk jÞpsÞ=ðpsðk jÞÞ ¼ wa cosððj lÞw þ pa=2Þdw; 2ps ps k¼1
j ¼ 1; 2; . . . ; n;
The integral for bk has to be evaluated numerically if a is non-integer. Numerical quadrature, such as the Romberg method and Gaussian quadrature, produces accurate result. Since the matrix A in (8) is a symmetric Toeplitz matrix, the system can be solved iteratively by using the well-known Levinson’s method [15]. This method is particularly effective when the matrix A is ill-conditioned. It is observed from (7) that El ! 0 as s ! 0 (or T ! 0Þ. This implies that the discrete fractional differentiator (4) yields more accurate approximation of a fractional derivative for smaller sampling interval T or a combined effect of the time-bandwidth product. This in turn also makes the matrix A ill-conditioned. Going through the derivations, it is readily seen that the filter coefficients in (4) are independent of function values at the sampling points. They only depend on the time-bandwidth product and the point of approximation, i,e, l, which are usually known a priori. The coefficients are fixed throughout the convolution process as the fractional derivative is computed, a standard feature of a digital filter such as the one in (4). It is clear that the discrete fractional differentiator (4) is based on uniform sampling. For the case of non-uniform sampling, the derivation is similar but with some modifications. First off, the FIR filter based on non-uniform sampling is given by
~ a f ðt t e Þ ¼ D
n X ak f ðt tk Þ; k¼1
~ a where t i < tiþ1 ; i ¼ 1; 2; . . . ; n 1. It can be shown that the error of approximation is D D ¼ Ee , where
ð9Þ
166
Y. Wu / Applied Mathematics and Computation 247 (2014) 162–168
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ffi u n u 1 Z pr X a Ee ¼ t ak eixtk dx ðixÞ eixte 2p pr k¼1
ð10Þ
The filter coefficients ak are obtained from minimizing (10), and they satisfy the system of equations Ax ¼ b where
2
sððt 2 t1 ÞrÞ . . . sððtn t1 ÞrÞ
sð0Þ
6 sððt t ÞrÞ 1 2 6 A¼6 4 ...
2
3
. . . sððtn t2 ÞrÞ 7 7 7; 5 ...
sð0Þ
sððt n t 1 ÞrÞ sððt n t 2 ÞrÞ . . .
sð0Þ
b1
3
6b 7 6 27 7 b¼6 6 .. 7 4. 5
ð11Þ
bn
R pr
1 where bk ¼ 2pr xa cosððtk te Þx þ ap=2Þdx; k ¼ 1; 2; . . . ; n. It is noted that the matrix A in (11) is symmetric but no pr longer Toeplitz. An interesting case of non-uniform sampling is the so-called interlaced sampling. The non-uniform sampling itself can be cumbersome as one has to keep track with the sampling points; however, if such sampling data are available, the proposed model (9) is ready to generate the fractional derivatives at the designated points, whereas most methods in the literature only work with uniform sampling. It is worth noting that the above derivations still hold for the case where a is a positive integer, corresponding to the conventional differentiator. Therefore, the proposed model for discretization of the fractional differentiator is also applicable to the integral order differentiator. Before presenting some numerical results, we analyze the performance of the proposed models (4) and (9). We derive an alternative expression on the error of approximation E2l , followed by an explicit formula on the minimized E2l . As mentioned earlier, the integral on E2l in (7) is in a quadratic form. To this end, we rewrite (7) into
2
El ¼ ð2pT ¼ ð2pT
2aþ1 1
Þ
Z ps " ps
2aþ1 1
Þ
Z ps
"
a ilw
ðiwÞ e
n X ak eikw
#"
a il w
ðiwÞ e
k¼1
# n X ikw dw ak e k¼1
# X XX a a a iðklÞw a iðklÞw iðjkÞw ðiwÞ ðiwÞ ak ððiwÞ e þ ðiwÞ e Þþ ak aj e dw
ps
k
k
j
After evaluating the integral over each term, we have
X XX E2l ¼ 2ðpsÞ2aþ1 =ð2a þ 1Þ 4ps ak bk þ 2ps ak aj Akj k
k
j
where bk and Akj , entries of the matrix A, are defined in Theorem 1. Consequently, the error E2l is written as
E2l ¼
r T 2a
where u ¼ ½ 1
" H¼
uT H u
ð12Þ
T
xT and x is the vector of filter coefficients ak , and the symmetric matrix H
h b b
T
#
A
ð13Þ
where h ¼ ðpsÞ2a =ð2a þ 1Þ . We observe that the matrix H is positive definite because the quadratic form (7) is positive definite. It is well known that, given a positive definite matrix H and a vector x, the following inequality holds [16], xT Hx P k0 kxk22 , and the equality is achieved if x ¼ v , where k0 is the smallest eigenvalue of H, andv is the corresponding eigenvector. Equipped with this result, T let v ¼ ½ v 0 v 1 . . . v n and kv k2 ¼ 1, the error of approximation (12) is minimized by choosing u ¼ v =v 0 and the minimal value is
min E2l ¼
k0 r T 2a v 20
ð14Þ
The implication of (14) is twofold: it first gives a quantitative measure of the error of approximation on the fractional derivative via (4); secondly, the optimal set of filter coefficients can be obtained from the eigenvector v above with ak ¼ v k =v 0 ; k ¼ 1; 2; . . . ; n. Hence, it provides an alternative approach for computing the filter coefficients. It is known that the matrix H is ill-conditioned since the Toeplitz matrix A is ill-conditioned with small time-bandwidth product. This implies that the smallest eigenvalue of H is close to zero, hence indicating accurate approximation of the fractional derivative due to (14). Numerical experiments are performed to demonstrate the accuracy of the proposed discrete fractional differentiator design and the results are consistent. We construct the Paley–Wiener functions and their fractional derivatives via Definition 1 and (5), respectively. For instance, the bandlimited function used in Fig. 1 is generated from the frequency domain with R pr FðwÞ ¼ w2 , i.e. f ðtÞ ¼ 1=ð2pÞ pr FðwÞeiwt dw, where r ¼ 6. We compute the precise fractional derivative, in this case the order is a ¼ 5=3, through (5). The sampling interval is T ¼ 0:05 and the length of the filter is n ¼ 6. The functional samples are con-
167
Y. Wu / Applied Mathematics and Computation 247 (2014) 162–168 6 Exact Approximation
4
2
0
-2
-4
-6
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
Fig. 1. Numerical approximation of fractional derivatives of a bandlimited function.
3 Exact Approximation
2.5
2
1.5
1
0.5
0
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5 1=2
Fig. 2. Numerical approximation of half-derivative of the function x with D
5
pffiffiffi pffiffiffiffi x ¼ 2 x= p.
voluted with the filter coefficients according to (4) to extrapolate the fractional derivative, i.e. l ¼ 0 in (4). The result is shown in Fig. 1. The point-wise error of approximation of the fractional derivative shown in Fig. 1 is between 104 and 103 . We also tested an interesting case with polynomial functions. Since polynomial functions have zero content in the frequency domain, we let the bandwidth parameter r approach zero when computing the filter coefficients of (4). Meanwhile, it is known that Da xn ¼ Cðn þ 1Þxna =Cðn a þ 1Þ as the exact fractional derivative of xn . When approximating the half-derivative of x, see Fig. 2, we choose T ¼ 0:1 and n ¼ 6 in (4). The approximation matches the exact half-derivative function well as shown in Fig. 2. Meanwhile, the upper bound on the error of approximation in this example is 0.082. The Grünwald–Letnikov methods for approximating fractional derivatives are frequently cited in the literature, see [10,11]. We performed comparison tests between the proposed approximation method and the Grünwald–Letnikov methods. The results show that our method yielded more accurate approximations than the Grünwald–Letnikov methods, and our method is more efficient as well, when we apply those methods to Paley–Wiener functions. A possible explanation is that the Grünwald–Letnikov methods are obtained from a generalization of the conventional numerical differentiation formulas, such as backward difference and central difference, which are conveniently derived from the Taylor series. Meanwhile, the truncation error of the Taylor series, which depends on the magnitude of the derivative, dictates the error of approximation. On the contrary, the proposed method in this paper is derived from minimizing the error of approximation directly. Hence, it is understood as the most accurate method in the sense of least squares optimization. The following table demonstrates the efficiency of our method versus Grünwald–Letnikov methods. The testing function is again obtained from the Fourier transform with FðwÞ ¼ w3 and bandwidth r ¼ 10, and its fractional derivatives are computed from (5) with
Table 1 Comparison between the proposed method and Grünwald–Letnikov methods. Error/methods
G–L backward difference
G–L central difference
Proposed
102, T = 0.1 103, T = 0.05 104, T = 0.02
16 21 28
13 16 22
7 8 10
168
Y. Wu / Applied Mathematics and Computation 247 (2014) 162–168
a ¼ 1=3. We fix the error of approximation at different levels and compare the number of samples needed to meet the error (see Table 1). In this paper, we present numerical methods for computing the fractional derivatives of functions from the Paley–Wiener space. These functions are also known as bandlimited signals in signal processing. Therefore, the proposed digital filters are applicable for generating fractional order sequences from signal samples, which can be used in signal encryption in secure communication. The filter is efficient as its coefficients are invariant with respect to functional samples and the computational process is a simple discrete convolution. The accuracy of the discrete model is demonstrated with a closed form formula on the minimized error of approximation as well as numerical experiments. All these make it possible to reconstruct a smooth fractional derivative of a rather complicated function, whereas its fractional derivative is impossible to obtain analytically. Acknowledgment The author would like to thank the anonymous reviewers for their valuable suggestions that helped improve the presentation of the original manuscript. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.amc.2014.08.089. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16]
K.B. Oldham, J. Spanier, The Fractional Calculus, Academic Press, New York and London, 1974. I. Podlubny, Geometric and physical interpretation of fractional integration and fractional differentiation, Fract. Calc. Appl. Anal. 5 (2002) 367–386. L. Debnath, Recent applications of fractional calculus to science and engineering, Int. J. Math. Math. Sci. 54 (2003) 3413–3442. L. Wang, Y. Ma, Z. Meng, Haar wavelet method for solving fractional partial differential equations numerically, Appl. Math. Comput. 227 (2014) 66–76. S. Wheatcraft, M. Meerschaert, Fractional conservation of mass, Adv. Water Resour. 31 (2008) 1377–1381. C. Tseng, S. Pei, S. Hsia, Computation of fractional derivatives using Fourier transform and digital FIR differentiator, Signal Process. 80 (2000) 151–159. A. Razminia, D.F.M. Torres, Control of a novel chaotic fractional order system using a state feedback technique, Mechatronics 23 (2013) 755–763. A. Bueno-Orovio, D. Kay, V. Grau, B. Rodriguez, K. Burrage, Fraction diffusion models cardiac electrical propagation reveal structural heterogeneity effects on dispersion of repolarization, J. R. Soc. Interface 11 (2014). M. Caputo, Linear model of dissipation whose Q is almost frequency independent-II, Geophys. J. R. Astr. Soc. 13 (1967) 529–539. E. Sousa, How to approximate the fractional derivative of order 1 < a < 2, Int. J. Bifur. Chaos 22 (2012). D. Baleanu, O. Defterli, O.P. Agrawal, A central difference numerical scheme for fractional optimal control problems, J. Vibr. Control 15 (2009) 583–597. K. Diethelm, N.J. Ford, A.D. Freed, Y. Luchko, Algorithms for the fractional calculus: a selection of numerical methods, Comput. Methods Appl. Mech. Eng. 194 (2005) 743–773. C. Lubich, Discretized fractional calculus, SIAM J. Math. Anal. 17 (1986) 704–719. J.R. Higgins, Sampling Theory in Fourier and Signal Analysis Foundations, Cearendon Press, Oxford, 1996. W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes, The Art of Scientific Computing, Cambridge University Press, 1986. G.H. Golub, C.F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore and London, 1989.