# Conservative polynomial approximations and applications to Fokker-Planck equations

Tino Laidin, Lorenzo Pareschi

► To cite this version:

Tino Laidin, Lorenzo Pareschi. Conservative polynomial approximations and applications to Fokker-Planck equations. SIAM Journal on Scientific Computing, In press, 47 (2), pp.A1279-A1303. 10.1137/24M1643086. hal-04449712v2

HAL Id: hal-04449712

https://hal.science/hal-04449712v2

Submitted on 3 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
CONSERVATIVE POLYNOMIAL APPROXIMATIONS AND
APPLICATIONS TO FOKKER-PLANCK EQUATIONS

TINO LAIDIN* AND LORENZO PARESCHI†

**Abstract.** We address the problem of constructing approximations based on orthogonal polynomials that preserve an arbitrary set of moments of a given function without losing the spectral convergence property. To this aim, we compute the constrained polynomial of best approximation for a generic basis of orthogonal polynomials. The construction is entirely general and allows us to derive structure preserving numerical methods for partial differential equations that require the conservation of some moments of the solution, typically representing relevant physical quantities of the problem. These properties are essential to capture with high accuracy the long-time behavior of the solution. We illustrate with the aid of several numerical applications to Fokker-Planck equations the generality and the performances of the present approach.

**KEYWORDS:** Fokker-Planck equation, Galerkin spectral method, conservative methods, spectral accuracy.

2020 MATHEMATICS SUBJECT CLASSIFICATION: 65N35 (Primary), 82C40 (Secondary).

**1. Introduction.** Computational techniques that maintain certain fundamental properties or structures of the underlying mathematical problem in their discrete approximations play a crucial role in the study and analysis of ODEs and PDEs, offering insights into complex physical phenomena that are often inaccessible through analytical means alone [15, 18, 21]. These structure-preserving properties could include symmetries, conservation laws, or other structural characteristics that are crucial for accurately representing the behavior of the system being modeled.

Capturing the long-time behavior of solutions to PDEs is also intricately linked to the structure-preserving properties of numerical methods. In fact, preserving key structural properties at a discrete level, such as conservation laws and invariant quantities, ensures that the numerical approximation retains the essential characteristics of the original problem over extended time intervals.

In this regard, the preservation of moments of the solution stands out as a fundamental aspect, especially when considering kinetic and mean-field equations, where the moments represent the macroscopic observable physical quantities [3, 7, 23, 34]. In kinetic equations, for example, the long time behavior of the systems is fully determined by the knowledge of some invariant moments. Achieving such properties, however, is particularly challenging due to the inherent complexity of many physical systems and the necessity to accurately capture these dynamics.

A substantial body of literature has been devoted to the development of structure-preserving methods for various types of PDEs, with a focus on preserving specific physical properties. For instance, in the context of Fokker-Planck equations, numerous studies have explored different numerical strategies aimed at maintaining conservation of mass, momentum and energy with the aim to describe accurately the steady state solution of the problem [2, 4, 6, 24, 26, 34].

Spectral methods, which rely on expansions in terms of orthogonal polynomials, have garnered significant attention for their effectiveness in solving collisional kinetic equations of Boltzmann-type. The spectral accuracy and efficiency of these methods

<footnote>*Univ. Lille, CNRS, Inria, UMR 8524 - Laboratoire Paul Painlevé, F-59000 Lille, France (tino.laidin@univ-lille.fr).</footnote>

<footnote>†Maxwell Institute for Mathematical Sciences and Department of Mathematics, Heriot-Watt University, Edinburgh, United Kingdom (l.pareschi@hw.ac.uk) and Department of Mathematics and Computer Science, University of Ferrara, Italy (lorenzo.pareschi@unife.it).</footnote>
make them particularly well-suited for capturing the non-equilibrium dynamics of
such systems [20, 27, 28, 31–33].

However, due to the lack of conservations, the long time behavior of such methods
may lead to accumulation of errors and their extension to a conservative setting
has represented a challenging task in numerical simulations. Recent advancements
in conservative spectral methods based on $L^2$-minimization frameworks have shown
promise in addressing these challenges [1, 12, 13, 30].

Building upon these developments, our main goal in this manuscript is to extend
the $L^2$-minimization framework to derive families of orthogonal polynomials capable
of preserving the moments of the solution, thereby enabling the construction of
spectrally accurate and efficient structure-preserving numerical methods. This exten-
sion leverages the well-established theory of orthogonal polynomials [5, 10, 14, 16, 22],
offering a robust framework for achieving high-accuracy approximations of PDEs
while preserving key physical properties. In particular, we have developed a general
framework for constructing constrained orthogonal polynomial approximations which
maintain the accuracy properties of the unconstrained polynomials for smooth solu-
tions. This framework enabled us to construct spectrally accurate moment-preserving
Galerkin-type approximations for several Fokker-Planck equations originating from
various fields, ranging from classical physics to social sciences and operations re-
search. Consequently, these equations are defined in different domains, both bounded
and unbounded, and possess different steady states, requiring the adoption of suitable
orthogonal polynomial bases.

Let us also remark, that the connection between orthogonal polynomials and prob-
ability theory to design uncertainty quantification methods further underscores the
importance of our approach, as classical families of continuous and discrete orthogonal
polynomials are intimately linked to probability distributions and our considerations
naturally generalize to such contexts [36].

In the following sections, we provide a detailed exposition of our methodology
and present numerical results demonstrating the efficacy of our approach. The rest of
the manuscript is organized as follows. In Section 2, we introduce some notations and
present our moment-preserving approach based on a constrained $L_2$-minimization set-
ting. We also study the convergence properties and prove a general result on spectral
accuracy. Section 3 is then devoted to testing the novel polynomial approximation,
first by approximating probability densities with given moments and subsequently by
considering different Fokker-Planck equations on bounded and unbounded domains.
Some concluding remarks are provided in the last section.

**2. Conservative approximations by orthogonal polynomials.** In order to
get conservative approximations, we construct a conservative projection on the space
of orthogonal polynomials by generalizing the constrained formulation approach in-
troduced in [30] for trigonometric polynomials. In particular, we will give an explicit
formulation of the orthogonal polynomial of best approximation in the weighted least
square sense, constrained by preservation of moments, and show that it preserves a
spectral accuracy property for smooth solutions like the classical orthogonal polyno-
mial approximation.

**2.1. Orthogonal polynomials.** Let us first set up the mathematical framework
of our analysis. To simplify our treatment we will consider the one-dimensional case.
Given a real function $f(x) \in L_\omega^2(\Omega)$, with $\Omega \subseteq \mathbb{R}$ and $\omega(x) > 0$ a positive function,
we define

$$ (2.1) \qquad \|f\|_{L^2_\omega} = \left( \int_{\Omega} |f(x)|^2 \omega(x) dx \right)^{1/2}, $$

which has the associated inner product

$$ (2.2) \qquad \langle f, g \rangle_{\omega} = \int_{\Omega} f(x)g(x)\omega(x) dx. $$

We assume that $\omega(x) > 0$ is such that

$$ (2.3) \qquad \int_{\Omega} |x|^q \omega(x) dx < \infty, \quad q = 0, 1, 2, \dots $$

We consider the problem of approximating functions in the weighted $L^2_\omega(\Omega)$ space defined by the norm (2.1). In the case $\omega(x) = 1$ we will use standard notations $\|\cdot\|_{L^2}$ and $\langle \cdot, \cdot \rangle$ to denote the $L^2$ norm and the associated inner product.

A polynomial basis $\{p_k\}$, $k = 0, 1, \dots$ for $L^2_\omega(\Omega)$ is a set of functions such that any $f$ in the space can be expressed uniquely as

$$ (2.4) \qquad f(x) = \sum_{k=0}^{\infty} \hat{f}_k p_k(x). $$

A set of polynomials $\{p_k\}$ is called orthogonal if $\langle p_h, p_k \rangle_\omega = 0$ for $h \neq k$. We refer to Table 1 for some examples of standard polynomial families

Orthogonal bases are particularly nice, both for theory and numerical approximation. In fact, from

$$ \langle f, p_k \rangle_\omega = \sum_{h=0}^{\infty} \hat{f}_h \langle p_h, p_k \rangle_\omega = \hat{f}_k \langle p_k, p_k \rangle_\omega, $$

the coefficients in (2.4) are easily represented as

$$ (2.5) \qquad \hat{f}_k = \frac{\langle f, p_k \rangle_\omega}{\langle p_k, p_k \rangle_\omega} = \frac{\langle f, p_k \rangle_\omega}{\|p_k\|_{L^2_\omega}^2}. $$

A first result that we recall is that the truncated approximation

$$ (2.6) \qquad f_N(x) = \sum_{k=0}^{N} \hat{f}_k p_k(x) $$

is the best approximation of $f$ in the subspace $S_N = \text{span}\{p_k | k = 0, 1, \dots, N\}$. More precisely, let $\mathcal{P}_{\omega,N}: L^2_\omega(\Omega) \to S_N$ be the orthogonal projection upon $S_N$ in the inner product of $L^2_\omega(\Omega)$

$$ (2.7) \qquad \langle f - \mathcal{P}_{\omega,N} f, p_k \rangle_\omega = 0, \quad \forall p_k \in S_N. $$

With these definitions $\mathcal{P}_{\omega,N}f = f_N$, is the solution of the following minimization problem

$$ f_N = \underset{\|g_N - f\|_{L^2_\omega} : g_N \in S_N}{\operatorname{argmin}} \{ \|g_N - f\|_{L^2_\omega} : g_N \in S_N \}. $$
By orthogonality, one also obtains the Parseval’s identities

$$
\|f\|_{L_{\omega}^2}^2 = \sum_{k=0}^{\infty} |\hat{f}_k|^2 \|p_k\|_{L_{\omega}^2}^2, \quad \|f_N\|_{L_{\omega}^2}^2 = \sum_{k=0}^{N} |\hat{f}_k|^2 \|p_k\|_{L_{\omega}^2}^2,
$$

which gives the error estimate

$$
\|f - f_N\|_{L_{\omega}^2}^2 = \sum_{k=N+1}^{\infty} |\hat{f}_k|^2 \|p_k\|_{L_{\omega}^2}^2.
$$

An important feature of the orthogonal polynomial approximations on $S_N$ is
related to their spectral convergence properties for smooth solutions. We report here
a result for Jacobi type approximations [10]

THEOREM 2.1. If $f \in H_{\omega}^{r}(\Omega)$, where $r \ge 0$ is an integer, then there exists a constant $C > 0$ dependent on $\alpha, \beta$ and $r$ such that

$$
(2.8) \quad \|f - f_N\|_{L_{\omega}^2} \le \frac{C}{N^r} \left\| (1-x^2)^{r/2} \frac{d^r f}{dx^r} \right\|_{L_{\omega}^2}, \quad N > r.
$$

<table><caption>TABLE 1
Families of orthogonal polynomials</caption>
  <thead>
    <tr>
      <th>Name</th>
      <th>&omega;(x)</th>
      <th>&Omega;</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>Jacobi</td>
      <td>(1 &minus; x)<sup>&alpha;</sup>(1 + x)<sup>&beta;</sup>, &alpha;, &beta; &gt; &minus;1</td>
      <td>&minus;1, 1</td>
    </tr>
    <tr>
      <td>&minus; Legendre</td>
      <td>1, &alpha; = &beta; = 0</td>
      <td>&minus;1, 1</td>
    </tr>
    <tr>
      <td>&minus; Chebyshev 1st kind</td>
      <td>
        <div style="display: inline-block; vertical-align: middle; text-align: center; line-height: 1.2em;">
          <div style="border-bottom: 1px solid black; padding: 0 0.2em;">1</div>
          <div style="padding: 0 0.2em;">&radic;(1 &minus; x<sup>2</sup>)</div>
        </div>, &alpha; = &beta; = &minus;
        <div style="display: inline-block; vertical-align: middle; text-align: center; line-height: 1.2em;">
          <div style="border-bottom: 1px solid black; padding: 0 0.2em;">1</div>
          <div style="padding: 0 0.2em;">2</div>
        </div>
      </td>
      <td>&minus;1, 1</td>
    </tr>
    <tr>
      <td>&minus; Chebyshev 2nd kind</td>
      <td>&radic;(1 &minus; x<sup>2</sup>), &alpha; = &beta; =
        <div style="display: inline-block; vertical-align: middle; text-align: center; line-height: 1.2em;">
          <div style="border-bottom: 1px solid black; padding: 0 0.2em;">1</div>
          <div style="padding: 0 0.2em;">2</div>
        </div>
      </td>
      <td>&minus;1, 1</td>
    </tr>
    <tr>
      <td>Laguerre</td>
      <td>x<sup>&alpha;</sup>e<sup>&minus;x</sup>, &alpha; &gt; &minus;1</td>
      <td>&#8477;<sup>+</sup></td>
    </tr>
    <tr>
      <td>Hermite</td>
      <td>|x|<sup>2&alpha;</sup>e<sup>&minus;x<sup>2</sup></sup>, &alpha; &gt; &minus;1</td>
      <td>&#8477;</td>
    </tr>
  </tbody>
</table>

**2.2. Approximation of moments.** Next, we discuss some properties of the projection operator $\mathcal{P}_{\omega,N}$, in particular those concerning approximation of moments of nonnegative functions. We will denote for $q = 0, 1, 2, \dots$

$$
(2.9) \quad m_q = \int_{\Omega} f(x) x^q dx = \langle f, x^q \rangle, \quad m_{q,N} = \int_{\Omega} f_N(x) x^q dx = \langle f_N, x^q \rangle.
$$

Note that

$$
m_q = \sum_{k=0}^{\infty} \hat{f}_k \mu_{q,k}, \quad m_{q,N} = \sum_{k=0}^{N} \hat{f}_k \mu_{q,k},
$$

where for $k, q \ge 0$

$$
(2.10) \qquad \mu_{q,k} = \langle p_k, x^q \rangle,
$$

is the *q*-th moment of the *k*-th order polynomial.
It is well-known that moments of the weight function $\omega(x)$ permit to generate explicitly the orthogonal polynomials [22]. In our case, however, we are interested in the behavior of the moments of the orthogonal polynomials themselves defined in (2.10). An important characterization of monic orthogonal polynomials is the classical three-term recurrence relation [10, 14, 16, 22]

$$ (2.11) \qquad p_k(x) = (a_k x + b_k) p_{k-1}(x) + c_k p_{k-2}(x), \quad k \ge 1, $$

where we assumed $p_0(x) = 1$, $p_{-1}(x) = 0$ and the sequences $\{a_k\}_{k \ge 1}$, $\{b_k\}_{k \ge 1}$ and $\{c_k\}_{k \ge 1}$ depend on the particular polynomial family under consideration.

Equation (2.11) permits to recursively compute the value of $p_k(x)$ starting from the values of $p_0(x)$ and $p_1(x)$.

From (2.11) one gets the following relations for the moments coefficients

$$ (2.12) \qquad \mu_{q,k} = a_k \mu_{q+1,k-1} + b_k \mu_{q,k-1} + c_k \mu_{q,k-2}, \quad k \ge 1, \quad q \ge 0, $$

with

$$ \mu_{q,0} = \frac{x^{q+1}}{q+1}\Big|_{\Omega}, \quad \mu_{q,-1} = 0, \quad q \ge 0, $$

which, therefore, can also be computed recursively from (2.12) up to the moment $q$ provided we know $q+h$ moments of $p_{k-h}(x)$ for $h = 1, \dots, k$.

**2.2.1. Unbounded domains.** Note, however, that for unbounded domains $\Omega$, like the case of Laguerre or Hermite polynomials (see Table 1), the moments $\mu_{q,0}$, $q \ge 0$, are not finite. In such cases, one resorts on suitable orthogonal function families with respect to the inner product $\langle \cdot, \cdot \rangle$ defined in the symmetrically weighted case [3, 19, 25] as

$$ P_k(x) = p_k(x) w^{1/2}(x), $$

so that $\langle P_h, P_k \rangle = \langle p_h, p_k \rangle_{\omega}, \forall h, k \ge 0$.

For functions that we write in the form $f(x) = h(x)w^{1/2}(x)$ we have

$$ (2.13) \qquad f(x) = \sum_{k=0}^{\infty} \hat{h}_k P_k(x), $$

with

$$ (2.14) \qquad \hat{h}_k = \frac{\langle f, P_k \rangle}{\langle P_k, P_k \rangle} = \frac{\langle h, p_k \rangle_{\omega}}{\|p_k\|_{L_{\omega}^2}^2}. $$

Clearly, the orthogonal functions $P_k$ satisfy an analogous three term relation as (2.11) with $P_0(x) = w^{1/2}(x)$, and similarly their moments satisfy (2.12) with

$$ \langle P_0, x^q \rangle = \frac{x^{q+1}}{q+1} w^{1/2}(x) \Big|_{\Omega}, \quad \mu_{q,-1} = 0, \quad q \ge 0. $$

The above quantity is bounded for Laguerre and Hermite polynomials thanks to (2.3).

We can then define $f_N(x)$ as the truncated approximation to the first $N+1$ terms of (2.13). We have

$$ (2.15) \qquad \langle f - f_N, P_k \rangle = 0, \quad \forall P_k = p_k w^{1/2}, p_k \in S_N. $$
In the weighted function family representation Parseval identity reads

$$ \|f\|_2 = \sum_{k=0}^{\infty} |\hat{h}_k|^2 \|P_k\|_2^2. $$

We report here a result for Hermite approximations [9, 10]

**THEOREM 2.2.** If $f \in H_{\omega}^r(\Omega)$, where $r \ge 0$ is an integer, then there exist a constant $C > 0$ and $r$ such that

$$ (2.16) \qquad \|f - f_N\|_{L_{\omega}^2} \le \frac{C}{N^{r/2}} \|f\|_{H_{\omega}^r}, \quad N > r. $$

Let us remark that, in general, moments are not preserved by the projection. Neither in standard nor in the weighted case. In fact, we have

$$ (2.17) \qquad m_q - m_{q,N} = \langle f - f_N, x^q \rangle = \left\langle f - f_N, \frac{x^q}{\omega(x)} \right\rangle_{\omega} $$

which, by (2.7), is equal to zero only if $x^q/\omega(x)$ belongs to $S_N$. A special case is represented by the Legendre polynomials which by construction preserve all moments of the approximated function $f(x)$ up to the order $q \le N$.

Finally, we have for each $\varphi \in L_{\omega}^2(\Omega)$ by Hölder's inequality

$$ (2.18) \qquad |\langle f, \varphi \rangle_{\omega} - \langle f_N, \varphi \rangle_{\omega}| \le \|\varphi\|_{L_{\omega}^2} \|f - f_N\|_{L_{\omega}^2}. $$

Therefore, the projection error on the moments decays faster than algebraically, for the weighted norms, when the solution is infinitely smooth. Although this is in general a guarantee of the accuracy of the methods, in practice the loss of conservations can lead to accumulations of error in the approximation of PDEs over long time scales, resulting in the determination of inaccurate steady states [15, 34].

**2.3. Moment preserving orthogonal polynomials.** Although the error on moments is spectrally small for smooth functions, in many applications moments represent structural properties of the PDE and are related to conservation of physical properties such as mass, momentum and energy. In such cases, the design of a numerical method that is structure preserving requires the construction of a projection operator that does not affect the relevant moments. To this end, we will construct a different projection operator on the space of polynomials, $P_{\omega,N}^c : L_{\omega}^2(\Omega) \to S_N$ such that it preserves a finite number of moments. More precisely, we require the projection to satisfy for $q = 0, 1, \dots, Q$

$$ \langle P_{\omega,N}^c f, x^q \rangle = \langle f, x^q \rangle, $$

but maintaining the convergence properties of the original orthogonal projection.

To this aim, it is natural to consider the following constrained best approximation problem

$$ (2.19) \qquad f_N^c = \underset{g_N \in S_N}{\operatorname{argmin}} \left\{ \|g_N - f\|_{L_{\omega}^2}^2 : g_N \in S_N, \langle g_N, x^q \rangle = \langle f, x^q \rangle, q = 0, 1, \dots, Q \right\}. $$

Now, since $g_N \in S_N$ we can represent it in the form

$$ g_N(x) = \sum_{k=0}^{N} \hat{g}_k p_k(x) $$
and then by Parseval's identity

$$ \|g_N - f\|_{L_\omega^2}^2 = \sum_{k=0}^{\infty} |\hat{g}_k - \hat{f}_k|^2 \|p_k\|_{L_\omega^2}^2, $$

where we assumed $\hat{g}_k = 0, k > N$.

We require that

$$ (2.20) \qquad \langle g_N, x^q \rangle = \sum_{k=0}^{N} \hat{g}_k \mu_{q,k} = \langle f, x^q \rangle, \quad q = 0, 1, \dots, Q $$

or, after introducing the vector of moments $\Phi = (1, x, x^2, \dots, x^Q)^T$, in vector form

$$ \langle g_N, \Phi \rangle = \sum_{k=0}^{N} \hat{g}_k \hat{\Phi}_k = \langle f, \Phi \rangle, $$

where $\hat{\Phi}_k = (\mu_{0,k}, \mu_{1,k}, \dots, \mu_{Q,k})^T$.

Let us now solve the minimization problem (2.19) using the Lagrange multiplier method. Let $\lambda \in \mathbb{R}^{Q+1}$ be the vector of Lagrange multipliers, we consider the objective function

$$ \mathcal{L}(\hat{g}, \lambda) = \sum_{k=0}^{\infty} |\hat{g}_k - \hat{f}_k|^2 \|p_k\|_{L_\omega^2}^2 + \lambda^T \left( \sum_{k=0}^{N} \hat{g}_k \hat{\Phi}_k - \langle f, \Phi \rangle \right) = 0, $$

with $\hat{g} \in \mathbb{R}^{N+1}$ the vector of coefficients $\hat{g}_k, k = 0, \dots, N$.

Stationary points are found by imposing

$$ \frac{\partial \mathcal{L}(\hat{g}, \lambda)}{\partial \hat{g}_k} = 0, \quad k = 0, 1, \dots, N; \qquad \frac{\partial \mathcal{L}(\hat{g}, \lambda)}{\partial \lambda_q} = 0, \quad q = 0, 1, \dots, Q. $$

From the first condition, one gets

$$ (2.21) \qquad 2(\hat{g}_k - \hat{f}_k) \|p_k\|_{L_\omega^2}^2 + \lambda^T \hat{\Phi}_k = 0, $$

whereas the second condition corresponds to (2.20).

Multiplying the above equation by $\hat{\Phi}_k / \|p_k\|_{L_\omega^2}^2$ and summing up over $k$ we can write

$$ 2 \sum_{k=0}^{N} (\hat{g}_k - \hat{f}_k) \hat{\Phi}_k + \sum_{k=0}^{N} \frac{1}{\|p_k\|_{L_\omega^2}^2} \hat{\Phi}_k \hat{\Phi}_k^T \lambda = 0 $$

and, using (2.20) as well as the fact that $\hat{\Phi}_k \hat{\Phi}_k^T$ are symmetric and positive definite matrices of size $Q+1$ one obtains

$$ (2.22) \qquad \lambda = -2 \left( \sum_{k=0}^{N} \frac{1}{\|p_k\|_{L_\omega^2}^2} \hat{\Phi}_k \hat{\Phi}_k^T \right)^{-1} (\langle f, \Phi \rangle - \langle f_N, \Phi \rangle). $$

Now, rewriting the first condition (2.21) as

$$ \hat{g}_k = \hat{f}_k - \frac{1}{2 \|p_k\|_{L_\omega^2}^2} \hat{\Phi}_k^T \lambda, $$

and plugging the expression (2.22) of $\lambda$ in it, one obtains that the minimum is achieved for $\hat{g}_k = \hat{f}_k^c$, given by the following definition.
DEFINITION 2.3. For $f \in L_{\omega}^2(\Omega)$, we define the conservative orthogonal projection $\mathcal{P}_{N,\omega}^c f = f_N^c$ in $S_N$, where $f_N^c$ is given by

$$ (2.23) \qquad f_N^c(x) = \sum_{k=0}^{N} \hat{f}_k^c p_k(x), $$

and we define the moment constrained coefficients as

$$ (2.24) \qquad \hat{f}_k^c = \hat{f}_k + \hat{C}_k^T (\langle f, \Phi \rangle - \langle f_N, \Phi \rangle), \qquad \hat{C}_k^T = \frac{1}{\|p_k\|_{L_\omega^2}^2} \hat{\Phi}_k^T \left( \sum_{h=0}^{N} \frac{1}{\|p_h\|_{L_\omega^2}^2} \hat{\Phi}_h \hat{\Phi}_h^T \right)^{-1}. $$

The following result states the spectral accuracy of the conservative best approximation in least square sense (2.23), and generalizes Theorem 2.1. In order to prove this result, one needs the following hypothesis:

**Hypothesis 2.4.** The moments of the approximated solution are spectrally accurate in the classical $L^2$-norm. Namely, there exists $C > 0$ such that

$$ (2.25) \qquad \|U - U_N\| \le \frac{C}{N^r} \|f\|_{H_{\omega}^r} \|\Phi\|_{2, L_{\omega}^2} $$

where $U$ and $U_N$ denote the moment vectors of $f$ and $f_N$ respectively:

$$ (2.26) \qquad U = \langle f, \Phi \rangle, \quad U_N = \langle f, \Phi \rangle. $$

This assumption is different from the spectral accuracy on the moment, in weighted norm, one could obtain from (2.18). This assumption is a consequence of a (classical) spectral accuracy in $L^2$ of the function $f$. Therefore, it is a little bit weaker than asking further regularity on $f$ in $L^2$ in addition to the one in weighted $L^2$. This property depends heavily on the function under consideration and we will see in numerical simulations that it may only be a technical detail.

**THEOREM 2.5.** If $f \in H_{\omega}^r(\Omega)$, where $r \ge 0$ is an integer, we have

$$ (2.27) \qquad \|f - f_N^c\|_{L_\omega^2} \le \frac{C_\Phi}{N^r} \|f\|_{H_\omega^r} $$

where the constant $C_\Phi$ depends on the spectral radius of the matrix

$$ M = \sum_{h=0}^{N} \frac{1}{\|p_h\|_{L_\omega^2}^2} \hat{\Phi}_h \hat{\Phi}_h^T, $$

and on $\|\Phi\|_{2, L_\omega^2}^2 = \sum_{q=0}^Q \|\Phi_q\|_{L_\omega^2}^2$ where $\Phi_q, q = 0, 1, \dots, Q$ are the components of the vector $\Phi$.

*Proof.* We can split

$$ \|f - f_N^c\|_{L_\omega^2}^2 \le \|f - f_N\|_{L_\omega^2}^2 + \|f_N - f_N^c\|_{L_\omega^2}^2. $$

The first term is bounded by the spectral estimate of truncated approximation of the form (2.8), whereas for the second term, by Parseval's identity, we have

$$ \|f_N^c - f_N\|_{L_\omega^2}^2 = \sum_{k=0}^{N} |\hat{f}_k^c - \hat{f}_k|^2 \|p_k\|_{L_\omega^2}^2. $$
Now, using the definition (2.24) we get

$$ \sum_{k=0}^{N} |\hat{f}_k^c - \hat{f}_k|^2 \|p_k\|_{L^2_\omega}^2 = \sum_{k=0}^{N} |\hat{C}_k^T (U - U_N)|^2 \|p_k\|_{L^2_\omega}^2 \le \|U - U_N\|_{2}^2 \sum_{k=0}^{N} \|\hat{C}_k^T\|_{2}^2 \|p_k\|_{L^2_\omega}^2. $$

where $\|\cdot\|_2$ denotes the euclidean norm of the vector. The assumption 2.4 of spectral accuracy, in classical $L^2$ norm, on the moments implies

$$ \|U - U_N\|_2 \le \frac{C}{N^r} \|f\|_{H^r_\omega} \left( \sum_{q=0}^{Q} \|\Phi_q\|_{L^2_\omega}^2 \right)^{1/2} = \frac{C}{N^r} \|f\|_{H^r_\omega} \|\Phi\|_{2, L^2_\omega} $$

where for $\Phi = (\Phi_0, \dots, \Phi_Q)^T$ we defined $\|\Phi\|_{2,L^2_\omega}^2 = \sum_{q=0}^Q \|\Phi_q\|_{L^2_\omega}^2$. Finally, by definition of $\hat{C}_k^T = \frac{1}{\|p_k\|_{L^2_\omega}^2} \hat{\Phi}_k^T \left( \sum_{h=0}^N \frac{1}{\|p_h\|_{L^2_\omega}^2} \hat{\Phi}_h \hat{\Phi}_h^T \right)^{-1}$, one has

$$ \sum_{k=0}^{N} \|\hat{C}_k^T\|_2^2 \|p_k\|_{L^2_\omega}^2 \le \sum_{k=0}^{N} \|\hat{\Phi}_k\|_2^2 \left\| \left( \sum_{h=0}^{N} \frac{1}{\|p_h\|_{L^2_\omega}^2} \hat{\Phi}_h \hat{\Phi}_h^T \right)^{-1} \right\|_2^2. $$

Now let $M = \sum_{h=0}^N \hat{\Phi}_h \hat{\Phi}_h^T / \|p_h\|_{L^2_\omega}^2$. It follows that

$$ (2.28) \qquad \sum_{k=0}^{N} \|\hat{C}_k^T\|_2^2 \|p_k\|_{L^2_\omega}^2 \le \|\Phi\|_{2, L^2_\omega}^2 \|M^{-1}\|_2^2. $$

As a sum of positive definite symmetric matrices, $M$ enjoys the same properties. Consequently,

$$ (2.29) \qquad \sum_{k=0}^{N} \|\hat{C}_k^T\|_2^2 \|p_k\|_{L^2_\omega}^2 \le \|\Phi\|_{2, L^2_\omega}^2 \rho^2(M^{-1}), $$

where $\rho(M^{-1})$ denotes the spectral radius of the matrix $M^{-1}$. Combining everything, one finally obtains

$$ \|f_N^c - f_N\|_{L^2_\omega} \le \frac{C}{N^r} \|f\|_{H^r_\omega} \|\Phi\|_{2, L^2_\omega} \rho(M^{-1}) $$

which proves (2.27). □

Remark 2.6. The conservative best approximation in least square (2.24) can be represented in term of the standard projection as

$$ (2.30) \qquad \mathcal{P}_N^c f = \mathcal{P}_N f + \sum_{k=0}^{N} \hat{C}_k^T \langle f - \mathcal{P}_N f, \Phi \rangle. $$

Remark 2.7. The result of Theorem 2.5 guarantees the spectral accuracy of the constrained approximation. It is however worth mentioning that the constant $C_\Phi$ from the error estimate involves the spectral radius of the matrix $M^{-1}$, which itself depends on the number of modes. Nonetheless, this dependency doesn't seem to significantly affect the convergence rate of the constrained method as shown through the various numerical experiments in Section 3. Note that the spectral radius of $M^{-1}$, which appears in $C_\Phi$, doesn't grow with $N$. It is more affected by a high number of constraints. We present in Figure 1 this spectral radius as a function of $N$ and $Q$ in the case of the Chebyshev 2nd kind polynomials.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Spectral radius of the matrix $M^{-1}$ as a function of N (Left), and as a function of Q (Right) for Chebyshev 2nd kind polynomials.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2. Constrained Legendre approximation for (3.1).</figcaption></figure>

### 3. Numerical examples and applications.
In this section, we conduct numerical experiments to illustrate the properties of the proposed method. The initial test focuses on evaluating the accuracy of the moment-preserving approach. Subsequently, in Tests 2 to 4, we employ the conservative method within a PDE framework, considering scenarios involving both bounded and unbounded domains. The ensuing analysis includes the assessment of the numerical scheme's accuracy, its conservation properties, and an exploration of the long-time behavior of solutions. In the following, all integrals are approximated using Gaussian quadrature rules associated to the chosen polynomial basis. The number of quadrature points is set to 80.

#### 3.1. Test 1: Approximation of functions.

*Bounded domain:* $[-1, 1]$. Let us consider the approximation of the oscillating and non symetric function

$$ (3.1) \qquad f(x) = \sin(2\pi x) + x^2 \cos(2\pi x), \quad x \in [-1, 1]. $$

We first consider a standard approximation using Legendre polynomials. Figure 2 illustrates the approximate solutions and Table 2 shows that the standard method is indeed spectrally accurate. In addition, we observe that moments are indeed preserved without constraining the projection as predicted by (2.17). Note that errors are
<table><caption>TABLE 2</caption><thead><tr><th>N</th><th>||f - f<sub>N</sub>||<sub>2</sub></th><th>|m<sub>0</sub> - m<sub>0,N</sub>|</th><th>|m<sub>1</sub> - m<sub>1,N</sub>|</th><th>|m<sub>2</sub> - m<sub>2,N</sub>|</th><th>|m<sub>3</sub> - m<sub>3,N</sub>|</th></tr></thead><tbody><tr><td>8</td><td>9.004e-03</td><td>2.429e-15</td><td>2.387e-15</td><td>2.387e-15</td><td>2.331e-15</td></tr><tr><td>16</td><td>2.197e-05</td><td>1.388e-17</td><td>0</td><td>0</td><td>0</td></tr><tr><td>32</td><td>1.050e-13</td><td>8.327e-17</td><td>1.110e-16</td><td>0</td><td>0</td></tr></tbody></table>

Error Analysis: Legendre approximation for (3.1).

Standard

<table><thead><tr><th>N</th><th>||f - f<sub>N</sub>||<sub>2</sub></th><th>|m<sub>0</sub> - m<sub>0,N</sub>|</th><th>|m<sub>1</sub> - m<sub>1,N</sub>|</th><th>|m<sub>2</sub> - m<sub>2,N</sub>|</th><th>|m<sub>3</sub> - m<sub>3,N</sub>|</th></tr></thead><tbody><tr><td>8</td><td>1.034e-01</td><td>1.746e-03</td><td>1.404e-02</td><td>1.958e-02</td><td>1.538e-02</td></tr><tr><td>16</td><td>2.477e-05</td><td>1.813e-07</td><td>6.126e-09</td><td>1.858e-07</td><td>6.264e-09</td></tr><tr><td>32</td><td>2.094e-14</td><td>1.388e-17</td><td>5.551e-17</td><td>2.776e-17</td><td>5.551e-17</td></tr></tbody></table>

Constrained

<table><caption>TABLE 3</caption><thead><tr><th>N</th><th>||f - f<sub>N</sub><sup>c</sup>||<sub>2</sub></th><th>|m<sub>0</sub> - m<sub>0,N</sub><sup>c</sup>|</th><th>|m<sub>1</sub> - m<sub>1,N</sub><sup>c</sup>|</th><th>|m<sub>2</sub> - m<sub>2,N</sub><sup>c</sup>|</th><th>|m<sub>3</sub> - m<sub>3,N</sub><sup>c</sup>|</th></tr></thead><tbody><tr><td>8</td><td>1.039e-01</td><td>5.551e-17</td><td>5.551e-17</td><td>0.</td><td>5.551e-17</td></tr><tr><td>16</td><td>2.477e-05</td><td>0.</td><td>0.</td><td>2.776e-17</td><td>0.</td></tr><tr><td>32</td><td>2.093e-14</td><td>2.776e-17</td><td>5.551e-17</td><td>2.776e-17</td><td>5.551e-17</td></tr></tbody></table>

Error Analysis: Standard vs. Conservative with Chebyshev 1st kind approximation of (3.1).

sometimes extremely close to 0 but not exactly 0 because the Gaussian quadrature used may introduces small machine epsilons.

Let us now compare the standard and constrained approximations for three polynomial families: Chebyshev 1st kind, Chebyshev 2nd kind, and general Jacobi with parameters $\alpha = 1$ and $\beta = -1/2$. The associated outcomes are illustrated in Figure 3 and errors are detailed in Tables 3, 4, and 5. Within the Jacobi family, the deliberate choice of $\alpha$ and $\beta$ serves the purpose of disrupting the basis symmetry to evaluate the method's robustness in this context. Across all three cases, it is clear that standard approximations exhibit spectral accuracy, with moments converging in a spectral manner. Consequently, our analysis falls within the range of Theorem 2.5, affirming spectral accuracy for constrained approximations and exact conservation of the first $M = 4$ moments. This property shows clearly in the Tables mentioned earlier.

It is noteworthy that the constraining matrix $\hat{C}_k$ may exhibit pronounced ill-conditioning for large values of $Q$, potentially introducing numerical artifacts. This behaviour is illustrated in Figure 4 where we plot the condition number of the matrix $M$ appearing in the proof of Theorem 2.5 as well as the minimum over $k = 1, \dots, N$ of the condition number of the matrices $\hat{C}_k$. In both cases, we observe large condition numbers as well as their increase with the number of constrained moments. In addition, we want to point out that enough modes should be considered when applying the constrained method. Indeed, a significant increase in the condition number occurs once the number of constraints surpasses the number of modes. Preconditioners may therefore be required to improve numerical stability for large numbers of constrained moments.

*Unbounded domain:* $\mathbb{R}$. Continuing, we extend our approximation to $\mathbb{R}$ utilizing Hermite functions. Let $H_k$ denote the $k$th Hermite polynomial. We define the $k$th symmetrically weighted Hermite functions as

$$ (3.2) \qquad \psi_k(x) = H_k(x)e^{-x^2/2}. $$
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Constrained approximations for Chebyshev 1st kind (Top Left), Chebyshev 2nd kind (Top Right) and Jacobi $\alpha = 1, \beta = -1/2$ (Bottom) for (3.1).</figcaption></figure>

<table><caption>TABLE 4
Error Analysis: Standard vs. Conservative with Chebyshev 2nd kind approximation of (3.1).</caption><thead><tr><th colspan="6">Standard</th></tr></thead><tbody><tr><th>N</th><th>||f - f<sub>N</sub>||<sub>2</sub></th><th>|m<sub>0</sub> - m<sub>0,N</sub>|</th><th>|m<sub>1</sub> - m<sub>1,N</sub>|</th><th>|m<sub>2</sub> - m<sub>2,N</sub>|</th><th>|m<sub>3</sub> - m<sub>3,N</sub>|</th></tr><tr><td>8</td><td>1.081e-01</td><td>8.262e-03</td><td>5.914e-03</td><td>8.499e-03</td><td>6.049e-03</td></tr><tr><td>16</td><td>2.701e-05</td><td>1.388e-06</td><td>4.889e-08</td><td>1.398e-06</td><td>4.921e-08</td></tr><tr><td>32</td><td>1.132e-14</td><td>8.327e-17</td><td>0</td><td>2.776e-17</td><td>5.551e-17</td></tr></tbody><thead><tr><th colspan="6">Constrained</th></tr></thead><tbody><tr><th>N</th><th>||f - f<sub>N</sub><sup>c</sup>||<sub>2</sub></th><th>|m<sub>0</sub> - m<sub>0,N</sub><sup>c</sup>|</th><th>|m<sub>1</sub> - m<sub>1,N</sub><sup>c</sup>|</th><th>|m<sub>2</sub> - m<sub>2,N</sub><sup>c</sup>|</th><th>|m<sub>3</sub> - m<sub>3,N</sub><sup>c</sup>|</th></tr><tr><td>8</td><td>1.029e-01</td><td>1.388e-17</td><td>0</td><td>0</td><td>0</td></tr><tr><td>16</td><td>2.629e-05</td><td>8.327e-17</td><td>0</td><td>2.776e-17</td><td>5.551e-17</td></tr><tr><td>32</td><td>1.135e-14</td><td>1.388e-17</td><td>0</td><td>0</td><td>5.551e-17</td></tr></tbody></table>

<table><caption>TABLE 5
Error Analysis: Standard vs. Conservative with Jacobi ($\alpha = 1, \beta = -1/2$) approximation of (3.1).</caption><thead><tr><th colspan="6">Standard</th></tr></thead><tbody><tr><th>N</th><th>||f - f<sub>N</sub>||<sub>2</sub></th><th>|m<sub>0</sub> - m<sub>0,N</sub>|</th><th>|m<sub>1</sub> - m<sub>1,N</sub>|</th><th>|m<sub>2</sub> - m<sub>2,N</sub>|</th><th>|m<sub>3</sub> - m<sub>3,N</sub>|</th></tr><tr><td>8</td><td>2.738e-01</td><td>4.374e-02</td><td>4.507e-02</td><td>0.437e-02</td><td>4.514e-02</td></tr><tr><td>16</td><td>8.873e-05</td><td>7.374e-06</td><td>7.433e-06</td><td>7.373e-06</td><td>7.434e-06</td></tr><tr><td>32</td><td>2.237e-12</td><td>8.263e-14</td><td>9.082e-14</td><td>8.271e-14</td><td>9.093e-14</td></tr></tbody><thead><tr><th colspan="6">Constrained</th></tr></thead><tbody><tr><th>N</th><th>||f - f<sub>N</sub><sup>c</sup>||<sub>2</sub></th><th>|m<sub>0</sub> - m<sub>0,N</sub><sup>c</sup>|</th><th>|m<sub>1</sub> - m<sub>1,N</sub><sup>c</sup>|</th><th>|m<sub>2</sub> - m<sub>2,N</sub><sup>c</sup>|</th><th>|m<sub>3</sub> - m<sub>3,N</sub><sup>c</sup>|</th></tr><tr><td>8</td><td>1.787e-01</td><td>6.939e-17</td><td>0</td><td>2.776e-17</td><td>5.551e-17</td></tr><tr><td>16</td><td>6.695e-05</td><td>1.110e-16</td><td>0</td><td>1.110e-16</td><td>5.551e-17</td></tr><tr><td>32</td><td>1.868e-12</td><td>1.388e-16</td><td>5.551e-17</td><td>5.551e-17</td><td>0</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Condition number of the matrix M (Left), and minimum of the condition number of the matrices $\tilde{C}_k$ (Right) as functions of Q for Chebyshev 2nd kind polynomials</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Constrained Hermite approximation for (3.3).</figcaption></figure>

Specifically, we want to approximate the difference between two Gaussian functions deliberately crafting a non-symmetric outcome:

$$ (3.3) \qquad f(x) = \frac{3}{\sqrt{2\pi}} e^{(x+3)/2} - \frac{1}{\sqrt{\pi}} e^{(x-2)}, \quad x \in \mathbb{R}. $$

We illustrate the results in Figure 5 and, again, the results (See Table 6) on the standard approximation ensures that we are within the scope of Theorem 2.5. As predicted, the conservative method exhibits spectral accuracy, along with conservation of moments.

Unbounded domain: $\mathbb{R}^+$. We now look at the approximation on $\mathbb{R}^+$ using symmetric Laguerre functions. We define for a Laguerre polynomial $\mathcal{L}_k$ the $k$th Laguerre function as

$$ (3.4) \qquad \xi_k(x) = \mathcal{L}_k(x)e^{-x/2}. $$

Let us consider the function

$$ (3.5) \qquad f(x) = (x^3 - 2x + \sin(x))e^{-x}, \quad x \in \mathbb{R}^+. $$
Standard

<table><thead><tr><th>N</th><th>||f - f<sub>N</sub>||<sub>2</sub></th><th>|m<sub>0</sub> - m<sub>0,N</sub>|</th><th>|m<sub>1</sub> - m<sub>1,N</sub>|</th><th>|m<sub>2</sub> - m<sub>2,N</sub>|</th><th>|m<sub>3</sub> - m<sub>3,N</sub>|</th></tr></thead><tbody><tr><td>8</td><td>1.264e-02</td><td>0.583e-01</td><td>1.718e00</td><td>1.130e+01</td><td>3.686e+01</td></tr><tr><td>16</td><td>1.716e-03</td><td>6.692e-03</td><td>2.060e-02</td><td>2.293e-01</td><td>7.457e-01</td></tr><tr><td>32</td><td>4.995e-09</td><td>1.698e-09</td><td>5.175e-09</td><td>9.470e-08</td><td>3.552e-07</td></tr></tbody></table>

Constrained

<table><caption>TABLE 6</caption><thead><tr><th>N</th><th>||f - f<sub>N</sub><sup>c</sup>||<sub>2</sub></th><th>|m<sub>0</sub> - m<sub>0,N</sub><sup>c</sup>|</th><th>|m<sub>1</sub> - m<sub>1,N</sub><sup>c</sup>|</th><th>|m<sub>2</sub> - m<sub>2,N</sub><sup>c</sup>|</th><th>|m<sub>3</sub> - m<sub>3,N</sub><sup>c</sup>|</th></tr></thead><tbody><tr><td>8</td><td>4.324e-01</td><td>0</td><td>1.776e-15</td><td>3.553e-15</td><td>2.842e-14</td></tr><tr><td>16</td><td>2.879e-03</td><td>2.220e-16</td><td>1.776e-15</td><td>0</td><td>0</td></tr><tr><td>32</td><td>5.055e-09</td><td>2.220e-16</td><td>1.776e-15</td><td>0</td><td>0</td></tr></tbody></table>

*Error Analysis: Standard vs. Conservative Hermite function approximation for (3.3).*

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Constrained Laguerre approximation for (3.5).</figcaption></figure>

The approximation is illustrated in Figure 6 and, consistent with prior cases, the observations in Table 7 shows the spectral accuracy of the conservative method. Although errors on the moments exhibit a slight increase compared to previous cases, it is noteworthy that this discrepancy can be attributed to numerical artifacts arising from the ill-conditioning of the constraining matrix $\hat{C}_k$ and numerical integration intricacies.

To conclude this first test let us again consider Laguerre functions but for the approximation of

$$ (3.6) \qquad f(x) = \frac{1}{\sqrt{2\pi\sigma x}} \exp\left(-\frac{(\ln(x) - \mu)^2}{2\sigma^2}\right), \quad x \in \mathbb{R}^{+,*}. $$

where $\sigma = 0.2$ and $\mu = \ln(40) - 0.2$. This function is badly approximated by negative exponentials at infinity, yet it holds significance for Test 4. As depicted in Figure 7, spectral accuracy is not attained neither on the standard approximation of (3.6) nor on its moments. Consequently, this experiment falls outside the assumptions of Theorem 2.5, since one lacks both spectral accuracy properties needed. Nevertheless, it is noteworthy that the resultant conservative approximation successfully preserves the 0th order moment and achieves fifth-order accuracy.

In the next tests, we apply the conservative method in a PDE setting to various Fokker-Planck equations.
Standard

<table><thead><tr><th>N</th><th>||f - f<sub>N</sub>||<sub>2</sub></th><th>|m<sub>0</sub> - m<sub>0,N</sub>|</th><th>|m<sub>1</sub> - m<sub>1,N</sub>|</th><th>|m<sub>2</sub> - m<sub>2,N</sub>|</th><th>|m<sub>3</sub> - m<sub>3,N</sub>|</th></tr></thead><tbody><tr><td>8</td><td>3.761e-02</td><td>3.780e-02</td><td>1.225e+00</td><td>3.968e+01</td><td>1.295e+03</td></tr><tr><td>16</td><td>9.552e-05</td><td>1.227e-04</td><td>7.954e-03</td><td>5.157e-01</td><td>3.350e+01</td></tr><tr><td>32</td><td>2.337e-11</td><td>3.174e-11</td><td>4.092e-09</td><td>5.276e-07</td><td>6.807e-05</td></tr></tbody></table>

Constrained

<table><caption>TABLE 7</caption><thead><tr><th>N</th><th>||f - f<sub>N</sub><sup>c</sup>||<sub>2</sub></th><th>|m<sub>0</sub> - m<sub>0,N</sub><sup>c</sup>|</th><th>|m<sub>1</sub> - m<sub>1,N</sub><sup>c</sup>|</th><th>|m<sub>2</sub> - m<sub>2,N</sub><sup>c</sup>|</th><th>|m<sub>3</sub> - m<sub>3,N</sub><sup>c</sup>|</th></tr></thead><tbody><tr><td>8</td><td>4.834e-02</td><td>8.882e-16</td><td>3.553e-15</td><td>3.126e-13</td><td>3.638e-12</td></tr><tr><td>16</td><td>1.174e-04</td><td>8.882e-16</td><td>0</td><td>5.684e-14</td><td>9.095e-13</td></tr><tr><td>32</td><td>2.622e-11</td><td>0</td><td>3.553e-15</td><td>2.842e-14</td><td>2.274e-13</td></tr></tbody></table>

Error Analysis: Standard vs. Conservative using Laguerre Function approximation for (3.5).

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Error Analysis: Solution and first Moment using Laguerre functions for (3.6).</figcaption></figure>

## 3.2. Test 2: Kinetic Fokker-Planck equation.
Let us now consider the following kinetic Fokker-Planck equation on $\mathbb{R}$:

$$ (3.7) \qquad \begin{cases} \frac{\partial}{\partial t} f = \partial_v ((v-\mu)f + T\partial_v f) := L_{FP}(f), \\ f(0, v) = f^0(v). \end{cases} $$

In this case, the unknown $f(t, v)$ describes the amount of particles moving with velocity $v$ at time $t$. In this setting, the particle will be forced towards a mean velocity $\mu$ and an average temperature $T$. It is well known that this equation admits a steady state $f_\infty$ given by a Maxwellian distribution:

$$ (3.8) \qquad f_{\infty} = \frac{\rho}{\sqrt{2\pi T}} \exp\left(\frac{-(v-\mu)^2}{2T}\right), $$

where $\rho$ is the total mass of particles. In the following, let us consider the initial data

$$ (3.9) \qquad f^0(v) = \frac{1}{\sqrt{\pi}} \left( e^{-(v+2)^2} + e^{-(v-2)^2} \right) $$

This distribution admits a mass $\rho^0$, mean velocity $\mu^0$ and temperature $T^0$ given by

$$ (3.10) \qquad \rho^0 = \int_{\mathbb{R}} f^0(v) dv, \quad \mu^0 = \frac{1}{\rho^0} \int_{\mathbb{R}} v f^0(v) dv, \quad T^0 = \frac{1}{\rho^0} \int_{\mathbb{R}} |v - \mu^0|^2 f^0(v) dv. $$
In addition, we set the parameters of $L_{FP}$ in (3.7) as

$$ (3.11) \qquad \mu = \mu^0 \text{ and } T = T^0. $$

With these choice of parameters, one can in particular expect preservation of mass, mean velocity and temperature. We now consider the symmetrically weighted Hermite functions: $\psi_k(v) = H_k(v)e^{-v^2/2}$. By definition of the Hermite polynomials, we have the relations:

$$ (3.12) \qquad \begin{aligned} H_k &= 2vH_{k-1} - 2(k-1)H_{k-2} \\ H'_k &= 2kH_{k-1}. \end{aligned} $$

Therefore, the Hermite functions $\psi_k$ satisfy:

$$ (3.13) \qquad \begin{aligned} \psi_k &= 2v\psi_{k-1} - 2(k-1)\psi_{k-2}, \\ \psi'_k &= k\psi_{k-1} - \frac{1}{2}\psi_{k+1}, \\ \psi''_k &= k(k-1)\psi_{k-2} - \left(k + \frac{1}{2}\right)\psi_k + \frac{1}{4}\psi_{k+2}. \end{aligned} $$

Using these relations, one can explicitly compute

$$ (3.14) \qquad \begin{aligned} L_{FP}(\psi_k) &= \psi_k + (v - \mu)\psi'_k + T\psi''_k \\ &= k(k-1)(1+T)\psi_{k-2} - \frac{\mu}{2}\psi_{k-1} + \left(-kT - \frac{T-1}{2}\right)\psi_k \\ &\quad + \frac{\mu}{2}\psi_{k+1} + \frac{1}{4}(T-1)\psi_{k+2}. \end{aligned} $$

We can now proceed to use a Galerkin method to solve (3.7). We denote by $\mathcal{P}_N$ and $\mathcal{P}_N^c$ the standard and constrained Hermite function approximation. We want to find $f_N \in S_N$ solution to

$$ (3.15) \qquad \left\{ \begin{aligned} \frac{\partial}{\partial t} f_N &= \mathcal{P}_N (L_{FP}(f_N)), \\ f_N(0, v) &= \mathcal{P}_N(f^0)(v). \end{aligned} \right. $$

The moment constrained problem writes

$$ (3.16) \qquad \left\{ \begin{aligned} \frac{\partial}{\partial t} f_N^c &= L_{FP,N}^c(f_N^c), \\ f_N^c(0, v) &= \mathcal{P}_N^c(f^0)(v), \end{aligned} \right. $$

where $L_{FP,N}^c$ is solution to

$$ (3.17) \qquad L_{FP,N}^c(f) = \underset{g_N \in S_N, \langle g_N, v^q \rangle = 0, q = 0, 1, 2, 3}{\operatorname{argmin}} \left\{ \|g_N - L_{FP}(f)\|_{L^2}^2 : \right. $$

Let us emphasize that (3.17) corresponds to (2.19) where we ensure the first three moments of $L_{FP}$ to be zero. For the time stepping method, we utilize a classical fourth order Runge-Kutta (RK4) method with fixed time step $\Delta t = 10^{-4}$. In Figure 8 one can observe a very good agreement between the two overlapping solutions as well as a trend towards the steady state $f_\infty$. This observation suggests that the constrained approximation may be able to finely capture long-time properties of the solution.
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Fokker-Planck model: Snapshots of the distribution at t = 0, 0.1 and 5, N = 32.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9. Fokker-Planck model, non-zero mean velocity: Snapshots of the distribution at t = 0, 0.1 and 20, N = 32.</figcaption></figure>

We also present in Figure 9 a test case for the following initial data that has a non-zero mean velocity:

$$f^0(v) = \frac{1}{\sqrt{\pi}} (e^{-(v+3)^2} + e^{-v^2}).$$

This test becomes relevant for spatially inhomogeneous kinetic problem. We observe that our constrained method is still able to nicely capture the solution, even in this more challenging case.

**Spectral accuracy and conservations.** Let us now investigate the accuracy of the method. In Figure 10 we set a final time $T_f = 0.1$ and compare the solution to a reference one computed with $N = 32$ modes. We do the same comparison for the first three moments. A first observation is that the standard method does exhibit spectral accuracy both on the solution and on the first three moments. This therefore falls within the scope of Theorem (2.5) and the constrained method is indeed spectrally accurate. Note that in this particular case, the mean velocity is 0, hence the machine accuracy in the bottom left panel of Figure 10.

**Long-time behaviour.** Let us conclude this test by investigating the long time behaviour of the solution observed in Figure 8. In Figure 11 The top left panel shows the evolution of the norms $\|f_N(t^n) - f_\infty\|_2$ and $\|f_N^c(t^n) - f_\infty\|_2$ while the remaining panels show the variation of the moments $|\rho^n - \rho^0|$, $|\mu^n - \mu^0|$ and $|T^n - T^0|$. We observe that the constrained method indeed preserves the first three moments up to machine accuracy. As a consequence, the long-time behaviour should be better approximated. The standard approximation does seem to converge towards the steady state but ultimately, accumulation of errors on the moments overcomes the dynamic
<figure><img src="image_11.png" /><figcaption>FIGURE 10. Fokker-Planck model: Approximation error on the solution (Top Left) and absolute error on the mass (Top Right), mean velocity (Bottom Left) and temperature (Bottom Right) at $T_f = 0.1$.</figcaption></figure>

of the solution, pushing it away from the equilibrium. On the contrary, our proposed method does not exhibit this re-bounce and converges to a plateau that decreases as one increases the number of modes. One can mention that equilibrium preserving techniques [8, 29] have been developed to overcome this saturation but these are beyond the scope of this article.

**3.3. Test 3: A model of Opinion formation.** For this second test, we turn our attention to the following opinion model [11] where $v \in [-1, 1]$:

$$ (3.18) \qquad \begin{cases} \frac{\partial}{\partial t}g = \frac{\lambda}{2}\partial_{vv}((1-v^2)g) + \partial_v((v-m)g) := L_{Op}g, \\ g(0,v) = g^0(v). \end{cases} $$

In this case, the unknown $g(t, v)$ describes the distribution of individuals at time $t$ with an opinion trait $v$. The modelling parameters are chosen so that $|m| < 1$ and $\lambda < 1 + |m|$ to prevent blow-ups as $v \to \pm 1$. In addition, we supplement this equation with a no flux boundary condition to ensure the conservation of mass:

$$ (3.19) \qquad \frac{\lambda}{2} \partial_v ((1 - v^2)g) + (v - m)g \Big|_{v=\pm 1} = 0. $$

In particular, one can show that this condition in fact boils down to the homogeneous Dirichlet boundary condition $g(t, \pm 1) = 0$. This equation also admits a steady state given by

$$ (3.20) \qquad g_{\infty}(v) = c_{m,\lambda}(1+v)^{\frac{1+m}{\lambda}-1}(1-v)^{\frac{1-m}{\lambda}-1} $$
<figure><img src="image_12.png" /><figcaption>FIGURE 11. Fokker Planck model: Convergence towards the steady state (Top Left) and absolute error on the mass (Top Right), mean velocity (Bottom Left) and temperature (Bottom Right).</figcaption></figure>

where the constant $c_{m,\lambda}$ is for normalization. In order to approximate solutions to (3.18) we want to consider general Jacobi polynomials $\binom{p_k^{\alpha,\beta}}{k=0,...,N}$ associated to the weight $\omega(v) = (1-v)^{\alpha}(1+v)^{\beta}$. However, this polynomial basis does not satisfy the homogeneous Dirichlet boundary conditions. To solve this issue, we construct a new polynomial basis obtained from the original Jacobi polynomials (see e.g. [35]). For $\alpha, \beta > -1$, let

$$ (3.21) \qquad \zeta_k = p_k^{\alpha,\beta} + a_k p_{k+1}^{\alpha,\beta} + b_k p_{k+2}^{\alpha,\beta}, \quad k = 0, \dots, N-2 $$

where the constants $a_k$ and $b_k$ are solution to the linear system:

$$ (3.22) \qquad \begin{cases} \zeta_k(-1) = 0 \\ \zeta_k(1) = 0. \end{cases} $$

This new polynomial family now satisfies the boundary conditions but it is, by construction, no longer orthogonal. Since orthogonality is a key ingredient of the conservative spectral method, we need to apply a Gram-Schmidt algorithm to obtain an orthonormal basis. Note that in practice the process of orthonormalization may induce an accumulation of machine errors that can lead to a significant loss in orthogonality for large number of modes. We can now expand $g$ in the basis of polynomials $(\zeta_k)_{k=0,...,N-2}$:

$$ (3.23) \qquad g(t, v) = \sum_{k=0}^{N-2} \hat{g}_k(t) \zeta_k(v), $$
and then proceed as before to solve (3.18) using a Galerkin method. We define the standard and constrained problems as

$$ (3.24) \qquad \begin{cases} \frac{\partial}{\partial t} g_N = \mathcal{P}_N (L_{Op}(g_N)), \\ g_N(0, v) = \mathcal{P}_N(g^0)(v), \end{cases} $$

and

$$ (3.25) \qquad \begin{cases} \frac{\partial}{\partial t} g_N^c = L_{Op,N}^c(g_N^c), \\ g_N(0, v) = \mathcal{P}_N^c(g^0)(v). \end{cases} $$

Since (3.18) preserves only the first moment, we define $L_{Op,N}^c$ as the solution to

$$ (3.26) \qquad L_{Op,N}^c(g) = \operatorname*{argmin}\left\{\|g_N - L_{Op}(g)\|_{L^2_\omega}^2 : g_N \in S_N, \langle g_N, 1 \rangle = 0\right\}. $$

The time steeping is again achieved using an RK4 time integrator with time step $\Delta t = 10^{-4}$. Before presenting numerical results, let us mention that the space spanned by the newly constructed basis does not contain the monomials $v^q$, $q = 0, 1, 2, ...$ anymore. As a consequence, in regards of (2.17), even if one considers Legendre polynomials ($\omega = 1$) to construct the new basis, one does not expect exact conservation of moments. However, as will be shown through numerical experiments, one still recovers spectral accuracy on the moments for the standard method which ensures spectral accuracy on the constrained solution. In regards of this discussion, we now present only the choice of Legendre polynomials for the construction of the new basis ($\zeta_k)_{k=0,...,N-2}$. Note however that classical families such as Chebyshev first and second kind provide similar results. As an initial data, we take

$$ (3.27) \qquad g^0(v) = c_0 \left( (1+v)^{12} (1-v)^6 + (1+v)^{13} (1-v)^{25} \right), $$

where $c_0$ is for normalizing the first moment. Then we set $m = 0$ and $\lambda = 0.1$ in the definition of $L_{Op}$ in (3.18), therefore expecting an equilibrium of the form

$$ g_\infty(v) = c_\infty(1 - v^2)^9. $$

Figure 12 shows the evolution of the approximate distributions $g_N$ and $g_N^c$. As in the previous case, we observe a very good agreement between the standard and constrained method with the two curves overlapping. Moreover, the solutions also appear to converge towards the steady state.

*Spectral accuracy and conservations.* Let us now investigate the accuracy of the method and its mass conservation properties. Indeed, at the continuous level, equation (3.18) along with no flux boundary conditions preserves only the 0th order moment. In the next study, we set a final time $T_f = 0.1$ and compare the approximation error between the standard and constrained method. We observe in Figure 13 that we again fall into the framework of the convergence theorem 2.5 and the constrained method behaves as predicted.

*Long-time behaviour.* Let us conclude this test by looking at the long time behaviour of the approximated solution to (3.18). In order to quantify the observation made on Figure 12, we show in Figure 14 the evolution of the norms $\|g_N - g_\infty\|_2$ and $\|g_N^c - g_\infty\|_2$ as well as the evolution of the mass variation $|m_0^0 - m_0(t^n)|$. Similarly as before, we observe that the conservation of moments allows to stabilize the long time dynamic of the solution.
<figure><img src="image_13.png" /><figcaption>FIGURE 12. Opinion model: Snapshots of the distribution at times t = 0, 0.4 and 5, N = 24.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 13. Opinion model: Approximation error on the solution (Left) and on the 0th order moment (Right) at T<sub>f</sub> = 0.1.</figcaption></figure>

**3.4. Test 4: A Call center service time model.** This final test deals with a
call center service time model introduced in [17]:

$$
(3.28) \qquad \begin{cases} \frac{\partial}{\partial t}h = \frac{\lambda}{2}\partial_{vv}(v^2h) + \frac{\gamma}{2}\partial_v\left(v\ln\left(\frac{v}{v_L}\right)h\right) := L_{CC}h, \\ h(0, v) = h^0(v). \end{cases}
$$

In this context, the distribution density $h(t, v)$ describes the service time of agents in a call center. The constants $0 < \gamma < 1$, $\lambda > 0$ and $v_L > 0$ are modelling parameters. In particular, the latter corresponds to an ideal time for agents to complete their task. To supplement this equation, we consider a no-flux boundary condition

$$
(3.29) \qquad \frac{\lambda}{2} \partial_v (v^2 h) + \frac{\gamma}{2} \ln \left( \frac{v}{v_L} \right) h \Big|_{v=0} = 0.
$$

It is clear that this condition is automatically satisfied by (3.28). It can also be shown
that (3.28) admits an equilibrium distribution of the form:

$$
(3.30) \qquad h_{\infty}(v) = \frac{1}{\sqrt{2\pi\sigma v}} \exp\left(-\frac{(\ln(v) - \mu)^2}{2\sigma}\right)
$$

where $\sigma = \frac{\lambda}{\gamma}$ and $\mu = \ln(v_L) - \sigma$. As in Test 1, in order to deal with integrability at infinity, we want to consider the symmetrically weighted Laguerre functions: $\xi_k = L_k \omega^{\frac{1}{2}}$. In addition, to deal with the poor approximation of the steady state observed
<figure><img src="image_15.png" /><figcaption>FIGURE 14. Opinion model: Convergence towards the steady state (Left) and mass variation (Right).</figcaption></figure>

in Figure 7 we consider a micro-macro approach. Since the equilibrium is known, we can decompose the unknown as

$$ (3.31) \qquad h = h_{\infty} + \tilde{h}, $$

where $\tilde{h}$ is the perturbation that can take negative values. Since (3.28) is linear in $h$, and by definition of the steady state, the perturbation also solves (3.28). The only difference is the initial data.

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial t} \tilde{h} = \frac{\lambda}{2} \partial_{vv} (v^2 \tilde{h}) + \frac{\gamma}{2} \partial_v \left( v \ln \left( \frac{v}{v_L} \right) \tilde{h} \right), \\[1.5ex] \displaystyle \tilde{h}(0, v) = h^0(v) - h_{\infty}(v). \end{array} \right. $$

From this procedure, we expect to observe that $\tilde{h}$ converges towards 0 for long time. Consequently, the long time behaviour of the solution should be better approximated since 0 now belongs to the approximation space. However, one cannot say anything about the short time behaviour as the initial data could still be poorly approximated.

As previously, we can then define the standard and constrained problems on the perturbation as

$$ (3.32) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial t} \tilde{h}_N = \mathcal{P}_N (L_{CC}(\tilde{h}_N)), \\[1.5ex] \displaystyle \tilde{h}_N(0, v) = \mathcal{P}_N(\tilde{h}^0 - h_{\infty})(v), \end{array} \right. $$

and

$$ (3.33) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial t} \tilde{h}_N^c = L_{CC,N}^c(\tilde{h}_N^c), \\[1.5ex] \displaystyle \tilde{h}_N(0, v) = \mathcal{P}_N^c(\tilde{h}^0 - h_{\infty})(v). \end{array} \right. $$

The original distribution is then reconstructed using (3.31) for visualisations. Since (3.28) preserves only the first moment, we define $L_{CC,N}^c$ as the solution to

$$ (3.34) \qquad L_{CC,N}^c(h) = \operatorname*{argmin}\left\{\|g_N - L_{CC}(h)\|_{L^2}^2 : g_N \in S_N, \langle g_N, 1 \rangle = 0\right\}. $$

The time stepping is yet again achieved through a RK4 method with time step $\Delta t = 10^{-4}$. Let us now fix the parameters $\lambda = 0.5$, $\gamma = 0.9$ and $v_L = 40$ in (3.28). As an
<figure><img src="image_16.png" /><figcaption>FIGURE 15. Call Center model: Snapshots of the distributions at times t = 0, 2 and 5, N = 32.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 16. Call Center model: Approximation error on the solution (Left) and on the 0th order moment (Right) at T<sub>f</sub> = 0.1.</figcaption></figure>

initial data, we consider

$$ (3.35) \qquad h^0(v) = (x^3 - 2x + \sin(x))e^{-x}. $$

This function is actually well approximated by Laguerre functions. We show on Figure 15 the evolution of the standard and constrained distributions. A first observation is that there is again a good agreement between the two methods. In addition, as expected from the micro-macro approach, the steady state is well approximated.

**Accuracy and conservation.** To assess the spectral accuracy of the method, we consider a final time $T_f = 0.1$ along with a time step $\Delta t = 10^{-4}$ for the time integration. We present on Figure 16 the study of the convergence of the method. As expected, in short times, the distribution is not very well approximated in our chosen basis. As a consequence, we are only able to recover a convergence that is quite slow. However, the up-side is that the mass is still exactly preserved by the constrained approximation. Now, contrary to previous tests, it is important to note that the standard method is not spectrally accurate on the mass.

**Long-time behaviour.** To conclude this section, and as in previous tests, we investigate the longtime behaviour of the solution through the norms $\|h_n - h_\infty\|_2$ and $\|h_n^c - h_\infty\|_2$. We observe in Figure 17 that the micro-macro approach allows to approximate really well the steady state with both methods. In addition, the mass variation induced by the standard approach decays exponentially fast towards 0 therefore ensuring that it also converges towards $h_\infty$ albeit more slowly than the conservative approach.
<figure><img src="image_18.png" /><figcaption>FIGURE 17. Call Center model: Convergence towards the steady state (Left), and mass variation (Right).</figcaption></figure>

**4. Conclusions.** We extended and analyzed a class of Galerkin spectral methods for PDEs based on general families of orthogonal functions that can preserve the moments of the solution. Leading examples of such problems arise in kinetic and mean-field theory where preservation of the moments of the distribution function is of paramount importance to describe correctly the long-time behavior. The methods were designed using an $L_2$ best constrained approximation formalism in the space of orthogonal polynomials. Thanks to this setting, the resulting methods allow the spectral accuracy property of the solution to be maintained for the conservative polynomial series. The new approximation was then used to derive moment preserving methods for several Fokker-Planck equations in distinct physical domains and by adopting different standard Galerkin projections, including the use of micro-macro decompositions.

Given its generality, the method admits numerous extensions. Among the most interesting are certainly the construction of spectrally accurate and conservative methods for the Vlasov equation. Another interesting direction is the application to kinetic models in the socio-economic domain where equilibrium states are often unknown and thanks to the present approach can be computed with spectral accuracy.

*Acknowledgment.* This work has been written within the activities of GNCS and GNFM groups of INdAM (Italian National Institute of High Mathematics). LP has been supported by the Royal Society under the Wolfson Fellowship "Uncertainty quantification, data-driven simulations and learning of multiscale complex systems governed by PDEs". The partial support by ICSC – Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union – NextGenerationEU and by MUR PRIN 2022 Project No. 2022KKJP4X "Advanced numerical methods for time dependent parametric partial differential equations with applications" is also acknowledged. T.L. would like to thank Thomas Rey for fruitful discussions on the implementation of the method. T.L. acknowledges support from Labex CEMPI (ANR-11-LABX-0007-01) and the MSCA DN-2022 program DATAHYKING.

## REFERENCES

[1] ALONSO, R. J., GAMBA, I. M., AND THARKABHUSHANAM, S. H. Convergence and error estimates for the Lagrangian-based conservative spectral method for Boltzmann equations. *SIAM J. Numer. Anal.* 56, 6 (2018), 3534–3579.

[2] BAILO, R., CARRILLO, J. A., AND HU, J. Fully discrete positivity-preserving and energy-
dissipating schemes for aggregation-diffusion equations with a gradient-flow structure. Commun. Math. Sci. 18, 5 (2020), 1259–1303.

[3] BESSEMOULIN-CHATARD, M., AND FILBET, F. On the stability of conservative discontinuous Galerkin/Hermite spectral methods for the Vlasov-Poisson system. ArXiv abs/2106.07468 (2021).

[4] BUET, C., AND DELLACHERIE, S. On the Chang and Cooper scheme applied to a linear Fokker-Planck equation. Commun. Math. Sci. 8, 4 (2010), 1079–1090.

[5] CANUTO, C., HUSSAINI, M., QUARTERONI, A., AND ZANG, T. Spectral Methods in Fluid Dynamics. Springer Series in Computational Physics. Springer-Verlag, New York, 1988.

[6] CHANG, J., AND COOPER, G. A practical difference scheme for Fokker-Planck equations. Journal of Computational Physics 6, 1 (1970), 1–16.

[7] DIMARCO, G., LI, Q., PARESCHI, L., AND YAN, B. Numerical methods for plasma physics in collisional regimes. J. Plasma Phys. 81, 1 (2015), 305810106.

[8] FILBET, F., PARESCHI, L., AND REY, T. On steady-state preserving spectral methods for homogeneous Boltzmann equations. Comptes Rendus Mathematique 353 (2014), 309–314.

[9] FOK, J. C. M., YU GUO, B., AND TANG, T. Combined Hermite spectral-finite difference method for the Fokker-Planck equation. Math. Comput. 71 (2002), 1497–1528.

[10] FUNARO, D. Polynomial approximations of differential equations. Lecture Notes in Physics. Springer-Verlag, New York, 1992.

[11] FURIOLI, G., PULVIRENTI, A., TERRANEO, E., AND TOSCANI, G. Fokker-planck equations in the modeling of socio-economic phenomena. Mathematical Models and Methods in Applied Sciences 27 (2017), 115–158.

[12] GAMBA, I. M., HAACK, J. R., HAUCK, C. D., AND HU, J. A fast spectral method for the Boltzmann collision operator with general collision kernels. SIAM J. Sci. Comput. 39, 4 (2017), B658–B674.

[13] GAMBA, I. M., AND THARKABHUSHANAM, S. H. Shock and Boundary Structure formation by Spectral-Lagrangian methods for the Inhomogeneous Boltzmann Transport Equation. J. Comput. Math. 28, 4 (2010), 430–460.

[14] GAUTSCHI, W. Orthogonal Polynomials: Computation and Approximation. Oxford University Press, Oxford, 2004.

[15] GOSSE, L. Computing Qualitatively Correct Approximations of Balance Laws. Exponential-Fit, Well-Balanced and Asymptotic-Preserving. SEMA SIMAI Springer Series. Springer, 2013.

[16] GOTTLIEB, D., AND ORSZAG, S. A. Numerical Analysis of Spectral Methods: Theory and Applications. SIAM, Philadelphia, 1977.

[17] GUALANDI, S., AND TOSCANI, G. Call center service times are lognormal: A Fokker-Planck description. Mathematical Models and Methods in Applied Sciences (2018).

[18] HAIRER, E., WANNER, G., AND LUBICH, C. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, vol. 31 of Springer Series in Computational Mathematics. Springer, 2006.

[19] HOLLOWAY, J. P. Spectral velocity discretizations for the Vlasov-Maxwell equations. Transport Theory and Statistical Physics 25 (1996), 1–32.

[20] HU, J., AND MA, Z. A fast spectral method for the inelastic Boltzmann collision operator and application to heated granular gases. J. Comput. Phys. 385 (2019), 119–134.

[21] JIN, S. Asymptotic-preserving schemes for multiscale physical problems. Acta Numerica 31 (2022), 415–489.

[22] KOEKOEK, R., LESKY, P. A., AND SWARTTOUW, R. F. Hypergeometric Orthogonal Polynomials and their q-Analogues. Springer, Berlin, 2010.

[23] KRAUS, M., AND HIRVIJOKI, E. Metriplectic integrators for the Landau collision operator. Physics of Plasmas 24, 10 (2017), 102311.

[24] LARSEN, E. W., LEVERMORE, C. D., POMRANING, G. C., AND SANDERSON, J. G. Discretization methods for one-dimensional Fokker-Planck operators. J. Comput. Phys. 61, 3 (1985), 359–390.

[25] MANZINI, G., FUNARO, D., AND DELZANNO, G. L. Convergence of spectral discretizations of the vlasov-poisson system. SIAM J. Numer. Anal. 55 (2016), 2312–2335.

[26] MOHAMMADI, M., AND BORZÍ, A. Analysis of the Chang-Cooper discretization scheme for a class of Fokker-Planck equations. J. Numer. Math. 23, 3 (2015), 271–288.

[27] MOUHOT, C., AND PARESCHI, L. Fast algorithms for computing the Boltzmann collision operator. Math. Comp. 75, 256 (2006), 1833–1852 (electronic).

[28] MOUHOT, C., PARESCHI, L., AND REY, T. Convective Decomposition and Fast Summation Methods for Discrete-Velocity Approximations of the Boltzmann Equation. ESAIM Math. Model. Numer. Anal. 47, 5 (2013), 1515–1531.

[29] PARESCHI, L., AND REY, T. On the stability of equilibrium preserving spectral methods for the
homogeneous Boltzmann equation. *Applied Mathematics Letters abs/2011.05811* (2020).

[30] PARESCHI, L., AND REY, T. Moment preserving Fourier-Galerkin spectral methods and application to the Boltzmann equation. *SIAM J. Num. Anal.*, to appear (2022).

[31] PARESCHI, L., AND RUSSO, G. Numerical Solution of the Boltzmann Equation I : Spectrally Accurate Approximation of the Collision Operator. *SIAM J. Numer. Anal.* **37**, 4 (2000), 1217–1245.

[32] PARESCHI, L., AND RUSSO, G. On the stability of spectral methods for the homogeneous Boltzmann equation. *Transport Theory Statist. Phys.* **29**, 3-5 (2000), 431–447.

[33] PARESCHI, L., RUSSO, G., AND TOSCANI, G. Fast spectral methods for the Fokker-Planck-Landau collision operator. *J. Comput. Phys.* **165**, 1 (2000), 216–236.

[34] PARESCHI, L., AND ZANELLA, M. Structure preserving schemes for nonlinear Fokker-Planck equations and applications. *J. Sci. Comput.* **74**, 3 (2018), 1575–1600.

[35] SHEN, J., TANG, T., AND WANG, L.-L. *Spectral Methods: Algorithms, Analysis and Applications*, 1st ed. Springer Publishing Company, Incorporated, 2011.

[36] XIU, D. *Numerical Methods for Stochastic Computations: A Spectral Method Approach*. Princeton University Press, Princeton, 2010.