Numerical solution of non-linear Fokker–Planck equation using finite differences method and the cubic spline functions

Numerical solution of non-linear Fokker–Planck equation using finite differences method and the cubic spline functions

Applied Mathematics and Computation 262 (2015) 187–190 Contents lists available at ScienceDirect Applied Mathematics and Computation journal homepag...

227KB Sizes 0 Downloads 24 Views

Applied Mathematics and Computation 262 (2015) 187–190

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Numerical solution of non-linear Fokker–Planck equation using finite differences method and the cubic spline functions Behnam Sepehrian a, Marzieh Karimi Radpoor b,∗ a b

Department of Mathematics, Faculty of Science, Arak University, Arak 38156-8-8349, Iran Young Researchers & Elite Club, Hamedan Branch, Islamic Azad University, Hamedan, Iran

a r t i c l e

i n f o

Keywords: Collocation technique Cubic spline Finite differences method Nonlinear Fokker–Planck equation Partial differential equations

a b s t r a c t In this paper we proposed a finite difference scheme for solving the nonlinear Fokker–Planck equation. We apply a finite difference approximation for discretizing spatial derivatives. Then we use the cubic C1 -spline collocation method which is an A-stable method for the time integration of the resulting nonlinear system of ordinary differential equations. The proposed method has second-order accuracy in space and fourth-order accuracy in time variables. The numerical results confirm the validity of the method. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Fokker–Planck equation arises in a number of different fields in natural science, including solid-state physics, chemical physics, quantum optics, theoretical biology, and circuit theory. A Fokker–Planck equation describes the change of probability of a random function in space and time; hence it is naturally used to describe solute transport. The Fokker–Planck equation was first used by Fokker and Planck (for instance, see [11]) to describe the Brownian motion of particles. There is a more general form of Fokker–Planck equation. Nonlinear Fokker–Planck equation has important applications in various areas such as plasma physics, surface physics, population dynamics, biophysics, engineering, neurosciences, nonlinear hydrodynamics, polymer physics, laser physics, pattern formation, psychology and marketing (see [1] and references therein). In one variable case the nonlinear Fokker–Planck equation is written in the following form





∂ u(x, t) ∂ ∂2 = − A(x, t, u) + 2 B(x, t, u) u(x, t), ∂t ∂x ∂x (x, t) ∈ [a, b] × [0, T], with initial condition

u(x, 0) = ϕ(x), and the boundary conditions

u(a, t) = ψ1 (t), ∗

u(b, t) = ψ2 (t),

t ≥ 0.

Corresponding author. Tel.: +98 8134494143. E-mail address: [email protected], [email protected], [email protected] (M.K. Radpoor).

http://dx.doi.org/10.1016/j.amc.2015.03.062 0096-3003/© 2015 Elsevier Inc. All rights reserved.

(1)

188

B. Sepehrian, M.K. Radpoor / Applied Mathematics and Computation 262 (2015) 187–190

Where A(x, t, u) is the drift vector and B(x, t, u) is the diffusion tensor. We assume that ψ 1 and ψ 2 are smooth functions. Note that when A(x, t, u) = A(x, t) and B(x, t, u) = B(x, t) the nonlinear Fokker–Planck equation (1) reduces to the linear Fokker–Planck equation. It is worth noting that some semi-analytic techniques are employed to solve the Fokker–Planck equation. For example this equation is investigated in [4] using the Adomian decomposition method. Also the variational iteration method is developed in [15] to solve this equation. For some other investigations on this model or some other similar models the interested readers can see references [5,6,9,10,16]. Authors of [18] developed a finite difference technique to solve the type of Fokker–Planck equations describing the stochastic dynamics of a particle in a storage ring. In [17] a finite difference procedure is given for solving the Fokker–Planck equation in two dimensions. Lakestani and Dehghan in [7] proposed a numerical scheme for Fokker–Planck equation using the cubic B-spline scaling functions. For more applications of the model studied in this work the interested reader can see [2,3,8]. The organization of this paper is as follows. In Section 2, we present our method for solving nonlinear Fokker–Planck equations. Validation of this method is shown in Section 3 through some examples. A conclusion is drawn in Section 4. 2. Method of solution In this section we will combine second-order central difference in space with cubic C1 -spline collocation method to obtain a high order method for solving the non-linear Fokker–Planck equation (1). At first we discretize partial differential equation (1) in space with central difference to obtain a system of ordinary differential equations with unknown function at each spatial grid point. Then we will apply the cubic C1 -spline collocation method for solving the resulting nonlinear system of ordinary differential equations (see [14]). Also this method can give a closed form approximation for the solution. For positive integers n and T, let h = b−a n denotes the step size of spatial derivatives and k denotes the step size of temporal derivative. So we define

xi = ih,

i = 0, 1, . . . , n,

tj = jk,

j = 0, 1, . . . .

We first rewrite the Eq. (1) as follows

∂u ∂ ∂2 = − [u(x, t)A(x, t, u)] + [u(x, t)B(x, t, u)]. ∂t ∂x ∂ x2

(2)

We have the following relation

∂ A(x, t, u) ∂ A ∂ u ∂ A = + , ∂x ∂u ∂x ∂x

(3)

∂ 2 A(x, t, u) ∂ 2 A ∂ 2A ∂ u ∂ A ∂ 2u + = + . ∂ x2 ∂ x2 ∂ x∂ u ∂ x ∂ u ∂ x2

(4)

and

Using Eqs. (3) and (4) in Eq. (2) we get:







∂u ∂A ∂u ∂A ∂u ∂ 2B ∂ 2B ∂ u ∂ B ∂ 2u = −u + A+u + + − ∂t ∂u ∂x ∂x ∂x ∂ x2 ∂ x∂ u ∂ x ∂ u ∂ x2   ∂u ∂B ∂u ∂B ∂ 2u +2 + + B. ∂x ∂u ∂x ∂x ∂ x2



(5)

Now we define

Q (x) =

∂u ∂A ∂ 2B +u −u 2, ∂t ∂x ∂x

(6)

Then Eqs. (5) and (6) give

Q (x) =











∂u ∂A ∂ 2B ∂B ∂ 2u ∂ B ∂u −u −A+u +2 + +B + u ∂x ∂u ∂ x∂ u ∂x ∂x ∂ x2 ∂ u

2 



2

∂B . ∂u

(7)

If we discretize the above equation with second-order central difference in space at each grid point, we obtain the following relation:

       1 −2 ∂A ∂ 2B ∂B ∂B ∂B 1 −A+u +2 +B +B −u + ui 2 u + 2 u 2h ∂u ∂ x∂ u ∂x ∂u ∂u h h i i      2 1 ∂A ∂ B ∂B ∂B 1 +A+u +2 +B + ui−1 − −u + 2 u 2h ∂u ∂ x∂ u ∂x ∂u h i   u2i+1 − 2ui+1 ui−1 + u2i−1 ∂ B + . 2 ∂u i 4h2 

Qi = ui+1

(8)

B. Sepehrian, M.K. Radpoor / Applied Mathematics and Computation 262 (2015) 187–190

189

By substituting Eq. (6) in Eq. (8), where ui = u(xi , t), we get

ui = a(i 1)ui+1 ui + a(i 2)ui+1 + a(i 3)ui−1 ui + a(i 4)ui−1 + a(i 5)u2i + a(i 6)ui   + a(i 7) u2i+1 − 2ui−1 ui+1 + u2i−1 , where







1 ∂ 2 B  1 ∂ B  ∂ A  + + ,  ∂ u i 2h ∂ x∂ u i h2 ∂ u i    1 ∂ A  1 ∂ 2 B  1 ∂ B  ai (3) = − + 2 ,   2h ∂ u i 2h ∂ x∂ u i h ∂ u i   −2 ∂ B  −2Bi ∂ A  ai (5) = 2 , ai (6) = − +  2 ∂ x i h ∂u i h ai (1) = −

1 2h

ai (2) = −

1 Ai + 2h h

(9)



∂ B  Bi + ∂ x i h2



∂ B  Bi + ∂ x i h2   1 ∂ B  ∂ 2 B  , ai (7) =  2 2 ∂x i 2h ∂ u i ai (4) =

1 Ai − 2h h

If we write Eq. (9) for each grid point we obtain a system of ordinary differential equations which is as follows:

u (t) = F (u(t), t).

(10)

Now we apply the cubic C1 -spline collocation method [12] for solving the system of ordinary differential equations (9). The cubic C1 -spline collocation method is an A-stable method for solving the first-order ordinary differential equations and has fourth-order accuracy. Let U(t) be a vector that approximates u(t) such that each of its component is a cubic spline function and satisfies in (10) at collocation points tj − 1 , tj and tj− 1 in the time interval [tj − 1 , tj ] i.e. U (tl ) = F(U(tl ), tl ), l = j − 1, j − 12 , j and U(t) = [u1 (t), . . . , 2

un − 1 (t)]T . From [12] we have the following relations:

U (t) = U j−1 + kT1 (m)U 

j−1

+ kT2 (m)U 

j− 12

+ kT3 (m)U  , j

(11)

where

3 2 2 3 4 m + m , T2 (m) = 2m2 − m3 , 2 3 3 1 2 T3 (m) = − m2 + m3 , t = tj−1 + km, m ∈ [0, 1] 2 3

T1 (m) = m −

and

U j = U j−1 +

k 1 [F (U j−1 ) + 4F (U j− 2 ) + F (U j )], 6

(12)

and 1

U j− 2 = U j−1 +

k 1 [5F (U j−1 ) + 8F (U j− 2 ) − F (U j )], 24

(13)

in which Uj = U(tj ), Uj = U (tj ). Eqs. (13) and (12) give a nonlinear system of 2n − 2 equations. By solving this system via QuasiNewton method approximate solutions (9) in [tj − 1 , tj ] are obtained. According to order of convergence of the cubic C1 - spline collocation method of [13] the presented method is fourth-order accurate in time component. Note that substituting obtained 1

values U j−1 , U j− 2 and Uj in Eq. (9) give an analytical solution as a polynomial of degree 3 in the interval [tj − 1 , tj ] which is a cubic spline function. 3. Numerical experiments In this section we present the numerical results of the new method on several test problems. We performed our computations using the well-known symbolic software Maple 13. 3.1. Test problem 1 We consider [7] Eq. (1) with a = 0, b = 1 and A(x, t, u) = x, B(x, t, u) =

x2 2 .

The exact solution is given with

u(x, t) = x exp(t). The boundary and initial conditions can be obtained easily from the exact solution. In Table 1 we show the errors obtained for Test problem 1 with different values of k for T = 1. We compared the numerical results, with the results of method presented in [7]. The numerical results are shown in Table 1.

190

B. Sepehrian, M.K. Radpoor / Applied Mathematics and Computation 262 (2015) 187–190 Table 1 Error obtained for Test problem 1 at T = 1 with k = h. Grid point

Method of [7]

h = 0.05

h = 0.025

h = 0.0125

0.2 0.4 0.6 0.8

8.0 × 10−10 3.2 × 10−9 1.3 × 10−8 1.3 × 10−8

3.6010 × 10−10 5.0079 × 10−9 1.2281 × 10−8 2.0546 × 10−8

2.1879 × 10−11 3.1322 × 10−10 7.6753 × 10−10 1.2839 × 10−9

1.3576 × 10−12 1.9580 × 10−11 4.7970 × 10−11 8.0241 × 10−11

Table 2 Error obtained with k = h at T = 1 for Test problem 2. Grid point 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

h = 0.1

h = 0.05 −9

3.4228 × 10 1.5611 × 10−8 3.7232 × 10−8 6.8629 × 10−8 1.1003 × 10−7 1.6143 × 10−7 2.2187 × 10−7 2.8778 × 10−7 3.1989 × 10−7

h = 0.025 −10

2.1295 × 10 9.7022 × 10−10 2.3165 × 10−9 4.2741 × 10−9 6.8582 × 10−9 1.0079 × 10−8 1.3947 × 10−8 1.8433 × 10−8 2.3037 × 10−8

1.3259 × 10−11 6.0452 × 10−11 1.4446 × 10−10 2.6669 × 10−10 4.2809 × 10−10 6.2939 × 10−10 8.7115 × 10−10 1.1538 × 10−9 1.4744 × 10−9

3.2. Test problem 2 We consider Eq. (1) with a = 0, b = 1 and A(x, t, u) =

4u x

− 3x , B(x, t, u) = u. The exact solution is given with

u(x, t) = x2 exp(t). The boundary and initial conditions can be obtained easily from the exact solution. In Table 2 we show the errors obtained for Test problem 2 with different values of k for T = 1. 4. Conclusion In this article, we presented a numerical scheme for solving the nonlinear Fokker–Planck equation. We applied a finite difference approximation for discretizing spatial derivatives, and we used the cubic C1 -spline collocation technique for the time integration of the resulting nonlinear system of ordinary differential equations. One of the advantages of this method is that one can obtain pointwise and analytical approximations for the unknown function. The numerical results confirm the validity of this method and the method is second-order accurate in space and fourth-order in time variable. References [1] M. Dehghan, A computational study of the one-dimensional parabolic equation subject to nonclassical boundary specifications, Numer. Methods Partial Diff. Equations 22 (2006) 220–257. [2] M. Dehghan, M. Lakestani, Numerical solution of nonlinear system of second-order boundary value problems using cubic B-spline scaling functions, Int. J. Comput. Math. 9 (2008) 1455–1461. [3] M. Dehghan, M. Lakestani, The use of cubic B-spline scaling functions for solving the one-dimensional hyperbolic equation with a nonlocal conservation condition, Numer. Methods Partial Diff. Equations 23 (2007) 1277–1289. [4] M. Dehghan, M. Tatari, The use of He’s variational iteration method for solving a Fokker–Planck equation, Phys. Scr. 74 (2006) 310–316. [5] S.A. El-Wakil, M.A. Abdou, A. Elhanbaly, Adomian decomposition method for solving the diffusion–convection–reaction equations, Appl. Math. Comput. 177 (2006) 729–736. [6] G. Harrison, Numerical solution of the Fokker–Planck equation using moving finite elements, Numer. Methods Partial Diff. Equations 4 (1988) 219–232. [7] M. Lakestani, M. Dehghan, Numerical solution of Fokker–Planck equation using the cubic B-spline scaling functions, Numer. Methods Partial Diff. Equations 25 (2009) 418–429. [8] M. Lakestani, M. Dehghan, S. Irandoust-pakchin, The construction of operational matrix of fractional derivatives using B-spline functions, Commun. Nonlinear Sci. Numer. Simul. 17 (2012) 1149–1162. [9] V. Palleschi, M. de Rosa, Numerical solution of the Fokker–Planck equation. II. Multidimensional case, Phys. Lett. A 163 (1992) 381–391. [10] V. Palleschi, F. Sarri, G. Marcozzi, M.R. Torquati, Numerical solution of the Fokker–Planck equation: a fast and accurate algorithm, Phys. Lett. A 146 (1990) 378–386. [11] H. Risken, The Fokker–Planck Equation: Method of Solution and Applications, Springer Verlag, Berlin, Heidelberg, 1989. [12] S. Sallam, M. Naim Anwar, Stabilized cubic C1 spline collocation method for solving first-order ordinary initial value problems, Int. J. Comput. Math. 74 (2000) 87–96. [13] S. Sallam, M. Naim Anwar, M.R. Abdel-Aziz, Unconditionally stable C1-cubic spline collocation method for solving parabolic equations, Int. J. Comput. Math. 81 (2004) 813–821. [14] B. Sepehrian, M. Karimi Radpoor, Numerical solution of Schrodinger equation using compact finite differences method and the cubic spline functions, Int. J. Appl. Math. Res. 3 (2014) 572–578. [15] M. Tatari, M. Dehghan, M. Razzaghi, Application of the Adomian decomposition method for the Fokker–Planck equation, Math. Comput. Model. 45 (2007) 639–650. [16] V. Vanaja, Numerical solution of simple Fokker–Planck equation, Appl. Numer. Math. 9 (1992) 533–540. [17] M.P. Zorzano, H. Mais, L. Vazquez, Numerical solution of two-dimensional Fokker–Planck equations, Appl. Math. Comput. 98 (1999) 109–117. [18] M.P. Zorzano, H. Mais, L. Vazquez, Numerical solution for Fokker–Planck equations in accelerators, Physica D 113 (1998) 379–381.