Yingda Cheng

Wei Guo, Yingda Cheng

In this paper, we develop a sparse grid discontinuous Galerkin (DG) scheme for transport equations and applied it to kinetic simulations. The method uses the weak formulations of traditional Runge-Kutta DG (RKDG) schemes for hyperbolic problems and is proven to be $L^2$ stable and convergent. A major advantage of the scheme lies in its low computational and storage cost due to the employed sparse finite element approximation space. This attractive feature is explored in simulating Vlasov and Boltzmann transport equations. Good performance in accuracy and conservation is verified by numerical tests in up to four dimensions.

Yingda Cheng, Irene M. Gamba, Philip J. Morrison

In this paper we consider Runge-Kutta discontinuous Galerkin (RKDG) schemes for Vlasov-Poisson systems that model collisionless plasmas. One-dimensional systems are emphasized. The RKDG method, originally devised to solve conservation laws, is seen to have excellent conservation properties, be readily designed for arbitrary order of accuracy, and capable of being used with a positivity-preserving limiter that guarantees positivity of the distribution functions. The RKDG solver for the Vlasov equation is the main focus, while the electric field is obtained through the classical representation by Green's function for the Poisson equation. A rigorous study of recurrence of the DG methods is presented by Fourier analysis, and the impact of different polynomial spaces and the positivity-preserving limiters on the quality of the solutions is ascertained. Several benchmark test problems, such as Landau damping, two-stream instability and the KEEN (Kinetic Electrostatic Electron Nonlinear) wave, are given.

Olivier Bokanowski, Yingda Cheng, Chi-Wang Shu

We study semi-Lagrangian discontinuous Galerkin (SLDG) and Runge-Kutta discontinuous Galerkin (RKDG) schemes for some front propagation problems in the presence of an obstacle term, modeled by a nonlinear Hamilton-Jacobi equation of the form $\min(u_t + c u_x, u - g(x))=0$, in one space dimension. New convergence results and error bounds are obtained for Lipschitz regular data. These "low regularity" assumptions are the natural ones for the solutions of the studied equations.

Zhanjing Tao, Wei Guo, Yingda Cheng

In this paper, we develop sparse grid discontinuous Galerkin (DG) schemes for the Vlasov-Maxwell (VM) equations. The VM system is a fundamental kinetic model in plasma physics, and its numerical computations are quite demanding, due to its intrinsic high-dimensionality and the need to retain many properties of the physical solutions. To break the curse of dimensionality, we consider the sparse grid DG methods that were recently developed in \cite{guo2016sparse,guo2017adaptive} for transport equations. Such methods are based on multiwavelets on tensorized nested grids and can significantly reduce the numbers of degrees of freedom. We formulate two versions of the schemes: sparse grid DG and adaptive sparse grid DG methods for the VM system. Their key properties and implementation details are discussed. Accuracy and robustness are demonstrated by numerical tests, with emphasis on comparison of the performance of the two methods, as well as with their full grid counterparts.

Yingda Cheng, Irene M. Gamba, Fengyan Li et al.

Discontinuous Galerkin methods are developed for solving the Vlasov-Maxwell system, methods that are designed to be systematically as accurate as one wants with provable conservation of mass and possibly total energy. Such properties in general are hard to achieve within other numerical method frameworks for simulating the Vlasov-Maxwell system. The proposed scheme employs discontinuous Galerkin discretizations for both the Vlasov and the Maxwell equations, resulting in a consistent description of the distribution function and electromagnetic fields. It is proven, up to some boundary effects, that charge is conserved and the total energy can be preserved with suitable choices of the numerical flux for the Maxwell equations and the underlying approximation spaces. Error estimates are established for several flux choices. The scheme is tested on the streaming Weibel instability: the order of accuracy and conservation properties of the proposed method are verified.

Anqi Chen, Fengyan Li, Yingda Cheng

In this paper, we develop an ultra-weak discontinuous Galerkin (DG) method to solve the one-dimensional nonlinear Schrödinger equation. Stability conditions and error estimates are derived for the scheme with a general class of numerical fluxes. The error estimates are based on detailed analysis of the projection operator associated with each individual flux choice. Depending on the parameters, we find out that in some cases, the projection can be defined element-wise, facilitating analysis. In most cases, the projection is global, and its analysis depends on the resulting $2\times2$ block-circulant matrix structures. For a large class of parameter choices, optimal $\textit{a priori}$ $L^2$ error estimates can be obtained. Numerical examples are provided verifying theoretical results.

Yingda Cheng, Andrew J. Christlieb, Wei Guo et al.

In plasma simulations, where the speed of light divided by a characteristic length is at a much higher frequency than other relevant parameters in the underlying system, such as the plasma frequency, implicit methods begin to play an important role in generating efficient solutions in these multi-scale problems. Under conditions of scale separation, one can rescale Maxwell's equations in such a way as to give a magneto static limit known as the Darwin approximation of electromagnetics. In this work, we present a new approach to solve Maxwell's equations based on a Method of Lines Transpose (MOL$^T$) formulation, combined with a fast summation method with computational complexity $O(N\log{N})$, where $N$ is the number of grid points (particles). Under appropriate scaling, we show that the proposed schemes result in asymptotic preserving methods that can recover the Darwin limit of electrodynamics.

Anqi Chen, Yingda Cheng, Yong Liu et al.

We analyze the superconvergence properties of ultra-weak discontinuous Galerkin (UWDG) methods with various choices of flux parameters for one-dimensional linear Schrödinger equation. In our previous work [10], stability and optimal convergence rate are established for a large class of flux parameters. Depending on the flux choices and if the polynomial degree $k$ is even or odd, in this paper, we prove $2k$ or $(2k-1)$-th order superconvergence rate for cell averages and numerical flux of the function, as well as $(2k-1)$ or $(2k-2)$-th order for numerical flux of the derivative. In addition, we prove superconvergence of $(k+2)$ or $(k+3)$-th order of the DG solution towards a special projection. At a class of special points, the function values and the first and second order derivatives of the DG solution are superconvergent with order $k+2, k+1, k$, respectively. The proof relies on the correction function techniques initiated in [8], and applied to [6] for direct DG (DDG) methods for diffusion problems. Compared with [6], Schrödinger equation poses unique challenges for superconvergence proof because of the lack of the dissipation mechanism from the equation. One major highlight of our proof is that we introduce specially chosen test functions in the error equation and show the superconvergence of the second derivative and jump across the cell interfaces of the difference between numerical solution and projected exact solution. This technique was originally proposed in [12] and is essential to elevate the convergence order for our analysis. Finally, by negative norm estimates, we apply the post-processing technique and show that the accuracy of our scheme can be enhanced to order $2k.$ Theoretical results are verified by numerical experiments.

Wei Guo, Yingda Cheng

In this paper, we develop an adaptive multiresolution discontinuous Galerkin (DG) scheme for time-dependent transport equations in multi-dimensions. The method is constructed using multiwavlelets on tensorized nested grids. Adaptivity is realized by error thresholding based on the hierarchical surplus, and the Runge-Kutta DG (RKDG) scheme is employed as the reference time evolution algorithm. We show that the scheme performs similarly to a sparse grid DG method when the solution is smooth, reducing computational cost in multi-dimensions. When the solution is no longer smooth, the adaptive algorithm can automatically capture fine local structures. The method is therefore very suitable for deterministic kinetic simulations. Numerical results including several benchmark tests, the Vlasov-Poisson (VP) and oscillatory VP systems are provided.

Zhanjing Tao, Anqi Chen, Mengping Zhang et al.

In this paper, we develop sparse grid central discontinuous Galerkin (CDG) scheme for linear hyperbolic systems with variable coefficients in high dimensions. The scheme combines the CDG framework with the sparse grid approach, with the aim of breaking the curse of dimensionality. A new hierarchical representation of piecewise polynomials on the dual mesh is introduced and analyzed, resulting in a sparse finite element space that can be used for non-periodic problems. Theoretical results, such as $L^2$ stability and error estimates are obtained for scalar problems. CFL conditions are studied numerically comparing discontinuous Galerkin (DG), CDG, sparse grid DG and sparse grid CDG methods. Numerical results including scalar linear equations, acoustic and elastic waves are provided.

Peter DelMastro, Daniel Appelö, Yingda Cheng

We propose a family of low-rank, completely positive and trace preserving schemes for the Lindblad equation, a common model for open quantum systems. Low-rank representation is employed at two levels: the density matrix is factorized into the product of tall-skinny matrices, and the columns of these matrices are further represented using the tensor train (TT) format, also know as matrix product states (MPS). This two-level low-rank format fits naturally into our existing Kraus is King scheme (arXiv:2409.08898v2 [math.NA]) for the Lindblad equation, whose underlying operations are arithmetic on the columns of the tall-skinny matrices. We show how these operations can be performed efficiently in the TT/MPS format, with particular emphasis on density matrix rank-truncation. We conclude with extensive numerical experiments demonstrating the convergence of this scheme and its efficiency in simulating systems with up to $10^{19}$ degrees of freedom using only modest compute resources.

Jiuhua Hu, Daniel Appelo, Yingda Cheng

Numerical simulation of quantum computing hardware and open quantum systems governed by the Lindblad equation is challenging due to the high dimensionality of the density matrix and the need to preserve fundamental physical properties. In our previous work, we developed an arbitrary-order, low-rank, completely positive and trace preserving (CPTP) method for the Lindblad equation with time-dependent Hamiltonians by nested Picard iteration (NPI). In this work, we develop Gregory NPI schemes, which are CPTP schemes constructed by Gregory-type quadrature on equispaced nodes. The methods, which are of order up to nine, substantially reduce the computational cost compared to our previously proposed NPI schemes with Gaussian quadrature rules, while retaining high-order accuracy and structure preservation. We analyze the stability of the resulting scheme for a physics-based test equation. Numerical experiments verify the convergence of the method and demonstrate the effectiveness of the low-rank approximation. We study the performance of a previously constructed CNOT gate for both closed and open quantum systems.

Kun Huang, Yingda Cheng, Irene M. Gamba

We propose a fast deterministic scheme for the space-homogeneous Boltzmann equation that exploits the low-rank structure of the velocity distribution. This paper consists of two independent contributions. The first is a \emph{lifting-projection (LP) scheme}, inspired by the approach in the recent theoretical breakthroughs \cite{guillen2025landau, imbert2026monotonicity, guillen2025landau2} on the well-posedness of the Landau and Boltzmann equations. In particular, the approach lifts the nonlinear 3D Boltzmann equation to the 6D linear Kac master equation, advanced over a single time step, and projected back to its marginal in 3D. The second contribution is a \emph{low-rank tensor method} for evaluating the collision operator, in which the lifted solution is represented in tensor train (TT) format and computed via a TT cross approximation algorithm with interpolation, complemented by a TT-friendly conservation correction that enforces conservation of mass, momentum, and energy. When the solution is low-rank in velocity, the method scales linearly in $n$ when cubic interpolation is used (and quadratic in $n$ when spectral interpolation is used), where $n$ is the number of grid points in each velocity direction. Therefore, our methods offer significant computational savings over existing deterministic solvers in such cases. Numerical experiments on 2D and 3D benchmarks, including the BKW exact solution and anisotropic initial data, confirm the computational scaling, the expected order of accuracy and verify the effectiveness of the conservation correction.

Juntao Huang, Yingda Cheng, Andrew J. Christlieb et al.

This is the third paper in a series in which we develop machine learning (ML) moment closure models for the radiative transfer equation (RTE). In our previous work \cite{huang2021gradient}, we proposed an approach to learn the gradient of the unclosed high order moment, which performs much better than learning the moment itself and the conventional $P_N$ closure. However, while the ML moment closure has better accuracy, it is not able to guarantee hyperbolicity and has issues with long time stability. In our second paper \cite{huang2021hyperbolic}, we identified a symmetrizer which leads to conditions that enforce that the gradient based ML closure is symmetrizable hyperbolic and stable over long time. The limitation of this approach is that in practice the highest moment can only be related to four, or fewer, lower moments. In this paper, we propose a new method to enforce the hyperbolicity of the ML closure model. Motivated by the observation that the coefficient matrix of the closure system is a lower Hessenberg matrix, we relate its eigenvalues to the roots of an associated polynomial. We design two new neural network architectures based on this relation. The ML closure model resulting from the first neural network is weakly hyperbolic and guarantees the physical characteristic speeds, i.e., the eigenvalues are bounded by the speed of light. The second model is strictly hyperbolic and does not guarantee the boundedness of the eigenvalues. Several benchmark tests including the Gaussian source problem and the two-material problem show the good accuracy, stability and generalizability of our hyperbolic ML closure model.

Juntao Huang, Yingda Cheng, Andrew J. Christlieb et al.

This is the second paper in a series in which we develop machine learning (ML) moment closure models for the radiative transfer equation (RTE). In our previous work \cite{huang2021gradient}, we proposed an approach to directly learn the gradient of the unclosed high order moment, which performs much better than learning the moment itself and the conventional $P_N$ closure. However, the ML moment closure model in \cite{huang2021gradient} is not able to guarantee hyperbolicity and long time stability. We propose in this paper a method to enforce the global hyperbolicity of the ML closure model. The main idea is to seek a symmetrizer (a symmetric positive definite matrix) for the closure system, and derive constraints such that the system is globally symmetrizable hyperbolic. It is shown that the new ML closure system inherits the dissipativeness of the RTE and preserves the correct diffusion limit as the Knunsden number goes to zero. Several benchmark tests including the Gaussian source problem and the two-material problem show the good accuracy, long time stability and generalizability of our globally hyperbolic ML closure model.

Juntao Huang, Yingda Cheng, Andrew J. Christlieb et al.

In this paper, we take a data-driven approach and apply machine learning to the moment closure problem for radiative transfer equation in slab geometry. Instead of learning the unclosed high order moment, we propose to directly learn the gradient of the high order moment using neural networks. This new approach is consistent with the exact closure we derive for the free streaming limit and also provides a natural output normalization. A variety of benchmark tests, including the variable scattering problem, the Gaussian source problem with both periodic and reflecting boundaries, and the two-material problem, show both good accuracy and generalizability of our machine learning closure model.