Abstract

Metzler et al. introduced a fractional Fokker-Planck equation (FFPE) describing a subdiffusive behavior of a particle under the combined influence of external nonlinear force field and a Boltzmann thermal heat bath. In this paper, we present an interval Shannon wavelet numerical method for the FFPE. In this method, a new concept named “dynamic interval wavelet” is proposed to solve the problem that the numerical solution of the fractional PDE is usually sensitive to boundary conditions. Comparing with the traditional wavelet defined in the interval, the Newton interpolator is employed instead of the Lagrange interpolation operator, so, the extrapolation points in the interval wavelet can be chosen dynamically to restrict the boundary effect without increase of the calculation amount. In order to avoid unlimited increasing of the extrapolation points, both the error tolerance and the condition number are taken as indicators for the dynamic choice of the extrapolation points. Then, combining with the finite difference technology, a new numerical method for the time fractional partial differential equation is constructed. A simple Fokker-Planck equation is taken as an example to illustrate the effectiveness by comparing with the Grunwald-Letnikov central difference approximation (GL-CDA).

1. Introduction

Due to the fact that 1/ signal gains the increasing interests in the field of biomedical signal processing and engineering systems [1], the differential equations of fractional order appear more and more frequently in various research areas and engineering applications [2, 3]. As a matter of fact, the applications of fractional differential equations and their corresponding time series have been developed in various fields of sciences and technologies [4, 5] in recent years, ranging from computer science to physics [6, 7]. An effective and easy-to-use method for solving such equations is needed. However, known methods have certain disadvantages. Methods, described in detail in [3] for fractional differential equations of rational order, do not work in the case of arbitrary real order. On the other hand, there is an iteration method described in [8], which allows solution of fractional differential equations of arbitrary real order but it works effectively only for relatively simple equations, in addition to the series method. Up to now, most studies on the numerical methods for the fractional PDEs concentrate on the finite difference methods. Li [9] proposed an analytical method taking the fractal time series as the solution to a differential equation of fractional order or a response of a fractional system or a fractional filter driven with a white noise in the domain of stochastic processes and gave the exact solution of impulse response to a class of fractional oscillators [10]. According to this idea, Li and his coresearchers solved many problems in science and technology [1114]. In addition, Wavelet numerical method is another way to get the solution of the fractional PDEs. In fact, the wavelet transform theory has been widely used in numerical analysis such as PDEs-based image processing [1517], option pricing model [18], integrodifferential operators [1923], and other nonlinear PDEs [2428]. The wavelet functions possess many excellent numerical properties, such as orthogonality, interpolation, smoothness, and compact support, which are helpful in improving numerical accuracy and efficiency. In recent decades, many wavelets which have compact support, smoothness, and other properties have been constructed. Among these wavelets, Shannon wavelet is paid little attention in applications as it does not possess compact support property although it possesses orthogonality, smooth continuity, and analytical expression. Cattani studied the properties of the Shannon wavelet function, which possesses many advantages such as orthogonality, continuous and differentiable. It also has the advantage over the Hermite DAF in that it is an interpolating function, producing matrix equations that have the potential to be relatively sparse. In addition, the second order approximation of a C2-function, based on Shannon wavelet functions, is given [29]. The approximation is compared with the wavelet reconstruction formula and the error of approximation is explicitly computed [30].

A perceived disadvantage of the Shannon scaling function is that it tends to zero quite slowly as. A direct consequence of this is that when calculating the derivatives, a large number of the nodal values will contribute significantly. It is for this reason that Hoffman et al. [31] have suggested using the Shannon-Gabor wavelet, which introduces the Gaussian window function to improve the compact support property of Shannon wavelet function in required precision range. However, the presence of the Gaussian window destroys the orthogonal properties possessed by the Shannon wavelet, effectively worsening the approximation to a Dirac delta function.

Comparing with the common PDEs, the solutions of the fractional PDEs are more sensitive to the boundary condition. Using the wavelet transform defined in infinite domain to solve the engineering problems in finite interval, the wavelet transform coefficients at the boundary are usually very large. It will bring server boundary effect which affects the calculation accuracy and efficiency. Vasilyev and Paolucci [32] construct an interval wavelet using external wavelets, which can decrease the boundary effect to some extent. Based on the same principle, a more general construction method for the interval interpolation wavelet [33, 34] was given in the framework of generalized variational principle and has been widely used in many areas [3537]. But the choice of parameter(that is the amount of the external collocation points) was not discussed in detail. It just points out that the value ofshould be taken between 1 and 3 based on experience. In fact, the value ofdepends on the smoothness and derivative of the approximated function at boundary points. That is, if the approximated function is the solution of the diffusion PDEs with respect to the time parameter, the value ofshould be taken dynamically. In addition, we should take into account that the impact of the external collocation points to the condition number of the system of the discretized algebraic equations. So, it is necessary to construct a dynamic interval wavelet in solving the PDEs with dynamic boundary conditions such as the fractional PDEs.

In this paper, a dynamic interval Shannon wavelet collocation method for the fractional FPDs is proposed. In this method, the relation between the parameterand the wavelet approximation error was discussed based on the interpolation error theory, and an adaptive choice procedure onwas constructed. Therefore, the so-called dynamic interval Shannon wavelet is constructed. Next, based on the Grünwald-Letnikov definition of the fractional order derivative, we construct a Shannon wavelet numerical method for the fraction Fokker-Planck equation.

2. Fractional Fokker-Planck Equation

The fractional Fokker-Planck equation has been used in many physical transport problems which take place under the influence of an external force field [2, 38].

In the presence of an external force field, the evolution of a test particle is usually described in terms of the Fokker-Planck equation (FPE)

which defines the probabilityof finding the particle at a certain position at a given time.denotes the mass of the diffusing particle,denotes the generalized diffusion coefficient with dimension, andis the generalized friction coefficient with dimension. The corresponding initial condition is and the boundary conditions are Equation (1) uses the Riemann-Liouville fractional derivative of order, defined by whereis the gamma function.

According to the properties of the Riemann-Liouville fractional derivative, it is easy to know that, if, (1) can be rewritten as follows: Metzler et al. [2] proposed three implicit approximations for solving (5) as follows.

 The Grünwald-Letnikovexpansion and the backward Euler implicit approximation (GL-BDIA) where,,, and andare positive integers. The local truncation error is.

-approximation and the central difference implicit approximation (-CDIA) The local truncation error is.

-approximation and the backward difference implicit approximation (-BDIA) The local truncation error is.

In fact, (6)–(8) are not perfect approximation as the boundary effect is not taken into account. So, it will introduce boundary effect in solving the PDEs with the Nuemann boundary conditions. It is well known that the finite difference method is equivalent to the Faber-Schauder wavelet collocation method, so the construction method of the dynamic interval wavelet introduced in this paper can also be used to deal with the boundary problem in the finite difference method.

According to the Shannon sample theory, it can improve the calculation precision that combining the Grünwald-Letnikov expansion or-approximation of the fractional derivative in (5) with the Shannon scaling function as the weight function instead of the various difference operators as follows:

3. Construction of the Interval Interpolation Wavelet

3.1. Shannon Wavelet and Shannon-Gabor Wavelet

The representation of Shannon wavelet is based upon approximating the Dirac delta function as a band-limited function and is given by and the Shannon-Gabor scaling function is defined as [17] where is the window size.

Consider a one-dimensional function,. A discrete point sequence of the variableis defined as and the corresponding discrete point sequence of the scaling function andcan be defined, respectively, as where.

The first and second order derivatives ofat the discrete point are

And the first and second order derivatives ofat the discrete point are In fact, there is no difference between the construction method of the Interval Shannon wavelet and the interval Shannon-Gabor wavelet. So, we just take one uniform symbolas the representation of the Shannon wavelet and the Shannon-Gabor wavelet in the following.

3.2. Interval Interpolation Wavelet

According to the definition of the interval wavelet, the interval interpolation basis functions can be expressed as:

where, whereis the amount of the external collocation points, the amount of discrete points in the definition domain isis the definition domain of the approximated function.

Equations (16) and (17) show that the interval wavelet is derived from the domain extension. The supplementary discrete points in the extended domain are called external points. The value of the approximated function at the external points can be obtained by Lagrange extrapolation method. Using the interval wavelet to approximate a function, the boundary effect can be left in the supplementary domain; that is, the boundary effect is eliminated in the definition domain.

According to (16) and (17), the interval wavelet approximant of the functioncan be expressed as is the given value at the discrete point. At the external points, can be obtained by extrapolation; that isSo the interval wavelet approximant ofcan be rewritten as Let then andcorrespond to the left and the right external points, respectively. They are obtained by Lagrange extrapolation using the internal collocation points near the boundary. So, the interval wavelet’s influence on the boundary effect can be attributed to Lagrange extrapolation. It should be pointed out that we did not care about the reliability of the extrapolation. The only function of the extrapolation is enlarging the definition domain of the given function which can avoid the boundary effect occurred in the domain. Therefore, we can discuss the choice ofby means of Lagrange inner-and extrapolation error polynomial as follows: Equation (23) indicates that the approximation error is both related to the smoothness and the gradient of the original function near the boundary. Setting differentcan satisfy the error tolerance.

3.3. Adaptive Interval Interpolation Wavelet

The interval interpolation wavelet is often used to solve the diffusion PDEs with Neumann boundary conditions. The smoothness and gradient of the PDE’s solution usually vary with the time parameter. If the parameteris a constant, we have to take a bigger value in order to obtain results with higher calculation precision. But the biggerusually introduces the famous Gibbs phenomena into the numerical solution, which usually results in the algorithm becoming invalid. In addition, the biggerwill bring much more calculation. To keep higher numerical precision and save calculation, the best way is to design a procedure thatcan vary with the curve’s smoothness and gradient dynamically.

In this dynamic procedure, the error estimation equation (23) can be taken as the criterion about. But in most cases, we cannot know the smoothness and the derivative’s order of the original function. This can be solved by substituting the difference coefficient for the derivative. This is coincident with the Newton interpolation equation which is equivalent with Lagrange interpolation equation. In addition, the Lagrange interpolation algorithm has no inheritance which is the key feature of Newton interpolation. So, the basis function has to be calculated repeatedly as interpolation points are added into the calculation, which increases the computation complexity greatly. In contracst to the Lagrange method, the advantage of Newton interpolation method is that the Newton divided difference form is employed, which can produce a mathematically equivalent result by using recurrence relations, which reduces the number of compute operation, especially the multiplication. So it is convenient using the Newton interpolation method to construct the dynamic procedure.

3.3.1. Newton Interpolation

The expression of Newton interpolation can be written as Substituting the Newton interpolation instead of the Lagrange interpolation into (22) can be rewritten as where

3.3.2. Relation between the Newton Interpolation Error and the Choice of

It is well known that the Newton interpolation is equivalent to the Lagrange interpolation. The corresponding error estimation can be expressed as And the simplest criterion to terminate the dynamic choice onis (is the absolute error tolerance). Obviously, it is difficult to definewhich should meet with the precision requirement of all approximated curves. In fact, the difference coefficient ) can be used directly as the criterion; that is As mentioned above, once the curves with lower order smoothness are approximated by higher order polynomial expression, the errors will become bigger on the contrary. In fact, even if theis infinite, the computational precision cannot be satisfied except by increasing computational complexity. To avoid this, we design the termination procedure of dynamic choice aboutas follows:If , then else if or , then else if or , then

3.3.3. and the Condition Number of the System of Algebraic Equations

In the field of numerical analysis, the condition number of a function with respect to an argument measures how much the output value of the function can change for a small change in the input argument. This is used to measure how sensitive a function is to changes or errors in the input and how much error in the output results from an error in the input. It is no doubt that the choice ofcan change the condition number of the system of algebraic equations discretized by the wavelet interpolation operator or the finite difference method. Therefore, the choice ofshould take the condition number into account. In fact, if the condition number, then you may lose up todigits of accuracy on top of what would be lost to the numerical method due to loss of precision from arithmetic methods [24]. According to the general rule of thumb, the choice should follow the rule as follows:

3.3.4. Relation between and Computation Complexity

The computational complexity of interpolation calculation is not proportional to the increasing points. The former is mainly up to the computation amount of and the derivative operations. Obviously, according to (5), the increase in computational complexity iswhen the number of extension pointsincreases by 1. But the computational complexity of adaptively increasing collocation points is related to the different wavelet functions. For the wavelet with compact support property such as Daubechies wavelet and Shannon wavelet, the value ofis impossible to be infinite. For Haar wavelet which has no smoothness property,can be taken as 0 at most since it need not to be extended. For Faber-Schauder wavelet,can be taken as 1 at most. For Daubechies wavelet,can be taken as different values according to the order of vanishing moments, but it must be finite. For the wavelets without compact support property,can take value dynamically, such as Shannon wavelet. The computational complexity of increasing points is mainly up to the wavelet function of itself.

4. Numerical Results and Discussion

Fractional Fokker-Planck equation is a typical fractional PDE, which is often used to describe a subdiffusive behavior of a particle under the combined influence of external nonlinear force field, and a Boltzmann thermal heat bath. This section considers the accuracy and efficiency of the proposed method for a fractional Fokker-Planck equation. Comparisons are made with results obtained with Chen’s finite difference approximations and the exact analytic solution.

It has been pointed out that the finite difference approximation formats proposed in [2] are not perfect as they do not take the boundary problems into account. In this section, we take the Grünwald-Letnikovexpansion and the central difference implicit approximation (GL-CDIA) to solve the example. That is, According to the wavelet collocation method, the fractional Fokker-Planck equation can be approximately represented as where . Let Then the system of (31) can be expressed as the matrix format: Next, we will discuss the precision of the method proposed in this paper with numerical experience. Consider the Fokker-Planck equation as follows: with the initial condition and the boundary conditions The exact analytic solution is

All the comparisons in this section are made qualitatively by comparing the calculation precision in the same time step and space mesh grid size. The first measure of error is given by which provides a measure of the accuracy of the solution near the boundary. The second measure of erroris given by which provides a general measure of the accuracy of the solution over the main body of the distribution and was often used to investigate the accuracy of the FEM.

The comparisons between the static interval Shannon-Gabor wavelets withandare showen in Figure 1. The boundary effect of the interval wavelet with(Figure 1(a)) is almost eliminated compared to (Figure 1(b)). FFPE is a 2-order PDEs with respect to, sois the necessary condition for the interval wavelet satisfying the requirement of FFPE. We also noticed that the condition number of FFPE from the Table 1 that the condition number ofincreases more rapid thanwith the increase ofand the decrease of . It has been mentioned in Section 2 that the larger condition number can decrease the calculation precision greatly. This also can be illustrated in Figure 2. The condition number in Figure 2(a) is greatly larger than in Figure 2(b), although the approximation ofis better than. The former has failed to solve FFPE obviously. In fact, this explained the reason why we construct the dynamic interval wavelet.

The numerical errors comparisons among the dynamic, static interval wavelet method and the interval finite difference method are showen in Figure 3. The result also can be illustrated in Table 2.

The robustness of the dynamic interval wavelet collocation method (DIWCM) is the best compared to the interval FDM and the static interval WCM, as it avoids both of the larger condition number and the error of the approximation simultaneity. The varied process of is showen in Table 3. It shows that the value of is fixed atafter a short time of vibration. This reflects the properties of the FFPE to some extent.

In addition, it also has to be noticed that we can get the higher precision solution with the interval finite difference method (FDM) as the amount of the collocation points decreases (Figure 4). It is well known that increasing the collocation points can impove the approximation although it can increase the condition number in FFPE. In fact, it profits from the smoothness of the solution, which would not work in solving the nonlinear problems.

All above numerical experiments are done with the Shannon-Gabor wavelet. It is well known that the presence of the Gaussian window destroys the orthogonal properties possessed by the Shannon wavelet, effectively worsening the approximation efficiency to a Dirac delta function. Comparing with the Shannon wavelet collocation method (Figure 5), the Shannon-Gabor wavelet numerical method has higher precision and more complicated calculation amount. But it is showen in Figure 5 that dynamic interpolation wavelet construction scheme can be applied to both of the Shannon-Gabor wavelet and the Shannon wavelet. As a matter of fact, the dynamic scheme is designed for the interpolation wavelet, which has no connection with certain concrete wavelet function.

5. Conclusions

In solving the fractional Fokker-Planck equations, there are two factors related to the choice of. The first factor is the condition number, which relates to the parameters , and the time step. The largercan decrease the calculation precision greatly. Another factor is the approximation of the function and its derivatives, especially near the boundary. Using the interval wavelet with constant to solve the fraction Fokker-Planck equations cannot eliminate the boundary effect completely as the condition number is sensitive to the parameter. With regard to the accuracy and time complexity of the solution in comparison with those obtained with other algorithms, the dynamic interval wavelet onconstructed in this paper is more reasonable. The numerical experiments illustrate that it is necessary to construct the dynamic interval wavelet collocation method for the fractional PDEs.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (no. 41171337) and the National High Technology Research and Development Program of China (no. 2006AA10Z235).