Accurate spectral numerical schemes for kinetic equations with energy diffusion

Abstract

We examine the merits of using a family of polynomials that are orthogonal with respect to a non-classical weight function to discretize the speed variable in continuum kinetic calculations. We consider a model one-dimensional partial differential equation describing energy diffusion in velocity space due to Fokker-Planck collisions. This relatively simple case allows us to compare the results of the projected dynamics with an expensive but highly accurate spectral transform approach. It also allows us to integrate in time exactly, and to focus entirely on the effectiveness of the discretization of the speed variable. We show that for a fixed number of modes or grid points, the non-classical polynomials can be many orders of magnitude more accurate than classical Hermite polynomials or finite-difference solvers for kinetic equations in plasma physics. We provide a detailed analysis of the difference in behavior and accuracy of the two families of polynomials. For the non-classical polynomials, if the initial condition is not smooth at the origin when interpreted as a three-dimensional radial function, the exact solution leaves the polynomial subspace for a time, but returns (up to roundoff accuracy) to the same point evolved to by the projected dynamics in that time. By contrast, using classical polynomials, the exact solution differs significantly from the projected dynamics solution when it returns to the subspace. We also explore the connection between eigenfunctions of the projected evolution operator and (non-normalizable) eigenfunctions of the full evolution operator, as well as the effect of truncating the computational domain.

A Spectral Transform Method for Singular Sturm--Liouville Problems with Applications to Energy Diffusion in Plasma Physics

“…Thus, with limited computational resources, the solution of u t = −Lu often cannot be resolved to the desired level of accuracy until t surpasses a critical value, t * , where the decay rate ofû(λ, t * ) becomes fast enough. Remarkably, the same is true of the projected dynamics in some spaces of orthogonal polynomials [49]. For singular initial conditions, the projected dynamics is a poor approximation of the true solution initially, regardless of which space of polynomials is used to represent the solution.…”

“…In subsequent work [49], jointly with Landreman, we will study the projected dynamics of this equation in finite-dimensional spaces of orthogonal polynomials. Roughly speaking, we show in this paper how to efficiently evaluate the exact solution by discretizing a continuous transform, while in [49] we discretize the PDE before evolving the solution.…”

“…In subsequent work [49], jointly with Landreman, we will study the projected dynamics of this equation in finite-dimensional spaces of orthogonal polynomials. Roughly speaking, we show in this paper how to efficiently evaluate the exact solution by discretizing a continuous transform, while in [49] we discretize the PDE before evolving the solution. The latter approach is faster and better suited to large scale computations of the full Fokker-Planck equation, while the current approach clarifies the role of the continuous spectrum in the dynamics and provides an independent means of validating the orthogonal polynomial approach.…”

“…If each term of the right-hand side of (3.28) is perturbed by O(δ), the accumulated effect remains O(nδ). In practice, the errors really are this small, as shown in the follow-up paper [48] by comparing doubleprecision results to "exact" solutions computed in quadruple-precision, and indirectly in [49] by comparing this method of solving (2.4) to a projected dynamics approach using orthogonal polynomials.…”

“…Such agreement provides strong evidence that monitoring Chebyshev and Fourier mode decay rates provides accurate a posteriori error estimates. The current algorithm may be viewed as a method of approximating exact integral formulas for the solution, while that in [49] may be regarded as a (nearly) exact evolution of a finite-dimensional approximation of the PDE. To gauge the performance of the new algorithm, in follow-up work [48], we compare our method of computing the spectral density function (steps 2 and 3) to the popular software package SLEDGE [36,19,18] and to the algorithm of Fulton, Pearson, and Pruess (FPP) [15,16].…”