A linearly implicit conservative difference scheme is applied to discretize the attractive coupled nonlinear Schr\"odinger equations with fractional Laplacian. Complex symmetric linear systems can be obtained, and the system matrices are indefinite and Toeplitz-plus-diagonal. Neither efficient preconditioned iteration method nor fast direct method is available to deal with these systems. In this paper, we propose a novel matrix splitting iteration method based on a normal splitting of an equivalent real block form of the complex linear systems. This new iteration method converges unconditionally, and the quasi-optimal iteration parameter is deducted. The corresponding new preconditioner is obtained naturally, which can be constructed easily and implemented efficiently by fast Fourier transform. Theoretical analysis indicates that the eigenvalues of the preconditioned system matrix are tightly clustered. Numerical experiments show that the new preconditioner can significantly accelerate the convergence rate of the Krylov subspace iteration methods. Specifically, the convergence behavior of the related preconditioned GMRES iteration method is spacial mesh-size-independent, and almost fractional order insensitive. Moreover, the linearly implicit conservative difference scheme in conjunction with the preconditioned GMRES iteration method conserves the discrete mass and energy in terms of a given precision.
Quadratization of polynomial and nonpolynomial systems of ordinary differential equations is advantageous in a variety of disciplines, such as systems theory, fluid mechanics, chemical reaction modeling and mathematical analysis. A quadratization reveals new variables and structures of a model, which may be easier to analyze, simulate, control, and provides a convenient parametrization for learning. This paper presents novel theory, algorithms and software capabilities for quadratization of non-autonomous ODEs. We provide existence results, depending on the regularity of the input function, for cases when a quadratic-bilinear system can be obtained through quadratization. We further develop existence results and an algorithm that generalizes the process of quadratization for systems with arbitrary dimension that retain the nonlinear structure when the dimension grows. For such systems, we provide dimension-agnostic quadratization. An example is semi-discretized PDEs, where the nonlinear terms remain symbolically identical when the discretization size increases. As an important aspect for practical adoption of this research, we extended the capabilities of the QBee software towards both non-autonomous systems of ODEs and ODEs with arbitrary dimension. We present several examples of ODEs that were previously reported in the literature, and where our new algorithms find quadratized ODE systems with lower dimension than the previously reported lifting transformations. We further highlight an important area of quadratization: reduced-order model learning. This area can benefit significantly from working in the optimal lifting variables, where quadratic models provide a direct parametrization of the model that also avoids additional hyperreduction for the nonlinear terms. A solar wind example highlights these advantages.
We consider the low-rank alternating directions implicit (ADI) iteration for approximately solving large-scale algebraic Sylvester equations. Inside every iteration step of this iterative process a pair of linear systems of equations has to be solved. We investigate the situation when those inner linear systems are solved inexactly by an iterative methods such as, for example, preconditioned Krylov subspace methods. The main contribution of this work are thresholds for the required accuracies regarding the inner linear systems which dictate when the employed inner Krylov subspace methods can be safely terminated. The goal is to save computational effort by solving the inner linear system as inaccurate as possible without endangering the functionality of the low-rank Sylvester-ADI method. Ideally, the inexact ADI method mimics the convergence behaviour of the more expensive exact ADI method, where the linear systems are solved directly. Alongside the theoretical results, also strategies for an actual practical implementation of the stopping criteria are developed. Numerical experiments confirm the effectiveness of the proposed strategies.
Regularization is a long-standing challenge for ill-posed linear inverse problems, and a prototype is the Fredholm integral equation of the first kind with additive Gaussian measurement noise. We introduce a new RKHS regularization adaptive to measurement data and the underlying linear operator. This RKHS arises naturally in a variational approach, and its closure is the function space in which we can identify the true solution. Also, we introduce a small noise analysis to compare regularization norms by sharp convergence rates in the small noise limit. Our analysis shows that the RKHS- and $L^2$-regularizers yield the same convergence rate when their optimal hyper-parameters are selected using the true solution, and the RKHS-regularizer has a smaller multiplicative constant. However, in computational practice, the RKHS regularizer significantly outperforms the $L^2$-and $l^2$-regularizers in producing consistently converging estimators when the noise level decays or the observation mesh refines.
This work proposes novel techniques for the efficient numerical simulation of parameterized, unsteady partial differential equations. Projection-based reduced order models (ROMs) such as the reduced basis method employ a (Petrov-)Galerkin projection onto a linear low-dimensional subspace. In unsteady applications, space-time reduced basis (ST-RB) methods have been developed to achieve a dimension reduction both in space and time, eliminating the computational burden of time marching schemes. However, nonaffine parameterizations dilute any computational speedup achievable by traditional ROMs. Computational efficiency can be recovered by linearizing the nonaffine operators via hyper-reduction, such as the empirical interpolation method in matrix form. In this work, we implement new hyper-reduction techniques explicitly tailored to deal with unsteady problems and embed them in a ST-RB framework. For each of the proposed methods, we develop a posteriori error bounds. We run numerical tests to compare the performance of the proposed ROMs against high-fidelity simulations, in which we combine the finite element method for space discretization on 3D geometries and the Backward Euler time integrator. In particular, we consider a heat equation and an unsteady Stokes equation. The numerical experiments demonstrate the accuracy and computational efficiency our methods retain with respect to the high-fidelity simulations.
We present an efficient matrix-free geometric multigrid method for the elastic Helmholtz equation, and a suitable discretization. Many discretization methods had been considered in the literature for the Helmholtz equations, as well as many solvers and preconditioners, some of which are adapted for the elastic version of the equation. However, there is very little work considering the reciprocity of discretization and a solver. In this work, we aim to bridge this gap. By choosing an appropriate stencil for re-discretization of the equation on the coarse grid, we develop a multigrid method that can be easily implemented as matrix-free, relying on stencils rather than sparse matrices. This is crucial for efficient implementation on modern hardware. Using two-grid local Fourier analysis, we validate the compatibility of our discretization with our solver, and tune a choice of weights for the stencil for which the convergence rate of the multigrid cycle is optimal. It results in a scalable multigrid preconditioner that can tackle large real-world 3D scenarios.
We propose a method for computing the Lyapunov exponents of renewal equations (delay equations of Volterra type) and of coupled systems of renewal and delay differential equations. The method consists in the reformulation of the delay equation as an abstract differential equation, the reduction of the latter to a system of ordinary differential equations via pseudospectral collocation, and the application of the standard discrete QR method. The effectiveness of the method is shown experimentally and a MATLAB implementation is provided.
We analyze a goal-oriented adaptive algorithm that aims to efficiently compute the quantity of interest $G(u^\star)$ with a linear goal functional $G$ and the solution $u^\star$ to a general second-order nonsymmetric linear elliptic partial differential equation. The current state of the analysis of iterative algebraic solvers for nonsymmetric systems lacks the contraction property in the norms that are prescribed by the functional analytic setting. This seemingly prevents their application in the optimality analysis of goal-oriented adaptivity. As a remedy, this paper proposes a goal-oriented adaptive iteratively symmetrized finite element method (GOAISFEM). It employs a nested loop with a contractive symmetrization procedure, e.g., the Zarantonello iteration, and a contractive algebraic solver, e.g., an optimal multigrid solver. The various iterative procedures require well-designed stopping criteria such that the adaptive algorithm can effectively steer the local mesh refinement and the computation of the inexact discrete approximations. The main results consist of full linear convergence of the proposed adaptive algorithm and the proof of optimal convergence rates with respect to both degrees of freedom and total computational cost (i.e., optimal complexity). Numerical experiments confirm the theoretical results and investigate the selection of the parameters.
Of all the possible projection methods for solving large-scale Lyapunov matrix equations, Galerkin approaches remain much more popular than Petrov-Galerkin ones. This is mainly due to the different nature of the projected problems stemming from these two families of methods. While a Galerkin approach leads to the solution of a low-dimensional matrix equation per iteration, a matrix least-squares problem needs to be solved per iteration in a Petrov-Galerkin setting. The significant computational cost of these least-squares problems has steered researchers towards Galerkin methods in spite of the appealing minimization properties of Petrov-Galerkin schemes. In this paper we introduce a framework that allows for modifying the Galerkin approach by low-rank, additive corrections to the projected matrix equation problem with the two-fold goal of attaining monotonic convergence rates similar to those of Petrov-Galerkin schemes while maintaining essentially the same computational cost of the original Galerkin method. We analyze the well-posedness of our framework and determine possible scenarios where we expect the residual norm attained by two low-rank-modified variants to behave similarly to the one computed by a Petrov-Galerkin technique. A panel of diverse numerical examples shows the behavior and potential of our new approach.
We develop a novel discontinuous Galerkin method for solving the rotating thermal shallow water equations (TRSW) on a curvilinear mesh. Our method is provably entropy stable, conserves mass, buoyancy and vorticity, while also semi-discretely conserving energy. This is achieved by using novel numerical fluxes and splitting the pressure and convection operators. We implement our method on a cubed sphere mesh and numerically verify our theoretical results. Our experiments demonstrate the robustness of the method for a regime of well developed turbulence, where it can be run stably without any dissipation. The entropy stable fluxes are sufficient to control the grid scale noise generated by geostrophic turbulence, eliminating the need for artificial stabilization.
This manuscript deals with the analysis of numerical methods for the full discretization (in time and space) of the linear heat equation with Neumann boundary conditions, and it provides the reader with error estimates that are uniform in time. First, we consider the homogeneous equation with homogeneous Neumann boundary conditions over a finite interval. Using finite differences in space and the Euler method in time, we prove that our method is of order 1 in space, uniformly in time, under a classical CFL condition, and despite its lack of consistency at the boundaries. Second, we consider the nonhomogeneous equation with nonhomogeneous Neumann boundary conditions over a finite interval. Using a tailored similar scheme, we prove that our method is also of order 1 in space, uniformly in time, under a classical CFL condition. We indicate how this numerical method allows for a new way to compute steady states of such equations when they exist. We conclude by several numerical experiments to illustrate the sharpness and relevance of our theoretical results, as well as to examine situations that do not meet the hypotheses of our theoretical results, and to illustrate how our results extend to higher dimensions.