Search

Variable stepsize, variable order (VSVO) methods have limited impact in timestepping methods in complex applications due to their computational complexity and the difficulty in implementing them in...

In the limit of infinite stiffness, the differential equations of motion of stiff mechanical systems become differential algebraic equations whose solutions stay in a constraint submanifold $\widehat{\mathcal{P}}$ of the phase space. Even though solutions of the stiff differential equations are typically oscillatory with large frequency, there exists a slow manifold $\widetilde{\mathcal{P}}$ consisting of nonoscillatory solutions; $\widetilde{\mathcal{P}}$ has the same dimension as $\widehat{\mathcal{P}}$ and converges to it as the stiffness approaches infinity. We introduce an iterative projection algorithm, IPA, that projects points in the phase space of a stiff mechanical system onto the associated slow manifold $\widetilde{\mathcal{P}}$. The algorithm is based on ideas such as micro-integration and filtering coming from the field of multiscale simulation and is applicable to initializing integration algorithms for both stiff ODEs and DAEs, including the initialization of Lagrange multipliers. We also ...

In contrast to stiff deterministic systems of ordinary differential equations, in general, the implicit Euler method for stiff stochastic differential equations is not effective. This paper introduces a new numerical method for stiff differential equations which consists of interlacing large implicit Euler time steps with a sequence of small explicit Euler time steps. We emphasize that uniform convergence with respect to the time scale separation parameter e is a desirable property of a stiff solver. We prove that the means and variances of this interlaced method converge uniformly in e for a suitably chosen test problem. We also illustrate the effectiveness of this method via some numerical examples.

In this paper, we present error corrected Euler methods for solving stiff initial value problems, which not only avoid unnecessary iteration process that may be required in most implicit methods but also have such a good stability as all implicit methods possess. The proposed methods use a Chebyshev collocation technique as well as an asymptotical linear ordinary differential equation of first-order derived from the difference between the exact solution and the Euler's polygon. These methods with or without the Jacobian are analyzed in terms of convergence and stability. In particular, it is proved that the proposed methods have a convergence order up to 4 regardless of the usage of the Jacobian. Numerical tests are given to support the theoretical analysis as evidences.

This paper constructs extrapolated implicit-explicit time stepping methods that allow one to efficiently solve problems with both stiff and nonstiff components. The proposed methods are based on Euler steps and can provide very high order discretizations of ODEs, index-1 DAEs, and PDEs in the method-of-lines framework. Implicit-explicit schemes based on extrapolation are simple to construct, easy to implement, and straightforward to parallelize. This work establishes the existence of perturbed asymptotic expansions of global errors, explains the convergence orders of these methods, and studies their linear stability properties. Numerical results with stiff ODE, DAE, and PDE test problems confirm the theoretical findings and illustrate the potential of these methods to solve multiphysics multiscale problems.

We introduce a new class of exponential integrators for the numerical integration of large-scale systems of stiff differential equations. These so-called Rosenbrock-type methods linearize the flow in each time step and make use of the matrix exponential and related functions of the Jacobian. In contrast to standard integrators, the methods are fully explicit and do not require the numerical solution of linear systems. We analyze the convergence properties of these integrators in a semigroup framework of semilinear evolution equations in Banach spaces. In particular, we derive an abstract stability and convergence result for variable step sizes. This analysis further provides the required order conditions and thus allows us to construct pairs of embedded methods. We present a third-order method with two stages, and a fourth-order method with three stages, respectively. The application of the required matrix functions to vectors are computed by Krylov subspace approximations. We briefly discuss these implementation issues, and we give numerical examples that demonstrate the efficiency of the new integrators.

We show that there exist classes of explicit numerical integration methods that can handle very stiff problems if the eigenvalues are separated into two clusters, one containing the "stiff," or fast, components, and one containing the slow components. These methods have large average step sizes relative to the fast components. Conventional implicit methods involve the solution of nonlinear equations at each step, which for large problems requires significant communication between processors on a multiprocessor machine. For such problems the methods proposed here have significant potential for speed improvement.

A second-order, L-stable Rosenbrock method from the field of stiff ordinary differential equations is studied for application to atmospheric dispersion problems describing photochemistry, advective, and turbulent diffusive transport. Partial differential equation problems of this type occur in the field of air pollution modeling. The focal point of the paper is to examine the Rosenbrock method for reliable and efficient use as an atmospheric chemical kinetics box-model solver within Strang-type operator splitting. In addition, two W-method versions of the Rosenbrock method are discussed. These versions use an inexact Jacobian matrix and are meant to provide alternatives for Strang-splitting. Another alternative for Strang-splitting is a technique based on so-called source-splitting. This technique is briefly discussed.

The Fourier--Gegenbauer (FG) method, introduced by [Gottlieb, Shu, Solomonoff, and Vandeven, ICASE Report 92-4, Hampton, VA, 1992] is aimed at removing the Gibbs phenomenon; that is, recovering the point values of a nonperiodic function from its Fourier coefficients. In this paper, we discuss some numerical aspects of the FG method related to its {\em pseudospectral\/} implementation. In particular, we analyze the behavior of the Gegenbauer series with a moderate (several hundred) number of terms suitable for computations. We also demonstrate the ability of the FG method to get a spectrally accurate approximation on small subintervals for rapidly oscillating functions or functions having steep profiles. Bearing on the previous analysis, we suggest a high-order spectral Fourier method for the solution of nonperiodic differential equations. It includes a polynomial subtraction technique to accelerate the convergence of the Fourier series and the FG algorithm to evaluate derivatives on the boundaries of nonperiodic functions. The present hybrid Fourier--Gegenbauer (HFG) method possesses better resolution properties than the original FG method. The precision of this method is demonstrated by solving stiff elliptic problems with steep solutions.

A general class of two-step Runge–Kutta methods that depend on stage values at two consecutive steps is studied. These methods are special cases of general linear methods introduced by Butcher and are quite efficient with respect to the number of function evaluations required for a given order. General order conditions are derived using the approach proposed recently by Albrecht, and examples of methods are given up to the order 5. These methods can be divided into four classes that are appropriate for the numerical solution of nonstiff or stiff differential equations in sequential or parallel computing environments.

Runge–Kutta methods applied to stiff systems in singular perturbation form are shown to give accurate approximations of phase portraits near hyperbolic stationary points. Over arbitrarily long time intervals, Runge–Kutta solutions shadow solutions of the differential equation and vice versa. Precise error bounds are derived. The proof uses attractive invariant manifolds to reduce the problem to the nonstiff case, which was previously studied by Beyn.

This work analyzes the integration of initial value problems for stiff systems of ordinary differential equations by Runge–Kutta methods. The author uses the characterization of stiff initial value problems due to Kreiss; the Jacobian matrix is essentially negative dominant and satisfies a relative Lipschitz condition. The existence and regularity of the numerical solution are established, and conditions under which the Runge–Kutta formula is stable are given.

This paper deals with the relation between differential/algebraic equations (DAEs) and certain stiff ODEs and their respective discretizations by implicit Runge–Kutta methods. For that purpose for any DAE a singular perturbed ODE is constructed such that the DAE is its reduced problem and the solution of the ODE converges in some sense to that of the DAE. Thus the DAE can be interpreted as an infinitely stiff ODE. An analysis of the discretization error of this singular perturbed system gives insight into the relationship of order-reduction phenomena observed for stiff ODEs to that for DAEs. Analysis of a general class of singularly perturbed problems and their discretizations is not attempted; however, the technique of treating singularly perturbed problems and DAEs in a unified way is new and can possibly be applied to other systems and their discretizations as well. Since asymptotic expansions are not used, but an approach similar to the ones used in B-convergence theory is applied, one can derive erro...

In a recent publication the authors showed how to reformulate the N-solition solution of the Korteweg-deVries equation in such a way that the relationship between the N elements of the solution and the solitary wave solution of the equation is clearly apparent. Here, this new formulation of the N-soliton solution is used to analyse some of the properties during interaction of a three-soliton solution of the equation.

Explicit Runge–Kutta methods are a popular way to solve the initial value problem for a system of ordinary differential equations. Although very effective for nonstiff problems, they are impractical for stiff problems. This paper is concerned with how to diagnose stiffness as the reason for unsatisfactory performance by a code based on explicit Runge–Kutta formulas.

A special class of implicit Runge–Kutta methods is developed for the numerical solution of stiff initial value problems. These formulae are derived from known singly implicit methods by adding one or more extra diagonally implicit stages. It is hoped that this modification of the original method will lead to an overall gain in efficiency, and an analysis of the advantages of making this enhancement is presented.

In this short note, the order reduction results of Petzold [SIAM J. Numer. Anal. 23(1986), pp. 837–852] for implicit Runge–Kutta methods applied to index 1 differential/algebraic systems are extended to include a larger class of methods. The relationship between the order reduction results for differential/algebraic systems and recent results for stiff systems of ordinary differential equations is explained.

The study of singularly perturbed initial value problems for nonlinear systems of ordinary differential equations parallels the analysis underlying the development of numerical algorithms for obtaining solutions to systems of stiff differential equations. This paper seeks to emphasize the advantages of combining these two substantial research efforts. It develops insight and intuition based on a sequence of solvable model problems, and it relates a variety of literature scattered throughout asymptotic and numerical analyses, stability and control theory and specific topics in applied mathematical modeling.

Certain differential equations that arise when solving chemical kinetics problems which have widely differing time constants are analyzed by a method that implicitly separates the fast reacting components from the remaining components of the system. A method for determining equilibrium of the fast reacting components is provided, and these values provide the initial conditions for an associated system of differential equations that yields a first and second order solution of the original system. Computational experiments have shown that a suitably chosen method for solving the modified system may be superior to a direct solution with a stiff differential equation solver.

The paper analyzes a number of known methods for the solution of two-point boundary value problems for ordinary differential equations; the methods considered have in common that they lead to a seq...

To be A-stable, and possibly useful for stiff systems, a Runge–Kutta formula must be implicit. There is a significant computational advantage in diagonally implicit formulae, whose coefficient matrix is lower triangular with all diagonal elements equal. We derive new, strongly S-stable diagonally implicit Runge–Kutta formulae of order 2 in 2 stages and of order 3 in 3 stages, and show that it is impossible for a strongly S-stable diagonally implicit method to attain order 4 in 4 stages. Merely A-stable diagonally implicit formulae, of order 3 in 2 stages and of order 4 in 3 stages, were previously known; we prove that no 4-stage method of this type has order 5. We describe a computer program for stiff differential equations which uses these methods, and compare them to each other and to the GEAR package.

This paper aims to assist the person who needs to solve stiff ordinary differential equations.First we identify the problem area and the basic difficulty by responding to some fundamental questions...

We present new and general numerical methods for dealing with problems whose solutions develop sharp transition layers or “near-shocks”. These methods allow many nodes automatically to concentrate in the critical regions and move with them. For clarity of exposition we concentrate on the space of piecewise linear functions with movable nodes, with Burgers’ equation as our test equation; but the generalization to much more general spaces and equations (including even certain previous “moving vorticity blobs” of the first author and S. Doss for the Navier–Stokes equations) becomes clear. In this paper we present the theoretical and computational details of our scheme, along with computational trial runs for Burgers’ equation showing that the nodes do concentrate sharply and move with the shocks as desired. The conclusiveness of these preliminary numerical trials is marred somewhat by the fact that we never successfully debugged a Newton’s method for our implicit stiff ODE solver and were thus limited to ver...

The paper describes and evaluates DSTIFF, a set of subroutines for solving stiff ordinary differential equations. The code is somewhat similar to the well-known packages LSODE, GEAR and DIFSUB but the present set of subroutines are based on least squares multistep formulas rather than the BDF. The paper describes the formulas used in the code, the structure of the code and the heuristics used, and evaluates its performance. The code seems to be much more efficient than LSODE in solving stiff equations which have Jacobians with eigenvalues having large imaginary parts.On other problems, DSTIFF is as efficient as LSODE on larger tolerances and somewhat less efficient than LSODE on stringent tolerances.

The existence of asymptotic expansions of the global discretization error for a general class of nonlinear stiff differential equations \[ y'(t) = A(t)y(t) + \varphi (t,y(t)), \] where $A(t)$ has a “stiff spectrum” characterized by a small parameter $\varepsilon $ and where $\varphi (t,y)$ is smooth, is discussed. The following methods are considered: implicit Euler, implicit midpoint, and trapezoidal rules. In strongly stiff situations ($\varepsilon $ significantly smaller than the stepsize h) the implicit Euler scheme admits a full asymptotic expansion; the same is true for the midpoint rule and for the trapezoidal rule under certain coupling conditions. In those strongly stiff cases where a full expansion does not exist for the midpoint or trapezoidal rule, the remainder term is of a reduced order but shows a regular, oscillating behavior that is described in detail. In mildly stiff situations, order reductions of the remainder term inevitably occur in any case after the start or after the change of stepsize but—as can be shown by discrete singular perturbation techniques—these order reductions are rapidly damped out as the integration proceeds. Results are illustrated by various numerical examples; in particular, numerical experience with extrapolation and defect correction is reported.

A new family of predictor-corrector schemes is designed for the numerical solution of stiff differential systems. Based on split Adams–Moulton formulas through sixth order, members of the new family achieve higher order and possess smaller error constants than corresponding split backward differentiation formulas of the same stepnumber, while maintaining similar stability properties. Some confirmation of this is obtained using a variable step implementation on test problems from the literature.

Each integration step of a stiff equation involves the solution of a nonlinear equation, usually by a quasi-Newton method which leads to a set of linear problems involving the Jacobian, J, of the differential equation. Iterative methods for these linear equations are studied. Of particular interest are methods which do not require an explicit Jacobian but can work directly with differences of function values using $J\delta \cong f(x + \delta ) - f(x)$. Some numerical experiments using a modification of LSODE are reported

A fifth-order second derivative formula is developed for stiff systems of ordinary differential equations. The exponentially fitted formulas possess a large region of absolute stability, and numerical results demonstrate increased accuracy with the same computational effort when compared with similar fourth-order formulas.

This paper is concerned with methods which consist of sequential application of several different schemes for stiff ordinary differential systems, with a predetermined ratio of step-lengths. We show that, under some conditions, such composite methods have higher order and better stability than the constituent schemes.Compositions of Obrechkoff methods and of implicit Runge–Kutta processes are investigated.

The difficulty associated with the numerical solution of stiff ordinary differential equations is considered and the stability requirements of methods suitable for stiff equations are described. A ...

The amount of work needed for a numerical method to solve a given problem depends mainly on the stepsize sequence that is used. For a given method and a given problem consider all stepsize sequences that give a global error whose norm is everywhere less than a prescribed tolerance. By an optimal stepsize sequence we mean a sequence that requires fewer steps than any other sequence. A large part of this study contains characterizations of optimal stepsize sequences for some classes of problems. From the characterization of the optimal stepsize sequences for stiff systems one can conclude that the stepsize sequence obtained with a fixed bound on the local error per unit step is far from optimal. Existing numerical methods that use a variable stepsize choose their stepsizes on the basis of either a bound on the local error per step or on a bound on the local error per unit step. For systems of stiff differential equations the stepsize sequence obtained with the first criterion, in general, requires less work...

New stability concepts—$BS$-stability and $BSI$-stability (internal $BS$-stability)—are introduced. $BS$-stability is a modification and extension of B-stability and enables the derivation of order results for Runge–Kutta methods applied to general nonlinear stiff initial value problems. These order results (B-consistency, B-convergence) are not affected by stiffness. Several classes of implicit Runge–Kutta methods are shown to be $BS$-stable: Gauss, Radau IA and Radau IIA schemes.

A useful concept of convergence for nonlinear stiff ODE’s will be developed, which permits the derivation of uniform global error bounds independent of the stiffness of the considered problem. This concept will be discussed for three simple methods for stiff systems, the implicit Euler scheme, the implicit midpoint rule and the implicit trapezoidal rule.

We study here a matrix-free method for solving stiff systems of ordinary differential equations (ODE’s). In the numerical time integration of stiff ODE initial value problems by BDF methods, the resulting nonlinear algebraic system is usually solved by a modified Newton method and an appropriate linear system algorithm. In place of that, we substitute Newton’s method (unmodified) coupled with an iterative linear system method. The latter is a projection method called the Incomplete Orthogonalization Method (IOM), developed mainly by Y. Saad. A form of IOM, with scaling included to enhance robustness, is studied in the setting of Inexact Newton Methods. The implementation requires no Jacobian matrix storage whatever. Tests on several stiff problems, of sizes up to 16,000, show the method to be quite effective and much more economical, in both computational cost and storage, than standard solution methods, at least when the problem has a certain amount of clustering in its spectrum.

Linear one step methods of a novel design are given for the numerical solution of stiff systems of ordinary differential equations. These methods permit fast integration with increments h of the independent variable adjusted to the slowly varying, dominant solutions. The integration formulas owe their efficiency to combining a reasonable order of accuracy for $h \to 0$ with a precise simulation of specific rapidly varying solutions, and with A-stability in the sense of Dahlquist.