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.     

Numerical solution of fractional order differential equations by extrapolation

We present an extrapolation type algorithm for the numerical solution of fractional order differential equations. It is based on the new result that the sequence of approximate solutions of these equations, computed by means of a recently published algorithm by Diethelm [6], possesses an asymptotic expansion with respect to the stepsize. From this we conclude that the application of extrapolation is justified, and we obtain a very efficient differential equation solver with practically no additional numerical costs. This is also illustrated by a number of numerical examples. This revised version was published online in August 2006 with corrections to the Cover Date.     

Super-fast validated solution of linear systems

Validated solution of a problem means to compute error bounds for a solution in finite precision. This includes the proof of existence of a solution. The computed error bounds are to be correct including all possible effects of rounding errors. The fastest known validation algorithm for the solution of a system of linear equations requires twice the computing time of a standard (purely) numerical algorithm. In this paper we present a super-fast validation algorithm for linear systems with symmetric positive definite matrix. This means that the entire computing time for the validation algorithm including computation of an approximated solution is the same as for a standard numerical algorithm. Numerical results are presented.     

An algorithm for the finite difference approximation of derivatives with arbitrary degree and order of accuracy

In this paper, we introduce an algorithm and a computer code for numerical differentiation of discrete functions. The algorithm presented is suitable for calculating derivatives of any degree with any arbitrary order of accuracy over all the known function sampling points. The algorithm introduced avoids the labour of preliminary differencing and is in fact more convenient than using the tabulated finite difference formulas, in particular when the derivatives are required with high approximation accuracy. Moreover, the given Matlab computer code can be implemented to solve boundary-value ordinary and partial differential equations with high numerical accuracy. The numerical technique is based on the undetermined coefficient method in conjunction with Taylor’s expansion. To avoid the difficulty of solving a system of linear equations, an explicit closed form equation for the weighting coefficients is derived in terms of the elementary symmetric functions. This is done by using an explicit closed formula for the Vandermonde matrix inverse. Moreover, the code is designed to give a unified approximation order throughout the given domain. A numerical differentiation example is used to investigate the validity and feasibility of the algorithm and the code. It is found that the method and the code work properly for any degree of derivative and any order of accuracy.     

Accurate computation of the zeros of the generalized Bessel polynomials

Summary. A general method for approximating polynomial solutions of second-order linear homogeneous differential equations with polynomial coefficients is applied to the case of the families of differential equations defining the generalized Bessel polynomials, and an algorithm is derived for simultaneously finding their zeros. Then a comparison with several alternative algorithms is carried out. It shows that the computational problem of approximating the zeros of the generalized Bessel polynomials is not an easy matter at all and that the only algorithm able to give an accurate solution seems to be the one presented in this paper. Received July 25, 1997 / Revised version received May 19, 1999 / Published online June 8, 2000     

Error and extrapolation of a compact LOD method for parabolic differential equations

This paper is concerned with a compact locally one-dimensional (LOD) finite difference method for solving two-dimensional nonhomogeneous parabolic differential equations. An explicit error estimate for the finite difference solution is given in the discrete infinity norm. It is shown that the method has the accuracy of the second-order in time and the fourth-order in space with respect to the discrete infinity norm. A Richardson extrapolation algorithm is developed to make the final computed solution fourth-order accurate in both time and space when the time step equals the spatial mesh size. Numerical results demonstrate the accuracy and the high efficiency of the extrapolation algorithm.     

Multilevel Schwarz methods

Summary We consider the solution of the algebraic system of equations which result from the discretization of second order elliptic equations. A class of multilevel algorithms are studied using the additive Schwarz framework. We establish that the condition number of the iteration operators are bounded independent of mesh sizes and the number of levels. This is an improvement on Dryja and Widlund's result on a multilevel additive Schwarz algorithm, as well as Bramble, Pasciak and Xu's result on the BPX algorithm. Some multiplicative variants of the multilevel methods are also considered. We establish that the energy norms of the corresponding iteration operators are bounded by a constant less than one, which is independent of the number of levels. For a proper ordering, the iteration operators correspond to the error propagation operators of certain V-cycle multigrid methods, using Gauss-Seidel and damped Jacobi methods as smoothers, respectively.This work was supported in part by the National Science Foundation under Grants NSF-CCR-8903003 at Courant Institute of Mathematical Sciences, New York University and NSF-ASC-8958544 at Department of Computer Science, University of Maryland.     

Multigrid methods for boundary integral equations

Summary Multigrid methods are applied for solving algebraic systems of equations that occur to the numerical treatment of boundary integral equations of the first and second kind. These methods, originally formulated for partial differential equations of elliptic type, combine relaxation schemes and coarse grid corrections. The choice of the relaxation scheme is found to be essential to attain a fast convergent iterative process. Theoretical investigations show that the presented relaxation scheme provides a multigrid algorithm of which the rate of convergence increases with the dimension of the finest grid. This is illustrated for the calculation of potential flow around an aerofoil.     

A modified accelerated monotone iterative method for finite difference reaction-diffusion-convection equations

This paper is concerned with monotone algorithms for the finite difference solutions of a class of nonlinear reaction-diffusion-convection equations with nonlinear boundary conditions. A modified accelerated monotone iterative method is presented to solve the finite difference systems for both the time-dependent problem and its corresponding steady-state problem. This method leads to a simple and yet efficient linear iterative algorithm. It yields two sequences of iterations that converge monotonically from above and below, respectively, to a unique solution of the system. The monotone property of the iterations gives concurrently improving upper and lower bounds for the solution. It is shown that the rate of convergence for the sum of the two sequences is quadratic. Under an additional requirement, quadratic convergence is attained for one of these two sequences. In contrast with the existing accelerated monotone iterative methods, our new method avoids computing local maxima in the construction of these sequences. An application using a model problem gives numerical results that illustrate the effectiveness of the proposed method.     

Parallel domain decomposition procedures of improved D-D type for parabolic problems

Two parallel domain decomposition procedures for solving initial-boundary value problems of parabolic partial differential equations are proposed. One is the extended D-D type algorithm, which extends the explicit/implicit conservative Galerkin domain decomposition procedures, given in [5], from a rectangle domain and its decomposition that consisted of a stripe of sub-rectangles into a general domain and its general decomposition with a net-like structure. An almost optimal error estimate, without the factor H^−1/2 given in Dawson-Dupont’s error estimate, is proved. Another is the parallel domain decomposition algorithm of improved D-D type, in which an additional term is introduced to produce an approximation of an optimal error accuracy in L²-norm.     

Frozen divided difference scheme for solving systems of nonlinear equations

The development of an inverse first-order divided difference operator for functions of several variables, as well as a direct computation of the local order of convergence of an iterative method is presented. A generalized algorithm of the secant method for solving a system of nonlinear equations is studied and the maximum computational efficiency is computed. Furthermore, a sequence that approximates the order of convergence is generated for the examples and it confirms in a numerical way that the order of the methods is well deduced.     

Minimum residual methods for augmented systems

For large systems of linear equations, iterative methods provide attractive solution techniques. We describe the applicability and convergence of iterative methods of Krylov subspace type for an important class of symmetric and indefinite matrix problems, namely augmented (or KKT) systems. Specifically, we consider preconditioned minimum residual methods and discuss indefinite versus positive definite preconditioning. For a natural choice of starting vector we prove that when the definite and indenfinite preconditioners are related in the obvious way, MINRES (which is applicable in the case of positive definite preconditioning) and full GMRES (which is applicable in the case of indefinite preconditioning) give residual vectors with identical Euclidean norm at each iteration. Moreover, we show that the convergence of both methods is related to a system of normal equations for which the LSQR algorithm can be employed. As a side result, we give a rare example of a non-trivial normal(1) matrix where the corresponding inner product is explicitly known: a conjugate gradient method therefore exists and can be employed in this case. This work was supported by British Council/German Academic Exchange Service Research Collaboration Project 465 and NATO Collaborative Research Grant CRG 960782     

The fast quartic solver

A fast and highly accurate algorithm for solving quartic equations is introduced. This new algorithm is more than six times as fast and several times more accurate than the quasi-standard Companion matrix eigenvalue quartic solver. Moreover, the method is exceptionally robust in cases of extreme root spread. The new algorithm is based on a factorization of the quartic in two quadratics, which are solved in closed form. The performance key at this point is a fixed-point iteration based fitting algorithm for backward optimization of the underlying quartic-to-quadratic polynomial decomposition. Detailed experimental results confirm our claims.     

A Crank-Nicolson scheme for the Landau-Lifshitz equation without damping

An accurate and efficient numerical approach, based on a finite difference method with Crank-Nicolson time stepping, is proposed for the Landau-Lifshitz equation without damping. The phenomenological Landau-Lifshitz equation describes the dynamics of ferromagnetism. The Crank-Nicolson method is very popular in the numerical schemes for parabolic equations since it is second-order accurate in time. Although widely used, the method does not always produce accurate results when it is applied to the Landau-Lifshitz equation. The objective of this article is to enumerate the problems and then to propose an accurate and robust numerical solution algorithm. A discrete scheme and a numerical solution algorithm for the Landau-Lifshitz equation are described. A nonlinear multigrid method is used for handling the nonlinearities of the resulting discrete system of equations at each time step. We show numerically that the proposed scheme has a second-order convergence in space and time.     

An ODE-based trust region method for unconstrained optimization problems

In this paper, a new trust region algorithm is proposed for solving unconstrained optimization problems. This method can be regarded as a combination of trust region technique, fixed step-length and ODE-based methods. A feature of this proposed method is that at each iteration, only a system of linear equations is solved to obtain a trial step. Another is that when a trial step is not accepted, the method generates an iterative point whose step-length is defined by a formula. Under some standard assumptions, it is proven that the algorithm is globally convergent and locally superlinear convergent. Preliminary numerical results are reported.     

BiCGStab,VPAStab and an adaptation to mildly nonlinear systems

The key equations of BiCGStab are summarised to show its connections with Padé and vector-Padé approximation. These considerations lead naturally to stabilised vector-Padé approximation of a vector-valued function (VPAStab), and an algorithm for the acceleration of convergence of a linearly generated sequence of vectors. A generalisation of this algorithm for the acceleration of convergence of a nonlinearly generated system is proposed here, and comparative numerical results are given.     

A fourth-order derivative-free algorithm for nonlinear equations

A two-step derivative-free iterative algorithm is presented for solving nonlinear equations. Error analysis shows that the algorithm is fourth-order with efficiency index equal to 1.5874. A lot of numerical results show that the algorithm is effective and is preferable to some existing derivative-free methods in terms of computation cost.     

On computation of solution curves for semilinear elliptic problems

We consider computation of solution curves for semilinear elliptic equations. In case solution is stable, we present an algorithm with monotone convergence, which is a considerable improvement of the corresponding schemes in [4] and [5]. For the unstable solutions, we show how to construct a fourth-order evolution equation, for which the same solution will be stable.     

A communication-less parallel algorithm for tridiagonal Toeplitz systems

Diagonally dominant tridiagonal Toeplitz systems of linear equations arise in many application areas and have been well studied in the past. Modern interest in numerical linear algebra is often focusing on solving classic problems in parallel. In McNally [Fast parallel algorithms for tri-diagonal symmetric Toeplitz systems, MCS Thesis, University of New Brunswick, Saint John, 1999], an m processor Split & Correct algorithm was presented for approximating the solution to a symmetric tridiagonal Toeplitz linear system of equations. Nemani [Perturbation methods for circulant-banded systems and their parallel implementation, Ph.D. Thesis, University of New Brunswick, Saint John, 2001] and McNally (2003) adapted the works of Rojo [A new method for solving symmetric circulant tri-diagonal system of linear equations, Comput. Math. Appl. 20 (1990) 61–67], Yan and Chung [A fast algorithm for solving special tri-diagonal systems, Computing 52 (1994) 203–211] and McNally et al. [A split-correct parallel algorithm for solving tri-diagonal symmetric Toeplitz systems, Internat. J. Comput. Math. 75 (2000) 303–313] to the non-symmetric case. In this paper we present relevant background from these methods and then introduce an m processor scalable communication-less approximation algorithm for solving a diagonally dominant tridiagonal Toeplitz system of linear equations.     

Reliable parallel solution of bidiagonal systems

Summary. This paper presents a new efficient algorithm for solving bidiagonal systems of linear equations on massively parallel machines. We use a divide and conquer approach to compute a representative subset of the solution components after which we solve the complete system in parallel with no communication overhead. We address the numerical properties of the algorithm in two ways: we show how to verify the à posteriori backward stability at virtually no additional cost, and prove that the algorithm is à priori forward stable. We then show how we can use the algorithm in order to bound the possible perturbations in the solution components. Received March 13, 1998 / Revised version received December 21, 1999 / Published online June 20, 2001