 Applied Mathematics, 2011, 2, 914-917 doi:10.4236/am.2011.27124 Published Online July 2011 (http://www.SciRP.org/journal/am) Copyright © 2011 SciRes. AM Fourier Truncation Method for Fractional Numerical Differentiation Ailin Qian, Jianfeng Mao Department of Mathematics and Statistics, Xianning University, Xianning, China E-mail: junren751113@126.com Received February 13, 2011; revised May 31, 2011; accepted June 3, 2011 Abstract We consider an ill-posed problem-fractional numerical differentiation with a new method. We propose Fou-rier truncation method to compute fractional numerical derivatives. A Hölder-type stability estimate is ob-tained. A numerical implementation is described. Numerical examples show that the proposed method is ef-fective and stable. Keywords: Inverse Problems, Ill-Posed Problem, Fractional Numerical Differentiation, Fourier Regularization, Error Estimate1. Introduction In this paper we shall consider the problem of stably calculating the fractional derivative of a function f given in , pLR  d1d1dxfttDfx xxt , (1.1) for . Such problem is frequently encountered in many practical contexts [1-3]. It is well known that is a formal solution of the Abel integral equa-tion. 0,1Dfx  11d,.xutIuxtf xxtx   (1.2) That is very important in various areas of mechanics, spectroscopy, computerized tomography . The process of numerical fractional differentiation is well known to be an ill-posed problem [1-3], and it has been discussed by many authors, and a large number of different solu-tion methods have been proposed. For references we refer the reader to [1-5] and references therein. Finite difference approaches for numerical differentiation have been used, for example, in [6-9]. However, these ap-proaches require knowledge of a bound of the second or third derivatives of the function under consideration that are not always available. Furthermore, there exist infi-nitely many functions that do not have bounded sec- ond derivatives at all. The same situation occurs also in numerical fractional differentiation ; one requires a high smoothness of the functions under consideration that does not always exist. In the present paper, as an alterna-tive way of dealing with fractional numerical differentia-tion, we introduce a new regularization method, i.e., Fourier truncation. We simply analyze the cause of ill- posedness of fractional numerical differentiation and then propose the method. The idea of Fourier truncation method is very simple and natural: since the ill-posed-ness of fractional numerical differentiation is caused by the high frequency components, we cut off them. Actu-ally, such a similar idea of solving numerical differentia-tion has occurred in [10,11]. We can easily find this fact by studying [10,11] in frequency space. However, Fou-rier truncation method is more direct, natural and simple. 2. Regularization in the Fourier Domain In this section we simply analyze the ill-posedness of fractional numerical differentiation and discuss how to stabilize the numerical derivaties. We set a function pfxHR which is a Soblev space, p. Let f be the Fourier transform of f, i.e.,  1ˆd.2πixffxe x  (2.1) Now we consider the -th order derivative of func-tionf. Taking the Fourier transform, we have A. L. QIAN ET AL.915  ^ˆ,fx if (2.2) i.e.,   1ˆd.2πixfx ife (2.3) From the right hand side of (2.2) or (2.3), we know that i can be seen as an amplification factor of ˆf. Therefore, when we consider our problem in LR2, the exact data function ˆf must decay rap-idly as . But in the applications the input data fx can only be measured and never be exact. We assume, say that, the measured data function  2fxLR satisfies ,ff (2.4) where  denotes - norm, the constant 2L0 represents a noise level. Thus, if we try to obtain frac-tional numerical derivatives, high frequency components in the error are magnified and can destroy the solution. In this sense it is impossible to solve the problem using classical numerical methods and requires special tech-niques to be employed. In the following, we will propose our regularization strategy to deal with the ill-posed problem. However, before doing that, we impose a priori bound on the input data (this is necessary in solving ill-posed problems), i.e., ,pfEp, (2.5) where is a constant, 0Epdenotes the norm in Soblev space pHR defined by 1222ˆ:1 dppff  (2.6) Note that (2.2) or (2.3), since the ill-posedness of the problem is caused by the high frequency components, a natural way to stabilize the problem is to eliminate all high frequencies and instead consider (2.3) only for max, where max will be seen that it exists as a regularization parameter. Then we get a regularized frac-tional derivative  max1ˆ:d.2πixRfxi fe  (2.7) where max is the characteristic function of the interval max max,, i.e., max maxmaxmax.1, ,0,  In the following sections we will derive an error esti-mate for the approximate derivative (2.7) and discuss how to compute it numerically. 3. Error Estimate In this section, we derive a bound on the difference be-tween the derivatives (2.3) and (2.7). We assume that we have an a priori bound on the exact input data. pfE (See (2.5)). The relation between any two regularized solutions (2.7) is given by the following lemma. Lemma 3.1. Suppose that we have two regularized derivatives 1Rf x1 and defined by (2.7) with data 2Rf xf and 2f, satisfying 12 .ff If we select 1max ,pE 1p, then we get the error bound 1()12 .ppRf RfE  (3.1) Proof. From the Parseval relation we get maxmaxmaxmax2() ()1221221222122222max2max 12dd1Rf RfRf Rfiffffff ff  Using 1maxpE and12ff , we obtain 112 ,.ppRfRfE p  From Lemma 3.1 we see that the derivative defined by (2.7) depends continuously on the input data f. Next, we will investigate the difference between the derivatives (2.3) and (2.7) with the same exact f. Lemma 3.2. Let fx and be the de-rivatives (2.3) and (2.7) with the same exact data Rf xf, and let 1max ,1pEp . Suppose that pfE. Then  1.ppfRf E  (3.2) Proof. As in Lemma 3.1 we start with the Parseval re-lation, and using the fact that the derivatives coincide for max max, we get Copyright © 2011 SciRes. AM A. L. QIAN ET AL. 916   maxmaxmaxmax2222222222222maxˆdˆdˆ1d1ˆˆ1dppppp.pfRfifffff  Now we use the bound pfE,and as before we have 1maxpEwhich leads to the error bound  1.ppfRf E  Now we are ready to formulate the main result of this section: Theorem 3.3. Suppose that fx is given by (2.3) with exact data f and that is given by Rfx(2.7) with measured data fx. If we have a bound pfE, and the measured function fx satisfies ff, and if we choose 1maxpE where 1p, then we get the error bound 12ppfRf E.  (3.3) Proof. Let be the derivative defined by (2.7) with exact data Rf xf. Then by using the triangle inequality and the two previous Lemmas we get 12.ppfRf fRfRf RfE   From Theorem 3.3, we find that (2.7) is an approxima-tion of the exact derivative xf. The approximation error depends continuously on the measurement error. Remark 3.4. In our application pf is usually un-known, therefore we have no the exact a priori bound E. However, if we select 1maxpCC, where is a positive constant, we can also obtain 110, when 0,,pfRf Cp  where the constant C depends on ,,ppf. This choice is helpful in our realistic computation. 4. Numerical Implement (2.7) numerically. Given a vector F containing samples from f on an equidistant grid1,nthen the unique 0jjxtrigon etric polynomial interpofom lating  on the grid can be written 1221ˆˆ(), 2π,2πjnixjjnjfxfe k (4.1) where the sequence122ˆnnjjfare the discrete Fourier coefficients, and we have assumed thatis even. The n pudiscrete Fourier coefficients can be comted by taking the FFT of the vector F. Thus we can approximate the derivative operator (2.7s follow: ) aHRFL L,F where is the Fourier matrix , and is a diago-L tri nal max corresponding to differentiation of trigono-metric interpolant, but where the frequency components with maxjare explicitly set to zero. Thus the di-agonas of l element are maxmax, ,0, ,jjijji (4.2) where j anare defined as in (4.1). The product of L and a vector c be computed using the Fast Fourier Transform (FFT) which leads to an efficient way to compute the derivative (2.7). When using the FFT algorithm we im-plicitly assume that the vector F represents a periodic function. This is not realistic in our application, and thus we need to modify the algorithm. A discussion on how to make the function “periodic” can be found in . 5. Numerical Examples For verifying the effect of the proposed algorithm, a smooth function and a non-smooth function will be tested in various cases. In the numerical experiment, we always fix 01x. For an exact data function fx, its discrete noiion is sy versff randn size f, (5.1) where  21211,,,1,2,11.NjNjjljjf,,fxfxxjNNff fx fxN  The function “randn” generates arrays of random numbers whose eare normally distributed with lements Here we discuss how to compute the regularized solution Copyright © 2011 SciRes. AM A. L. QIAN ET AL. Copyright © 2011 SciRes. AM 917mean 0, variance 21, and standard deviation 1, “(())randn sizef s an array of random s thize as “return entrieat is the same sf. In our computations, we al-ways take 4097N (Ife take N = 100,513,1025,, we can also satisfactory result.)The derivative errors are measured by the weighted 2L-norms defined as follows:  in aw obta A. K. Louis, “Inverse und Schlecht Gesteellte Problem,” Teubner Verlag, Wiesbaden, 1989.  D. A. Murio, “Automatic Numerical Differentiation by Discrete Molification,” Computers and Mathematics with Applications, Vol. 13, No. 4, 1987, pp. 381-386. doi:10.1016/0898-1221(87)90006-X  211NjNtion pa. (5.2) D. A. Murio and L. Guo, “Discrete Stability Analysis of Molification Method for Numerical Differentiation,” Journal of Computational Applied Mathematics, Vol. 19, No. 6, 1990, pp. 15-25. jjx RxeterERfThe regularizaframf  I. Daubechies, “Ten Lectures on Wavelets (CBMS-NSF Regional Conference Series in Applied Mathematics),” SIAM: Society for Industrial and Applied Mathematics, Philadelphia, 1992. max was selected ac-4, , cording to Remark 3. i.e.1maxpc. The erro norms in r es- timates of Section 3 contain2LR,and the  J. N. Lyness, “Has Numerical Differentiation a Future?” Utilitas Mathematics, Winnipeg, 1978. numerical experiments are performed witite inter-val. However, the method for selecting maxh a fin works well,also in the case where the norms are comed for a finite interval. put J. Oliver, “An Algorithm for Numerical Differentiation of a Funcion of One Real Variable,” Journal of Computa-tional Applied Mathematics, Vol. 6, No. 2, 1980, pp. 145-160. doi:10.1016/0771-050X(80)90008-X  T. Strom and J. N. Lyness, “On Numerical Differentia-tion,” BIT Numerical Mathematics, Vol. 15, No. 3, 1975, pp. 314-322. doi:10.1007/BF01933664 6. AcknowleThe author wants to7. References  J. Baumeistdgements exer, “Spress table Sohis thanks to the refelution of Inverse Problems,” F. ral Equree for ations,” D. N. Hao, “A Molification Method for Ill-Posed Prob-lems,” Numerische Mathematik, Vol. 68. No. 4, 1994, pp. 469-506. doi:10.1007/s002110050073 many valuable comments. This work is supported by Educational Commission of Hubei Province of China (Q20102804,T201009).  D. A. Murio, C. E. Mejia and S. Zhan, “Discrete Molifi-cation and Automatic Numerical Differentiation,” Com-puters and Mathematics Application, Vol. 35, No. 5, 1998, pp. 1-16. doi:10.1016/S0898-1221(98)00001-7  L. Elden, F. Berntsson and T. Reginska, “Wavelet and Fourier Methods For Solving the Sideways Heat Equa-tion,” SIAM Journal on Scientific Computing, Vol. 21, No. 6, 2000, pp. 2187-2205. Vieweg and Sohn, Braunschweig, 1987.  R. Gorenflo and S. Vessela, “Abel Integ Springer-Verlag, Berlin, Heidelberg, New York, 1991.