This work presents an adaptive time-stepping method for the incompressible Navier-Stokes equations (NSE) based on the deferred correction scheme. It is motivated by the capability of increasing the order by one at each iteration. This capability not only enhances accuracy but also naturally provides a framework for developing adaptive methods. Criterion for adaptively adjusting the time step is derived from approximations of local error estimation that utilizes the difference between consecutive solutions of the deferred correction scheme. By employing this criterion, we can effectively capture the dynamics of the flow. In addition, we prove unconditional nonlinear energy stability for the predicted step in the variable timestep scheme and also provide an error analysis. Finally, several challenging numerical tests are conducted to demonstrate the accuracy, effectiveness and robustness of our proposed adaptive scheme.
Let G=(V,E) be a connected graph. A dominating set is a subset D ⊂ V such that every vertex of G is either in D or adjacent to a vertex in D. A connected dominating set is a dominating set in G which induces a connected subgraph. The connected domatic number of G is the maximum number of pairwise disjoint, connected dominating sets in V(G). Finding the connected domatic number of general graphs is NP-hard. In this paper, we study the connected domatic number for a well-known family of graphs - the generalized Petersen graphs GP(n, k). Determining the connected domatic number of GP(n, k) is equivalent to determining whether it has two disjoint connected dominating sets. We identify several classes of generalized Petersen graphs with two disjoint connected dominating sets. Moreover, for small k (2 ≤ k ≤ 5), we provide necessary and sufficient conditions for a generalized Petersen graph GP(n, k) to contain two disjoint connected dominating cycles.
This paper presents a structure-preserving space-time POD reduction method for solving variable-coefficient parabolic equations. A space-time coupled full-order model with Kroneckerproduct structure is first derived through space-time finite element discretization, leading to a space-time reduced-order model. This formulation inherently possesses spatiotemporal separation properties and enables simultaneous dimensionality reduction in both space and time. Moreover, an error estimate between the numerical solution and the reduced solution is provided. Numerical experiments demonstrate the effectiveness of the proposed approach. Compared with the existing space reduction techniques, the advantages of the proposed method lie in the derivation of a space-time coupled model with Kronecker-product form, which yields a natural spatiotemporal separable structure that is inherently well-suited for space-time model order reduction.
In this paper, we propose a Riemannian conjugate gradient model order reduction algorithm for linear switched systems. An analytical expression for the error between the original switched system and its reduced-order system in the $\mathcal {H}_{2}$-norm is derived, where involved matrices exhibit a block diagonal structure. We reformulate the model order reduction task as an optimization problem with orthogonality constraints, and transform it into an unconstrained Riemannian optimization framework by leveraging the Stiefel manifold. By computing the Riemannian gradient of the objective function, a conjugate gradient direction on the Stiefel manifold is constructed. The proposed algorithm avoids redundant vector transport computations. Theoretically, it allows parallel implementation by exploiting the block diagonal structure (mode-independence) of the matrices. Global convergence of the algorithm is analyzed, and two numerical examples are conducted to validate the efficiency of the proposed algorithm.
This paper presents a novel interpolation-based model order reduction technique designed specifically for high-order polynomial systems. By building upon the existing two-sided projection interpolation approach, a new subspace construction method is proposed. In this method, the basis of the matrix successfully retains the crucial property of matching the generalized moments of the polynomial system. To further exploit the characteristics of high-order tensors, a tensor decomposition method, namely t-eigendecomposition, is introduced based on tensor t-product and t-SVD. This allows for the derivation of low-rank approximations for the Hessian tensor H and cross term tensor N, effectively capitalizing on the high-dimensional linear algebraic structure of the polynomial system and significantly reducing the complexity and storage requirements during the offline stage. The feasibility and efficiency of the proposed method are verified through a series of numerical experiments.
In this paper, we are concerned with two-scale integrators for the non-relativistic Klein-Gordon (NRKG) equation with a dimensionless parameter 0 < e << 1, which is inversely proportional to the speed of light. The highly oscillatory property in time of this model corresponds to the parameter e and the equation in the form of partial derivative(tt)u-Delta/epsilon(2)u + 1/epsilon(4)u + lambda/epsilon(2)f(u) = 0 has a factor 1/epsilon(2) in front of the nonlinearity which means that this part becomes strong when e is small. These two aspects bring significantly numerical burdens in designing numerical methods. We propose a class of two-scale integrators which is constructed based on some reformulations to the system, Fourier pseudo-spectral method and exponential integrators. Two practical integrators up to order three and four are constructed by using some symmetric conditions and the stiff order conditions of implicit exponential integrators. The convergence of the obtained integrators is rigorously studied, and it is shown that the uniform accuracy in time is O(h(3)) and O(h(4)) for the time stepsize h. The near energy conservation over long times is also established for the multi-stage integrators by using modulated Fourier expansions. Numerical results on a NRKG equation show that the proposed integrators have high accuracy, excellent long time energy conservation and competitive efficiency.
In this work, we study the discrete approximation of the continuous data assimilation (CDA) algorithm based on Runge-Kutta method and nudging approach for solving reaction-diffusion equations with unknown or incomplete initial data. For the spatial discretization, this paper only investigate the finite element method, although the investigation can be also applied to other spatial discretization methods. Under suitable conditions on the nudging parameter, the stability of the Runge-Kutta semi-discrete and fully discrete methods for data assimilation equations are obtained by exploring the algebraical stability of the methods. The uniform error estimates are then derived for the Runge-Kutta time semi-discrete and fully discrete CDA algorithms. The error estimates demonstrate that the time semi-discrete and fully discrete approximations converge to the true solution exponentially over time. A numerical study is also provided to support the theoretical findings.
Numerical simulation of time-periodic problems is a special area of research, since the time periodicity modifies the problem structure, and then it is desirable to use parallel methods to solve such problems. The classical parareal algorithm for time-periodic problems, which is parallel in time, solving an initial value coarse problem, called the periodic parareal algorithm with initial value coarse problem (PP-IC), usually converges very slowly, and even diverges for wave propagation problems. In this paper, we first present a new PP-IC algorithm based on a diagonalization technique proposed recently. In this new algorithm, we approximate the coarse propagator G in the classical PP-IC algorithm with a head-tail coupled condition such that G can be parallelized using diagonalization in time. We analyze the convergence factors of the diagonalization-based PP-IC algorithm for both the linear and nonlinear cases. Then, we further design and analyze a new parallel-intime algorithm for time-periodic problems by combining the Krylov subspace method with the diagonalization-based PP-IC algorithm to accelerate the convergence. Finally, we also determine an appropriate choice of the parameter alpha in the head-tail coupling condition, and illustrate our theoretical results with several numerical experiments, both for model problems and the realistic application of Maxwell's equations.
A model order reduction method based on shifted Legendre polynomials for solving convection-diffusion equations with variable coefficients is presented in this paper. The ordinary differential system of the convection-diffusion equation is obtained by finite element discretization procedure. Then, approximating the system state via shifted Legendre polynomials, the reduced-order system is produced, which can be solved efficiently to obtain the numerical solution. Error analysis is presented, and numerical examples are used to verify the feasibility of the presented method.
Small-signal stability analysis focuses on evaluating the ability of the power system to maintain synchronous operation and return to steady state under small disturbances. To address the escalating computational complexity and storage requirements caused by exponential growth in system dimensionality, this paper proposes a structure-preserving low-rank balanced truncation method via ε -embedding procedure and Hermite polynomials for small-signal stability analysis of power systems. This method employs Hermite polynomials for low-rank approximation of the controllability and observability Gramians to avoid solving two large-scale Lyapunov equations to compute the Gramians, significantly reducing the computational complexity. This approach reduces system dimensionality from thousands to tens of dimensions while preserving dominant dynamic characteristics and capturing the response of the system to small disturbances. The efficiency of the proposed algorithm is demonstrated through a numerical example, with the reduced-order model preserving the important properties of the original high-order model, including time and frequency domain responses as well as eigenvalues.
In the field of financial derivative pricing, option pricing under stochastic market environments has long been a core focus of both academic research and practical application. However, traditional numerical methods often face severe challenges in balancing computational accuracy, efficiency, and stability when dealing with the partial integro-differential equations (PIDEs) derived from jump-diffusion models, especially under non-smooth initial conditions. In this paper, we consider the stochastic Merton’s jump-diffusion model for European option pricing, which can be transformed into a PIDE. In view of the non-smoothness of the initial data, a variable-step Crank-Nicolson (CN) method is proposed for temporal discretization, and the quadratic spline collocation (QSC) method is used for spatial discretization. The prior estimate and convergence of the variable-step QSC-CN method are presented. Since a spatial non-local integration is handled implicitly, which leads to dense algebraic equations. Then a matrix-free fast Krylov subspace iterative solver is proposed to increase the computational efficiency, and preconditioning technique is applied to further accelerate the convergence. We show that the preconditioned fast QSC-CN method significantly reduces the computational cost from O((M+2)3) to O((M+2)log(M+2)), and the memory requirement from O((M+2)2) to O((M+2)), where M is the number of spatial discretization grids. This numerical method is further generalized for American option pricing. Numerical experiments are attached to verify the theoretical convergence order and the effectiveness of the fast computing technique.
This paper presents projection-based model order reduction methods for time-delay systems with initial history functions, addressing both time and frequency-domain formulations. In the time domain, a shifted Legendre polynomial expansion is employed to resolve the coupling between delay states and non-zero initial history functions. The projection matrix is derived from coefficients obtained via a Sylvester-type matrix equation, ensuring that the reduced model’s expansion coefficients match those of the original system. In the frequency domain, a unified MOR framework is proposed through generalized transfer function theory, integrating low-rank spatiotemporal decomposition of initial history functions. Moreover, a hybrid projection strategy combining Taylor and Laguerre series expansions is introduced to jointly capture transient and steady-state behaviors. Numerical validation underscores the accuracy and efficiency of both approaches.
This paper presents a novel parametric model order reduction (MOR) method specifically developed for discrete-time parameter-varying systems. The main contribution lies in extending non-parametric reduction techniques, previously established for continuous-time systems, to the discrete-time framework. By employing a bilinearization strategy, the parameter-varying system is transformed into an equivalent discrete-time bilinear form through the introduction of auxiliary input functions derived from the Taylor expansion of the parameter-dependent matrices. To efficiently approximate the Gramians of the resulting bilinear system, a Laguerre function-based low-rank approximation is proposed, enabling accurate and computationally efficient evaluation of the system’s dominant dynamics. The reduced-order model is then constructed using projection matrices obtained via singular value decomposition. The proposed approach effectively addresses non-affine parameter dependencies, providing a more precise and scalable reduction framework for large-scale systems. Numerical examples demonstrate the accuracy and computational efficiency of the proposed method in approximating complex discrete-time parameter-varying systems.
Dirichlet-Neumann and Neumann-Neumann methods are not only the parallel strategies in the spatial domain, but also can be used as a class of parallel methods in time. In this paper, we propose the Dirichlet-Neumann and Dirichlet-Neumann algorithms for time-periodic parabolic optimal control problems. By the Lagrange multiplier approach, a coupled system is obtained with a special time-periodic condition. For this coupled system, the Dirichlet-Neumann and Neumann-Dirichlet algorithms and their three variants are derived. We present the convergence analysis for all proposed algorithms. The numerical performance of the convergence factors is shown to illustrate our theoretical analysis. Based on our analysis, there is a class of algorithms with the better convergence compared with the natural Dirichlet-Neumann algorithm. Finally, numerical experiments are provided to illustrate the theoretical results.
We present and analyze in this paper a new space-time parallel method for solving evolution equations, i.e., the Parareal optimized Schwarz waveform relaxation (POSWR) algorithm. Since the classical Dirichlet transmission conditions inhibit the information exchange between subdomains to slow down the convergence speed of the Parareal Schwarz waveform relaxation (PSWR) algorithm, we introduce a class of optimized transmission conditions of Robin type to improve the convergence performance. We provide a convergence factor estimate based on the Laplace transform when the POSWR algorithm both with and without overlap is applied to the representative one-dimensional heat equation. Furthermore, we also analyze the optimized choice for the free parameter in the Robin transmission conditions to optimize the convergence behavior of the POSWR algorithm, which is different for the overlapping and nonoverlapping cases. Finally, we illustrate our theoretical analysis with several numerical experiments. We show that the new optimized algorithm even if without overlap converges much faster than the classical one, and the POSWR algorithm has different convergence behaviors in the overlapping and nonoverlapping cases.
Discrete time-delay systems arise widely in engineering and scientific applications. These systems are often modelled as high-dimensional dynamical systems, and the presence of delays increases the computational complexity and cost of both simulation and analysis. To alleviate this challenge, this paper proposes a new model reduction method in a parallel manner tailored for discrete time-delay systems. The proposed method begins with the explicit difference relations of Krawtchouk polynomials and a structural analysis of the shift-transformation matrices. We expand the original system over a Krawtchouk polynomial basis, which yields a system of linear algebraic equations characterised by a block lower triangular Toeplitz matrix. Further, applying the block discrete Fourier transform to the involved block alpha-circulant matrices, we design a parallel strategy to efficiently construct the projection basis used for reducing discrete time-delay systems. This is the main contribution of this paper. We provide rigorous theoretical results on the invertibility of the block alpha-circulant matrices and present error estimation bounds to ensure the validity and stability of the reduced model. Numerical examples demonstrate that our method achieves accurate approximation while offering substantial computational speed-ups, especially for large-scale time-delay systems.
This paper studies the data-driven balanced truncation (BT) method for second-order systems based on the measurements in the frequency domain. The basic idea is to approximate Gramians via the numerical quadrature rules, and establish the relationship between the main quantities in the procedure of BT and the sample data, which paves the way for the execution of BT in a nonintrusive manner. We construct the structure-preserving reduced models approximately based on the sample data of second-order systems with proportional damping, and provide a detailed algorithm in real-valued arithmetic to establish the data-driven counterpart of BT. In order to address the issue of large amount of sample data, we exploit the fact that the main quantities satisfy a couple of Sylvester matrix equations. The low-rank approximation to the solution of Sylvester equations is employed to avoid the explicit calculation of the main quantities, leading to an acceleration of the process of the data-driven BT. The performance of our approach is illustrated in detail via two numerical examples.
In this work, the penalty method is studied for the mixed Stokes-Darcy problem, motivated by the penalty method applied to Stokes equation. This work first proposes the penalty Stokes- Darcy model at the continuous level. Then we prove that the solution of the penalty model converges strongly to the original solution as O(epsilon) in which the penalty parameter is epsilon -> 0 . What is more, the finite element method is used to solve the penalty model and the optimal error estimates are presented. Finally, several numerical tests are carried out to verify our theoretical results.
This paper considers the H2 optimal model reduction problem of linear dynamical systems with quadratic output on the Riemannian manifolds. A one-sided projection is used to reduce the state equation, while a suitable symmetric matrix is chosen to determine the output equation of the reduced system. Since the projection matrix is an orthonormal matrix, it can be seen as a point on the Stiefel manifold. Because symmetric matrices of the same dimension allow a manifold structure, it is used to define a product manifold combined with the Stiefel manifold. The H2 error between the original system and the reduced system is treated as a function defined on the product manifold. Then, the H2 optimal model reduction problem is formulated as an unconstrained optimization problem defined on the product manifold. Concerning the symmetric matrix, the H2 error is proved to be convex. In terms of the orthonormal matrix and the symmetric matrix, the gradients of the H2 error are derived respectively. Then, the Riemannian BFGS method is used to obtain the orthonormal matrix, and the symmetric matrix is calculated by the convexity and the related gradient. By introducing the Riemannian manifolds to the H2 optimal model reduction problem, the constrained optimization problem in the Euclidean space is transformed into an unconstrained optimization problem on the manifolds, and the gradients of the H2 error are equipped with relatively concise formulas. Finally, numerical results illustrate the performance of the proposed model reduction method.
Research on nonlinear model order reduction has revealed that as nonlinearity increases, the subspaces capturing dominant information require more complex bases. The complexity is influenced by two main factors: coefficients and approximation criteria. On one hand, it is affected by the characteristics of all coefficients. Therefore, we begin by introducing a generalized Gramian-based method for estimating eigenvalue decay, which demonstrates the factors on the reduced order. On the other hand, existing interpolation conditions based on transfer functions must match multiple features with different bases. This article confirms that kernels, as new criteria, can match all features of a large class of nonlinear systems using the same basis as the linear part. We propose a $\mathit{k}$-dimensional refined space for affine input/output nonlinear systems, in contrast to existing methods that require at least an $\mathit{O(nk)}$-dimensional space to match $\mathit{n}$ transfer functions at $\mathit{k}$ interpolating points, even for $\mathbf{0}$th moment matching. To expand the applicability, we present a rank-$\mathbf{R}$ quadratic approximation to transform general systems into normalized affine input/output nonlinear systems. In terms of computational efficiency, we propose a parallel partial columnwise least-squares method to further reduce the rank. Finally, we provide two numerical examples to illustrate the effectiveness of our method.