In this paper, we present an abstract framework which describes algebraically the derivation of order conditions independently of the nature of differential equations considered or the type of integrators used to solve them. Our structure includes a Hopf algebra of functions, whose properties are used to answer several questions of prime interest in numerical analysis. In particular, we show that, under some mild assumptions, there exist integrators of arbitrarily high orders for arbitrary (modified) vector fields.

In long-time numerical integration of Hamiltonian systems, and especially in molecular dynamics simulation, it is important that the energy is well conserved. For symplectic integrators applied with sufficiently small step size, this is guaranteed by the existence of a modified Hamiltonian that is exactly conserved up to exponentially small terms. This article is concerned with the simplified Takahashi-Imada method, which is a modification of the Störmer-Verlet method that is as easy to implement but has improved accuracy. This integrator is symmetric and volume-preserving, but no longer symplectic. We study its long-time energy conservation and give theoretical arguments, supported by numerical experiments, which show the possibility of a drift in the energy (linear or like a random walk). With respect to energy conservation, this article provides empirical and theoretical data concerning the importance of using a symplectic integrator.

We show that while Runge-Kutta methods cannot preserve polynomial invariants in general, they can preserve polynomials that are the energy invariant of canonical Hamiltonian systems.

In this paper, we study the linear Schrödinger equation over the d-dimensional torus, with small values of the perturbing potential. We consider numerical approximations of the associated solutions obtained by a symplectic splitting method (to discretize the time variable) in combination with the Fast Fourier Transform algorithm (to discretize the space variable). In this fully discrete setting, we prove that the regularity of the initial datum is preserved over long times, i.e. times that are exponentially long with the time discretization parameter. We here refer to Gevrey regularity, and our estimates turn out to be uniform in the space discretization parameter. This paper extends [G. Dujardin and E. Faou, Numer. Math.97 (2004) 493–535], where a similar result has been obtained in the semi-discrete situation, i.e. when the mere time variable is discretized and space is kept a continuous variable.

The plane wave stability properties of the conservative schemes of Besse [SIAM J. Numer. Anal.42 (2004) 934–952] and Fei et al. [Appl. Math. Comput.71 (1995) 165–177] for the cubic Schrödinger equation are analysed. Although the two methods possess many of the same conservation properties, we show that their stability behaviour is very different. An energy preserving generalisation of the Fei method with improved stability is presented.

The integration to steady state of many initial value ODEs and PDEs using the forward Euler method can alternatively be considered as gradient descent for an associated minimization problem. Greedy algorithms such as steepest descent for determining the step size are as slow to reach steady state as is forward Euler integration with the best uniform step size. But other, much faster methods using bolder step size selection exist. Various alternatives are investigated from both theoretical and practical points of view. The steepest descent method is also known for the regularizing or smoothing effect that the first few steps have for certain inverse problems, amounting to a finite time regularization. We further investigate the retention of this property using the faster gradient descent variants in the context of two applications. When the combination of regularization and accuracy demands more than a dozen or so steepest descent steps, the alternatives offer an advantage, even though (indeed because) the absolute stability limit of forward Euler is carefully yet severely violated.

The numerical solution of ill-posed problems requires suitable regularization techniques. One possible option is to consider time integration methods to solve the Showalter differential equation numerically. The stopping time of the numerical integrator corresponds to the regularization parameter. A number of well-known regularization methods such as the Landweber iteration or the Levenberg-Marquardt method can be interpreted as variants of the Euler method for solving the Showalter differential equation. Motivated by an analysis of the regularization properties of the exact solution of this equation presented by [U. Tautenhahn, Inverse Problems10 (1994) 1405–1418], we consider a variant of the exponential Euler method for solving the Showalter ordinary differential equation. We discuss a suitable discrepancy principle for selecting the step sizes within the numerical method and we review the convergence properties of [U. Tautenhahn, Inverse Problems10 (1994) 1405–1418], and of our discrete version [M. Hochbruck et al., Technical Report (2008)]. Finally, we present numerical experiments which show that this method can be efficiently implemented by using Krylov subspace methods to approximate the product of a matrix function with a vector.

Many complex systems occurring in various application share the property that the underlying Markov process remains in certain regions of the state space for long times, and that transitions between such metastable sets occur only rarely. Often the dynamics within each metastable set is of minor importance, but the transitions between these sets are crucial for the behavior and the understanding of the system. Since simulations of the original process are usually prohibitively expensive, the effective dynamics of the system, i.e. the switching between metastable sets, has to be approximated in a reliable way. This is usually done by computing the dominant eigenvectors and eigenvalues of the transfer operator associated to the Markov process. In many real applications, however, the matrix representing the spatially discretized transfer operator can be extremely large, such that approximating eigenvectors and eigenvalues is a computationally critical problem. In this article we present a novel method to determine the effective dynamics via the transfer operator without computing its dominant spectral elements. The main idea is that a time series of the process allows to approximate the sampling kernel of the process, which is an integral kernel closely related to the transition function of the transfer operator. Metastability is taken into account by representing the approximative sampling kernel by a linear combination of kernels each of which represents the process on one of the metastable sets. The effect of the approximation error on the dynamics of the system is discussed, and the potential of the new approach is illustrated by numerical examples.

We present a Monte Carlo technique for sampling from the canonical distribution in molecular dynamics. The method is built upon the Nosé-Hoover constant temperature formulation and the generalized hybrid Monte Carlo method. In contrast to standard hybrid Monte Carlo methods only the thermostat degree of freedom is stochastically resampled during a Monte Carlo step.

In [C.W. Gear, T.J. Kaper, I.G. Kevrekidis and A. Zagaris, SIAM J. Appl. Dyn. Syst. 4 (2005) 711–732], we developed a class of iterative algorithms within the context of equation-free methods to approximate low-dimensional, attracting, slow manifolds in systems of differential equations with multiple time scales. For user-specified values of a finite number of the observables, the mth member of the class of algorithms ($m = 0, 1, \ldots$) finds iteratively an approximation of the appropriate zero of the (m+1)st time derivative of the remaining variables and uses this root to approximate the location of the point on the slow manifold corresponding to these values of the observables. This article is the first of two articles in which the accuracy and convergence of the iterative algorithms are analyzed. Here, we work directly with fast-slow systems, in which there is an explicit small parameter, ε, measuring the separation of time scales. We show that, for each $m = 0, 1, \ldots$, the fixed point of the iterative algorithm approximates the slow manifold up to and including terms of ${\mathcal O}(\varepsilon^m)$. Moreover, for each m, we identify explicitly the conditions under which the mth iterative algorithm converges to this fixed point. Finally, we show that when the iteration is unstable (or converges slowly) it may be stabilized (or its convergence may be accelerated) by application of the Recursive Projection Method. Alternatively, the Newton-Krylov Generalized Minimal Residual Method may be used. In the subsequent article, we will consider the accuracy and convergence of the iterative algorithms for a broader class of systems – in which there need not be an explicit small parameter – to which the algorithms also apply.

In this paper, we consider linear ordinary differential equations originating in electronic engineering, which exhibit exceedingly rapid oscillation. Moreover, the oscillation model is completely different from the familiar framework of asymptotic analysis of highly oscillatory integrals. Using a Bessel-function identity, we expand the oscillator into asymptotic series, and this allows us to extend Filon-type approach to this setting. The outcome is a time-stepping method that guarantees high accuracy regardless of the rate of oscillation.

This paper deals with modeling the passive behavior of skeletal muscle tissue including certain microvibrations at the cell level. Our approach combines a continuum mechanics model with large deformation and incompressibility at the macroscale with chains of coupled nonlinear oscillators. The model verifies that an externally applied vibration at the appropriate frequency is able to synchronize microvibrations in skeletal muscle cells. From the numerical analysis point of view, one faces here a partial differential-algebraic equation (PDAE) that after semi-discretization in space by finite elements possesses an index up to three, depending on certain physical parameters. In this context, the consequences for the time integration as well as possible remedies are discussed.