Numerical simulation of Feller's diffusion equation

This article is devoted to Feller's diffusion equation which arises naturally in probabilities and physics (e.g. wave turbulence theory). If discretized naively, this equation may represent serious numerical difficulties since the diffusion coefficient is practically unbounded and most of its solutions are weakly divergent at the origin. In order to overcome these difficulties we reformulate this equation using some ideas from the Lagrangian fluid mechanics. This allows us to obtain a numerical scheme with a rather generous stability condition. Finally, the algorithm admits an elegant implementation and the corresponding Matlab code is provided with this article under an open source license.


Introduction
The celebrated Feller equation was introduced in two seminal papers published by William Feller (1951/1952 in Annals of Mathematics [10,11]. These publications studied mathematically (and, henceforth, gave the name) to the following equation * : where the subscripts t , x denote the partial derivatives, i.e. (·) t def := ∂(·) ∂t , (·) x def := ∂(·) ∂x . Two parameters γ and η > 0 can be time dependent in some physical applications, even if in this study we assume they are constants, for the sake of simplicity † . Equation (1.1) can be seen as the Fokker-Planck (or the forward Kolmogorov) equation with γ x being the drift and η x being the diffusion coefficients. One can notice also that Equation (1.1) becomes singular at x = 0 and x = + ∞ . We remind that practically important solutions to Feller's equation might be unbounded near x = 0 . In order to attempt at solving Equation (1.1), one has to prescribe an initial condition p (x, 0) = p 0 (x) presumably with a boundary condition at x = 0 . A popular choice is to prescribe the homogeneous boundary condition p (0, t) ≡ 0 . For this choice of the boundary condition it is not difficult to show that the Feller equation dynamics would preserve solution positivity provided that p 0 (x) 0 (see Appendix A for a proof). The solution norm is also preserved (see Appendix B). Moreover, Feller using the Laplace transform techniques has shown in [10] that the initial condition p 0 (x) determines uniquely the solution. In other words, no boundary condition at x = 0 should be prescribed. This conclusion might appear, perhaps, to be counter-intuitive.
The great interest in Feller's equation can be explained by its connection to Feller's processes, which can be described by the following stochastic differential Langevin equation ‡ : where W t is the standard Wiener process, i.e. ξ (t) def := dW t dt is zero-mean Gaussian white noise, i.e. ξ (t) = 0 , ξ (t) ξ (s) = δ (t − s) , * To be more accurate, W. Feller studied the following equation [10]: where a > 0 and 0 < x < + ∞ . † The numerical method we are going to propose can be straightforwardly generalized for this case when γ = γ (t) and η = η (t) . Moreover, Feller's processes with time-varying coefficients were studied recently in [19]. ‡ The stochastic differential equations is understood in the sense of Itō.
where the brackets · denote an ensemble averaging operator. Then, the Probability Density Function (PDF) p (x, t; x 0 ) of the process X (t) , i.e.
satisfies Equation (1.1) with the following initial condition [13]: The point x = 0 is a singular boundary that the process X (t) cannot cross. The Feller process is a continuous representation of branching and birth-death processes, which never attains negative values. This property makes it an ideal model not only in physical, but also in biological and social sciences [12,19,20]. As a general comprehensive reference on generalized Feller's equations we can mention the book [18]. Since at least a couple of years there is again a growing interest for studying equation (1.1). Some singular solutions to Feller's equation with constant coefficients were constructed in [12] via spectral decompositions. Feller's equation and Feller's processes with time-varying coefficients were studied analytically (always using the Laplace transform) and asymptotically in [19].
Recently, the Feller equation was derived in the context of the weakly interacting random waves dominated by four-wave interactions [6]. Wave Turbulence * (WT) is a common name for such processes [23]. In WT the Feller equation governs the PDF of squared Fourier wave amplitudes, i.e. x ∼ | a | 2 . In [6] some steady solutions to this equation with finite flux in the amplitude space were constructed † . See also [21,Chapter 11] for a detailed discussion and interpretations.
The behaviour of solutions p (x, t) for large x describes the appearance probability of extreme waves. In the context of ocean waves, these extreme events are known as rogue (or freak ) waves [9]. In the WT literature, any noticeable deviation from the Rayleigh distribution for x ≫ 1 is referred to as the anomalous probability distribution of large amplitude waves [6]. For Gaussian wave fields all statistical properties can be derived from the spectrum. However, the PDFs and other higher order moments are compulsory tools to study such deviations.
The present study focuses on the numerical discretization and simulation of Feller equation. The naive approach to solve this equation numerically encounters notorious difficulties. The first question, which arises is what is the (numerical) boundary condition to be imposed at x = 0 ? Moreover, one can notice that Equation (1.1) is posed on a semi-infinite domain. There are three main strategies to tackle this difficulty: (1) Map R + on a finite interval [ 0, ℓ ] (2) Use spectral expansions on R + (e.g. Laguerre or associated Laguerre polynomials) * We could define the Wave Turbulence (WT) as a physical and mathematical study of systems where random and coherent waves co-exist and interact [23]. † There is a misprint in [6, p. 366]. To obtain mathematically correct solutions one has to define n k def := η γ on the line below Equation (14).
In most studies the latter option is retained by imposing some appropriate boundary conditions at the artificial boundary x = L . In our study we shall propose a method, which is able to handle the semi-infinite domain R + without any truncations or simplifications. Finally, the diffusion coefficient in the Feller equation (1.1) is unbounded. If the domain is truncated at x = L , then the diffusion coefficient takes the maximal value ν max def := η L ≫ 1 , which depends on the truncation limit L and it can become very large in practice. We remind also that explicit schemes for diffusion equations are subject to the so-called Courant-Friedrichs-Lewy (CFL) stability conditions [7]: Taking into account the fact that ν max can be arbitrarily large, no explicit scheme can be usable with Feller equation in practice. Moreover, the dynamics of the Feller equation spreads over the space R + even localized initial conditions. In general, one can show that the support of p (x, t) , t > 0 is strictly larger * than the one of p (x, 0) . It is the so-called retention property. Thus, longer simulation times require larger domains. For all these reasons, it becomes clear that numerical discretization of the Feller equation requires special care.
In this study we demonstrate how to overcome this assertion as well. The main idea behind our study is to bring together PDEs and Fluid Mechanics. First, we observe that the classical Eulerian description is not suitable for this equation, even if the problem is initially formulated in the Eulerian setting. Consequently, the Feller equation will be recast in special material or the so-called Lagrangian variables † , which make the resolution easier and naturally adaptive [14,Chapter 7].
The present manuscript is organized as follows. The symmetry analysis of Equation (1.1) is performed in Section 2. Then, the governing equation is reformulated in Lagrangian variables in Section 3. The numerical results are presented in Section 4. Finally, the main conclusions and perspectives are outlined in Section 5.

Symmetry analysis
In general, a linear PDE admits an infinity of conservation laws, with integrating multipliers being solutions to the adjoint PDE [2]. Here we provide an interesting conservation law, which was found using the GeM Maple package [5]: , * Using modern analytical techniques it is possible to show even sharper results on the solution support, see e.g. [4].
† It is known that both Eulerian and Lagrangian descriptions were proposed by the same person -Leonhard Euler.

/ 21
where dt is the so-called exponential integral function [1] and the flux G is defined as The symmetry group of point transformations can be computed using GeM package as well. The infinitesimal generators are given below: where c ∈ R , M (a, b, z) and U (a, b, z) are Kummer special functions [1,17] (see also Appendix C). The corresponding point transformations, which map solutions of (1.1) into other solutions, can be readily obtained by integrating several ODE systems (we do not provide integration details here): The first symmetry is the time translation. The second one is the scaling of the dependent variable (the governing equation is linear). Symmetries 3 and 4 are exponential scalings. Two last symmetries express the fact that we can always add to the solution a particular solution to the homogeneous equation to obtain another solution. For instance, the solutions invariant under time translations (ξ 1 ) are steady states and their general form is the following: where C 1, 2 are 'arbitrary' constants, which have to be determined from imposed conditions. Of course, they should be chosen so that the resulting steady solution is a valid probability distribution. It is not difficult to check that the imposed flux F on the steady state solution is equal to C 1 η . Please, notice also an important property of the exponential integral function, which is useful in manipulating its values for negative arguments: where Ei (z) is the following exponential integral: The last relation can be also extended to the entire complex plain: We provide here also the general solutions invariant under the symmetry (ξ 3 ) : These solutions might be used, for example, to validate numerical codes.

Reformulation
By following the lines of [14,Chapter 7], we are going to rewrite Feller's Equation (1.1) with the so-called Lagrangian or material variables. The main advantage of this formulation consists in the fact that we can handle infinite domains without any truncations, transformations, etc. It becomes possible to carry computations in infinite domains. Our domain is semi-infinite (x ∈ R + ) with the left boundary x = 0 being a reflection point.
As the first step, we introduce the distribution function associated to the probability density p (x, t) : The same can be done for the initial condition as well: We notice also two obvious properties of the function P (x, t) : Due to the positivity preservation property (see Appendix A) the function P (x, t) is non-decreasing in variable x . Thus, we can define its pseudo-inverse * : which is defined as Similarly, the initial condition does possess a pseudo-inverse as well: such that X (P, 0) ≡ X 0 (P) . If Feller equation (1.1) holds in the sense of distributions, then the following equation holds as well: along with the initial condition

Equation (3.3) can be readily obtained by exploiting the obvious property
Thus, the implicit function theorem [24,25] guarantees the existence of derivatives of the inverse mapping X (P, t) . Let us compute them by differentiating with respect toP and t the following obvious identity: P (X (P, t), t) ≡P .
Thus, one can easily show that Using these expressions of partial derivatives, we derive the following evolution equation for the inverse mapping X (P, ·) : The last equation can be rewritten also by introducing a new dynamic variable: It is not difficult to see that Equation (3.4) becomes: The last equation will be solved numerically in the following Section. * This mapping is sometimes called the reciprocal mapping [14] or an order preserving string [3].

Numerical discretization
Earlier we derived Equation (3.4), which governs the dynamics of the pseudo-inverse mapping X (P, ·) . The initial condition for Equation (3.4) is given by the pseudo-inverse (3.2) of the initial condition P 0 (x) . We discretize Equation (3.4) with an explicit discretization in time since it yields the most straightforward implementation.
The first step in our algorithm consists in choosing the initial sampling interval. We make this choice depending on the provided initial condition. Typically, we want to sample only where it is needed. Thus, it seems reasonable to choose the initial segment [ 0, ℓ 0 ] with ℓ 0 being the leftmost location such that In simulations presented below we chose tol ∼ O (10 −5 ) . Then, we choose the initial sampling X 0 It is desirable that the initial sampling be adapted to the initial condition, since errors made initially cannot be corrected later. We define also P k = P 0 (X 0 k ) . We stress out that P k N k = 0 stand for a discrete cumulative mass variable and, thus, they are time independent. More generally, we introduce the following notation: with t n def := n ∆t , n ∈ N and ∆t > 0 is a chosen time step * . We introduce also similar notation for the dynamic variable: Now we can state the fully discrete scheme for Equation (3.6): with n 0 , k = 0, 1, . . . , N − 1 and The quantity ∆P k can be defined as the arithmetic or geometric mean of two neighbouring discretization steps ∆P k ± 1 2 : . * We present our algorithm with a constant time step for the sake of simplicity. However, in realistic simulations presented in Section 4 the time step will be chosen adaptively and automatically to meet the stability and accuracy requirements prescribed by the user.
To be specific, in our code we implemented the arithmetic mean. The fully discrete scheme can be easily rewritten under the form of a discrete dynamical system: Remark 1. We would like to say a few words about the implementation of boundary conditions. First of all, no boundary condition is required on the left side, where X n 0 = Y n 0 ≡ 0 . On the right boundary we prefer to impose the homogeneous Neumann-type boundary condition, which yields the exact 'mass' conservation at the discrete level as well. Namely, at the rightmost cell we have the following fully discrete scheme: . As a result, we obtain the exact conservation of 'mass' at the discrete level: To summarize, our numerical strategy consists in: (1) We compute the pseudo-inverse of the initial data p 0 (x) to obtain X (P, 0) ≡ Y (P, 0) . (2) This initial condition Y (P, 0) is evolved in (discrete) time using an explicit marching scheme in order to obtain numerical approximation to Y (P, t) , t > 0 . (3) The variable X (P, t) is recovered by inverting (3.5), i.e. X (P, t) = e −γ t Y (P, t) .
(4) Thanks to (3.1) we can deduce the values of p X (P, t), t ∈ [ 0, 1 ] by applying a favourite finite difference formula * .
Working with the pseudo-inverse allows to overcome the issue of the retention phenomenon, which manifests as the expanding support of p (x, t) for positive (and possibly large) times t > 0 , t ≫ 1 , since the computational domain was transformed to [ 0, 1 ] . This method is the Lagrangian counterpart of the moving mesh technique in the Eulerian setting [15,16]. A simple Matlab code, which implements the scheme we described hereinabove, is freely available for reader's convenience at this URL address: https://github.com/dutykh/Feller/ * In our code we employed the simplest forward finite differences and it lead satisfactory results. This point can be easily improved when necessary.

Numerical results
In this Section we validate and illustrate the application of the proposed algorithm on several examples. However, first we begin with a straightforward validation test. The only difference with the proposed above algorithm is that we are using a higher order adaptive time stepping for our practical simulations. The explicit first order scheme was used to simplify the presentation. In practice, much more sophisticated time steppers can be used. For instance, we shall employ the explicit embedded Dormand-Prince Runge-Kutta pair (4, 5) [8] implemented in Matlab under the ode45 routine [22]. Conceptually, this method is similar to the explicit Euler scheme presented above. It conserved all good properties we showed, but it provides additionally the higher accuracy order and totally automatic adaptivity of the time step, which matches very well with adaptive features of the Lagrangian scheme in space. The values of absolute and relative tolerances used in the time step choice are systematically reported below.

Steady state preservation
In order to validate the numerical algorithm we have at our disposal a family of steady state solutions (2.1). Hence, if we take such a solution as an initial condition, normally the algorithm has to keep it intact under the discretized dynamics. The parameters η , γ of the equation, those of the steady solution and numerical parameters used in our computation are reported in Table 1. The initial condition at t = 0 along with the final state at t = T are shown in Figure 1. Up to graphical resolution they coincide completely. We can easily check that during the whole simulation the points barely moved, i.e.
We can check also other quantities. For instance, P (ξ, t) is preserved up to the machine precision. If we reconstruct the probability distribution p (x, t) , we obtain: The last error comes essentially from the fact that we apply simple first order finite difference to reconstruct the variable p (x, t) from its primitive P (x, t) . We can improve this point, but even this simple reconstruction seems to be consistent with the overall scheme accuracy. Thus, this example shows that our implementation of the proposed algorithm is also practically well-balanced [14], since steady state solutions are well preserved.

Transient computations
In this Section we present a couple of extra truly unsteady computations in order to illustrate the capabilities of our method. Namely, we shall simulate the probability distributions emerging from a family of initial conditions (normalized to have the probability

Parameter Value
Drift coefficient, γ Final simulation time, T 10.0 Number of discretization points, N 500 Absolute tolerance, tol a 10 −5 Relative tolerance, tol r 10 −5 Table 1. Numerical parameters used in the steady state computations. distribution): The primitive of the last distribution can be easily computed as well:  We design two different experiments in silico to show two completely different behaviours of solutions to Feller equation (1.1) depending on the sign of the drift coefficient γ . These will constitute additional tests for the proposed numerical method. In both cases, the initial positions of particles are chosen according to the logarithmic distribution (logspace function in Matlab) on the segment [ 0, 20 ] . This choice is made to represent more accurately the exponentially decaying initial condition since the errors made in the initial condition cannot be corrected in the dynamics.

Expanding drift
If the drift coefficient γ < 0 , this term will cause Feller's equation (1.1) to transport information in the rightward direction. This situation is quite uncomfortable from the numerical point of view since the initial condition expands quickly towards the (positive) infinity. We submitted our method to this case. All numerical and physical parameters are provided in Table 2 (the middle column). The initial condition along with the probability distribution at the end of our simulation are shown simultaneously in Figure 2(a,b) in Cartesian and semi-logarithmic coordinates correspondingly. As expected, the support of the initial condition more than triples from t = 0 to t = T = 3.0 . We remind that the diffusion and drift coefficients are proportional to x and the scheme handles the growth of these terms automatically. The smooth decay of the solution on the semi-logarithmic plot (see Figure 2(b)) shows the absence of any numerical instabilities. The trajectory of grid nodes is shown in Figure 2(c). One can see that points follow the expansion of the solution. Nevertheless, they concentrate in the areas where the probability takes significant values.  Table 2 (middle column).

Confining drift
In the case of the positive drift coefficient γ > 0 , the Feller dynamics gets even more interesting since we have two competing effects: (1) Positive drift moving information towards the origin x = 0 , (2) Diffusion spreading solution towards the positive infinity.
It might happen that both effects balance each other and result in a steady solution. Such a simulation is presented in this Section. The numerical and physical parameters are given in Table 2 (the right column). The initial condition along with the probability distribution at the end of our simulation are shown simultaneously in Figure 3(a,b) in Cartesian and semi-logarithmic coordinates respectively. The trajectories of grid nodes are shown in Figure 2(c). One can see the rapid initial expansion stage, which is followed by a stabilization process, when we almost achieved the expected steady state. Again, the grid nodes trajectories show excellent adaptive features of the numerical scheme: at the end of the simulation the relative density of nodes turns out to be preserved. The semilogarithmic plot shown in Figure 3(b) shows that the numerical solution is free of numerical instabilities. The obtained steady solution will be preserved by discrete dynamics thanks to the well-balanced property demonstrated in Section 4.1.

Discussion
Above we presented some rationale about the discretization, existence and uniqueness theory for Feller's equation. The main conclusions and perspectives are outlined below.

Conclusions
The celebrated Feller equation was studied mathematically in two seminal papers published by William Feller (1951/1952 in Annals of Mathematics [10,11]. In particular, the uniqueness was shown in [11] without any boundary condition required at x = 0 . This result is notable and we use it in our study. The main goal of our work was to present an efficient numerical scheme, which were able to handle the unbounded diffusion inherent to Equation (1.1). Moreover, the retention phenomenon causes the solution support to expand all the time. Thus, if we wait sufficiently long time, it will reach the computational domain boundaries (since R + ∋ x was truncated to a finite segment to make the simulation in silico possible). To overcome this difficulty, we proposed a simple and explicit Lagrangian scheme using the notion of the pseudo-inverse. In this way, we obtained a stable numerical scheme (under not so stringent stability conditions), which turns out to be naturally adaptive as well, since particles move together with the growing support (the rightmost particle position defines the support) and they tend to concentrate where it is really needed. The performance of our scheme was illustrated on several examples. We share also the Matlab code so that anybody can check the claims we made hereinabove and use it to solve numerically the Feller equation in their applications: https://github.com/dutykh/Feller/   Table 2 (right column).

Perspectives
All the numerical schemes and results presented in this paper were in one 'spatial' dimension. The Feller equation considered here is 1−D as well. However, it seems interesting * to consider generalized Feller equations in higher space dimensions and to extend the proposed numerical strategy to the multi-dimensional case as well. Another possible extension direction consists in proposing higher order schemes to capture numerical solutions with better accuracy with the same number of degrees of freedom. On the more theoretical side, we would like to obtain an alternative well-posedness theory for Feller equation by taking a continuous limit in our numerical scheme.
Taking into account the fact that the point x ⋆ is a local minimum (p x | x ⋆ = 0), where the solution takes zero value (p | x ⋆ = 0), the last equation greatly simplifies at this point: since in the minimum p x x | x ⋆ > 0 . Thus, zero values of the solution are repulsive and for (at least small) times t > t ⋆ the function t → p (x ⋆ , t) will be increasing.

B. 'Mass' conservation
It is straightforward to show that the L 1 norm of the solution to Equations (1.1) is preserved. Indeed, taking into account that the solution p (x, t) > 0 is positive for all times t > 0 provided that the initial condition p 0 (x) 0 , ∀x ∈ R + , we have | p (x, t) | ≡ p (x, t) . By integrating Equation (1.1), we have Taking into account that F | x = 0 ≡ 0 and the solution p (x, t) is decaying sufficiently fast as x → + ∞ together with its derivative, we obtain ∂ tˆR + p (x, t) dx ≡ 0 .
The last constant can be in general taken equal to one after the appropriate rescaling (provided that the initial condition is non-trivial). It is this scaling, which is assumed throughout the whole text above. This renormalization is consistent with the 'physical sense' of the variable p being the density of a probability distribution.

C. Kummer functions
The Kummer functions M (a, b, z) and U (a, b, z) are two linearly independent solutions of the following ordinary differential equation: This equation admits two singular points: z = 0 (regular) and z = +∞ (irregular).
There exist a connection between Kummer and hyper-geometric functions [1].