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

Abstract

Abstract. We develop a spectrally accurate numerical method to compute solutions of a model PDE used in plasma physics to describe diffusion in velocity space due to Fokker-Planck collisions. The solution is represented as a discrete and continuous superposition of normalizable and nonnormalizable eigenfunctions via the spectral transform associated with a singular Sturm-Liouville operator. We present a new algorithm for computing the spectral density function of the operator that uses Chebyshev polynomials to extrapolate the value of the Titchmarsh-Weyl m-function from the complex upper half-plane to the real axis. The eigenfunctions and density function are rescaled, and a new formula for the limiting value of the m-function is derived to avoid amplification of roundoff errors when the solution is reconstructed. The complexity of the algorithm is also analyzed, showing that the cost of computing the spectral density function at a point grows less rapidly than any fractional inverse power of the desired accuracy. A WKB analysis is used to prove that the spectral density function is real analytic. Using this new algorithm, we highlight key properties of the PDE and its solution that have strong implications on the optimal choice of discretization method in large-scale plasma physics computations. 1. Introduction. Partial differential equations involving singular Sturm-Liouville operators with continuous spectra arise frequently in computational physics. Common approaches to solving them include domain truncation, which often regularizes the operator and makes the spectrum discrete, or projection onto finite dimensional orthogonal polynomial or finite element subspaces, which also leads to discrete spectra. Here we develop an alternative approach in which the continuous spectrum is treated analytically via a spectral transform, and the numerical challenge is in accurately representing and evaluating the integrals giving the exact solution.While the methods developed in this paper to diagonalize singular Sturm-Liouville operators are quite general, we will describe them in the context of velocity-space diffusion in one dimension,

Accurate spectral numerical schemes for kinetic equations with energy diffusion

“…The point spectrum of L consists of λ = 0 with eigenfunction u ≡ 1, and the continuous spectrum is (0, ∞). A spectral transform algorithm is developed in [14] that diagonalizes the evolution operator and expresses the solution as a discrete and continuous superposition of normalizable and non-normalizable eigenfunctions of L. We will use this computationally expensive method to assess the accuracy of the projected dynamics below. In this work, we approximate solutions of the PDE (3.12) by projecting onto finite dimensional subspaces of orthogonal polynomials.…”

“…For comparison, we also plot the spectral transform of the solution in (C), which was computed at 512 equally spaced points between σ = −4 and σ = 8 using the algorithm described in [14]. Here σ = ln λ is the spectral parameter used in [14] to represent the solution as a discrete and continuous superposition of eigenfunctions. In more detail, the solution of the PDE in the infinite dimensional space H may be written…”

“…where v(x, σ) = u 1 (x, e σ )/Y (e σ ), u 1 (x, λ) is a solution of Lu = λu with appropriate boundary conditions at x = 0, Y (λ) is a scale factor defined in [14],…”

“…The purpose of this paper is to explore the suitability of the non-classical polynomials for initial-value calculations of turbulent plasma transport in the presence of collisions [4,13]. To do so, we consider a model one-dimensional problem describing energy diffusion due to Fokker-Planck collisions [14]: Here, Ψ differs from the usual Chandrasekhar function by an additional factor of 1/x. Our choice to focus on this particular model problem is motivated by the following characteristic features.…”

“…In both cases we find that the non-standard polynomials are effective at representing the solution, and that they are much more accurate than classical polynomials defined by orthogonality conditions on the interval (−∞, ∞). The improved behavior is explained by comparing the mode amplitudes in the eigenbasis of the projected evolution operator to the spectral transform of the solution, computed as described in [14]. This transform method is used in §4.2- §4.4 as an independent means of validating the accuracy of the Galerkin approach.…”