# Interval Shannon Wavelet Collocation Method for Fractional Fokker-Planck Equation

Authors: Shu-Li Mei; De-Hai Zhu
Journal: Advances in Mathematical Physics
 (2013)

---

## 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).


---

## Body

## 1. Introduction

Due to the fact that 1/f 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 [11–14]. 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 [15–17], option pricing model [18], integrodifferential operators [19–23], and other nonlinear PDEs [24–28]. 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|x|→∞. 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 [35–37]. But the choice of parameterL(that is the amount of the external collocation points) was not discussed in detail. It just points out that the value ofLshould be taken between 1 and 3 based on experience. In fact, the value ofLdepends 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 ofLshould 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 parameterLand the wavelet approximation error was discussed based on the interpolation error theory, and an adaptive choice procedure onLwas 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 fieldF(x)=-ν′(x), the evolution of a test particle is usually described in terms of the Fokker-Planck equation (FPE)
(1)∂u(x,t)∂t=Dt1-α0[∂∂xν′(x)mηα+Kα∂2∂x2]u(x,t),hhhhhhhhhhhhhha≤x≤b,0≤t≤T,which defines the probabilityu(x,t)of finding the particle at a certain position xat a given timet.mdenotes the mass of the diffusing particle,Kα>0denotes the generalized diffusion coefficient with dimension[Kα]=cm2sec-α, andηαis the generalized friction coefficient with dimension[ηα]=secα-2. The corresponding initial condition is
(2)u(x,0)=φ(x),a≤x≤b,
and the boundary conditions are
(3)u(a,t)=p1(t),u(b,t)=p2(t),0<t≤T.
Equation (1) uses the Riemann-Liouville fractional derivative of order1-α, defined by
(4)Dt1-α0u(x,t)=1Γ(α)∂∂t∫01u(x,η)(t-η)1-αdη,mmmmmmmmmmmmmmll0≤α<1,
whereΓ(α)is the gamma function.According to the properties of the Riemann-Liouville fractional derivative, it is easy to know that, if(x,t)∈Cx,t2,1([a,b]×[0,T]), (1) can be rewritten as follows:
(5)Dtαu(x,t)-u(x,0)t-αΓ(1-α)=[∂∂xν′(x)mηα+Kα∂2∂x2]u(x,t),mmmmmmmmmmmmmmlma≤x≤b,0≤t≤T.
Metzler et al. [2] proposed three implicit approximations for solving (5) as follows.(
1
) The Grünwald-Letnikovexpansion and the backward Euler implicit approximation (GL-BDIA)
(6)τ-α[uin+∑k=1n-1gkuin-k-∑k=0n-1gkui0]=fiuin-fi-1ui-1nh+Kαui+1n-2uin+ui-1nh2,hhhhhhhhhhhhhhhhhhhii=1,2,…,M-1ui0=φ(xi),1≤i≤Mu0n=p1(tn),uMn=p2(tn),n≥1,gk=(1-1+αk)gk-1,g0=1,fi=f(xi)=ν′(xi)mηα,
whereh=(b-a)/M,τ=T/N,xi=a+ih, and tn=nτ.MandNare positive integers. The local truncation error isO(τ+h).(
2
) L1-approximation and the central difference implicit approximation (L1-CDIA)
(7)τ-αΓ(2-α)[uin+∑k=1n-1(an-k-1-an-k)uik-an-1ui0]=fi+1ui+1n-fi-1ui-1n2h+Kαui+1n-2uin+ui-1nh2hhhhhhhhhhhhhhhhhhhhhhi=1,2,…,M-1ui0=φ(xi),1≤i≤M,u0n=p1(tn),uMn=p2(tn),n≥1ak=(k+1)1-α-k1-α,fi=f(xi)=ν′(xi)mηα.
The local truncation error isO(τ2-α+h2).(
3
) L1-approximation and the backward difference implicit approximation (L1-BDIA)
(8)τ-αΓ(2-α)[uin+∑k=1n-1(an-k-1-an-k)uik-an-1ui0]=fiuin-fi-1ui-1nh+Kαui+1n-2uin+ui-1nh2hhhhhhhhhhhhhhhhhhhhhi=1,2,…,M-1ui0=φ(xi),1≤i≤Mu0n=p1(tn),uMn=p2(tn),n≥1fi=f(xi)=ν′(xi)mηα,ak=(k+1)1-α-k1-α.
The local truncation error isO(τ2-α+h).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 orL1-approximation of the fractional derivative in (5) with the Shannon scaling function as the weight function instead of the various difference operators as follows:
(9)τ-αΓ(2-α)[uin+∑k=1n-1(an-k-1-an-k)uik-an-1ui0]=∑i=02J[fiuiw′(x)+Kαuiw′′(x)]i=1,2,…,2J,considerJisthepositiveinteger.ui0=φ(xi),1≤i≤Mu0n=p1(tn),uMn=p2(tn),n≥1fi=f(xi)=ν′(xi)mηαak=(k+1)1-α-k1-α.


## 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(10)ϕ(x)=sin(πx)πx
and the Shannon-Gabor scaling function is defined as [17]
(11)G(x)=sin(πx)πxexp(-x22σ2),σ>0,
where σ is the window size.Consider a one-dimensional functionf(x),x∈[a,b]. A discrete point sequence of the variablexis defined as
(12)xn=a+b-a2j·n,j∈Z,n=0,1,2,…,2j,
and the corresponding discrete point sequence of the scaling functionϕ(x) andG(x)can be defined, respectively, as
(13)ϕj,n(x)=ϕj(x-xn)=sin(2jπ/(b-a))(x-xn)(2jπ/(b-a))(x-xn),Gj,n(x)=Gj(x-xn)=sin(2jπ/(b-a))(x-xn)(2jπ/(b-a))(x-xn)×exp(-22j-1(x-xn)2r2(b-a)2),
wherer=2jσ/(b-a).The first and second order derivatives ofϕj(x-xn)at the discrete pointxk are
(14)ϕj′(xk-xn)={0,k=n,2jcos[π(k-n)](k-n)(b-a),k≠n,ϕj′′(xk-xn)={-π23((b-a)/2j)2,k=n,-2cos[π(k-n)]((b-a)/2j)2(k-n)2,k≠n.And the first and second order derivatives ofGj(x-xn)at the discrete pointxk are(15)Gj′(xk-xn)={0,k=n2jcos[π(k-n)]exp[-(k-n)2/2r2](k-n)(b-a),k≠nGj′′(xk-xn)={-3+π2r23r2((b-a)/2j)2,k=n-2cos[π(k-n)]exp[-(k-n)2/2r2]((b-a)/2j)2[1(k-n)2+1r2],k≠n. 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 symbolw(x)as 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:(16)wjk(x)={ϕ(2jx-k)+∑n=-L+1-1ankϕ(2jx-n),k=0,…,Lϕ(2jx-k),k=L+1,…,2j-L-1ϕ(2jx-k)+∑n=2j+12j+L-1bnkϕ(2jx-n),k=2j-L,…,2j,where,(17)ank=∏i=L-1,i≠k-1xj,n-xj,ixj,k-xj,i,bnk=∏i=2j+1,i≠k2j+1+Lxj,n-xj,ixj,k-xj,ixj,k=kxmax-xmin2j,k∈ℤ,
whereLis the amount of the external collocation points, the amount of discrete points in the definition domain is2j+1(j∈ℤ),and[xmin,xmax]is 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 functionf(x)x∈[xmin,xmax]can be expressed as
(18)fj(x)=∑fj(xn)wj(2jx-n),xn=xmin+nxmax-xmin2j.fj(xn) is the given value at the discrete pointxn. At the external points,fj(xn) can be obtained by extrapolation; that is(19)fj(xn)={∑k=0L-1(fj(xk)∏i=0,i≠kL-1xn-xixk-xi),n=-1,…,-L∑k=2j-L+12j(fj(xk)∏i=2j-L+1,k≠i2jxn-xixk-xi),n=2j+1,…,2j+L.So the interval wavelet approximant off(x)can be rewritten as
(20)fj(x)=∑n=-L-1(∑k=0L-1fj(xk)∏i=0L-1xn-xixk-xi)ω(2jx-n)+∑n=02jfj(xk)ω(2jx-n)+∑n=2j+12j+L(∑k=2j-L2jfj(xk)∏i=2j-L2jxn-xixk-xi)ω(2jx-n).
Let
(21)LSL(xn)=∑k=0L-1fj(xk)∏i=0L-1xn-xixk-xi,LEL(xn)=∑k=2j-L2jfj(xk)∏i=2j-L2jxn-xixk-xi,
then
(22)fj(x)=∑n=-L-1LSL(xn)ω(2jx-n)+∑n=02jfj(xk)ω(2jx-n)+∑n=2j+12j+LLEL(xn)ω(2jx-n).LSL(xn)andLEL(xn)correspond 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 ofLby means of Lagrange inner-and extrapolation error polynomial as follows:
(23)RL(x)=f(L+1)(ξ)(L+1)!∏i=0L(x-xi),forsomeξbetweenhhhhhhhhhhhhhhhhhhhhhhhhhhx,x0,…,xL.
Equation (23) indicates that the approximation error is both related to the smoothness and the gradient of the original function near the boundary. Setting differentLcan 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 parameterLis a constant, we have to take a bigger value in order to obtain results with higher calculation precision. But the biggerLusually introduces the famous Gibbs phenomena into the numerical solution, which usually results in the algorithm becoming invalid. In addition, the biggerLwill bring much more calculation. To keep higher numerical precision and save calculation, the best way is to design a procedure thatLcan vary with the curve’s smoothness and gradient dynamically.In this dynamic procedure, the error estimation equation (23) can be taken as the criterion aboutL. 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(24)Nn(x)=f(x0)+(x-x0)f(x0,x1)+(x-x0)(x-x1)f(x0,x1,x2)+⋯+(x-x0)(x-x1)⋯(x-xn-1)×f(x0,x1,…,xn).
Substituting the Newton interpolation instead of the Lagrange interpolation into (22) can be rewritten as
(25)fj(x)=∑n=-L-1(NSL(xn))ω(2jx-n)+∑n=02jfj(xn)ω(2jx-n)+∑n=2j+12j+L(NEL(xn))ω(2jx-n),
where
(26)NSL(xn)=f(x0)+(xn-x0)f(x0,x1)+(xn-x0)(xn-x1)f(x0,x1,x2)+⋯+(xn-x0)(xn-x1)⋯(xn-xL-1)×f(x0,x1,…,xL),NSR(xn)=f(x2j)+(xn-x2j)f(x2j,x2j-1)+(xn-x2j)(xn-x2j-1)×f(x2j,x2j-1,x2j-2)+⋯+(xn-x2j)(xn-x2j-1)⋯(xn-x2j-L)×f(x2j,x2j-1,…,x2j-L).


#### 3.3.2. Relation between the Newton Interpolation Error and the Choice ofL

It is well known that the Newton interpolation is equivalent to the Lagrange interpolation. The corresponding error estimation can be expressed as(27)Rn(x)=(x-x0)(x-x1)⋯(x-xn)f(x,x0,…,xn).
And the simplest criterion to terminate the dynamic choice onLis|Rn(x)|≤Ta (Tais the absolute error tolerance). Obviously, it is difficult to defineTawhich should meet with the precision requirement of all approximated curves. In fact, the difference coefficient f(x,x0,…,xn) can be used directly as the criterion; that is
(28)|f(x,x0,…,xn)|<ε.
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 theLis infinite, the computational precision cannot be satisfied except by increasing computational complexity. To avoid this, we design the termination procedure of dynamic choice aboutLas follows:Iff(x0,x1)<Ta, then L=1else iff(x0,x1,x2)<Ta or f(x0,x1,x2)<f(x0,x1), then L=2else iff(x0,x1,x2,x3)<Taor f(x0,x1,x2,x3)<f(x0,x1,x2), then L=3…


#### 3.3.3.L 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 ofLcan change the condition number of the system of algebraic equations discretized by the wavelet interpolation operator or the finite difference method. Therefore, the choice ofLshould take the condition number into account. In fact, if the condition numbercond(A)=10k, then you may lose up tokdigits 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:
(29)Cond(AL+1)Cond(AL)<10.


#### 3.3.4. Relation betweenL 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(x-x0)(x-x1)⋯(x-xn) and the derivative operations. Obviously, according to (5), the increase in computational complexity isO(L3)when the number of extension pointsLincreases 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 ofLis impossible to be infinite. For Haar wavelet which has no smoothness property,Lcan be taken as 0 at most since it need not to be extended. For Faber-Schauder wavelet,Lcan be taken as 1 at most. For Daubechies wavelet,Lcan be taken as different values according to the order of vanishing moments, but it must be finite. For the wavelets without compact support property,Lcan take value dynamically, such as Shannon wavelet. The computational complexity of increasing points is mainly up to the wavelet function of itself.


## 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(10)ϕ(x)=sin(πx)πx
and the Shannon-Gabor scaling function is defined as [17]
(11)G(x)=sin(πx)πxexp(-x22σ2),σ>0,
where σ is the window size.Consider a one-dimensional functionf(x),x∈[a,b]. A discrete point sequence of the variablexis defined as
(12)xn=a+b-a2j·n,j∈Z,n=0,1,2,…,2j,
and the corresponding discrete point sequence of the scaling functionϕ(x) andG(x)can be defined, respectively, as
(13)ϕj,n(x)=ϕj(x-xn)=sin(2jπ/(b-a))(x-xn)(2jπ/(b-a))(x-xn),Gj,n(x)=Gj(x-xn)=sin(2jπ/(b-a))(x-xn)(2jπ/(b-a))(x-xn)×exp(-22j-1(x-xn)2r2(b-a)2),
wherer=2jσ/(b-a).The first and second order derivatives ofϕj(x-xn)at the discrete pointxk are
(14)ϕj′(xk-xn)={0,k=n,2jcos[π(k-n)](k-n)(b-a),k≠n,ϕj′′(xk-xn)={-π23((b-a)/2j)2,k=n,-2cos[π(k-n)]((b-a)/2j)2(k-n)2,k≠n.And the first and second order derivatives ofGj(x-xn)at the discrete pointxk are(15)Gj′(xk-xn)={0,k=n2jcos[π(k-n)]exp[-(k-n)2/2r2](k-n)(b-a),k≠nGj′′(xk-xn)={-3+π2r23r2((b-a)/2j)2,k=n-2cos[π(k-n)]exp[-(k-n)2/2r2]((b-a)/2j)2[1(k-n)2+1r2],k≠n. 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 symbolw(x)as 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:(16)wjk(x)={ϕ(2jx-k)+∑n=-L+1-1ankϕ(2jx-n),k=0,…,Lϕ(2jx-k),k=L+1,…,2j-L-1ϕ(2jx-k)+∑n=2j+12j+L-1bnkϕ(2jx-n),k=2j-L,…,2j,where,(17)ank=∏i=L-1,i≠k-1xj,n-xj,ixj,k-xj,i,bnk=∏i=2j+1,i≠k2j+1+Lxj,n-xj,ixj,k-xj,ixj,k=kxmax-xmin2j,k∈ℤ,
whereLis the amount of the external collocation points, the amount of discrete points in the definition domain is2j+1(j∈ℤ),and[xmin,xmax]is 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 functionf(x)x∈[xmin,xmax]can be expressed as
(18)fj(x)=∑fj(xn)wj(2jx-n),xn=xmin+nxmax-xmin2j.fj(xn) is the given value at the discrete pointxn. At the external points,fj(xn) can be obtained by extrapolation; that is(19)fj(xn)={∑k=0L-1(fj(xk)∏i=0,i≠kL-1xn-xixk-xi),n=-1,…,-L∑k=2j-L+12j(fj(xk)∏i=2j-L+1,k≠i2jxn-xixk-xi),n=2j+1,…,2j+L.So the interval wavelet approximant off(x)can be rewritten as
(20)fj(x)=∑n=-L-1(∑k=0L-1fj(xk)∏i=0L-1xn-xixk-xi)ω(2jx-n)+∑n=02jfj(xk)ω(2jx-n)+∑n=2j+12j+L(∑k=2j-L2jfj(xk)∏i=2j-L2jxn-xixk-xi)ω(2jx-n).
Let
(21)LSL(xn)=∑k=0L-1fj(xk)∏i=0L-1xn-xixk-xi,LEL(xn)=∑k=2j-L2jfj(xk)∏i=2j-L2jxn-xixk-xi,
then
(22)fj(x)=∑n=-L-1LSL(xn)ω(2jx-n)+∑n=02jfj(xk)ω(2jx-n)+∑n=2j+12j+LLEL(xn)ω(2jx-n).LSL(xn)andLEL(xn)correspond 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 ofLby means of Lagrange inner-and extrapolation error polynomial as follows:
(23)RL(x)=f(L+1)(ξ)(L+1)!∏i=0L(x-xi),forsomeξbetweenhhhhhhhhhhhhhhhhhhhhhhhhhhx,x0,…,xL.
Equation (23) indicates that the approximation error is both related to the smoothness and the gradient of the original function near the boundary. Setting differentLcan 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 parameterLis a constant, we have to take a bigger value in order to obtain results with higher calculation precision. But the biggerLusually introduces the famous Gibbs phenomena into the numerical solution, which usually results in the algorithm becoming invalid. In addition, the biggerLwill bring much more calculation. To keep higher numerical precision and save calculation, the best way is to design a procedure thatLcan vary with the curve’s smoothness and gradient dynamically.In this dynamic procedure, the error estimation equation (23) can be taken as the criterion aboutL. 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(24)Nn(x)=f(x0)+(x-x0)f(x0,x1)+(x-x0)(x-x1)f(x0,x1,x2)+⋯+(x-x0)(x-x1)⋯(x-xn-1)×f(x0,x1,…,xn).
Substituting the Newton interpolation instead of the Lagrange interpolation into (22) can be rewritten as
(25)fj(x)=∑n=-L-1(NSL(xn))ω(2jx-n)+∑n=02jfj(xn)ω(2jx-n)+∑n=2j+12j+L(NEL(xn))ω(2jx-n),
where
(26)NSL(xn)=f(x0)+(xn-x0)f(x0,x1)+(xn-x0)(xn-x1)f(x0,x1,x2)+⋯+(xn-x0)(xn-x1)⋯(xn-xL-1)×f(x0,x1,…,xL),NSR(xn)=f(x2j)+(xn-x2j)f(x2j,x2j-1)+(xn-x2j)(xn-x2j-1)×f(x2j,x2j-1,x2j-2)+⋯+(xn-x2j)(xn-x2j-1)⋯(xn-x2j-L)×f(x2j,x2j-1,…,x2j-L).


### 3.3.2. Relation between the Newton Interpolation Error and the Choice ofL

It is well known that the Newton interpolation is equivalent to the Lagrange interpolation. The corresponding error estimation can be expressed as(27)Rn(x)=(x-x0)(x-x1)⋯(x-xn)f(x,x0,…,xn).
And the simplest criterion to terminate the dynamic choice onLis|Rn(x)|≤Ta (Tais the absolute error tolerance). Obviously, it is difficult to defineTawhich should meet with the precision requirement of all approximated curves. In fact, the difference coefficient f(x,x0,…,xn) can be used directly as the criterion; that is
(28)|f(x,x0,…,xn)|<ε.
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 theLis infinite, the computational precision cannot be satisfied except by increasing computational complexity. To avoid this, we design the termination procedure of dynamic choice aboutLas follows:Iff(x0,x1)<Ta, then L=1else iff(x0,x1,x2)<Ta or f(x0,x1,x2)<f(x0,x1), then L=2else iff(x0,x1,x2,x3)<Taor f(x0,x1,x2,x3)<f(x0,x1,x2), then L=3…


### 3.3.3.L 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 ofLcan change the condition number of the system of algebraic equations discretized by the wavelet interpolation operator or the finite difference method. Therefore, the choice ofLshould take the condition number into account. In fact, if the condition numbercond(A)=10k, then you may lose up tokdigits 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:
(29)Cond(AL+1)Cond(AL)<10.


### 3.3.4. Relation betweenL 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(x-x0)(x-x1)⋯(x-xn) and the derivative operations. Obviously, according to (5), the increase in computational complexity isO(L3)when the number of extension pointsLincreases 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 ofLis impossible to be infinite. For Haar wavelet which has no smoothness property,Lcan be taken as 0 at most since it need not to be extended. For Faber-Schauder wavelet,Lcan be taken as 1 at most. For Daubechies wavelet,Lcan be taken as different values according to the order of vanishing moments, but it must be finite. For the wavelets without compact support property,Lcan take value dynamically, such as Shannon wavelet. The computational complexity of increasing points is mainly up to the wavelet function of itself.


## 3.3.1. Newton Interpolation

The expression of Newton interpolation can be written as(24)Nn(x)=f(x0)+(x-x0)f(x0,x1)+(x-x0)(x-x1)f(x0,x1,x2)+⋯+(x-x0)(x-x1)⋯(x-xn-1)×f(x0,x1,…,xn).
Substituting the Newton interpolation instead of the Lagrange interpolation into (22) can be rewritten as
(25)fj(x)=∑n=-L-1(NSL(xn))ω(2jx-n)+∑n=02jfj(xn)ω(2jx-n)+∑n=2j+12j+L(NEL(xn))ω(2jx-n),
where
(26)NSL(xn)=f(x0)+(xn-x0)f(x0,x1)+(xn-x0)(xn-x1)f(x0,x1,x2)+⋯+(xn-x0)(xn-x1)⋯(xn-xL-1)×f(x0,x1,…,xL),NSR(xn)=f(x2j)+(xn-x2j)f(x2j,x2j-1)+(xn-x2j)(xn-x2j-1)×f(x2j,x2j-1,x2j-2)+⋯+(xn-x2j)(xn-x2j-1)⋯(xn-x2j-L)×f(x2j,x2j-1,…,x2j-L).


## 3.3.2. Relation between the Newton Interpolation Error and the Choice ofL

It is well known that the Newton interpolation is equivalent to the Lagrange interpolation. The corresponding error estimation can be expressed as(27)Rn(x)=(x-x0)(x-x1)⋯(x-xn)f(x,x0,…,xn).
And the simplest criterion to terminate the dynamic choice onLis|Rn(x)|≤Ta (Tais the absolute error tolerance). Obviously, it is difficult to defineTawhich should meet with the precision requirement of all approximated curves. In fact, the difference coefficient f(x,x0,…,xn) can be used directly as the criterion; that is
(28)|f(x,x0,…,xn)|<ε.
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 theLis infinite, the computational precision cannot be satisfied except by increasing computational complexity. To avoid this, we design the termination procedure of dynamic choice aboutLas follows:Iff(x0,x1)<Ta, then L=1else iff(x0,x1,x2)<Ta or f(x0,x1,x2)<f(x0,x1), then L=2else iff(x0,x1,x2,x3)<Taor f(x0,x1,x2,x3)<f(x0,x1,x2), then L=3…


## 3.3.3.L 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 ofLcan change the condition number of the system of algebraic equations discretized by the wavelet interpolation operator or the finite difference method. Therefore, the choice ofLshould take the condition number into account. In fact, if the condition numbercond(A)=10k, then you may lose up tokdigits 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:
(29)Cond(AL+1)Cond(AL)<10.


## 3.3.4. Relation betweenL 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(x-x0)(x-x1)⋯(x-xn) and the derivative operations. Obviously, according to (5), the increase in computational complexity isO(L3)when the number of extension pointsLincreases 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 ofLis impossible to be infinite. For Haar wavelet which has no smoothness property,Lcan be taken as 0 at most since it need not to be extended. For Faber-Schauder wavelet,Lcan be taken as 1 at most. For Daubechies wavelet,Lcan be taken as different values according to the order of vanishing moments, but it must be finite. For the wavelets without compact support property,Lcan 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,
(30)τ-α[uin+∑k=1n-1gkuin-k-∑k=0n-1gkui0]=fi+1ui+1n-fi-1ui-1n2h+Kαui+1n-2uin+ui-1nh2,hhhhhhhhhhhhhhhhhhhhhhii=1,2,…,M-1,τ-α[u0n+∑k=1n-1gku0n-k-∑k=0n-1gku00]=4f1u1n-3f0u0n-f2u2n2h+Kαu0n-2u1n+u2nh2,τ-α[uMn+∑k=1n-1gkuMn-k-∑k=0n-1gkuM0]=fM-2uM-2n-4fM-1uM-1n+3fMuMn2h+KαuM-2n-2uM-1n+uMnh2,ui0=φ(xi),1≤i≤M,u0n=p1(tn),uMn=p2(tn),n≥1,gk=(1-1+αk)gk-1,g0=1fi=f(xi)=ν′(xi)mηα.
According to the wavelet collocation method, the fractional Fokker-Planck equation can be approximately represented as
(31)τ-α[uj(xi,tn)+∑k=1n-1gkuj(xi,tn-k)-∑k=0n-1gkuj(xi,t0)]=f′(xi)uj(xi,tn)+∑m=02juj(xm,tn)×[w′(xi-xm)+Kαw′′(xi-xm)],
where i=0,1,2,…2j. Let
(32)Vjn=(uj(x0,tn),uj(x1,tn),…,uj(x2j,tn))T,F=diag(f′(x0),f′(x1),…,f′(x2j)),W1=[w′(x0-x0)w′(x0-x1)⋯w′(x0-x2j)w′(x1-x0)w′(x1-x1)⋯w′(x1-x2j)⋮⋮⋱⋮w′(x2j-x0)w′(x2j-x1)⋯w′(x2j-x2j)]W2=[w′′(x0-x0)w′′(x0-x1)⋯w′′(x0-x2j)w′′(x1-x0)w′′(x1-x1)⋯w′′(x1-x2j)⋮⋮⋱⋮w′′(x2j-x0)w′′(x2j-x1)⋯w′′(x2j-x2j)].
Then the system of (31) can be expressed as the matrix format:
(33)(W1+KαW2+F-τ-αI)Vjn=∑k=1n-1gkVjn-k-∑k=0n-1gkVj0.
Next, we will discuss the precision of the method proposed in this paper with numerical experience. Consider the Fokker-Planck equation as follows:
(34)∂u(x,t)∂t=Dt1-α0[∂∂x(-1)+∂2∂x2]u(x,t),hhhhhhhhhhhhhhhhhhhh0≤x≤1,t>0,
with the initial condition
(35)u(x,0)=x(1-x),0≤x≤1
and the boundary conditions
(36)u(0,t)=-3tαΓ(1+α)-2t2αΓ(1+2α),t>0,u(1,t)=-tαΓ(1+α)-2t2αΓ(1+2α),t>0.
The exact analytic solution is
(37)u(x,t)=x(1-x)+(2x-3)tαΓ(1+α)-2t2αΓ(1+2α).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 errore1  is given by
(38)e1=∥Vjn-Vexactn∥∞,
which provides a measure of the accuracy of the solution near the boundary. The second measure of errore2is given by
(39)e2=12j+1∑i=02j(u(xi)-uexact(xi))2,
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 withL=1andL=2are showen in Figure 1. The boundary effect of the interval wavelet withL=2(Figure 1(a)) is almost eliminated compared toL=1 (Figure 1(b)). FFPE is a 2-order PDEs with respect tox, soL≥2is 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 ofL=2increases more rapid thanL=1with the increase ofjand 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 ofL=2is better thanL=1. The former has failed to solve FFPE obviously. In fact, this explained the reason why we construct the dynamic interval wavelet.Table 1
Condition number of the Fokker-Planck equation.
j
α
τ
=
0.0001
τ
=
0.00001
Interval FDM
Interval wavelet
Interval FDM
Interval wavelet
L
=
1
L
=
2
L
=
1
L
=
2
4
0.8
1.9730
2.5810
2.9399
1.1364
1.2461
1.2738
0.6
11.3319
11.8703
20.8760
2.6645
3.5347
4.3019
0.4
198.8737
91.0582
365.5470
43.4632
31.2993
80.5876
5
0.8
6.2479
7.8009
11.6382
1.5798
2.0236
2.2010
0.6
83.0421
52.1757
160.0987
10.9216
12.0887
20.8948
0.4
1912.4
476.7221
3632.5
380.6050
152.3718
727.8708
6
0.8
39.1651
31.2950
76.7074
39.1651
5.2458
6.9255
0.6
764.3693
255.2801
1476.4
79.4793
51.7012
155.2663
0.4
19847.0000
2574.7
37997
3769.1
790.3041
7238.0
7
0.8
340.7877
145.7761
663.4654
19.8525
19.3224
38.9722
0.6
7757.6000
1333.9
14931
730.3466
249.7668
1416.8
0.4
214660.0000
14202
415100
39724
4266.4
76386Wavelet collocation method with constantL(α=0.8).
(a)
j
=
6
,
L
=
2, condition number is 76.7074,e1=6.2776×10-6, and e2=4.9985×10-6
(b)
j
=
6
,
L
=
1, condition number is 31.2950,e1=1.4×10-3, and e2=1.1×10-3Wavelet collocation method with constantL(α=0.6).
(a)
j
=
6
,
L
=
2, and τ=0.0001
(b)
j
=
6
,
L
=
1, and τ=0.0001The numerical errors comparisons among the dynamic, static interval wavelet method and the interval finite difference method are showen in Figure3. The result also can be illustrated in Table 2.Table 2
Influence ofα on the numerical precision (t=0.0001, T=0.1).
j
α
e
1
e
2
Interval FDM
Interval WCM (L=2)
Dynamic interval WCM
Interval FDM
Interval WCM (L=2)
Dynamic interval WCM
4
0.8
5.5367 × 10−6
8.1588 × 10−5
5.5920 × 10−6
4.1037 × 10−6
5.8298 × 10−5
4.1514 × 10−6
0.6
5.7907 × 10−6
4.1158 × 10−4
7.1813 × 10−6
4.2424 × 10−6
2.9468 × 10−4
5.5971 × 10−6
0.4
3.5309 × 10−6
9.1673 × 10−4
3.7967 × 10−5
2.6232 × 10−6
6.3692 × 10−4
2.1348 × 10−5
5
0.8
5.5551 × 10−6
8.2118 × 10−5
5.8424 × 10−6
4.1718 × 10−6
5.9148 × 10−5
4.4649 × 10−6
0.6
5.7760 × 10−6
4.0932 × 10−4
8.9707 × 10−6
4.2907 × 10−6
2.9753 × 10−4
7.5330 × 10−6
0.4
6.5910 × 10264
0.3971
0.0585
inf
0.1195
0.0491
6
0.8
5.5563 × 10−6
8.1965 × 10−5
1.3154 × 10−5
4.2041 × 10−6
5.9517 × 10−5
8.6913 × 10−6
0.6
3.7124 × 10265
0.0554
0.0096
inf
0.0105
0.0076
0.4
inf
1.2637 × 103
0.0588
inf
375.3305
0.0499
7
0.8
3.4932 × 10243
0.0031
0.0019
inf
3.8313 × 10−4
1.2937 × 10−4
0.6
inf
216.5596
23.7361
inf
40.0263
9.3964
0.4
inf
1.4462 × 106
327.6987
inf
4.2662 × 105
21.7694Numerical errors comparison among the dynamic, static interval wavelet method and the finite difference method (j=6,α=0.6).
(a)
Adapt interval wavelet collocation method
(b)
Interval-finite difference method
(c)
Interval wavelet collocation method with constantL=1The 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 ofL is showen in Table 3. It shows that the value ofL is fixed atL=2after a short time of vibration. This reflects the properties of the FFPE to some extent.Table 3
DynamicL and the iteration times at the same L value (j=6, T=0.1, and τ=0.0001).
L
3
1
2
3
2
Iteration times
11
14
4
3
968In 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 (Figure4). 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.Finite-difference method (α=0.8).
(a)
j
=
6,e1=5.5563×10-6,e2=4.2041×10-6, and condition number is 39.1651
(b)
j
=
3,e1=5.4691×10-6,e2=3.9526×10-6, and condition number is 1.2202All 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 (Figure5), 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.Comparison between the static and dynamic interval Shannon wavelet collocation method (α=0.8).
(a)
Static interval Shannon wavelet collocation method,j=5, and e1=9.95023×10-3,e2=5.23917×10-3
(b)
Dynamic interval Shannon wavelet collocation method,j=5,e1=1.1392×10-3, and e2=4.0893×10-4


## 5. Conclusions

In solving the fractional Fokker-Planck equations, there are two factors related to the choice ofL. The first factor is the condition number, which relates to the parameters α,j and the time stepτ. The largerLcan 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 constantL 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 onLconstructed 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.



---

*Source: 821820-2013-12-07.xml*