Blended extended linear multistep methods for the accurate numerical integration of stiff initial value problems

A class of blended extended linear multistep methods suitable for the approximate numerical integration of stiff systems of first order ordinary differential equations is described. These methods are formulated as a result of combining the second derivative extended backward differentiation formulae of Cash and the blended linear multistep methods of Skeel and Kong. The new methods combine a high order or accuracy with good stability properties and, as a direct consequence, they are often suitable for the numerical integration of stiff differential systems when high accuracy is requested. In the first part of the present paper we consider the derivation of these new blended methods and give the coefficients and stability regions for formulae of order up to and including 10. In the second half we consider their practical implementation. In particular we describe a variable order/variable step package based on these blended formulae and we evaluate the performance of this package on the well known DETEST test set. It is shown that the new code is reliable on this test set and is competitive with the well known second derivative method of Enright.     

Splitting methods for second‐order initial value problems

We consider implicit integration methods for the solution of stiff initial value problems for second-order differential equations of the special form y' = f(y). In implicit methods, we are faced with the problem of solving systems of implicit relations. This paper focuses on the construction and analysis of iterative solution methods which are effective in cases where the Jacobian of the right‐hand side of the differential equation can be split into a sum of matrices with a simple structure. These iterative methods consist of the modified Newton method and an iterative linear solver to deal with the linear Newton systems. The linear solver is based on the approximate factorization of the system matrix associated with the linear Newton systems. A number of convergence results are derived for the linear solver in the case where the Jacobian matrix can be split into commuting matrices. Such problems often arise in the spatial discretization of time‐dependent partial differential equations. Furthermore, the stability matrix and the order of accuracy of the integration process are derived in the case of a finite number of iterations. This revised version was published online in August 2006 with corrections to the Cover Date.     

Stein implicit Runge-Kutta methods with high stage order for large-scale ordinary differential equations

The Runge-Kutta method is one of the most popular implicit methods for the solution of stiff ordinary differential equations. For large problems, the main drawback of such methods is the cost required at each integration step for computing the solution of a nonlinear system of equations. In this paper, we propose to reduce the cost of the computation by transforming the linear systems arising in the application of Newton's method to Stein matrix equations. We propose an iterative projection method onto block Krylov subspaces for solving numerically such Stein matrix equations. Numerical examples are given to illustrate the performance of our proposed method.     

Approximation, integration and differentiation of time functions using a set of orthogonal hybrid functions (HF) and their application to solution of first order differential equations

Differential equations of different types and orders are of utmost importance for mathematical modeling of control system problems. State variable method uses the concept of expressing n number of first order differential equations in vector matrix form to model and analyze/synthesize control systems.The present work proposes a new set of orthogonal hybrid functions (HF) which evolved from synthesis of sample-and-hold functions (SHF) and triangular functions (TF). This HF set is used to approximate a time function in a piecewise linear manner with the mean integral square error (MISE) much less than block pulse function based approximation which always provides staircase solutions.The operational matrices for integration and differentiation in HF domain are also derived and employed for solving non-homogeneous and homogeneous differential equations of the first order as well as state equations. The results are compared with exact solutions, the 4th order Runge-Kutta method and its further improved versions proposed by Simos [6]. The presented HF domain theory is well supported by a few illustrations.     

Efficient transient noise analysis in circuit simulation

Stochastic differential algebraic equations (SDAEs) arise as a mathematical model for electrical network equations that are influenced by additional sources of Gaussian white noise. We discuss adaptive linear multi-step methods for their numerical integration, in particular stochastic analogues of the trapezoidal rule and the two-step backward differentiation formula, and we obtain conditions that ensure mean-square convergence of this methods. For the case of small noise we present a strategy for controlling the step-size in the numerical integration. It is based on estimating the mean-square local errors and leads to step-size sequences that are identical for all computed paths. Test results illustrate the performance of the presented methods. (© 2006 WILEY-VCH Verlag GmbH & Co. KGaA, Weinheim)     

On a class of cyclic methods for the numerical integration of stiff systems of O.D.E.s

A class of cyclic linear multistep methods suitable for the approximate numerical integration of stiff systems of first order ordinary differential equations is developed. Particular attention is paid to the problem of deriving schemes which are almostA-stable, self starting, have relatively high orders of accuracy and contain a built in error estimate. These requirements demand that the linear multistep methods which are used are solved iteratively rather than directly in the usual way and an efficient method for doing this is suggested. Finally the algorithms are illustrated by application to a particular test problem.     

Numerical solution of the finite horizon stochastic linear quadratic control problem

The treatment of the stochastic linear quadratic optimal control problem with finite time horizon requires the solution of stochastic differential Riccati equations. We propose efficient numerical methods, which exploit the particular structure and can be applied for large‐scale systems. They are based on numerical methods for ordinary differential equations such as Rosenbrock methods, backward differentiation formulas, and splitting methods. The performance of our approach is tested in numerical experiments.     

Adaptive partitioning techniques for ordinary differential equations

Many stiff systems of ordinary differential equations (ODEs) modeling practical problems can be partitioned into loosely coupled subsystems. In this paper the objective of the partitioning is to permit the numerical integration of one time step to be performed as the solution of a sequence of small subproblems. This reduces the computational complexity compared to solving one large system and permits efficient parallel execution under appropriate conditions. The subsystems are integrated using methods based on low order backward differentiation formulas.This paper presents an adaptive partitioning algorithm based on a classical graph algorithm and techniques for the efficient evaluation of the error introduced by the partitioning.The power of the adaptive partitioning algorithm is demonstrated by a real world example, a variable step-size integration algorithm which solves a system of ODEs originating from chemical reaction kinetics. The computational savings are substantial. In memory of Germund Dahlquist (1925–2005).AMS subject classification (2000) 65L06, 65Y05     

Splitting and composition methods for explicit time dependence in separable dynamical systems

We consider splitting methods for the numerical integration of separable non-autonomous differential equations. In recent years, splitting methods have been extensively used as geometric numerical integrators showing excellent performances (both qualitatively and quantitatively) when applied on many problems. They are designed for autonomous separable systems, and a substantial number of methods tailored for different structures of the equations have recently appeared. Splitting methods have also been used for separable non-autonomous problems either by solving each non-autonomous part separately or after each vector field is frozen properly. We show that both procedures correspond to introducing the time as two new coordinates. We generalize these results by considering the time as one or more further coordinates which can be integrated following either of the previous two techniques. We show that the performance as well as the order of the final method can strongly depend on the particular choice. We present a simple analysis which, in many relevant cases, allows one to choose the most appropriate split to retain the high performance the methods show on the autonomous problems. This technique is applied to different problems and its performance is illustrated for several numerical examples.     

High order approximations using spline-based differential quadrature method: Implementation to the multi-dimensional PDEs

A new differential quadrature method based on cubic B-spline is developed for the numerical solution of differential equations. In order to develop the new approach, the B-spline basis functions are used on the grid and midpoints of a uniform partition. Some error bounds are obtained by help of cubic spline collocation, which show that the method in its classic form is second order convergent. In order to derive higher accuracy, high order perturbations of the problem are generated and applied to construct the numerical algorithm. A new fourth order method is developed for the numerical solution of systems of second order ordinary differential equations. By solving some test problems, the performance of the proposed methods is examined. Also the implementation of the method for multi-dimensional time dependent partial differential equations is presented. The stability of the proposed methods is examined via matrix analysis. To demonstrate the applicability of the algorithms, we solve the 2D and 3D coupled Burgers’ equations and 2D sine-Gordon equation as test problems. Also the coefficient matrix of the methods for multi-dimensional problems is described to analyze the stability.     

Contractive methods for stiff differential equations part I

An integration method for ordinary differential equations is said to be contractive if all numerical solutions of the test equationx=x generated by that method are not only bounded (as required for stability) but non-increasing. We develop a theory of contractivity for methods applied to stiff and non-stiff, linear and nonlinear problems. This theory leads to the design of a collection of specific contractive Adams-type methods of different orders of accuracy which are optimal with respect to certain measures of accuracy and/or contractivity. Theoretical and numerical results indicate that some of these novel methods are more efficient for solving problems with a lack of smoothness than are the familiar backward differentiation methods. This lack of smoothness may be either inherent in the problem itself, or due to the use of strongly varying integration steps. In solving smooth problems, the efficiency of the low-order contractive methods we propose is approximately the same as that of the corresponding backward differentiation methods.This work was done during the first author's stay at the IBM Thomas J. Watson Research Center under his appointment as Senior Researcher of the Academy of Finland. It was sponsored in addition by the IBM Corporation and by the AirForce Office of Scientific Research (AFSC), United States Air Force, under contracts No. F44620-75-C-0058 and F49620-77-C-0088. The United States Government is authorized to reproduce and distribute reprints for governmental purposes notwithstanding any copyright notation hereon.     

Multiblock matrix methods for solution of integrodifferential equations of (0, k) type

Easily implementable numerical-analytical (so-called matrix) methods are developed and analyzed for approximate solution of one-dimensional linear integrodifferential equations with constant integration limits containing a differential expression only in the integrand. Zero-rank matrix algorithms of numerical differentiation and integration are used for construction. The original equation is reduced to a linear algebraic vector equation for the values of the sought solution at the grid points and to an interpolation formula.Translated from Vychislitel'naya i Prikladnaya Matematika, No. 55, pp. 69–75, 1985.     

Solving boundary value problems,integral, and integro-differential equations using Gegenbauer integration matrices

We introduce a hybrid Gegenbauer (ultraspherical) integration method (HGIM) for solving boundary value problems (BVPs), integral and integro-differential equations. The proposed approach recasts the original problems into their integral formulations, which are then discretized into linear systems of algebraic equations using Gegenbauer integration matrices (GIMs). The resulting linear systems are well-conditioned and can be easily solved using standard linear system solvers. A study on the error bounds of the proposed method is presented, and the spectral convergence is proven for two-point BVPs (TPBVPs). Comparisons with other competitive methods in the recent literature are included. The proposed method results in an efficient algorithm, and spectral accuracy is verified using eight test examples addressing the aforementioned classes of problems. The proposed method can be applied on a broad range of mathematical problems while producing highly accurate results. The developed numerical scheme provides a viable alternative to other solution methods when high-order approximations are required using only a relatively small number of solution nodes.     

Solving differential equations using modified Picard iteration

Many classes of differential equation are shown to be open to solution through a method involving a combination of a direct integration approach with suitably modified Picard iterative procedures. The classes of differential equations considered include typical initial value, boundary value and eigenvalue problems arising in physics and engineering and include non-linear as well as linear differential equations. Examples involving partial as well as ordinary differential equations are presented. The method is easy to implement on a computer and the solutions so obtained are essentially power series. With its conceptual clarity (differential equations are integrated directly), its uniform methodology (the overall approach is the same in all cases) and its straightforward computer implementation (the integration and iteration procedures require only standard commercial software), the modified Picard methods offer obvious benefits for the teaching of differential equations as well as presenting a basic but flexible tool-kit for the solution process itself.     

On a type of matrix splitting preconditioners for a class of block two-by-two linear systems

Recently, Bai et al. (2013) proposed an effective and efficient matrix splitting iterative method, called preconditioned modified Hermitian/skew-Hermitian splitting (PMHSS) iteration method, for two-by-two block linear systems of equations. The eigenvalue distribution of the iterative matrix suggests that the splitting matrix could be advantageously used as a preconditioner. In this study, the CGNR method is utilized for solving the PMHSS preconditioned linear systems, and the performance of the method is considered by estimating the condition number of the normal equations. Furthermore, the proposed method is compared with other PMHSS preconditioned Krylov subspace methods by solving linear systems arising in complex partial differential equations and a distributed control problem. The numerical results demonstrate the difference in the performance of the methods under consideration.     

Stable barrier-projection and barrier-Newton methods in linear programming

The present paper is devoted to the application of the space transformation techniques for solving linear programming problems. By using a surjective mapping the original constrained optimization problem is transformed to a problem in a new space with only equality constraints. For the numerical solution of the latter problem the stable version of the gradient-projection and Newton's methods are used. After an inverse transformation to the original space a family of numerical methods for solving optimization problems with equality and inequality constraints is obtained. The proposed algorithms are based on the numerical integration of the systems of ordinary differential equations. These algorithms do not require feasibility of the starting and current points, but they preserve feasibility. As a result of a space transformation the vector fields of differential equations are changed and additional terms are introduced which serve as a barrier preventing the trajectories from leaving the feasible set. A proof of a convergence is given.Dedicated to Professor George B. Dantzig on the occasion of his eightieth birthday.Research was supported by the grant N93-012-450 from Russian Scientific Fund.     

A-priori error estimates of Galerkin backward differentiation methods in time-inhomogeneous parabolic problems

Summary Backward differentiation methods up to orderk=5 are applied to solve linear ordinary and partial (parabolic) differential equations where in the second case the space variables are discretized by Galerkin procedures. Using a mean square norm over all considered time levels a-priori error estimates are derived. The emphasis of the results lies on the fact that the obtained error bounds do not depend on a Lipschitz constant and the dimension of the basic system of ordinary differential equations even though this system is allowed to have time-varying coefficients. It is therefore possible to use the bounds to estimate the error of systems with arbitrary varying dimension as they arise in the finite element regression of parabolic problems.     

A class of one-step one-stage methods for stiff systems of ordinary differential equations

A new class of one-step one-stage methods (ABC-schemes) designed for the numerical solution of stiff initial value problems for ordinary differential equations is proposed and studied. The Jacobian matrix of the underlying differential equation is used in ABC-schemes. They do not require iteration: a system of linear algebraic equations is once solved at each integration step. ABC-schemes are A- and L-stable methods of the second order, but there are ABC-schemes that have the fourth order for linear differential equations. Some aspects of the implementation of ABC-schemes are discussed. Numerical results are presented, and the schemes are compared with other numerical methods.     

A method for integration of unstable systems of ordinary differential equation subject to two-point boundary conditions

Instability problems in systems of differential equations are discussed. A matrix technique is given for producing numerical solutions to a system of ordinary differential equations with boundary conditions specified at each end of the interval when the system contains dominant solutions which give rise to numerical instability in conventional integration methods. A method of bringing up the initial conditions is described, whereby the two-point nature of the problem is made use of to stabilize the system. Three numerical examples are included.