We present a class of high-order Eulerian-Lagrangian Runge-Kutta finite volume methods that can numerically solve Burgers' equation with shock formations, which could be extended to general scalar conservation laws. Eulerian-Lagrangian (EL) and semi-Lagrangian (SL) methods have recently seen increased development and have become a staple for allowing large time-stepping sizes. Yet, maintaining relatively large time-stepping sizes post shock formation remains quite challenging. Our proposed scheme integrates the partial differential equation on a space-time region partitioned by linear approximations to the characteristics determined by the Rankine-Hugoniot jump condition. We trace the characteristics forward in time and present a merging procedure for the mesh cells to handle intersecting characteristics due to shocks. Following this partitioning, we write the equation in a time-differential form and evolve with Runge-Kutta methods in a method-of-lines fashion. High-resolution methods such as ENO and WENO-AO schemes are used for spatial reconstruction. Extension to higher dimensions is done via dimensional splitting. Numerical experiments demonstrate our scheme's high-order accuracy and ability to sharply capture post-shock solutions with large time-stepping sizes.
We extend the fourth order, two stage Multi-Derivative Runge Kutta (MDRK) scheme to the Flux Reconstruction (FR) framework by writing both stages in terms of a time averaged flux and then using the approximate Lax-Wendroff procedure to compute the time averaged flux. Numerical flux is carefully constructed to enhance Fourier CFL stability and accuracy. A subcell based blending limiter is developed for the MDRK scheme which ensures that the limited scheme is provably admissibility preserving. Along with being admissibility preserving, the blending scheme is constructed to minimize dissipation errors by using Gauss-Legendre solution points and performing MUSCL-Hancock reconstruction on subcells. The accuracy enhancement of the blending scheme is numerically verified on compressible Euler's equations, with test cases involving shocks and small-scale structures.
We propose a high order discontinuous Galerkin (DG) scheme with subcell finite volume (FV) limiter to solve a monolithic first--order hyperbolic BSSNOK formulation of the coupled Einstein--Euler equations. The numerical scheme runs with adaptive mesh refinement (AMR) in three space dimensions, is endowed with time-accurate local time stepping (LTS) and is able to deal with both conservative and non-conservative hyperbolic systems. The system of governing partial differential equations was shown to be strongly hyperbolic and is solved in a monolithic fashion with one numerical framework that can be simultaneously applied to both the conservative matter subsystem as well as the non-conservative subsystem for the spacetime. Since high order unlimited DG schemes are well-known to produce spurious oscillations in the presence of discontinuities and singularities, our subcell finite volume limiter is crucial for the robust discretization of shock waves arising in the matter as well as for the stable treatment of puncture black holes. We test the new method on a set of classical test problems of numerical general relativity, showing good agreement with available exact or numerical reference solutions. In particular, we perform the first long term evolution of the inspiralling merger of two puncture black holes with a high order ADER-DG scheme.
This work concerns the analysis of the discontinuous Galerkin spectral element method (DGSEM) with implicit time stepping for the numerical approximation of nonlinear scalar conservation laws in multiple space dimensions. We consider either the DGSEM with a backward Euler time stepping, or a space-time DGSEM discretization to remove the restriction on the time step. We design graph viscosities in space, and in time for the space-time DGSEM, to make the schemes maximum principle preserving and entropy stable for every admissible convex entropy. We also establish well-posedness of the discrete problems by showing existence and uniqueness of the solutions to the nonlinear implicit algebraic relations that need to be solved at each time step. Numerical experiments in one space dimension are presented to illustrate the properties of these schemes.
We consider weak convergence of one-step schemes for solving stochastic differential equations (SDEs) with one-sided Lipschitz conditions. It is known that the super-linear coefficients may lead to a blowup of moments of solutions and their numerical solutions. When solutions to SDEs have all finite moments, weak convergence of numerical schemes has been investigated in [Wang et al (2023), Weak error analysis for strong approximation schemes of SDEs with super-linear coefficients, IMA Journal numerical analysis]. Some modified Euler schemes have been analyzed for weak convergence. In this work, we present a family of explicit schemes of first and second-order weak convergence based on classical schemes for SDEs. We explore the effects of limited moments on these schemes. We provide a systematic but simple way to establish weak convergence orders for schemes based on approximations/modifications of drift and diffusion coefficients. We present several numerical examples of these schemes and show their weak convergence orders.
We derive the Alternating-Direction Implicit (ADI) method based on a commuting operator split and apply the results to the continuous time algebraic Lyapunov equation with low-rank constant term and approximate solution. Previously, it has been mandatory to start the low-rank ADI (LR-ADI) with an all-zero initial value. Our approach retains the known efficient iteration schemes of low-rank increments and residual to arbitrary low-rank initial values for the LR-ADI method. We further generalize some of the known properties of the LR-ADI for Lyapunov equations to larger classes of algorithms or problems. We investigate the performance of arbitrary initial values using two outer iterations in which LR-ADI is typically called. First, we solve an algebraic Riccati equation with the Newton method. Second, we solve a differential Riccati equation with a first-order Rosenbrock method. Numerical experiments confirm that the proposed new initial value of the alternating-directions implicit (ADI) can lead to a significant reduction in the total number of ADI steps, while also showing a 17% and 8x speed-up over the zero initial value for the two equation types, respectively.
In this paper, we consider the task of efficiently computing the numerical solution of evolutionary complex Ginzburg--Landau equations on Cartesian product domains with homogeneous Dirichlet/Neumann or periodic boundary conditions. To this aim, we employ for the time integration high-order exponential methods of splitting and Lawson type with constant time step size. These schemes enjoy favorable stability properties and, in particular, do not show restrictions on the time step size due to the underlying stiffness of the models. The needed actions of matrix exponentials are efficiently realized by using a tensor-oriented approach that suitably employs the so-called $\mu$-mode product (when the semidiscretization in space is performed with finite differences) or with pointwise operations in Fourier space (when the model is considered with periodic boundary conditions). The overall effectiveness of the approach is demonstrated by running simulations on a variety of two- and three-dimensional (systems of) complex Ginzburg--Landau equations with cubic or cubic-quintic nonlinearities, which are widely considered in literature to model relevant physical phenomena. In fact, we show that high-order exponential-type schemes may outperform standard techniques to integrate in time the models under consideration, i.e., the well-known second-order split-step method and the explicit fourth-order Runge--Kutta integrator, for stringent accuracies.
Two numerical schemes are proposed and investigated for the Yang--Mills equations, which can be seen as a nonlinear generalisation of the Maxwell equations set on Lie algebra-valued functions, with similarities to certain formulations of General Relativity. Both schemes are built on the Discrete de Rham (DDR) method, and inherit from its main features: an arbitrary order of accuracy, and applicability to generic polyhedral meshes. They make use of the complex property of the DDR, together with a Lagrange-multiplier approach, to preserve, at the discrete level, a nonlinear constraint associated with the Yang--Mills equations. We also show that the schemes satisfy a discrete energy dissipation (the dissipation coming solely from the implicit time stepping). Issues around the practical implementations of the schemes are discussed; in particular, the assembly of the local contributions in a way that minimises the price we pay in dealing with nonlinear terms, in conjunction with the tensorisation coming from the Lie algebra. Numerical tests are provided using a manufactured solution, and show that both schemes display a convergence in $L^2$-norm of the potential and electrical fields in $\mathcal O(h^{k+1})$ (provided that the time step is of that order), where $k$ is the polynomial degree chosen for the DDR complex. We also numerically demonstrate the preservation of the constraint.
We present a new, monolithic first--order (both in time and space) BSSNOK formulation of the coupled Einstein--Euler equations. The entire system of hyperbolic PDEs is solved in a completely unified manner via one single numerical scheme applied to both the conservative sector of the matter part and to the first--order strictly non--conservative sector of the spacetime evolution. The coupling between matter and space-time is achieved via algebraic source terms. The numerical scheme used for the solution of the new monolithic first order formulation is a path-conservative central WENO (CWENO) finite difference scheme, with suitable insertions to account for the presence of the non--conservative terms. By solving several crucial tests of numerical general relativity, including a stable neutron star, Riemann problems in relativistic matter with shock waves and the stable long-time evolution of single and binary puncture black holes up and beyond the binary merger, we show that our new CWENO scheme, introduced two decades ago for the compressible Euler equations of gas dynamics, can be successfully applied also to numerical general relativity, solving all equations at the same time with one single numerical method. In the future the new monolithic approach proposed in this paper may become an attractive alternative to traditional methods that couple central finite difference schemes with Kreiss-Oliger dissipation for the space-time part with totally different TVD schemes for the matter evolution and which are currently the state of the art in the field.
We studied the dynamical properties of Rabi oscillations driven by an alternating Rashba field applied to a two-dimensional (2D) harmonic confinement system. We solve the time-dependent (TD) Schr\"{o}dinger equation numerically and rewrite the resulting TD wavefunction onto the Bloch sphere (BS) using two BS parameters of the zenith ($\theta_B$) and azimuthal ($\phi_B$) angles, extracting the phase information $\phi_B$ as well as the mixing ratio $\theta_B$ between the two BS-pole states. We employed a two-state rotating wave (TSRW) approach and studied the fundamental features of $\theta_B$ and $\phi_B$ over time. The TSRW approach reveals a triangular wave formation in $\theta_B$. Moreover, at each apex of the triangular wave, the TD wavefunction passes through the BS pole, and the state is completely replaced by the opposite spin state. The TSRW approach also elucidates a linear change in $\phi_B$. The slope of $\phi_B$ vs. time is equal to the difference between the dynamical terms, leading to a confinement potential in the harmonic system. The TSRW approach further demonstrates a jump in the phase difference by $\pi$ when the wavefunction passes through the BS pole. The alternating Rashba field causes multiple successive Rabi transitions in the 2D harmonic system. We then introduce the effective BS (EBS) and transform these complicated transitions into an equivalent "single" Rabi one. Consequently, the EBS parameters $\theta_B^{\mathrm{eff}}$ and $\phi_B^{\mathrm{eff}}$ exhibit mixing and phase difference between two spin states $\alpha$ and $\beta$, leading to a deep understanding of the TD features of multi-Rabi oscillations. Furthermore, the combination of the BS representation with the TSRW approach successfully reveals the dynamical properties of the Rabi oscillation, even beyond the TSRW approximation.
In this paper, we innovatively develop uniform/variable-time-step weighted and shifted BDF2 (WSBDF2) methods for the anisotropic Cahn-Hilliard (CH) model, combining the scalar auxiliary variable (SAV) approach with two types of stabilized techniques. Using the concept of $G$-stability, the uniform-time-step WSBDF2 method is theoretically proved to be energy-stable. Due to the inapplicability of the relevant G-stability properties, another technique is adopted in this work to demonstrate the energy stability of the variable-time-step WSBDF2 method. In addition, the two numerical schemes are all mass-conservative.Finally, numerous numerical simulations are presented to demonstrate the stability and accuracy of these schemes.