The integration to steady state of many initial value ODEs and PDEs using the forward Euler methodcan alternatively be considered as gradient descent for an associated minimization problem.Greedy algorithms such as steepest descent for determining the step size are asslow 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 thefirst few steps have for certain inverse problems,amounting to a finite time regularization. We further investigate the retention of thisproperty 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 steepestdescent steps, the alternatives offer an advantage, even though (indeed because)the absolute stability limit of forward Euler is carefully yet severely violated.

The problem of modeling acoustic waves scattered by an object with Neumann boundary condition is considered. The boundary condition is taken into account by means of the fictitious domain method, yielding a first order in time mixed variational formulation for the problem. The resulting system is discretized with two families of mixed finite elements that are compatible with mass lumping. We present numerical results illustrating that the Neumann boundary condition on the object is not always correctly taken into account when the first family of mixed finite elements is used. We, therefore, introduce the second family of mixed finite elements for which a theoretical convergence analysis is presented and error estimates are obtained. A numerical study of the convergence is also considered for a particular object geometry which shows that our theoretical error estimates are optimal.

We derive a posteriori estimates for a discretization in space of the standardCahn–Hilliard equation with a double obstacle free energy.The derived estimates are robust and efficient, and in practice are combinedwith a heuristic time step adaptation. We present numerical experiments in two and three space dimensions and compareour method with an existing heuristic spatial mesh adaptation algorithm.

The atomistic to continuum interface for quasicontinuum energiesexhibits nonzero forces under uniform strain that have beencalled ghost forces.In this paper,we prove for a linearization of a one-dimensional quasicontinuum energyaround a uniform strainthat the effect of the ghost forces on the displacementnearly cancels and has a small effect on the error away from the interface.We give optimal order error estimatesthat show that the quasicontinuum displacementconverges to the atomistic displacement at the rate O(h)in the discrete $\ell^\infty$ andw 1,1 norms where h is the interatomic spacing.We also give a proof that the error in the displacement gradientdecays away from the interface to O(h) at distance O(h|logh|)in the atomistic region and distance O(h) in the continuum region.Our work gives an explicit and simplified form for the decay of the effect of theatomistic to continuum coupling error in terms of a general underlying interatomic potential and gives the estimates described above in the discrete $\ell^\infty$ and w1,p norms.

We prove convergence and quasi-optimal complexity of an adaptive finite element algorithm on triangular meshes with standard mesh refinement. Our algorithm is based on an adaptive marking strategy. In each iteration, a simple edge estimator is compared to an oscillation term and the marking of cells for refinement is done according to the dominant contribution only. In addition, we introduce an adaptive stopping criterion for iterative solution which compares an estimator for the iteration error with the estimator for the discretization error.

The numerical solution of ill-posed problems requires suitableregularization techniques. One possible option is to consider timeintegration methods to solve the Showalter differential equationnumerically. The stopping time of the numerical integrator correspondsto the regularization parameter. A number of well-knownregularization methods such as the Landweber iteration or theLevenberg-Marquardt method can be interpreted as variants of theEuler method for solving the Showalter differential equation. Motivated by an analysis of the regularization properties of theexact solution of this equation presented by [U. Tautenhahn, Inverse Problems10 (1994) 1405–1418], we consider a variant of the exponential Euler methodfor solving the Showalter ordinary differential equation. We discuss asuitable discrepancy principle for selecting the step sizes withinthe 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 thismethod can be efficiently implemented by using Krylov subspacemethods to approximate the product of a matrix function with avector.


We show that the entropy method, that has been used successfully in orderto prove exponential convergence towards equilibrium with explicit constants in many contexts,among which reaction-diffusion systems coming out of reversible chemistry, can also be usedwhen one considers a reaction-diffusion system corresponding to an irreversible mechanism ofdissociation/recombination, for which no natural entropy is available.


In this paper, we present numerical methodsfor the determination of solitons, that consist in spatially localizedstationary states of nonlinear scalar equations or coupled systemsarising in nonlinear optics.We first use the well-known shooting method in order to findexcited states (characterized by the number k of nodes) for theclassical nonlinear Schrödinger equation. Asymptotics can thenbe derived in the limits of either large k are large nonlinear exponents σ.In a second part, we compute solitons for a nonlinearsystem governing the propagation of two coupled waves in a quadratic media in anyspatial dimension, starting from one-dimensional states obtained with a shooting method and considering the dimension as a continuation parameter. Finally, we investigate the case of three wavemixing, for which the shooting method is not relevant.

Many complex systems occurring in various application share the property that theunderlying 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 transferoperator 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 eigenvaluesis 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 thecanonical distribution in molecular dynamics. The method is built uponthe Nosé-Hoover constant temperature formulation and the generalizedhybrid Monte Carlo method. In contrast to standard hybrid Monte Carlo methodsonly 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 developeda class of iterative algorithmswithin the contextof equation-free methodsto approximatelow-dimensional,attracting,slow manifoldsin systemsof differential equationswith multiple time scales.For user-specified valuesof a finite numberof the observables,the mth memberof the classof algorithms( $m = 0, 1, \ldots$ )finds iterativelyan approximationof the appropriate zeroof the (m+1)st time derivativeof the remaining variablesanduses this rootto approximate the locationof the pointon the slow manifoldcorresponding to these valuesof the observables.This articleis the firstof two articlesin whichthe accuracy and convergenceof the iterative algorithmsare analyzed.Here,we work directlywith fast-slow systems,in which there isan explicit small parameter, ε,measuring the separationof time scales.We show that,for each $m = 0, 1, \ldots$ ,the fixed pointof the iterative algorithmapproximates the slow manifoldup to and includingterms of ${\mathcal O}(\varepsilon^m)$ .Moreover,for each m,we identify explicitlythe conditionsunder whichthe mth iterative algorithmconverges to this fixed point.Finally,we show thatwhenthe iterationis unstable(orconverges slowly)it may be stabilized(orits convergencemay be accelerated)by applicationof the Recursive Projection Method.Alternatively,the Newton-KrylovGeneralized Minimal Residual Methodmay be used.In the subsequent article,we will considerthe accuracy and convergenceof the iterative algorithmsfor a broader classof systems – in whichthere need not bean explicitsmall parameter – to whichthe 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 passivebehavior of skeletal muscle tissue includingcertain 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 appliedvibration 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 possessesan index up to three, depending on certain physicalparameters. In this context, the consequences forthe time integration as well as possible remedies are discussed.