Definite integrals with parameters of holonomic functions satisfy holonomic systems of linear partial differential equations. When we restrict parameters to a one dimensional curve, the system becomes a linear ordinary differential equation (ODE) with respect to a curve in the parameter space. We can evaluate the integral by solving the linear ODE numerically. This approach to evaluate numerically definite integrals is called the holonomic gradient method (HGM) and it is useful to evaluate several normalizing constants in statistics. We will discuss and compare methods to solve linear ODE's to evaluate normalizing constants.
In this paper, we propose a variationally consistent technique for decreasing the maximum eigenfrequencies of structural dynamics related finite element formulations. Our approach is based on adding a symmetric positive-definite term to the mass matrix that follows from the integral of the traction jump across element boundaries. The added term is weighted by a small factor, for which we derive a suitable, and simple, element-local parameter choice. For linear problems, we show that our mass-scaling method produces no adverse effects in terms of spatial accuracy and orders of convergence. We illustrate these properties in one, two and three spatial dimension, for quadrilateral elements and triangular elements, and for up to fourth order polynomials basis functions. To extend the method to non-linear problems, we introduce a linear approximation and show that a sizeable increase in critical time-step size can be achieved while only causing minor (even beneficial) influences on the dynamic response.
The asymptotic stable region and long-time decay rate of solutions to linear homogeneous Caputo time fractional ordinary differential equations (F-ODEs) are known to be completely determined by the eigenvalues of the coefficient matrix. Very different from the exponential decay of solutions to classical ODEs, solutions of F-ODEs decay only polynomially, leading to the so-called Mittag-Leffler stability, which was already extended to semi-linear F-ODEs with small perturbations. This work is mainly devoted to the qualitative analysis of the long-time behavior of numerical solutions. By applying the singularity analysis of generating functions developed by Flajolet and Odlyzko (SIAM J. Disc. Math. 3 (1990), 216-240), we are able to prove that both $\mathcal{L}$1 scheme and strong $A$-stable fractional linear multistep methods (F-LMMs) can preserve the numerical Mittag-Leffler stability for linear homogeneous F-ODEs exactly as in the continuous case. Through an improved estimate of the discrete fractional resolvent operator, we show that strong $A$-stable F-LMMs are also Mittag-Leffler stable for semi-linear F-ODEs under small perturbations. For the numerical schemes based on $\alpha$-difference approximation to Caputo derivative, we establish the Mittag-Leffler stability for semi-linear problems by making use of properties of the Poisson transformation and the decay rate of the continuous fractional resolvent operator. Numerical experiments are presented for several typical time fractional evolutional equations, including time fractional sub-diffusion equations, fractional linear system and semi-linear F-ODEs. All the numerical results exhibit the typical long-time polynomial decay rate, which is fully consistent with our theoretical predictions.
Numerical solving differential equations with fractional derivatives requires elimination of the singularity which is inherent in the standard definition of fractional derivatives. The method of integration by parts to eliminate this singularity is well known. It allows to solve some equations but increases the order of the equation and sometimes leads to wrong numerical results or instability. We suggest another approach: the elimination of singularity by substitution. It does not increase the order of equation and its numerical implementation provides the opportunity to define fractional derivative as the limit of discretization. We present a sufficient condition for the substitution-generated difference approximation to be well-conditioned. We demonstrate how some equations can be solved using this method with full confidence that the solution is accurate with at least second order of approximation.
We present a method for solving linear and nonlinear PDEs based on the variable projection (VarPro) framework and artificial neural networks (ANN). For linear PDEs, enforcing the boundary/initial value problem on the collocation points leads to a separable nonlinear least squares problem about the network coefficients. We reformulate this problem by the VarPro approach to eliminate the linear output-layer coefficients, leading to a reduced problem about the hidden-layer coefficients only. The reduced problem is solved first by the nonlinear least squares method to determine the hidden-layer coefficients, and then the output-layer coefficients are computed by the linear least squares method. For nonlinear PDEs, enforcing the boundary/initial value problem on the collocation points leads to a nonlinear least squares problem that is not separable, which precludes the VarPro strategy for such problems. To enable the VarPro approach for nonlinear PDEs, we first linearize the problem with a Newton iteration, using a particular form of linearization. The linearized system is solved by the VarPro framework together with ANNs. Upon convergence of the Newton iteration, the network coefficients provide the representation of the solution field to the original nonlinear problem. We present ample numerical examples with linear and nonlinear PDEs to demonstrate the performance of the method herein. For smooth field solutions, the errors of the current method decrease exponentially as the number of collocation points or the number of output-layer coefficients increases. We compare the current method with the ELM method from a previous work. Under identical conditions and network configurations, the current method exhibits an accuracy significantly superior to the ELM method.
We demonstrate the effectiveness of an adaptive explicit Euler method for the approximate solution of the Cox-Ingersoll-Ross model. This relies on a class of path-bounded timestepping strategies which work by reducing the stepsize as solutions approach a neighbourhood of zero. The method is hybrid in the sense that a convergent backstop method is invoked if the timestep becomes too small, or to prevent solutions from overshooting zero and becoming negative. Under parameter constraints that imply Feller's condition, we prove that such a scheme is strongly convergent, of order at least 1/2. Control of the strong error is important for multi-level Monte Carlo techniques. Under Feller's condition we also prove that the probability of ever needing the backstop method to prevent a negative value can be made arbitrarily small. Numerically, we compare this adaptive method to fixed step implicit and explicit schemes, and a novel semi-implicit adaptive variant. We observe that the adaptive approach leads to methods that are competitive in a domain that extends beyond Feller's condition, indicating suitability for the modelling of stochastic volatility in Heston-type asset models.
The paper provides a novel framework to study the accuracy and stability of numerical integration schemes when employed for the time domain simulation of power systems. A matrix pencil-based approach is adopted to evaluate the error between the dynamic modes of the power system and the modes of the approximated discrete-time system arising from the application of the numerical method. The proposed approach can provide meaningful insights on how different methods compare to each other when applied to a power system, while being general enough to be systematically utilized for, in principle, any numerical method. The framework is illustrated for a handful of well-known explicit and implicit methods, while simulation results are presented based on the WSCC 9-bus system, as well as on a 1, 479-bus dynamic model of the All-Island Irish Transmission System.
This study proposes a novel coupled-mode theory for two-dimensional exterior Helmholtz problems. The proposed approach is based on the separation of the entire space R2 into a fictitious disk and its exterior. The disk is allocated in such a way that it comprises all the inhomogeneity; therefore, the exterior supports cylindrical waves with a continuous spectrum. For the interior, we expand an unknown wave field using normal modes that satisfy some auxiliary boundary conditions on the surface of the disk. For the interior expansion, we propose combining the Neumann and Dirichlet normal modes. We show that the proposed expansion sacrifices L2 orthogonality but significantly improve the convergence. Finally, we present some numerical verifications of the proposed coupled-mode theory.
The recently developed technique of DOC kernels has been a great success in the stability and convergence analysis for BDF2 scheme with variable time steps. However, such an analysis technique seems not directly applicable to problems with initial singularity. In the numerical simulations of solutions with initial singularity, variable time-steps schemes like the graded mesh are always adopted to achieve the optimal convergence, whose first adjacent time-step ratio may become pretty large so that the acquired restriction is not satisfied. In this paper, we revisit the variable time-step implicit-explicit two-step backward differentiation formula (IMEX BDF2) scheme presented in [W. Wang, Y. Chen and H. Fang, \emph{SIAM J. Numer. Anal.}, 57 (2019), pp. 1289-1317] to compute the partial integro-differential equations (PIDEs) with initial singularity. We obtain the sharp error estimate under a mild restriction condition of adjacent time-step ratios $r_{k}: =\tau_{k}/\tau_{k-1} \; (k\geq 3) < r_{\max} = 4.8645 $ and a much mild requirement on the first ratio, i.e., $r_2>0$. This leads to the validation of our analysis of the variable time-step IMEX BDF2 scheme when the initial singularity is dealt by a simple strategy, i.e., the graded mesh $t_k=T(k/N)^{\gamma}$. In this situation, the convergence of order $\mathcal{O}(N^{-\min\{2,\gamma \alpha\}})$ is achieved with $N$ and $\alpha$ respectively representing the total mesh points and indicating the regularity of the exact solution. This is, the optical convergence will be achieved by taking $\gamma_{\text{opt}}=2/\alpha$. Numerical examples are provided to demonstrate our theoretical analysis.
The classical Smagorinsky model's solution is an approximation to a (resolved) mean velocity. Since it is an eddy viscosity model, it cannot represent a flow of energy from unresolved fluctuations to the (resolved) mean velocity. This model has recently been modified to incorporate this flow and still be well-posed. Herein we first develop some basic properties of the modified model. Next, we perform a complete numerical analysis of two algorithms for its approximation. They are tested and proven to be effective.
We introduce a new family of deep neural network models. Instead of specifying a discrete sequence of hidden layers, we parameterize the derivative of the hidden state using a neural network. The output of the network is computed using a black-box differential equation solver. These continuous-depth models have constant memory cost, adapt their evaluation strategy to each input, and can explicitly trade numerical precision for speed. We demonstrate these properties in continuous-depth residual networks and continuous-time latent variable models. We also construct continuous normalizing flows, a generative model that can train by maximum likelihood, without partitioning or ordering the data dimensions. For training, we show how to scalably backpropagate through any ODE solver, without access to its internal operations. This allows end-to-end training of ODEs within larger models.