High-dimensional transport equations frequently occur in science and engineering. Computing their numerical solution, however, is challenging due to its high dimensionality. In this work we develop an algorithm to efficiently solve the transport equation in moderately complex geometrical domains using a Galerkin method stabilized by streamline diffusion. The ansatz spaces are a tensor product of a sparse grid in space and discontinuous piecewise polynomials in time. Here, the sparse grid is constructed upon nested multilevel finite element spaces to provide geometric flexibility. This results in an implicit time-stepping scheme which we prove to be stable and convergent. If the solution has additional mixed regularity, the convergence of a $2d$-dimensional problem equals that of a $d$-dimensional one up to logarithmic factors. For the implementation, we rely on the representation of sparse grids as a sum of anisotropic full grid spaces. This enables us to store the functions and to carry out the computations on a sequence regular full grids exploiting the tensor product structure of the ansatz spaces. In this way existing finite element libraries and GPU acceleration can be used. The combination technique is used as a preconditioner for an iterative scheme to solve the transport equation on the sequence of time strips. Numerical tests show that the method works well for problems in up to six dimensions. Finally, the method is also used as a building block to solve nonlinear Vlasov-Poisson equations.
The numerical solution of continuum damage mechanics (CDM) problems suffers from convergence-related challenges during the material softening stage, and consequently existing iterative solvers are subject to a trade-off between computational expense and solution accuracy. In this work, we present a novel unified arc-length (UAL) method, and we derive the formulation of the analytical tangent matrix and governing system of equations for both local and non-local gradient damage problems. Unlike existing versions of arc-length solvers that monolithically scale the external force vector, the proposed method treats the latter as an independent variable and determines the position of the system on the equilibrium path based on all the nodal variations of the external force vector. This approach renders the proposed solver substantially more efficient and robust than existing solvers used in CDM problems. We demonstrate the considerable advantages of the proposed algorithm through several benchmark 1D problems with sharp snap-backs and 2D examples under various boundary conditions and loading scenarios. The proposed UAL approach exhibits a superior ability of overcoming critical increments along the equilibrium path. Moreover, the proposed UAL method is 1-2 orders of magnitude faster than force-controlled arc-length and monolithic Newton-Raphson solvers.
Over the last two decades, the field of geometric curve evolutions has attracted significant attention from scientific computing. One of the most popular numerical methods for solving geometric flows is the so-called BGN scheme, which was proposed by Barrett, Garcke, and Nurnberg (J. Comput. Phys., 222 (2007), pp. 441{467), due to its favorable properties (e.g., its computational efficiency and the good mesh property). However, the BGN scheme is limited to first-order accuracy in time, and how to develop a higher-order numerical scheme is challenging. In this paper, we propose a fully discrete, temporal second-order parametric finite element method, which incorporates a mesh regularization technique when necessary, for solving geometric flows of curves. The scheme is constructed based on the BGN formulation and a semi-implicit Crank-Nicolson leap-frog time stepping discretization as well as a linear finite element approximation in space. More importantly, we point out that the shape metrics, such as manifold distance and Hausdorff distance, instead of function norms, should be employed to measure numerical errors. Extensive numerical experiments demonstrate that the proposed BGN-based scheme is second-order accurate in time in terms of shape metrics. Moreover, by employing the classical BGN scheme as a mesh regularization technique when necessary, our proposed second-order scheme exhibits good properties with respect to the mesh distribution.
We develop easily accessible quantities for bounding the almost sure exponential convergence rate of push-sum algorithms. We analyze the scenario of i.i.d. synchronous gossip, every agent communicating towards its single target at every step. Multiple bounding expressions are developed depending on the generality of the setup, all functions of the spectrum of the network. While the most general bound awaits further improvement, with more symmetries, close bounds can be established, as demonstrated by numerical simulations.
The unified gas-kinetic wave-particle method (UGKWP) has been developed for the multiscale gas, plasma, and multiphase flow transport processes for the past years. In this work, we propose an implicit unified gas-kinetic wave-particle (IUGKWP) method to remove the CFL time step constraint. Based on the local integral solution of the radiative transfer equation (RTE), the particle transport processes are categorized into the long-$\lambda$ streaming process and the short-$\lambda$ streaming process comparing to a local physical characteristic time $t_p$. In the construction of the IUGKWP method, the long-$\lambda$ streaming process is tracked by the implicit Monte Carlo (IMC) method; the short-$\lambda$ streaming process is evolved by solving the implicit moments equations; and the photon distribution is closed by a local integral solution of RTE. In the IUGKWP method, the multiscale flux of radiation energy and the multiscale closure of photon distribution are constructed based on the local integral solution. The IUGKWP method preserves the second-order asymptotic expansion of RTE in the optically thick regime and adapts its computational complexity to the flow regime. The numerical dissipation is well controlled, and the teleportation error is significantly reduced in the optically thick regime. The computational complexity of the IUGKWP method decreases exponentially as the Knudsen number approaches zero, and the computational efficiency is remarkably improved in the optically thick regime. The IUGKWP is formulated on a generalized unstructured mesh, and multidimensional 2D and 3D algorithms are developed. Numerical tests are presented to validate the capability of IUGKWP in capturing the multiscale photon transport process. The algorithm and code will apply in the engineering applications of inertial confinement fusion (ICF).
Sequential models, such as Recurrent Neural Networks and Neural Ordinary Differential Equations, have long suffered from slow training due to their inherent sequential nature. For many years this bottleneck has persisted, as many thought sequential models could not be parallelized. We challenge this long-held belief with our parallel algorithm that accelerates GPU evaluation of sequential models by up to 3 orders of magnitude faster without compromising output accuracy. The algorithm does not need any special structure in the sequential models' architecture, making it applicable to a wide range of architectures. Using our method, training sequential models can be more than 10 times faster than the common sequential method without any meaningful difference in the training results. Leveraging this accelerated training, we discovered the efficacy of the Gated Recurrent Unit in a long time series classification problem with 17k time samples. By overcoming the training bottleneck, our work serves as the first step to unlock the potential of non-linear sequential models for long sequence problems.
We evaluate the performance of novel numerical methods for solving one-dimensional nonlinear fractional dispersive and dissipative evolution equations. The methods are based on affine combinations of time-splitting integrators and pseudo-spectral discretizations using Hermite and Fourier expansions. We show the effectiveness of the proposed methods by numerically computing the dynamics of soliton solutions of the the standard and fractional variants of the nonlinear Schr{\"o}dinger equation (NLSE) and the complex Ginzburg-Landau equation (CGLE), and by comparing the results with those obtained by standard splitting integrators. An exhaustive numerical investigation shows that the new technique is competitive when compared to traditional composition-splitting schemes for the case of Hamiltonian problems both in terms accuracy and computational cost. Moreover, it is applicable straightforwardly to irreversible models, outperforming high-order symplectic integrators which could become unstable due to their need of negative time steps. Finally, we discuss potential improvements of the numerical methods aimed to increase their efficiency, and possible applications to the investigation of dissipative solitons that arise in nonlinear optical systems of contemporary interest. Overall, the method offers a promising alternative for solving a wide range of evolutionary partial differential equations.
We present here a new splitting method to solve Lyapunov equations in a Kronecker product form. Although this resulting matrix is of order $n^2$, each iteration demands two operations with the matrix $A$: a multiplication of the form $(A-\sigma I) \tilde{B}$ and a inversion of the form $(A-\sigma I)^{-1}\tilde{B}$. We see that for some choice of a parameter the iteration matrix is such that all their eigenvalues are in absolute value less than 1. Moreover we present a theorem that enables us to get a good starting vector for the method.
Solving linear systems is of great importance in numerous fields. In particular, circulant systems are especially valuable for efficiently finding numerical solutions to physics-related differential equations. Current quantum algorithms like HHL or variational methods are either resource-intensive or may fail to find a solution. We present an efficient algorithm based on convex optimization of combinations of quantum states to solve for banded circulant linear systems whose non-zero terms are within distance $K$ of the main diagonal. By decomposing banded circulant matrices into cyclic permutations, our approach produces approximate solutions to such systems with a combination of quantum states linear to $K$, significantly improving over previous convergence guarantees, which require quantum states exponential to $K$. We propose a hybrid quantum-classical algorithm using the Hadamard test and the quantum Fourier transform as subroutines and show its PromiseBQP-hardness. Additionally, we introduce a quantum-inspired algorithm with similar performance given sample and query access. We validate our methods with classical simulations and actual IBM quantum computer implementation, showcasing their applicability for solving physical problems such as heat transfer.
A method of numerically solving the Maxwell equations is considered for modeling harmonic electromagnetic fields. The vector finite element method makes it possible to obtain a physically consistent discretization of the differential equations. However, solving large systems of linear algebraic equations with indefinite ill-conditioned matrices is a challenge. The high order of the matrices limits the capabilities of the Gaussian method to solve such systems, since this requires large RAM and much calculation. To reduce these requirements, an iterative preconditioned algorithm based on integral Laguerre transform in time is used. This approach allows using multigrid algorithms and, as a result, needs less RAM compared to the direct methods of solving systems of linear algebraic equations.
We propose an approach to compute inner and outer-approximations of the sets of values satisfying constraints expressed as arbitrarily quantified formulas. Such formulas arise for instance when specifying important problems in control such as robustness, motion planning or controllers comparison. We propose an interval-based method which allows for tractable but tight approximations. We demonstrate its applicability through a series of examples and benchmarks using a prototype implementation.