2025-04-02 | | Total: 21
In this work, we propose an exponentially convergent numerical method for the Caputo fractional propagator Sα(t) and the associated mild solution of the Cauchy problem with time-independent sectorial operator coefficient A and Caputo fractional derivative of order α∈(0,2) in time. The proposed methods are constructed by generalizing the earlier developed approximation of Sα(t) with help of the subordination principle. Such technique permits us to eliminate the dependence of the main part of error estimate on α, while preserving other computationally relevant properties of the original approximation: native support for multilevel parallelism, the ability to handle initial data with minimal spatial smoothness, and stable exponential convergence for all t∈[0,T]. Ultimately, the use of subordination leads to a significant improvement of the method's convergence behavior, particularly for small α<0.5, and opens up further opportunities for efficient data reuse. To validate theoretical results, we consider applications of the developed methods to the direct problem of solution approximation, as well as to the inverse problem of fractional order identification.
In this work, we address the efficient computation of parameterized systems of linear equations, with possible nonlinear parameter dependence. When the matrix is highly sensitive to the parameters, mean-based preconditioning might not be enough. For this scenario, we explore an approach in which several preconditioners are placed in the parameter space during a precomputation step. To determine the optimal placement of a limited number of preconditioners, we estimate the expected number of iterations with respect to a given preconditioner a priori and use a location-allocation strategy to optimize the placement of the preconditioners. We elaborate on our methodology for the Helmholtz problem with exterior Dirichlet scattering at high frequencies, and we estimate the expected number of GMRES iterations via a gray-box Gaussian process regression approach. We illustrate our approach in two practical applications: scattering in a domain with a parametric refractive index and scattering from a scatterer with parameterized shape. Using these numerical examples, we show how our methods leads to runtime savings of about an order of magnitude. Moreover, we investigate the effect of the parameter dimension and the importance of dimension anisotropy on their efficacy.
In this paper, we compare the performance of different numerical schemes in approximating Pontryagin's Maximum Principle's necessary conditions for the optimal control of nonholonomic systems. Retraction maps are used as a seed to construct geometric integrators for the corresponding Hamilton equations. First, we obtain an intrinsic formulation of a discretization map on a distribution D. Then, we illustrate this construction on a particular example for which the performance of different symplectic integrators is examined and compared with that of non-symplectic integrators.
The Eberlein diagonalization method is an iterative Jacobi-type method for solving the eigenvalue problem of a general complex matrix. In this paper we develop the block version of the Eberlein method. We prove the global convergence of our block method and present several numerical examples.
We propose a fast fourth-order cut cell method for solving constant-coefficient elliptic equations in two-dimensional irregular domains. In our methodology, the key to dealing with irregular domains is the poised lattice generation (PLG) algorithm that generates finite-volume interpolation stencils near the irregular boundary. We are able to derive high-order discretization of the elliptic operators by least squares fitting over the generated stencils. We then design a new geometric multigrid scheme to efficiently solve the resulting linear system. Finally, we demonstrate the accuracy and efficiency of our method through various numerical tests in irregular domains.
The two-phase Navier-Stokes Cahn-Hilliard (NSCH) mixture model is a key framework for simulating multiphase flows with non-matching densities. Developing fully discrete, energy-stable schemes for this model remains challenging, due to the possible presence of negative densities. While various methods have been proposed, ensuring provable energy stability under phase-field modifications, like positive extensions of the density, remains an open problem. We propose a simple, fully discrete, energy-stable method for the NSCH mixture model that ensures stability with respect to the energy functional, where the density in the kinetic energy is positively extended. The method is based on an alternative but equivalent formulation using mass-averaged velocity and volume-fraction-based order parameters, simplifying implementation while preserving theoretical consistency. Numerical results demonstrate that the proposed scheme is robust, accurate, and stable for large density ratios, addressing key challenges in the discretization of NSCH models.
The purpose of this note is to investigate the coupling of Dirichlet and Neumann numerical boundary conditions for the transport equation set on an interval. When one starts with a stable finite difference scheme on the lattice Z and each numerical boundary condition is taken separately with the Neumann extrapolation condition at the outflow boundary, the corresponding numerical semigroup on a half-line is known to be bounded. It is also known that the coupling of such numerical boundary conditions on a compact interval yields a stable approximation, even though large time exponentially growing modes may occur. We review the different stability estimates associated with these numerical boundary conditions and give explicit examples of such exponential growth phenomena for finite difference schemes with ''small'' stencils. This provides numerical evidence for the optimality of some stability estimates on the interval.
We construct a new compact semi-explicit three-level in time fourth-order finite-difference scheme for numerical solving the general multidimensional acoustic wave equation, where both the speed of sound and density of a medium are variable. The scheme is three-point in each spatial direction, has the truncation order O(|h|4+h4t) and is easily implementable. It seems to be the first compact scheme with such properties for the equation under consideration. It generalizes a semi-explicit compact scheme developed and studied recently in the much simpler case of the variable speed of sound only. Numerical experiments confirm the high precision of the scheme and its fourth error order not only in the mesh C norm but in the mesh C1 norm as well.
We develop and analyse residual-based a posteriori error estimates for the virtual element discretisation of a nonlinear stress-assisted diffusion problem in two and three dimensions. The model problem involves a two-way coupling between elasticity and diffusion equations in perturbed saddle-point form. A robust global inf-sup condition and Helmholtz decomposition for H(div,Ω) lead to a reliable and efficient error estimator based on appropriately weighted norms that ensure parameter robustness. The a posteriori error analysis uses quasi-interpolation operators for Stokes and edge virtual element spaces, and we include the proofs of such operators with estimates in 3D for completeness. Finally, we present numerical experiments in both 2D and 3D to demonstrate the optimal performance of the proposed error estimator.
This paper proposes an adaptive hyper-reduction method to reduce the computational cost associated with the simulation of parametric particle-based kinetic plasma models, specifically focusing on the Vlasov-Poisson equation. Conventional model order reduction and hyper-reduction techniques are often ineffective for such models due to the non-sparse nature of the nonlinear operators arising from the interactions between particles. To tackle this issue, we propose an adaptive, structure-preserving hyper-reduction method that leverages a decomposition of the discrete reduced Hamiltonian into a linear combination of terms, each depending on a few components of the state. The proposed approximation strategy allows to: (i) preserve the Hamiltonian structure of the problem; (ii) evaluate nonlinear non-sparse operators in a computationally efficient way; (iii) overcome the Kolmogorov barrier of transport-dominated problems via evolution of the approximation space and adaptivity of the rank of the solution. The proposed method is validated on numerical benchmark simulations, demonstrating stable and accurate performance with substantial runtime reductions compared to the full order model.
We present a variational integrator based on the Lobatto quadrature for the time integration of dynamical systems issued from the least action principle. This numerical method uses a cubic interpolation of the states and the action is approximated at each time step by Lobatto's formula. Numerical analysis is performed on a harmonic oscillator. The scheme is conditionally stable, sixth-order accurate, and symplectic. It preserves an approximate energy quantity. Simulation results illustrate the performance of the proposed method.
When solving partial differential equations using classical schemes such as finite difference or finite volume methods, sufficiently fine meshes and carefully designed schemes are required to achieve high-order accuracy of numerical solutions, leading to a significant increase in computational costs, especially for three-dimensional (3D) time-dependent problems. Recently, machine learning-assisted numerical methods have been proposed to enhance accuracy or efficiency. In this paper, we propose a data-driven finite difference numerical method to solve the hyperbolic equations with smooth solutions on coarse grids, which can achieve higher accuracy than classical numerical schemes based on the same mesh size. In addition, the data-driven schemes have better spectrum properties than the classical schemes, although the spectrum properties are not explicitly optimized during the training process. Numerical examples are presented to demonstrate the accuracy and efficiency of the proposed method, as well as its good performance on dispersion and dissipation.
Pre-filtering and post-filtering steps can be added to many of the traditional numerical methods to generate new, higher order methods with strong stability properties. Presented in this paper are a variable step pre-filter and post-filter that allow adaptive time stepping for a filtered method based on Implicit Euler.
The linear Boltzmann model for proton beams is a six-dimensional partial differential equation (PDE). We propose a deterministic solver for the linear Boltzmann model based on scattering decomposition and depth-splitting methods. The main idea is to first divide the protons into primary protons and scattering protons, whose equations are derived using the source iteration method. We then treat depth as the time variable in classical time-evolutionary problems and apply the depth-splitting method. In the depth-splitting method, the full operator is decomposed into three parts, with each subsystem being easily parallelizable, which is crucial for efficient simulations. The resulting discretization exhibits second-order convergence in both the depth and energy variables. The dose distributions obtained from our solver are compared with those from Monte Carlo simulations for various materials and heterogeneous cases.
This paper is devoted to examining the stability of Runge-Kutta methods for solving stiff nonlinear Volterra delay-integro-differential-algebraic equations (DIDAEs) with constant delay. Hybrid numerical schemes combining Runge-Kutta methods and compound quadrature rules are analyzed for nonlinear stiff DIDAEs. Criteria for ensuring the global and asymptotic stability of the proposed schemes are established. Several numerical examples are provided to validate the theoretical findings.
A cross matrix X can have nonzero elements located only on the main diagonal and the anti-diagonal, so that the sparsity pattern has the shape of a cross. It is shown that X can be factorized into products of matrices that are at most rank-two perturbations to the identity matrix and can be symmetrically permuted to block diagonal form with 2×2 diagonal blocks and, if n is odd, a 1×1 diagonal block. The permutation similarity implies that any well-defined analytic function of X remains a cross matrix. By exploiting these properties, explicit formulae for the determinant, inverse, and characteristic polynomial are derived. It is also shown that the structure of cross matrix can be preserved under matrix factorizations, including the LU, QR, and SVD decompositions.
A new Chebyshev-type family of stabilized explicit methods for solving mildly stiff ODEs is presented. Besides conventional conditions of order and stability we impose an additional restriction on the methods: their stability function must be monotonically increasing and positive along the largest possible interval of negative real axis. Although stability intervals of the proposed methods are smaller than those of classic Chebyshev-type methods, their stability functions are more consistent with the exponent, they have more convex stability regions and smaller error constants. These properties allow the monotonic methods to be competitive with contemporary stabilized second-order methods, as the presented results of numerical experiments demonstrate.
We consider Fredholm integral equation of the first kind, present an efficient new iterated Tikhonov method to solve it. The new Tikhonov iteration method has been proved which can achieve the optimal order under a-priori assumption. In numerical experiments, the new iterated Tikhonov regularization method is compared with the classical iterated Tikhonov method, Landweber iteration method to solve the corresponding discrete problem, which indicates the validity and efficiency of the proposed method.
Gaussian processes (GPs) are crucial in machine learning for quantifying uncertainty in predictions. However, their associated covariance matrices, defined by kernel functions, are typically dense and large-scale, posing significant computational challenges. This paper introduces a matrix-free method that utilizes the Non-equispaced Fast Fourier Transform (NFFT) to achieve nearly linear complexity in the multiplication of kernel matrices and their derivatives with vectors for a predetermined accuracy level. To address high-dimensional problems, we propose an additive kernel approach. Each sub-kernel in this approach captures lower-order feature interactions, allowing for the efficient application of the NFFT method and potentially increasing accuracy across various real-world datasets. Additionally, we implement a preconditioning strategy that accelerates hyperparameter tuning, further improving the efficiency and effectiveness of GPs.
In this paper, we study numerical approximations for optimal control of a class of stochastic partial differential equations with partial observations. The system state evolves in a Hilbert space, whereas observations are given in finite-dimensional space \rrd. We begin by establishing stochastic maximum principles (SMPs) for such problems, where the system state is driven by a cylindrical Wiener process. The corresponding adjoint equations are characterized by backward stochastic partial differential equations. We then develop numerical algorithms to solve the partially observed optimal control. Our approach combines the stochastic gradient descent method, guided by the SMP, with a particle filtering algorithm to estimate the conditional distributions of the state of the system. Finally, we demonstrate the effectiveness of our proposed algorithm through numerical experiments.
The canonical polyadic tensor decomposition has a long history. But it becomes challenging to find a tensor decomposition when the rank is between the largest and the second-largest dimension. In such cases, traditional optimization methods, such as nonlinear least squares or alternative least squares methods, often fail to find a tensor decomposition. There are also direct methods, such as the normal form algorithm and the method by Domanov and De Lathauwer, that solve tensor decompositions algebraically. However, these methods can be computationally expensive and require significant memory, especially when the rank is high. This paper proposes a novel two-stage algorithm for the order-3 nonsymmetric tensor decomposition problem when the rank is not greater than the largest dimension. It transforms the tensor decomposition problem into two optimization problems. When the first-stage optimization is not fully solved, the partial solution will also be leveraged in the second-stage optimization problem. We prove the equivalence between tensor decompositions and the global minimizers of the two-stage optimization problems. Our numerical experiments demonstrate the proposed two-stage optimization method is very efficient and robust, capable of finding tensor decompositions where other commonly used state-of-the-art methods fail.