Composite Laguerre Pseudospectral Method for Fokker-Planck Equations

. A composite generalized Laguerre pseudospectral method for the nonlinear Fokker-Planck equations on the whole line is developed. Some composite generalized Laguerre interpolation approximation results are established. As an application, a composite Laguerre pseudospectral scheme is provided for the problems of the relaxation of fermion and boson gases. Convergence and stability of the scheme are proved. Numerical results show the efficiency of this approach and coincide well with theoretical analysis


Introduction
Fokker-Planck equations describe the evolution of stochastic systems, such as the erratic motions of small particles immersed in fluids, and the velocity distributions of fluid particles in turbulent flows.Several methods have been developed for the linear Fokker-Planck equations, such as combined Hermite spectral-finite difference method, domain decomposition spectral method, a semi-analytical iterative technique and etc, see [1,6,9,13,22].However, the nonlinear Fokker-Planck equations can better reflect nonlinear characteristics of the corresponding to physical problems.Applications could be found in various fields such as astrophysics, the physics of polymer fluids and particle beams, biophysics and population dynamics, see [7].So it is interesting to solve various nonlinear problems, such as the nonlinear Fokker Planck equations, see [2,4,10,11,16,20].
Let R = { v | −∞ < v < ∞} be the velocity of particles.Denote by W (v, t) the probability density.Moreover, W 0 (v) is the initial state.For simplicity, let ∂ z W = ∂W ∂z , etc.The nonlinear Fokker-Planck equations for the relaxation of fermion and boson gases are given as follows (cf.[7,14]): where k = 1 for bosons and k = −1 for fermions.These models have been introduced as a simplification with respect to Boltzmann-based models as in [5,15].The entropy method was applied for quantifying explicitly the exponential decay towards Fermi-Dirac and Bose-Einstein distributions in the one-dimensional case, see [3].Further more, some numerical methods were developed for the problems (1.1).For example, the full-discrete generalized Hermite spectral and pseudo-spectral schemes were proposed in [4,21].However, in the stability and convergence analysis of the fully discrete pseudospectral scheme, the aliasing error brings a certain difficulty.Also, a composite Laguerre-Legendre pseudospectral scheme is presented in [19].Yet it requires more basis functions, which makes the calculation more complex.Recently, Wang [20] considered the composite Laguerre spectral method for the problems (1.1), in which the nonlinear term ∂ v (vW (v, t)(1 + kW (v, t))) exacerbates the difficulty of calculating the quadratures over the whole line.So we prefer composite generalized Laguerre interpolation method that only predicates on estimated values of unknown functions on the interpolation nodes, and handles nonlinear terms easily [18,24,25,27], which preserves the continuity and possesses the global spectral accuracy on the whole line [8,9,12,17,22,26].This paper is focused on developing a semidiscrete composite generalized Laguerre pseudospectral method for problem (1.1).The method needs fewer basis functions and avoids the aliasing errors analysis caused by the second order difference term in [21].To cope with the stationary solution whose value decays exponentially as |v| → ∞ in [3], we take the generalized Laguerre functions as the base functions, which is a complete L 2 (R)-orthogonal system in [9].Also, we can greatly ameliorate the accuracy of numerical errors by selecting the scaling factor involved in the base functions.Moreover, the numerical analysis is simplified and the algorithm scheme has better stability by using the basis function of weight function χ(v) ≡ 1.
This paper is organized as follows.In Section 2, we establish some results on the composite Laguerre interpolation approximation, which are pivotal to the error analysis of pseudospectral methods for various differential equations on the whole line.Then we construct a semidiscrete composite Laguerre pseudospectral scheme for (1.1) and present some numerical results to show the high accuracy of the proposed algorithm in Section 3. We prove its convergence and stability in Section 4. The final section is for some concluding remarks.

Preliminaries
In this section, we establish some basic results on the composite Laguerre-Gauss-Radau interpolations.
In particular, H 0 χ (R + ) = L 2 χ (R + ), with the inner product (u, w) χ,R + and the norm ∥u∥ χ,R + .We omit the subscript χ in the notations when χ(v) ≡ 1. Let The generalized Laguerre polynomial of degree l is defined by We denote by P N (R + ) the set of all polynomials of degree at most N .Let R+ = R + ∪ {0} and ξ N +1 (v), which are arranged in ascending order.Denote by ω We introduce the following discrete inner product and norm (cf.[22]): By the exactness of (2.1), By (2.4) of [22], we have Let R+ be the same as before.For any u ∈ C( R+ ), the generalized Laguerre-Gauss-Radau interpolation I R,N,α,β,R + u ∈ P N ( R+ ) is determined by To design proper pseudospectral method for the Fokker-Planck equation and many other similar problems, we shall use the orthogonal system of generalized Laguerre functions, defined by We now consider the new generalized Laguerre-Gauss-Radau interpolation corresponding to the weight function R,N,R + ,j , 0 ≤ j ≤ N. We introduce the following discrete inner product and norm (cf.[22]): In particular, (cf.(2.9) of [22]) Moreover, for any ϕ ∈ Q N +1,β (R + ), we have (cf.(2.10) of [22]) For pseudospectral method of nonlinear problems with varying coefficient, we need the following result.Lemma 1.For any ϕ ∈ Q N,β (R + ), there holds (2.3) Proof.Following the same line as in the proof of Lemma 2.1 of [19], we can obtain the desired result.⊓ ⊔ For any u ∈ C( R+ ), the generalized Laguerre-Gauss-Radau interpolation ĨR,N,α,β,R Furthermore, by the definitions of I R,N,α,β,R + and ĨR,N,α,β,R + , we have that (2.4) In particular, if |α| < 1, then the above result holds for any integer r ≥ 1.
We now consider the interpolation on the subdomain which are arranged in descending order.Denote by ω(α,β) R,N,R − ,j , 0 ≤ j ≤ N, the corresponding Christoffel numbers.For v ∈ R − , we introduce the discrete inner product and norm as follows (cf.[22]): 3), we have that ) ( In particular, if |α| < 1, then the above result holds for any integer r ≥ 1.

Composite generalized Laguerre interpolation on the whole
We are now in position of studying the composite generalized Laguerre-Gauss-Radau interpolation on the whole line R = R+ ∪ R− .The space L 2 χ (R) is defined as usual, with the inner product (u, w) χ,R and the norm ∥u∥ χ,R .We omit the subscript χ in the notations when χ(v) ≡ 1.Further, let We introduce the discrete inner product and norm as By virtue of (2.2)-(2.3)and (2.5)-(2.6),we have that ) ) ) Math.Model.Anal., 28(4):542-560, 2023.
In numerical analysis of the composite pseudospectral method for the Fokker-Planck equation on the whole line, we need a non-standard projection , which is defined by By (2.11) of [20], we have that if We need the following embedding inequality (cf.[21]).
By Lemma 2 and (2.12), we can get the estimate on L ∞ (R)-norm of P 1 N,β,R u.
(2.14) In the end of this section, we need some inverse inequality which will be used in the sequel (cf.[20]).

Composite pseudospectral method
In this section, we propose the composite pseudospectral method for the nonlinear Fokker-Planck equations on the whole line.We also describe the implementation and present some numerical results.

Pseudospectral scheme
Now, we deduce the pseudospectral scheme of (1.1), whose weak formulation is to seek W which belongs to the space L ∞ (0, T ; We now design the composite generalized Laguerre pseudospectral scheme for (3.1).It is to find Thanks to (2.8), the above problem is equivalent to Next, we describe the implementations for pseudospectral scheme (3.2) with a nonhomogeneous term f (v, t), we use the Crank-Nicolson discretization in time t, with the mesh size τ .
For simplicity of statements, we use the notation The fully discrete scheme of (3.2) is as follows: Then, at each time step, we need to solve the following nonlinear equation: Math.Model.Anal., 28(4):542-560, 2023.
For notational convenience, let L(β) The functions In actual computation, we expand the numerical solution as

Numerical results
In this subsection, we present some numerical results confirming the theoretical analysis.We use scheme (3.3) to solve (1.1) with k = 1.The numerical errors are measured by the discrete norm Now, we take the test function which decays exponentially at infinity,   In Figure 1, we plot the errors log 10 E N,τ (t) with t = 10 and β = 5.85 vs. √ N .Clearly, the errors decay fast when N increases and τ decreases.The above facts coincide very well with theoretical analysis in Theorem 1 on page 16.In particular, they show the spectral accuracy in the space of scheme (3.3).
In Figure 2, we plot log 10 E N (t) at t = 10, τ = 0.01 and different values of parameter β vs. √ N .It seems that the errors with suitably bigger β are smaller than those with smaller β.However, how to choose the best parameter β is still an open problem.Roughly speaking, if the exact solution decays faster as v increases, then it is better to take suitably bigger β.

Convergence and stability analysis of pseudospectral scheme
In this section, we consider the convergence and stability of scheme (3.2).

Convergence analysis
We next deal with the convergence of scheme (3.2).Let W N = P 1 N,β,R W .We derive from (3.1) that where Taking W N = w N − W N and subtracting (4.1) from (3.2), we obtain that where 2), we deduce that for 0 < t ≤ T , Therefore, it suffices to estimate the terms |G j (t, W N )|.Firstly, we use the Cauchy inequality and (2.12) to verify that for integers r ≥ 1, Obviously, According to (2.11) with α = 1, it is easy to derive that ), (2.14) with r = 1, (2.11) with α = 2 and (2.12), we have that A combination of the above two estimates gives that By (2.8) and (2.9) with α = 0, we deduce that Using (2.15), (2.16) with q = ∞ and (2.13), we derive that ), For fully big N , and r > 1, we have where d(W (t), β) is a positive constant, depend on ∂ r v (e ), where Integrating (4.6) with respect to t, we deduce that where λ = 1 for w N (0) = I N,0,β,R W 0 , and λ = 0 for w N (0) = P 1 N,β,R W 0 .We shall use the following lemma (cf.Lemma 3.1 of [11]).
Theorem 1.For 0 ≤ t ≤ T, we have , provided that the norms appearing in the previous statements are finite.

Stability analysis
We now consider the stability of scheme (3.2), which might be of the generalized stability as described in [10].Suppose that W 0 has the errors W0 .They induce the error of w N denoted by wN .Then, we obtain from (3.2) that for all (4.9) Taking ϕ = 2 wN (t) in (4.9), we derive that for 0 < t ≤ T , Substituting above inequality into (4.10),we deduce that Let Z(t) be the same as before (see the page 16 of the paper).We take in Lemma 4,

Conclusions
In this paper, we developed the composite generalized Laguerre pseudospectral method for the nonlinear Fokker-Planck equation on the whole line, which is distinguished from the methods as mentioned in the references in Section 1.The numerical results demonstrated spectral accuracy in space, and well confirmed the theoretical analysis.
The main advantages of the proposed approach are as follows: By using different generalized Laguerre interpolation approximations on different subdomains, we could deal with non-standard types of PDEs on the whole line properly.This trick also simplifies actual computations, especially for large modes N .
With the aid of composite generalized Laguerre interpolation approximations coupled with domain decomposition, we could exactly match the numerical solutions on the common boundary v = 0 of adjacent subdomains R + and R − , and the singularities of coefficients appearing in the underlying differential equations.Consequently, we could deal with the nonlinear Fokker-Planck equation on the whole domain properly.