Approximate solution for Fokker-Planck equation

In this paper, an approximate solution to a specific class of the Fokker-Planck equation is proposed. The solution is based on the relationship between the Schr\"{o}dinger type equation with a partially confining and symmetrical potential. To estimate the accuracy of the solution, a function error obtained from the original Fokker-Planck equation is suggested. Two examples, a truncated harmonic potential and non-harmonic polynomial, are analyzed using the proposed method. For the truncated harmonic potential, the system behavior as a function of temperature is also discussed.


Introduction
Differential equations are used to model different natural phenomena such as diffusion that are of great importance in many physical, chemical and biological processes [1].
In general, diffusive processes can be treated via the Fokker-Planck equation [2][3][4][5][6]. This equation is obtained from the Langevin equation and gives the probability of finding a given particle in a state x at time t . The usual representation of the Fokker-Planck equation is given as follows: where t is the time variable, and x is the variable characteristic of the system (which may be identified, for example, as velocity). Q is the diffusion coefficient and P (x, t ) is the probability distribution. The f (x) function is known as the external force that acts on the system, although this designation is only suitable when x represents velocity. This function can be identified as the derivative of a potential function V (x): Different methods for treating the Fokker-Planck equation have been suggested as, for example, its association with a Schrödinger type equation [2,7]. However, there are only a few cases where the equation (1.1) has an exact analytical solution.
In a recent study [8], an approximate solution for partially confining potentials was proposed, where part of the solution is written in terms of functions that arise from the solution of the Schrödinger type equation for bound states and part is composed of functions originating from the free particle. This solution method is discussed in reference [9] for the particular case of the Rosen-Morse potential. This potential has an exact solution to the Schrödinger equation. As expected, it can be concluded that the results obtained using the functions derived from the original Schrödinger type equation provide better solutions than those obtained by the proposed approximate method. However, this approach proves to be very restrictive, since the number of potentials with exact solutions to the Schrödinger equation is very small.
In this paper, an improvement to the previously proposed solution [8] is suggested, producing a more accurate result for cases of symmetrical partially-confining potentials. To check the validity of the results, the error on the results is estimated through direct substitution of the approximate solution into the Fokker-Planck equation.
In the study of the truncated harmonic potential, the amount of particles that escape from the potential well are calculated and a phase transition for the systems being studied is identified.
In section 2, there is a brief discussion of the Fokker-Planck equation and the approximate solution method is presented, with the changes suggested in relation to the previously proposed solution [8]. In sections 3 and 4, the proposed method is applied to two partially-confining potentials. In the first case, the region of the well potential is described by a harmonic potential and, for the second case, the well is given by a non-harmonic polynomial potential. Finally, the conclusions are presented in section 5.

The Fokker-Planck equation
Given the importance and difficulty of solving the Fokker-Planck equation (FPE), different methods have been proposed to study this equation. Such methods include numerical treatments [10,11] and mapping the FPE onto a Schrödinger type equation [2,7,12]. In the latter case, the expression of P (x, t ) is given by a series of functions as follows: where the eigenfunctions φ n (x) and the eigenvalues λ n are the solution of the Schrödinger type equation Comparing equation (2.2) with the Schrödinger equation, it can be seen that the term in parentheses is equivalent to a potential function. Thus, this term is commonly called the effective potential, V ef (x): The probability distribution indicated by equation (2.1) assumes that the initial condition for the time t = 0 is the probability distribution expressed by a delta function P (x, 0) = δ(x).
In order to obtain the numerical solution shown in equation (2.1), it is necessary to define a criterion to the cut off in the sum to truncate the series. Another problem arises when one cannot get the solution of the corresponding Schrödinger equation (2.2). To get around this last problem, an approximate analytical solution composed of two parts was suggested in a previous paper [8]. This solution involves the discrete levels of the problem and the other one uses a Gaussian distribution for all continuous states. Although the solution proposed in [8] agrees with the numerical results, the adopted approach can be improved upon.
The attention here is focused on specific issues involving symmetrical and partially confined potentials, where it is not possible to get the complete solution of the Schrödinger type equation (2.2). Thus, for these types of potential, a similar treatment to that of reference [9], for example, is unfeasible. For the potentials studied, there are regions where the spectrum is continuous and a region where there may be bound states (i.e., a discrete spectrum). For x > d and x < −d, it is assumed that the potential is constant and the spectrum is continuous. In −d < x < d, there is a potential well and the spectrum becomes discrete. It is also assumed that the potential is continuous, especially at the points x = ±d that correspond to the intersections between the region where the potential is constant and the well potential region (figures 1, 6 and 8 exemplify the type of the studied potential). The probability distribution in such cases is calculated separately for each region of the potential.

43003-2
The function ρ(x, t ) in (2.4) is the probability distribution for each region; N I , N II and N III are related to the normalization in each piece of the probability distribution. For regions (I and III) of the continuous spectrum, the ρ functions are given by a Gaussian function, as suggested previously [8] and for region II, the function ρ is expressed by the series indicated in (2.1) with a limited number of eigenvalues ( j ) The number j corresponds to the number of discrete eigenvalues present in the potential well under analysis.
As the functions ρ I (x, t ) and ρ III (x, t ) are equal, because of the symmetry of the problem, one can assume that the normalization parameters N I and N III are the same. It can also be assumed that at the interface points between the potentials (±d ), the distribution ρ(x, t ), equation (2.4) should be continuous, i.e., ρ II (d, t ) = ρ III (d, t ) and ρ I (−d, t ) = ρ II (−d, t ). Thus, the condition of the continuity of the distribution and the symmetry of the problem imply that N I = N III and This relationship (2.7) shows that the normalization N II depends on N I and also, in general, depends on time. To simplify the notation, we will rewrite equation (2.7) as N II (t ) = N I g (t ), such that: Therefore, the overall probability distribution for a problem with the discussed features (a partially confining and symmetric potential) is obtained by equation (2.4), subject to the condition (2.7), i.e., x > d. (2.9) Applying the normalization condition to the probability distribution (2.9), one gets According to this result, the normalization parameter N I is dependent on time and the probability distribution (2.9) should be normalized for each time value. Looking at the approximate solution given by equation (2.9), one can see that when d is close to zero the system approximated to a free particle system and the probability distribution approximates to a Gaussian. On the other hand, when d is very large (d → ∞), i.e., the size of the system tends to an infinite well, the solution of the problem is given by the usual discrete series of functions, as shown in equation (2.1).
The difference between the solution given by (2.9) and that shown in reference [8] is that here, the Gaussian is used only in areas where the potential is constant, while in reference [8], it was suggested that it should be included in all regions of the space. Numerical results show that the new approach is more suitable, leading to more accurate results.
In general, when there is a partial confinement, a temporal dependency arises already in the coefficient N I . Under these conditions, the probability distribution for large times in the region of the potential

43003-3
well (region II) may go to zero, which indicates the escape of particles from the minimum region of the potential.
From the proposed solution it is suggested that the escape of particles from the potential well can be quantified by the value Y (t ,Q) defined by: The function Y (t ,Q) gives the number of particles within the confinement region for each time t and different values of the diffusion coefficient (Q). The functions N I (t ) and g (t ) are given by expressions (2.10) and (2.8), respectively. Equation (2.11) also allows the evaluation of the influence of temperature in the escape process of the particles in the well. Assuming that the temperature is proportional to the diffusion coefficient [2], the calculation of the population for different values of Q allows for the analysis of the evolution of the system in terms of temperature. This information allows, in principle, the study of the thermodynamic properties of the system, such as phase transitions.
To check the accuracy of the proposed method, we introduce the function ε(x, t ) based on the Fokker- This function provides a quantitative parameter to check if the solution approximates to the actual solution of the problem. If the solution is accurate, then ε(x, t ) = 0. In this way, the further this function ε(x, t ) is close to zero, the better is the proposed function P (x, t ) to describe the real solution for the system under study.
Substituting the solution presented in (2.9) into expression (2.12), it can be seen that for x < −d and Since the construction of ε(x, t ) involves derivatives of the probability distribution P (x, t ), it is worth noting that, close to the points where the derivative is discontinuous, the use of this criterion is impaired and should be used with caution. Thus, in this study, ε(x, t ) was not defined for values of x close to ±d. However, the use of the function ε(x, t ) to evaluate the solution avoids comparisons with solutions obtained by other methods which could, in itself, introduce an additional error.

Truncated harmonic oscillator
In this section, we apply the approximate solution of the FPE to a model with a truncated harmonic oscillator, whose strength is given by with (±d ) being the interface points between the harmonic potential and the constant potential. Substituting this expression of force into equation (2.3), it can be seen that the effective potential can be identified by where v 0 = k 2 d 2 /4Q − k/2 is a constant value chosen such that the potential function is continuous.   figure 1 shows that there are two distinct regions: one is a potential well region described by a harmonic potential and the other is described by a constant potential.

43003-4
The harmonic oscillator problem is a case in which equation (2.2) has a full analytical solution [2] and the probability distribution is given by: where H n are Hermite polynomials, α = k/2Q and k is a constant. The eigenvalues λ n are equal to nk (n = 0, 1, 2, . . . , ∞). Therefore, replacing the function f (x), equation (3.1), in the FPE (2.2), using the approach presented in the previous section, the proposed solution for this example is given by: with j being the maximum number of discrete states in the potential well region. The function g (t ) in this case is found by substituting the discrete functions of equation ( and N I (t ) is obtained by normalization, equation (2.10), It is assumed that, within the well, the truncation of the potential only slightly alters the original solutions of the harmonic oscillator. Thus, for the region between −d x d, the eigenfunctions [equation (3.3)] and eigenvalues (λ n = nk) are the same as for the harmonic potential. The only difference is that the number of terms of the series was limited taking into account the height of the potential well.

43003-5
(a) (b)  Observing the figure 2 (a), one can see that the probability distribution is greater in the region of minimum potential, even for extended periods of time. Since the given solution was constructed using an approximate method, one can see from figure 2 (b) that the probability distribution (3.4) does not completely satisfy the Fokker-Planck equation, in other words, ε(x, t ) 0. However, it is noted that the calculated errors are small and decrease as time increases. The larger relative errors appear when the probability distribution is calculated within the region of the well potential.
In the vicinity of the interface points (±d ), the error of the solution shows a discontinuity. This discontinuity is expected since the potential behavior studied is composed by joining different functions, and their profile (figure 1) is not smooth for the whole curve.
Through the probability distribution P (x, t ), equation (3.4), the variation of the number of particles in the region of minimum potential can be calculated, equation (2.11). The curves in figure 3 show the variation in the number of particles in the region of the potential well for different values of the diffusion coefficient.
In the definition of the potential used, equation (3.2), the depth of the well (related to v 0 ) depends on the diffusion coefficient (Q) and the interface point (d ).Thus, to maintain the fixed value of v 0 (equal to 0.5) for different values of Q, it is necessary to change the value of d. Maintaining a fixed value v 0 ensures that, within the potential well, the number of eigenvalues j that are solutions of the Schrödinger type equation (2.2) is fixed. In the example discussed here, v 0 = 0.5, there is only one eigenvalue of this kind.
Since it is assumed that the diffusion coefficient is proportional to temperature [2], lower values of Q represent lower system temperatures. Figure 3 shows that the decrease in the population of the region of the potential well depends on the value of Q, i.e., it depends on the temperature.
Initially the number of particles in the region of the well has the maximum value and with the increase in time this number of particles decreases. This drop in the number of particles is more pronounced for larger values of Q. This behavior is expected and consistent with the behavior of a system subject to a non-confining potential.   (t = 10 4 ). In this figure, it can be seen that for small values of Q (typically Q lower than 0.1), the population is confined to the minimum potential region. On the other hand, for larger values of Q, the number of particles within the well of potential decreases to zero. This behavior shows a phase transition in which the particles remain in the well of potential at low temperatures, and at high temperatures the potential well becomes emptied.
Another example can be got by increasing the depth of the potential well for v 0 = k 2 d 2 /4Q − k/2 = 1.5 with k = 1.4 and d = 2.1. The solution of (3.4) under these conditions gives two terms of the series (λ 0 = 0 and λ 1 = k) and the probability distribution is represented in figure 5, along with the associated error ε(x, t ) [equation (2.12)].  Comparison of figures 2 (a) and 5 (a) shows that the curve of the probability distribution is smooth and the peak is more pronounced when the depth of the well is increased [ figure 5 (a). As previously discussed, the calculation of the error ε(x, t ) shows a discontinuity at the points (x = ±d). It can be seen that, as time increases, there is a decrease in the calculated error indicating that the proposed solution is best for longer times.
In figure 6, the evolution of the number of particles as a function of the diffusion coefficient Q for a long time, t = 10 4 is shown. Again, one can see a phase transition in the system. However, this transition is less sudden and its effects are noticeable at values of Q higher than in the previous case ( figure 4). This effect is a result of the increased depth of the well, which makes the particle escape more difficult.

A non-harmonic polynomial potential
As a second example, the force associated with the system is assumed to be given by where the constants are chosen to ensure the continuity of the potential. For the region corresponding to In general, for non-harmonic polynomial potentials, the Schrödinger equation (2.2) has no exact/analytical solution. However, it is possible to determine part of the solution (partially soluble potential [13, Approximate solution for Fokker-Planck equation 14]). In such cases, the approach introduced in this work can be used, approximating the solution for the original potential to that of the truncated potential and building P (x, t ) from function (2.9), that is, . In this case, just the ground state is determined. However, depending on the well depth more eigenfunctions are necessary. Then, one possibility to get around this problem is to use other approximate methods, for example, the variational method [15].
Considering the potential characteristics studied and the parameters used (a = 0.45, b = 1.75, Q = 1, d = ±1.37 and v 0 = 1), there is only one state in the potential well region. Thus, only the first term of the series, should be considered in the region −d x d: P (x, t ) = N I g (t )φ 0 (x) 2 e −t λ 0 . (4.5) Then, in this example, when the potential is not truncated, just the ground state solution is analytically determined. In this case the eigenvalue (λ 0 ) is equal to 0 and the function φ 2 0 (x) is the same adopted solution of stationary ground state of the Schrödinger type equation when the potential well is infinite. Thus, the function to be used in equation (4.5) is: (4.6) Therefore, the probability distribution for the force given by equation (4.1) is represented as, where the function g (t ) is given by (4.9) Figure 8 shows the probability distribution (4.7) and the associated error ε(x, t ), equations (2.12), for different values of time. The potential (4.2) has the constants a and b equal to 0.45 and 1.75, respectively.
For numerical calculations, the depth of the well was fixed as v 0 = 1, which implies a unique eigenvalue λ 0 = 0. The interface points between the regions are d = ±1.37 and the value used for the diffusion coefficient to construct the curves shown in figure 8 is Q = 1.
It can be seen from figure 8 (a) that the probability distribution P (x, t ) has a peak in the region of the potential well even for very long times. Initially, there is a very distinct peak probability (t = 1) and, as time passes, this peak decreases and there is an increase in the width of the curve P (x, t ) at its base.
The increased width of the probability distribution indicates the escape of the particles from the central region to the region of constant potential.
In the same way as for the harmonic case, the suggested solution is substituted in the Fokker-Planck equation to evaluate the error of the proposed method by determining the function ε(x, t ) [figure 8 (b)].
One can see that the approximate solution is better for large values of time than for shorter times. It is also noted that the largest error is in the region within the potential well and is smaller in the side regions where the potential is constant.

Conclusion
This paper presents an approximate analytical solution to the Fokker-Planck equation for partially confining potentials. The suggested solution corresponds to an adaptation of a previous proposal [8] from the same authors. Here, there is suggested the removal of the Gaussian function from the region of potential well, which permits a greater numerical accuracy of solution. This can be noted by using the expression ε(x, t ), equation (2.12).
In all the cases studied, following the initial condition P (x, 0) = δ(x) and with the values of the constants as given in the examples, the probability distribution has a peak in the central region of the po-43003-10 tential well. As the time increases, the curves P (x, t ) show a widening and a reduction in height. The calculation of ε(x, t ), equation (2.12), as a way of assessing the accuracy of the approximate solutions in each example, indicates that the solutions have smaller errors for longer values of time than for shorter times. The approach outlined above allows for the study of a large number of problems whose solution proves difficult or impossible to obtain by other methods. For example, the truncated potentials discussed here could not be handled by the procedure given in reference [9], since it is not possible to get the exact/analytical solutions of the associated Schrödinger type equation.
The suggested solution method has the advantage of being extended to classes of partially confining systems that do not have an exact analytical solution. Thus, an analytical expression, albeit approximate, of the probability distribution provides important information, allowing the study of a much larger number of systems. In addition to this, the use of the test function ε(x, t ) [equation (2.12)] permits a quantitative measure of the accuracy of the result.
Analyzing the first example studied, i.e., the truncated harmonic potential with different depths, a transition phase can be identified involving the escape of particles from the well region of the potential. At low temperatures, the particles are trapped, while for higher temperatures the particles can escape. This escape leads to the emptying of the well. These results are very reliable, since they are obtained for long periods of time, a condition at which the proposed method turns out to be more accurate.
As a final remark, one observes that the approach introduced here can be addressed to the well known problem of the diffusion controlled escaping from a potential well [16,17]. In this kind of problem, the calculation of the rate coefficients has a central importance [18] and the calculation developed in the present work can be used to compute these quantities. Particularly, the escape rate problem is hard to analyze when the system is trapped in a potential well which correspond to the only point of minimum in the potential [19]. In this context, the proposed function Y (t ,Q), equation (2.11), can be useful.