# Interval Wavelet Numerical Method on Fokker-Planck Equations for Nonlinear Random System

Authors: Li-wei Liu
Journal: Advances in Mathematical Physics
 (2013)

---

## Abstract

The Fokker-Planck-Kolmogorov (FPK) equation governs the probability density function (p.d.f.) of the dynamic response of a particular class of linear or nonlinear system to random excitation. An interval wavelet numerical method (IWNM) for nonlinear random systems is proposed using interval Shannon-Gabor wavelet interpolation operator. An FPK equation for nonlinear oscillators and a time fractional Fokker-Planck equation are taken as examples to illustrate its effectiveness and efficiency. Compared with the common wavelet collocation methods, IWNM can decrease the boundary effect greatly. Compared with the finite difference method for the time fractional Fokker-Planck equation, IWNM can improve the calculation precision evidently.


---

## Body

## 1. Introduction

The Fokker-Planck equation describes the time evolution of the probability density function of the velocity of a particle, and can be generalized to other observables as well. It is named after Adriaan Fokker and Max Planck and is also known as the Kolmogorov forward equation (diffusion), named after Andrey Kolmogorov, who first introduced it in a 1931 paper. When applied to particle position distributions, it is better known as the Smoluchowski equation. The case with zero diffusion is known in statistical mechanics as Liouville equation. In order to describe subdiffusive behavior of a particle under the combined influence of external nonlinear force field and a Boltzmann thermal heat bath, Metzler et al. introduced a fractional Fokker-Planck equation (FFPE) which was shown to obey generalized Einstein relations, and its stationary solution is the Boltzmann distribution [1].Many methods for calculating nonlinear random response have been developed by numerous scholars over a long period of time. One type of these methods is the diffuseness process theory method, and the primarily one is Fokker-Planck equation method. In practice, the most difficult problem of using Fokker-Planck equation method is how to solve the equation. For general nonlinear system, it is very difficult to obtain the exact solution. Various numerical methods have been used to solve the Fokker-Planck equation, such as the weighted residual method [2], the finite element method [3], and the path integral method [4] and so on. The Galerkin method for the numerical solution of the stationary Fokker-Planck equation developed by Bhandari and Sherrer is based on taking multiple Hermite polynomial as joint probability density, but the rate of convergence of this method is slow for strong nonlinear system. Based on Galerkin method, a finite element method for Fokker-Planck equation was developed by Langley; this method is more efficient than Galerkin method in computation. The common problem of these three methods is that the error at the tails and peak of the response distribution is larger. As a numerical approach method, the moment equation method, including the non-Gaussian moment-cutting method for nonlinear random vibration, was developed [5].Recently the wavelet analysis is getting high attention to many authors for the nonlinear dynamic system, not only for the signal analysis but also for developing new numerical methods for calculating the partial difference equations. Bertoluzza and Naldi presented a wavelet collocation method for solving the partial differential equations [6] in 1996. In this method, the autocorrelation function of the Daubechies scaling functions was used as the weight function. The autocorrelation function of the Daubechies wavelet functions possesses interpolation properties due to the orthogonality properties of Daubechies wavelets. Therefore, the numerical solution of the discrete format of the partial differential equation, which was obtained by wavelet collocation method, is an approximation solution of the partial differential equation. This method is not similar to the wavelet Galerkin method, which needs going back and forth between the space of wavelet coefficients and physical space. On the other hand, the Daubechies scaling functions and the corresponding wavelet function have no explicit expressions, which limit its applications to some extents. Cattani studied the properties of the Shannon wavelet function, which possesses many advantages such as orthogonality, being continuous, and being 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 [7]. The approximation is compared with the wavelet reconstruction formula and the error of approximation is explicitly computed [8]. The advantages of the Shannon wavelet have been illustrated in solving PDEs in engineering [9–11]. McWilliam et al. [12] presented a Shannon wavelet collocation method for stationary Fokker-Planck equation. The Shannon scaling function was taken as the weight function in this method, which can avoid the shortcomings of Daubechies wavelet such as the interpolation property. Based on the theory of this method, the specific procedure for solving stationary Fokker-Planck equation, the analysis of the problems existed in this method, and the remedies of the problems were presented in this paper. But the Shannon wavelet has no compact support property which can decrease the efficiency and computational accuracy. Compared with Shannon wavelet function, Shannon-Gabor wavelet [13, 14] possesses better compact support property, which has been widely employed in various mechanical analysis fields [15] including solving the FPK equation [16]. Furthermore, Shi et al. [17] constructed the Shannon-Gabor wavelet defined in the interval based on the Lagrange interpolation theory, which can decrease the boundary effect effectively.Mei proposed a general construction method of the interval wavelet based on the general variational principle [18] and took the Shannon-Gabor wavelet as example to illustrate its correctness, which has been employed to eliminate the boundary effect in solving PDEs such as the image processing [19], stochastic analysis [20, 21], option pricing [22], hydrodynamics [23–26], and image segmentation [27]. The purpose of this paper is to construct an interval Shannon-Gabor wavelet collocation method on solving Fokker-Planck equation for nonlinear oscillators and a time fractional Fokker-Planck equation describing a subdiffusive behavior of a particle under the combined influence of external nonlinear force field and a Boltzmann thermal heat bath [28].


## 2. Shannon-Gabor Wavelet Collocation Method

### 2.1. Construction of the Basis Function

Consider a one dimension functionf(x). In order to take Shannon scaling function as basic function, in term of the multiresolution analysis theory, the function f(x) needs to be discretized uniformly in range of [a, b], which is the domain of definition of it. Let the number of the discrete points is 2j+1, (j∈Z), and the definition of the discrete point of variable x is
(1)xi=a+b-a2j·i.
Then we can obtain the basic function as follows:
(2)wj(x-xi)=sin(2jπ/(b-a))(x-xi)(2jπ/(b-a))(x-xi)×exp(-22j-1(x-xi)2r2(b-a)2).
Using the theory of Fourier transform, it can be proved that the basic function ω(x-xi) satisfied the conclusions as follows:
(3)wj(xi-xk)=δik,δik={1,i=k,0,i≠k,∫-∞∞wj(x-xi)wj(x-xk)=b-a2jδik,∫-∞∞wj(x-xi)dnwj(x-xk)dxn=b-a2jdnwj(x-xk)dxn.The two dimension basic function can be expressed by tensor product of the one dimension basic function. In two dimension scaling function space, to eachj∈Z, Vj can be expressed as Vj=vj⊗vj, and so the corresponding basic function is
(4)Wj,n1,n2(x,y)=wj(x-xn1)wj(x-xn2).
In a similar way, the basic function of n dimension space can be expressed as
(5)Wj,n1,n2,…nn(x1,x2,…,xn)=∏m=1nwj(xm-xmnm).
According to the interval interpolation basic functions, the Shannon-Gabor interval wavelet can be obtained as follows:
(6)wjk=w(2jx-k)+∑n=-N+1-1ankw(2jx-n),k=0,…,L,wjk=w(2jx-k),k=L+1,…,2j-L-1,wjk=w(2jx-k)+∑n=2j+12j+N-1bnkw(2jx-n),k=2j-L,…,2j.ank and bnk could be calculated, respectively, as
(7)ank=ljk1(xjn),bnk=ljk2(xjn),lj,k1=∏-1i=-L-1i≠kx-xj,ixj,k-xj,i,lj,k2=∏2j+1+Li=2j+1i≠kx-xj,ixj,k-xj,i,
where L is the number of the external points, and N is the support domain of the wavelet function, that is, supϕ=[-N,N].It is easy to know that the Shannon-Gabor interval wavelet is a linear combination of the Shannon-Gabor scaling functionwj,k(x). Therefore, the Shannon-Gabor interval wavelet function possesses all the properties of the Shannon-Gabor scaling function.It should be noted that the construction method of the interval wavelet function reveals the close relationship between the restricted variational principle and the interval interpolation wavelet.


### 2.2. Discrete Format of the Fokker-Planck Equation

Consider the following two-dimension stationary Fokker-Planck equation:(8)∂(a1p)∂x+∂(a2p)∂y-12∂2(b11p)∂x2-12∂2(b22p)∂y2-12∂2(b12p)∂x∂y-12∂2(b21p)∂x∂y=0,
where
(9)a1=limΔt→01ΔtE[z1-x],a2=limΔt→01ΔtE[z2-y],b12=b21=limΔt→01ΔtE[(z1-x)(z2-y)],b11=limΔt→01ΔtE[(z1-x)2],b22=limΔt→01ΔtE[(z2-y)2].The unknown quantityp denotes transitional probability density function p(x,y), which can be expressed as
(10)p(x,y)≈pj(x,y)=∑n1=02j∑n2=02jpj(xn1,yn2)wj(x-xn1)wj(y-yn2).In general,a1, a2, b11, b12, b21, b22 are the known function on the  variables   x and y. To avoid having to perform the complicated integration presented in the collocation method, those aforementioned known functions can be expressed as a sum of shape functions, such that
(11)am(x,y)·p(x,y)≈ajm(x,y)pj(x,y)=∑n1=02j∑n2=02jajm(xn1,yn2)pj(xn1,yn2)×wj(x-xn1)wj(y-yn2)bmn(x,y)·p(x,y)≈bjmn(x,y)pj(x,y)=∑n1=02j∑n2=02jbjmn(xn1,yn2)pj(xn1,yn2)×wj(x-xn1)wj(y-yn2),
where m=1, n=2.In terms of the theory of collocation method, substituting (11) into (8), we can obtain the system of equations as follows:
(12)∫Wj,k1,k2[∂(aj1pj)∂x+∂(aj2pj)∂y-12∂2(bj11pj)∂x2-12∂2(bj22pj)∂y2-12∂2(bj12pj)∂x∂y-12∂2(bj21pj)∂x∂y]dxdy=0,
where k1,k2=0,1,2,…,2j.Using (3) and integrating with respect to (12), the system of equations can be obtained as follows:
(13)∑n1=02j∑n2=02j‍pj(xn1,yn2)×[aj1(xn1,yn2)δk2n2∂w(xk1-xn1)∂x+aj2(xn1,yn2)δk1n1∂w(yk2-yn2)∂y-12bj11(xn1,yn2)δk2n2∂w2(xk1-xn1)∂x2-12bj22(xn1,yn2)δk1n1∂w2(yk2-yn2)∂y2-bj12(xn1,yn2)∂w(xk1-xn1)∂x∂w(yk2-yn2)∂y]=0,
where, k1,k2=0,1,2,…,2j.As the basic function is the Shannon scaling function which possess explicit expression, so it is easy to deduce the first and second derivative of the basic functionw as follows:
(14)∂(xk-xn)∂x={0,k=n,2jcos[π(k-n)](b-a),k≠n,(15)∂2(xk-xn)∂x2={-22j+1π26(b-a)2,k=n,-22j+1cos[π(k-n)](b-a)2(k-n)2,k≠n.Instituting (14) and (15) into (13), a system of homogenous algebraic equations can be deduced, such that
(16)AP=0.
This is the discrete format of stationary Fokker-Planck equation in equinoctial points of the domain of definition.


### 2.3. Solution of the Discrete Format of Fokker-Planck Equation

In the system of (16),  A  is a square matrix; vector  P  contains the values of the j.p.d.f. on the mesh of equinoctial points. Equation (16) can be solved by prescribing the value of p(x,y) at node 1 to be some constant C, whichs lead to
(17)A1P1=CV,
where A1 is the matrix A without the first row or column, P1 is the vector of nodal values from node 2 onwards, and V is the first column of A with the first entry omitted. This then yields the value of p(x,y) at each node in terms of C. Finally, C can be evaluated by imposing the normalization property
(18)∫a1b1∫a2b2p(x,y)dxdy=1,
where a1, b1 and a2, b2 are numeric area of the variables x and y. So we should firstly evaluate the value of a1, b1 and a2, b2. They can be obtained approximately from the solution of the random vibration differential equation by using the equivalent linear linearization method. Based on the equivalent linearization results, the numeric area of the variables x and y can be evaluated.


## 2.1. Construction of the Basis Function

Consider a one dimension functionf(x). In order to take Shannon scaling function as basic function, in term of the multiresolution analysis theory, the function f(x) needs to be discretized uniformly in range of [a, b], which is the domain of definition of it. Let the number of the discrete points is 2j+1, (j∈Z), and the definition of the discrete point of variable x is
(1)xi=a+b-a2j·i.
Then we can obtain the basic function as follows:
(2)wj(x-xi)=sin(2jπ/(b-a))(x-xi)(2jπ/(b-a))(x-xi)×exp(-22j-1(x-xi)2r2(b-a)2).
Using the theory of Fourier transform, it can be proved that the basic function ω(x-xi) satisfied the conclusions as follows:
(3)wj(xi-xk)=δik,δik={1,i=k,0,i≠k,∫-∞∞wj(x-xi)wj(x-xk)=b-a2jδik,∫-∞∞wj(x-xi)dnwj(x-xk)dxn=b-a2jdnwj(x-xk)dxn.The two dimension basic function can be expressed by tensor product of the one dimension basic function. In two dimension scaling function space, to eachj∈Z, Vj can be expressed as Vj=vj⊗vj, and so the corresponding basic function is
(4)Wj,n1,n2(x,y)=wj(x-xn1)wj(x-xn2).
In a similar way, the basic function of n dimension space can be expressed as
(5)Wj,n1,n2,…nn(x1,x2,…,xn)=∏m=1nwj(xm-xmnm).
According to the interval interpolation basic functions, the Shannon-Gabor interval wavelet can be obtained as follows:
(6)wjk=w(2jx-k)+∑n=-N+1-1ankw(2jx-n),k=0,…,L,wjk=w(2jx-k),k=L+1,…,2j-L-1,wjk=w(2jx-k)+∑n=2j+12j+N-1bnkw(2jx-n),k=2j-L,…,2j.ank and bnk could be calculated, respectively, as
(7)ank=ljk1(xjn),bnk=ljk2(xjn),lj,k1=∏-1i=-L-1i≠kx-xj,ixj,k-xj,i,lj,k2=∏2j+1+Li=2j+1i≠kx-xj,ixj,k-xj,i,
where L is the number of the external points, and N is the support domain of the wavelet function, that is, supϕ=[-N,N].It is easy to know that the Shannon-Gabor interval wavelet is a linear combination of the Shannon-Gabor scaling functionwj,k(x). Therefore, the Shannon-Gabor interval wavelet function possesses all the properties of the Shannon-Gabor scaling function.It should be noted that the construction method of the interval wavelet function reveals the close relationship between the restricted variational principle and the interval interpolation wavelet.


## 2.2. Discrete Format of the Fokker-Planck Equation

Consider the following two-dimension stationary Fokker-Planck equation:(8)∂(a1p)∂x+∂(a2p)∂y-12∂2(b11p)∂x2-12∂2(b22p)∂y2-12∂2(b12p)∂x∂y-12∂2(b21p)∂x∂y=0,
where
(9)a1=limΔt→01ΔtE[z1-x],a2=limΔt→01ΔtE[z2-y],b12=b21=limΔt→01ΔtE[(z1-x)(z2-y)],b11=limΔt→01ΔtE[(z1-x)2],b22=limΔt→01ΔtE[(z2-y)2].The unknown quantityp denotes transitional probability density function p(x,y), which can be expressed as
(10)p(x,y)≈pj(x,y)=∑n1=02j∑n2=02jpj(xn1,yn2)wj(x-xn1)wj(y-yn2).In general,a1, a2, b11, b12, b21, b22 are the known function on the  variables   x and y. To avoid having to perform the complicated integration presented in the collocation method, those aforementioned known functions can be expressed as a sum of shape functions, such that
(11)am(x,y)·p(x,y)≈ajm(x,y)pj(x,y)=∑n1=02j∑n2=02jajm(xn1,yn2)pj(xn1,yn2)×wj(x-xn1)wj(y-yn2)bmn(x,y)·p(x,y)≈bjmn(x,y)pj(x,y)=∑n1=02j∑n2=02jbjmn(xn1,yn2)pj(xn1,yn2)×wj(x-xn1)wj(y-yn2),
where m=1, n=2.In terms of the theory of collocation method, substituting (11) into (8), we can obtain the system of equations as follows:
(12)∫Wj,k1,k2[∂(aj1pj)∂x+∂(aj2pj)∂y-12∂2(bj11pj)∂x2-12∂2(bj22pj)∂y2-12∂2(bj12pj)∂x∂y-12∂2(bj21pj)∂x∂y]dxdy=0,
where k1,k2=0,1,2,…,2j.Using (3) and integrating with respect to (12), the system of equations can be obtained as follows:
(13)∑n1=02j∑n2=02j‍pj(xn1,yn2)×[aj1(xn1,yn2)δk2n2∂w(xk1-xn1)∂x+aj2(xn1,yn2)δk1n1∂w(yk2-yn2)∂y-12bj11(xn1,yn2)δk2n2∂w2(xk1-xn1)∂x2-12bj22(xn1,yn2)δk1n1∂w2(yk2-yn2)∂y2-bj12(xn1,yn2)∂w(xk1-xn1)∂x∂w(yk2-yn2)∂y]=0,
where, k1,k2=0,1,2,…,2j.As the basic function is the Shannon scaling function which possess explicit expression, so it is easy to deduce the first and second derivative of the basic functionw as follows:
(14)∂(xk-xn)∂x={0,k=n,2jcos[π(k-n)](b-a),k≠n,(15)∂2(xk-xn)∂x2={-22j+1π26(b-a)2,k=n,-22j+1cos[π(k-n)](b-a)2(k-n)2,k≠n.Instituting (14) and (15) into (13), a system of homogenous algebraic equations can be deduced, such that
(16)AP=0.
This is the discrete format of stationary Fokker-Planck equation in equinoctial points of the domain of definition.


## 2.3. Solution of the Discrete Format of Fokker-Planck Equation

In the system of (16),  A  is a square matrix; vector  P  contains the values of the j.p.d.f. on the mesh of equinoctial points. Equation (16) can be solved by prescribing the value of p(x,y) at node 1 to be some constant C, whichs lead to
(17)A1P1=CV,
where A1 is the matrix A without the first row or column, P1 is the vector of nodal values from node 2 onwards, and V is the first column of A with the first entry omitted. This then yields the value of p(x,y) at each node in terms of C. Finally, C can be evaluated by imposing the normalization property
(18)∫a1b1∫a2b2p(x,y)dxdy=1,
where a1, b1 and a2, b2 are numeric area of the variables x and y. So we should firstly evaluate the value of a1, b1 and a2, b2. They can be obtained approximately from the solution of the random vibration differential equation by using the equivalent linear linearization method. Based on the equivalent linearization results, the numeric area of the variables x and y can be evaluated.


## 3. Numerical Experiments

Example 1.
In this section the method was applied to the random vibration of a Duffing oscillator, the equation of motion for which is [2]
(19)x¨+2βωnx˙+ωn2(1+γx2)x=F(t),
where F(t) is Gaussian white noise, taken here to have a spectral density of (1/π). Putting z1=x, z2=x˙, (19) is deduced to
(20)(z˙1z˙2)=(z2-2βωnz2-ωn2(z1+γz13))+(000π1/2)(0W(t)),
and so we obtain the Fokker-Planck equation as follows:
(21)∂∂z1(z2p)-∂∂z2{[2βωnz2+ωn2(z1+γz13)]p}-∂2p∂z22=0.
The exact analytical solution of (21) is
(22)p(z1,z2)=Cexp[-2βωn(12z22+12ωn2z12+14ωn2γz14)],
where C is chosen to satisfy the normalization condition, (14).For the caseωn=1, β=0.5, γ=0.1, the r.m.s. values of the response of displacement and velocity are σx=0.7407 and σx˙=1.0, which were obtained from the equivalent linear linearization method. Based on this result, a finite region of -3≤x≤3, -6≤y≤6 was chosen for the wavelet collocation method.The numerical resultspj(x,0), pj(0,y) and exact result pje(x,0), pje(0,y) are plotting in Figure 1. It can be seen that the proposed method gives excellent agreement with the exact result, both over the main body of the curve and up to the tails.The comparison between numerical solution and exact solution.
(a)
p
j
(
x
,
0
)
(b)
p
j
(
0
,
y
)
(c)
p
j
e
(
x
,
0
)
-
p
j
(
x
,
0
)
(d)
p
j
e
(
0
,
y
)
-
p
j
(
0
,
y
)The comparison between the wavelet and interval wavelet collocation method for FPK equation is shown in Figure2. It is easy to see that the error of the interval wavelet collocation method is smaller evidently than the common wavelet method in solving the FPK equation. It is not difficult to understand that the boundary effect of the common wavelet numerical method enlarges the error in the whole definition domain. This illustrates the effect of the interval wavelet numerical method to some extent.Error comparison between the wavelet and the interval wavelet numerical methods for FPK equation.
(a)
Shannon-Gabor wavelet numerical method (the maximum of the absolute error is2.2470×10-3)
(b)
Interval Shannon-Gabor wavelet numerical method (the maximum of the absolute error is4.1769×10-6)Example 2.
Time fractional Fokker-Planck equations have been recently treated by a number of authors and are found to be a useful approach for the description of transport dynamics in complex systems that are governed by anomalous diffusion and nonexponential relaxation patterns. Fractional derivatives play a key role in modeling particle transport in anomalous diffusion. Wei [14] introduced a time fractional extension of the FPK equation, namely, the time-fractional Fokker-Planck equations as follows:
(23)∂u(x,t)∂t=D0t1-α[∂∂x(-1)+∂2∂x2]u(x,t),0≤x≤1,t>0
with initial condition
(24)u(x,0)=x(1-x),0≤x≤1,
and boundary conditions
(25)u(0,t)=-3tαΓ(1+α)-2t2αΓ(1+2α),t>0,u(1,t)=-tαΓ(1+α)-2t2αΓ(1+2α),t>0.
The exact analytical solution of the time-fractional Fokker-Planck equation can be expressed as
(26)u(x,t)=x(1-x)+(2x-3)tαΓ(1+α)-2t2αΓ(1+2α).
As mentioned above, the solution of (23) can be expressed approximately as
(27)u(x)≈uj(x)=∑n=02juj(xn)wj(x-xn).
Substituting (27) into (23), we can obtain
(28)D0t1-α(∑n=02juj(xn,t)[ω′′(xk-xn)-ω′(xk-xn)])=duj(xk,t)dt,
where k=0,1,2,…,2j. Let
(29)U(t)=(uj(x0,t),uj(x1,t),…,uj(x2j,t))T,W1=[ω′(x0-x0)ω′(x0-x1)⋯ω′(x0-x2j)ω′(x1-x0)ω′(x1-x1)⋯ω′(x1-x2j)⋮⋮⋱⋮ω′(x2j-x0)ω′(x2j-x1)⋯ω′(x2j-x2j)],W2=[ω′′(x0-x0)ω′′(x0-x1)⋯ω′′(x0-x2j)ω′′(x1-x0)ω′′(x1-x1)⋯ω′′(x1-x2j)⋮⋮⋱⋮ω′′(x2j-x0)ω′′(x2j-x1)⋯ω′′(x2j-x2j)],M=W2-W1.
We can obtain the discretized format to solve the time-fractional FPK equations as follows:
(30)MU(tn+1)=-∑k=1ngkU(tn-k+1)+∑k=0ngkU(t0),MU(t1)=U(t0).
The comparison between the finite difference method and the wavelet numerical method in solving the time fractional FPK equation is shown in Figure 3. Although the iteration times are only 1000, the error of the numerical solution obtained by the finite difference method becomes larger evidently than one obtained by the wavelet numerical method. It should be noticed that the effect of the interval wavelet is limited as the explicit boundary condition in this example.Comparison between the finite difference method and the wavelet numerical method in solving the time fractional FPK equation (α=0.6, t=0.1, τ=0.0001, j=4).
(a)
Finite difference method
(b)
Wavelet numerical method


## 4. Conclusions

Compared with the finite difference method, the Shannon-Gabor wavelet numerical method possesses more excellent numerical properties in solving the FPK equations, which has been illustrated in solving the time-fractional FPK equations. Based on the interpolation properties of the Shannon-Gabor wavelet, the interval interpolation wavelet collocation method based on the wavelet interpolation technique has been developed to solve the two dimension Fokker-Planck equation in a finite domain in this paper. This new method can decrease the boundary effect evidently and then decrease the numerical error in the whole in the definition domain greatly. We believe that the method can be easily generalized to rectangular higher dimension case. The comparisons with other numerical algorithms show that the method is competitive and efficient. Furthermore, it should be noted that the method can also be used to solve general partial differential equation.



---

*Source: 651357-2013-10-27.xml*