Energy-preserving methods for Poisson systems

We present and analyze energy-conserving methods for the numerical integration of IVPs of Poisson type that are able to preserve some Casimirs. Their derivation and analysis is done following the ideas of Hamiltonian BVMs (HBVMs) (see Brugnano et al. [10] and references therein). It is seen that the proposed approach allows us to obtain the methods recently derived in Cohen and Hairer (2011) [17], giving an alternative derivation of such methods and a new proof of their order. Sufficient conditions that ensure the existence of a unique solution of the implicit equations defining the formulae are given. A study of the implementation of the methods is provided. In particular, order and preservation properties when the involved integrals are approximated by means of a quadrature formula, are derived.     

Extrapolation with spline-collocation methods for two-point boundary-value problems I: Proposals and justifications

Asymptotic expansions for the error in some spline interpolation schemes are used to derive asymptotic expansions for the truncation errors in some spline-collocation methods for two-point boundary-value problems. This raises the possibility of using Richardson extrapolation or iterated deferred corrections to develop efficient high-order algorithms based on low-order collocation in analogy with similar codes based on low-order finite difference methods; some specific such procedures are proposed.This research was supported in part by the United States Office of Naval Research under Contract N00014-67-A-0126-0015.     

Multistep Methods for Conservative Problems

We discuss the use of linear multistep methods for the solution of conservative (in particular Hamiltonian) problems. Despite the lack of good results concerning their behaviour, linear multistep methods have been satisfactorily used by researchers in applicative fields. We try to elucidate the requirements that a linear multistep method should fulfill in order to be suitable for the integration of such problems, relating our analysis to their use as boundary value methods. We collect a number of results obtained for the linear autonomous case and also consider some numerical tests on the Kepler problem.Work supported by GNCS and MIUR.     

A note on error expressions for reflected and averaged implicit Runge-Kutta methods

In addition to their usefulness in the numerical solution of initial value ODE's, the implicit Runge-Kutta (IRK) methods are also important for the solution of two-point boundary value problems. Recently, several classes of modified IRK methods which improve significantly on the efficiency of the standard IRK methods in this application have been presented. One such class is the Averaged IRK methods; a member of the class is obtained by applying an averaging operation to a non-symmetric IRK method and its reflection. In this paper we investigate the forms of the error expressions for reflected and averaged IRK methods. Our first result relates the expression for the local error of the reflected method to that of the original method. The main result of this paper relates the error expression of an averaged method to that of the method upon which it is based. We apply these results to show that for each member of the class of the averaged methods, there exists an embedded lower order method which can be used for error estimation, in a formula-pair fashion.This work was supported by the Natural Science and Engineering Research Council of Canada.     

NP-stability of Runge-Kutta methods based on classical quadrature

Recently Bellen, Jackiewicz and Zennaro have studied stability properties of Runge-Kutta (RK) methods for neutral delay differential equations using a scalar test equation. In particular, they have shown that everyA-stable collocation method isNP-stable, i.e., the method has an analogous stability property toA-stability with respect to the test equation. Consequently, the Gauss, Radau IIA and Lobatto IIIA methods areNP-stable. In this paper, we examine the stability of RK methods based on classical quadrature by a slightly different approach from theirs. As a result, we prove that the Radau IA and Lobatto IIIC methods equipped with suitable continuous extensions are alsoNP-stable by virtue of fundamental notions related to those methods such as simplifying conditions, algebraic stability, and theW-transformation.     

Runge–Kutta convolution quadrature methods for well-posed equations with memory

Runge–Kutta based convolution quadrature methods for abstract, well-posed, linear, and homogeneous Volterra equations, non necessarily of sectorial type, are developed. A general representation of the numerical solution in terms of the continuous one is given. The error and stability analysis is based on this representation, which, for the particular case of the backward Euler method, also shows that the numerical solution inherits some interesting qualitative properties, such as positivity, of the exact solution. Numerical illustrations are provided.     

An investigation of the internal structure of shock profiles for shock capturing schemes

The theoretical understanding of discrete shock transitions obtained by shock capturing schemes is very incomplete. Previous experimental studies indicate that discrete shock transitions obtained by shock capturing schemes can be modeled by continuous functions, so called continuum shock profiles. However, the previous papers have focused on linear methods. We have experimentally studied the trajectories of discrete shock profiles in phase space for a range of different high resolution shock capturing schemes, including Riemann solver based flux limiter methods, high resolution central schemes and ENO type methods. In some cases, no continuum profiles exists. However, in these cases the point values in the shock transitions remain bounded and appear to converge toward a stable limit cycle. The possibility of such behavior was anticipated in Bultelle, Grassin and Serre, 1998, but no specific examples, or other evidence, of this behavior have previously been given. In other cases, our results indicate that continuum shock profiles exist, but are very complicated. We also study phase space orbits with regard to post shock oscillations.     

Numerical solution of stochastic differential problems in the biosciences

Stochastic differential equations (SDEs) models play a prominent role in many application areas including biology, epidemiology and population dynamics, mostly because they can offer a more sophisticated insight through physical phenomena than their deterministic counterparts do. So, suitable numerical methods must be introduced to simulate the solutions of the resulting stochastic differential systems. In this work we take into account both Euler–Taylor expansion and Runge–Kutta-type methods for stochastic ordinary differential equations (SODEs) and the Euler–Maruyama method for stochastic delay differential equations (SDDEs), focusing on the most relevant implementation issues. The corresponding Matlab codes for both SODEs and SDDEs problems are tested on mathematical models arising in the biosciences.     

Structure preservation of exponentially fitted Runge–Kutta methods

The preservation of some structure properties of the flow of differential systems by numerical exponentially fitted Runge–Kutta (EFRK) methods is considered. A complete characterisation of EFRK methods that preserve linear or quadratic invariants is given and, following the approach of Bochev and Scovel [On quadratic invariants and symplectic structure, BIT 34 (1994) 337–345], the sufficient conditions on symplecticity of EFRK methods derived by Van de Vyver [A fourth-order symplectic exponentially fitted integrator, Comput. Phys. Comm. 174 (2006) 255–262] are obtained. Further, a family of symplectic EFRK two-stage methods with order four has been derived. It includes the symplectic EFRK method proposed by Van de Vyver as well as a collocation method at variable nodes that can be considered as the natural collocation extension of the classical RK Gauss method. Finally, the results of some numerical experiments are presented to compare the relative merits of several fitted and nonfitted fourth-order methods in the integration of oscillatory systems.     

Construction of a multirate RODAS method for stiff ODEs

Multirate time stepping is a numerical technique for efficiently solving large-scale ordinary differential equations (ODEs) with widely different time scales localized over the components. This technique enables one to use large time steps for slowly varying components, and small steps for rapidly varying ones. Multirate methods found in the literature are normally of low order, one or two. Focusing on stiff ODEs, in this paper we discuss the construction of a multirate method based on the fourth-order RODAS method. Special attention is paid to the treatment of the refinement interfaces with regard to the choice of the interpolant and the occurrence of order reduction. For stiff, linear systems containing a stiff source term, we propose modifications for the treatment of the source term which overcome order reduction originating from such terms and which we can implement in our multirate method.     

Stability of Runge-Kutta methods for the generalized pantograph equation

Summary. This paper deals with stability properties of Runge-Kutta (RK) methods applied to a non-autonomous delay differential equation (DDE) with a constant delay which is obtained from the so-called generalized pantograph equation, an autonomous DDE with a variable delay by a change of the independent variable. It is shown that in the case where the RK matrix is regular stability properties of the RK method for the DDE are derived from those for a difference equation, which are examined by similar techniques to those in the case of autonomous DDEs with a constant delay. As a result, it is shown that some RK methods based on classical quadrature have a superior stability property with respect to the generalized pantograph equation. Stability of algebraically stable natural RK methods is also considered. Received May 5, 1998 / Revised version received November 17, 1998 / Published online September 24, 1999     

A criterion forP-stability properties of Runge-Kutta methods

P-stability is an analogous stability property toA-stability with respect to delay differential equations. It is defined by using a scalar test equation similar to the usual test equation ofA-stability. EveryP-stable method isA-stable, but anA-stable method is not necessarilyP-stable. We considerP-stability of Runge-Kutta (RK) methods and its variation which was originally introduced for multistep methods by Bickart, and derive a sufficient condition for an RK method to have the stability properties on the basis of an algebraic characterization ofA-stable RK methods recently obtained by Schere and Müller. By making use of the condition we clarify stability properties of some SIRK and SDIRK methods, which are easier to implement than fully implicit methods, applied to delay differential equations.     

Twostep-by-twostep PIRK-type PC methods with continuous output formulas

This paper deals with parallel predictor–corrector (PC) iteration methods based on collocation Runge–Kutta (RK) corrector methods with continuous output formulas for solving nonstiff initial-value problems (IVPs) for systems of first-order differential equations. At n$n$th step, the continuous output formulas are used not only for predicting the stage values in the PC iteration methods but also for calculating the step values at (n+2)$(n+2)$th step. In this case, the integration processes can be proceeded twostep-by-twostep. The resulting twostep-by-twostep (TBT) parallel-iterated RK-type (PIRK-type) methods with continuous output formulas (twostep-by-twostep PIRKC methods or TBTPIRKC methods) give us a faster integration process. Fixed stepsize applications of these TBTPIRKC methods to a few widely-used test problems reveal that the new PC methods are much more efficient when compared with the well-known parallel-iterated RK methods (PIRK methods), parallel-iterated RK-type PC methods with continuous output formulas (PIRKC methods) and sequential explicit RK codes DOPRI5 and DOP853 available from the literature.     

Doubly quasi-consistent parallel explicit peer methods with built-in global error estimation

Recently, Kulikov presented the idea of double quasi-consistency, which facilitates global error estimation and control, considerably. More precisely, a local error control implemented in such methods plays a part of global error control at the same time. However, Kulikov studied only Nordsieck formulas and proved that there exists no doubly quasi-consistent scheme among those methods.Here, we prove that the class of doubly quasi-consistent formulas is not empty and present the first example of such sort. This scheme belongs to the family of superconvergent explicit two-step peer methods constructed by Weiner, Schmitt, Podhaisky and Jebens. We present a sample of s-stage doubly quasi-consistent parallel explicit peer methods of order s−1 when s=3. The notion of embedded formulas is utilized to evaluate efficiently the local error of the constructed doubly quasi-consistent peer method and, hence, its global error at the same time. Numerical examples of this paper confirm clearly that the usual local error control implemented in doubly quasi-consistent numerical integration techniques is capable of producing numerical solutions for user-supplied accuracy conditions in automatic mode.     

Beyond conventional Runge–Kutta methods in numerical integration of ODEs and DAEs by use of structures and local models

There are two parts in this paper. In the first part we consider an overdetermined system of differential-algebraic equations (DAEs). We are particularly concerned with Hamiltonian and Lagrangian systems with holonomic constraints. The main motivation is in finding methods based on Gauss coefficients, preserving not only the constraints, symmetry, symplecticness, and variational nature of trajectories of holonomically constrained Hamiltonian and Lagrangian systems, but also having optimal order of convergence. The new class of (s,s)$(s,s)$-Gauss–Lobatto specialized partitioned additive Runge–Kutta (SPARK) methods uses greatly the structure of the DAEs and possesses all desired properties. In the second part we propose a unified approach for the solution of ordinary differential equations (ODEs) mixing analytical solutions and numerical approximations. The basic idea is to consider local models which can be solved efficiently, for example analytically, and to incorporate their solution into a global procedure based on standard numerical integration methods for the correction. In order to preserve also symmetry we define the new class of symmetrized Runge–Kutta methods with local model (SRKLM).     

A stability property ofA-stable natural Runge-Kutta methods for systems of delay differential equations

A natural Runge-Kutta method is a special type of Runge-Kutta method for delay differential equations (DDEs); it is known that any one-step collocation method is equivalent to one of such methods. In this paper, we consider a linear constant-coefficient system of DDEs with a constant delay, and discuss the application of natural Runge-Kutta methods to the system. We show that anA-stable method preserves the asymptotic stability property of the analytical solutions of the system.     

A robust trigonometrically fitted embedded pair for perturbed oscillators

A new kind of trigonometrically fitted embedded pair of explicit ARKN methods for the numerical integration of perturbed oscillators is presented in this paper. This new pair is based on the trigonometrically fitted ARKN method of order five derived by Yang and Wu in [H.L. Yang, X.Y. Wu, Trigonometrically-fitted ARKN methods for perturbed oscillators, Appl. Numer. Math. 9 (2008) 1375–1395]. We analyze the stability properties, phase-lag (dispersion) and dissipation of the higher-order method of the new pair. Numerical experiments carried out show that our new embedded pair is very competitive in comparison with the embedded pairs proposed in the scientific literature.     

Probabilistic and deterministic convergence proofs for software for initial value problems

The numerical solution of initial value problems for ordinary differential equations is frequently performed by means of adaptive algorithms with user-input tolerance τ. The time-step is then chosen according to an estimate, based on small time-step heuristics, designed to try and ensure that an approximation to the local error commited is bounded by τ. A question of natural interest is to determine how the global error behaves with respect to the tolerance τ. This has obvious practical interest and also leads to an interesting problem in mathematical analysis. The primary difficulties arising in the analysis are that: (i) the time-step selection mechanisms used in practice are discontinuous as functions of the specified data; (ii) the small time-step heuristics underlying the control of the local error can break down in some cases. In this paper an analysis is presented which incorporates these two difficulties. For a mathematical model of an error per unit step or error per step adaptive Runge–Kutta algorithm, it may be shown that in a certain probabilistic sense, with respect to a measure on the space of initial data, the small time-step heuristics are valid with probability one, leading to a probabilistic convergence result for the global error as τ→0. The probabilistic approach is only valid in dimension m>1 this observation is consistent with recent analysis concerning the existence of spurious steady solutions of software codes which highlights the difference between the cases m=1 and m>1. The breakdown of the small time-step heuristics can be circumvented by making minor modifications to the algorithm, leading to a deterministic convergence proof for the global error of such algorithms as τ→0. An underlying theory is developed and the deterministic and probabilistic convergence results proved as particular applications of this theory. This revised version was published online in June 2006 with corrections to the Cover Date.     

Global error estimation and control in linearly-implicit parallel two-step peer W-methods

The class of linearly-implicit parallel two-step peer W-methods has been designed recently for efficient numerical solutions of stiff ordinary differential equations. Those schemes allow for parallelism across the method, that is an important feature for implementation on modern computational devices. Most importantly, all stage values of those methods possess the same properties in terms of stability and accuracy of numerical integration. This property results in the fact that no order reduction occurs when they are applied to very stiff problems. In this paper, we develop parallel local and global error estimation schemes that allow the numerical solution to be computed for a user-supplied accuracy requirement in automatic mode. An algorithm of such global error control and other technical particulars are also discussed here. Numerical examples confirm efficiency of the presented error estimation and stepsize control algorithm on a number of test problems with known exact solutions, including nonstiff, stiff, very stiff and large-scale differential equations. A comparison with the well-known stiff solver RODAS is also shown.