Showing papers in "Numerical Algorithms in 1997"
TL;DR: It is shown that apparently innocuous algorithmic modifications to the Padé iteration can lead to instability, and a perturbation analysis is given to provide some explanation.
Abstract: Any matrix with no nonpositive real eigenvalues has a unique square root for which every eigenvalue lies in the open right half-plane. A link between the matrix sign function and this square root is exploited to derive both old and new iterations for the square root from iterations for the sign function. One new iteration is a quadratically convergent Schulz iteration based entirely on matrix multiplication; it converges only locally, but can be used to compute the square root of any nonsingular M-matrix. A new Pade iteration well suited to parallel implementation is also derived and its properties explained. Iterative methods for the matrix square root are notorious for suffering from numerical instability. It is shown that apparently innocuous algorithmic modifications to the Pade iteration can lead to instability, and a perturbation analysis is given to provide some explanation. Numerical experiments are included and advice is offered on the choice of iterative method for computing the matrix square root.
207 citations
TL;DR: It is concluded that the application of extrapolation is justified, and the algorithm is obtained a very efficient differential equation solver with practically no additional numerical costs.
Abstract: 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.
182 citations
TL;DR: The cocycle description of a nonaut autonomous system and the concept of a cocycle attractor are reviewed in the context of nonautonomous ordinary differential equations and variable time-step numerical schemes for autonomous ordinary differential equation.
Abstract: The concept of an attractor for autonomous systems is generally too restrictive in the nonautonomous context An appropriate generalization is the cocycle attractor which consists of a family of equivariant sets Here the cocycle description of a nonautonomous system and the concept of a cocycle attractor are reviewed in the context of nonautonomous ordinary differential equations and variable time-step numerical schemes for autonomous ordinary differential equations In the latter case, theorems are stated for the existence and convergence of numerical cocycle attractors to an assumed attractor of an autonomous ordinary differential equations
143 citations
TL;DR: The results of numerical experiments suggest that sparse approximate inverse preconditioning is a viable approach for the solution of large-scale dense linear systems on parallel computers.
Abstract: We investigate the use of sparse approximate inverse preconditioners for the iterative solution of linear systems with dense complex coefficient matrices arising in industrial electromagnetic problems. An approximate inverse is computed via a Frobenius norm approach with a prescribed nonzero pattern. Some strategies for determining the nonzero pattern of an approximate inverse are described. The results of numerical experiments suggest that sparse approximate inverse preconditioning is a viable approach for the solution of large-scale dense linear systems on parallel computers.
134 citations
TL;DR: The so-called Linear-time Legendre Transform (LLT) improves all previously known Fast Legendre transform algorithms by reducing their log-linear worst-case time complexity to linear.
Abstract: A new algorithm to compute the Legendre–Fenchel transform is proposed and investigated. The so-called Linear-time Legendre Transform (LLT) improves all previously known Fast Legendre Transform algorithms by reducing their log-linear worst-case time complexity to linear. Since the algorithm amounts to computing several convex hulls and sorting, any convex hull algorithm well-suited for a particular problem gives a corresponding LLT algorithm. After justifying the convergence of the Discrete Legendre Transform to the Legendre–Fenchel transform, an extended computation time complexity analysis is given and confirmed by numerical tests. Finally, the LLT is illustrated with several examples and a LLT MATLAB package is described.
99 citations
TL;DR: The Crank-Nicolson difference scheme with modifications near the interface is used to solve for the solution u(x,t) and the interface α (t) simultaneously and second order accuracy on uniform grids is confirmed both for the solutions and the position of the interface.
Abstract: A second order difference method is developed for the nonlinear moving interface problem of the form
$$u_t + \lambda uu_x = \left( {\beta u_x } \right)_x - f\left( {x,t} \right),x \in \left[ {\left. {0,\alpha } \right) \cup \left( {\left. {\alpha ,1} \right]} \right.,} \right.$$
$$\frac{{d\alpha }}{{dt}} = w\left( {t,\alpha ;u,u_x } \right),$$
, where α (t) is the moving interface. The coefficient β(x,t) and the source term f(x,t) can be discontinuous across α (t) and moreover, f(x,t) may have a delta or/and delta-prime function singularity there. As a result, although the equation is parabolic, the solution u and its derivatives may be discontinuous across α (t). Two typical interface conditions are considered. One condition occurs in Stefan-like problems in which the solution is known on the interface. A new stable interpolation strategy is proposed. The other type occurs in a one-dimensional model of Peskin’s immersed boundary method in which only jump conditions are given across the interface. The Crank-Nicolson difference scheme with modifications near the interface is used to solve for the solution u(x,t) and the interface α (t) simultaneously. Several numerical examples, including models of ice-melting and glaciation, are presented. Second order accuracy on uniform grids is confirmed both for the solution and the position of the interface.
83 citations
TL;DR: A two-tage AOR method is presented, which particularly uses the AOR iteration as the inner iteration and is substantially a relaxed variant of the afore-presented method.
Abstract: The discretizations of many differential equations by the finite difference or the finite element methods can often result in a class of system of weakly nonlinear equations. In this paper, by applying the two-tage iteration technique and in accordance with the special properties of this weakly nonlinear system, we first propose a general two-tage iterative method through the two-tage splitting of the system matrix. Then, by applying the accelerated overrelaxation (AOR) technique of the linear iterative methods, we present a two-tage AOR method, which particularly uses the AOR iteration as the inner iteration and is substantially a relaxed variant of the afore-presented method. For these two classes of methods, we establish their local convergence theories, and precisely estimate their asymptotic convergence factors under some suitable assumptions when the involved nonlinear mapping is only B-differentiable. When the system matrix is either a monotone matrix or an H-matrix, and the nonlinear mapping is a P-bounded mapping, we thoroughly set up the global convergence theories of these new methods. Moreover, under the assumptions that the system matrix is monotone and the nonlinear mapping is isotone, we discuss the monotone convergence properties of the new two-tage iteration methods, and investigate the influence of the matrix splittings as well as the relaxation parameters on the convergence behaviours of these methods. Numerical computations show that our new methods are feasible and efficient for solving the system of weakly nonlinear equations.
78 citations
TL;DR: A method is described, DDVERK, which implements this approach and justifies the strategies and heuristics that have been adopted and shows how the assumptions related to error control, stepsize control, and discontinuity detection can be efficiently realized for a particular sixth-order numerical method.
Abstract: We have recently developed a generic approach for solving neutral delay differential equations based on the use of a continuous Runge–Kutta formula with defect control and investigated its convergence properties. In this paper, we describe a method, DDVERK, which implements this approach and justify the strategies and heuristics that have been adopted. In particular we show how the assumptions related to error control, stepsize control, and discontinuity detection (required for convergence) can be efficiently realized for a particular sixth-order numerical method. Summaries of extensive testing are also reported.
74 citations
TL;DR: A further improvement of the above method is proposed, based on a point-wise evaluation/interpolation at a suitable set of Fourier points, of the functional relations defining each step of cyclic reduction (Bini and Meini, 1996).
Abstract: The cyclic reduction technique (Buzbee et al., 1970), rephrased in functional form (Bini and Meini, 1996), provides a numerically stable, quadratically convergent method for solving the matrix equation X = ∑+ ∞
i=0 Xi Ai, where the Ai's are nonnegative k × k matrices such that ∑+ ∞
i=0 Ai is column stochastic. In this paper we propose a further improvement of the above method, based on a point-wise evaluation/interpolation at a suitable set of Fourier points, of the functional relations defining each step of cyclic reduction (Bini and Meini,1996). This new technique allows us to devise an algorithm based on FFT having a lower computational cost and a higher numerical stability. Numerical results and comparisons are provided.
59 citations
TL;DR: For nonlinear problems asynchronous multisplitting algorithms including both the basic situation of O'Leary and White and the discrete analogue of Schwarz's alternating method and its multisubdomain extensions and moreover their two-stage counterparts are presented.
Abstract: Our aim is to present for nonlinear problems asynchronous multisplitting algorithms including both the basic situation of O'Leary and White and the discrete analogue of Schwarz's alternating method and its multisubdomain extensions and moreover their two-stage counterparts The analysis of these methods is based on El Tarazi’s convergence theorem for asynchronous iterations and leads to a good level of asynchronism in each of the considered situations
56 citations
TL;DR: Methods for stabilizing such an invariant manifold defined explicitly by algebraic equations are reviewed and theoretical convergence results for these methods are summarized.
Abstract: Many problems of practical interest can be modeled by differential systems where the solution lies on an invariant manifold defined explicitly by algebraic equations. In computer simulations, it is often important to take into account the invariant's information, either in order to improve upon the stability of the discretization (which is especially important in cases of long time integration) or because a more precise conservation of the invariant is needed for the given application. In this paper we review and discuss methods for stabilizing such an invariant. Particular attention is paid to post-stabilization techniques, where the stabilization steps are applied to the discretized differential system. We summarize theoretical convergence results for these methods and describe the application of this technique to multibody systems with holonomic constraints. We then briefly consider collocation methods which automatically satisfy certain, relatively simple invariants. Finally, we consider an example of a very long time integration and the effect of the loss of symplecticity and time-reversibility by the stabilization techniques.
TL;DR: Algorithms for computing updatable rank-revealing UTV decompositions that are efficient whenever the numerical rank of the matrix is much less than its dimensions are presented.
Abstract: A UTV decomposition of an m × n matrix is a product of an orthogonal matrix, a middle triangular matrix, and another orthogonal matrix In this paper we present and analyze algorithms for computing updatable rank-revealing UTV decompositions that are efficient whenever the numerical rank of the matrix is much less than its dimensions
TL;DR: A Detection Algorithm for the localisation of unknown fault lines of a surface from scattered data using thin plate splines is given, and it is shown that this yields approximation of second order accuracy instead of first order as in the global case.
Abstract: A Detection Algorithm for the localisation of unknown fault lines of a surface from scattered data is given. The method is based on a local approximation scheme using thin plate splines, and we show that this yields approximation of second order accuracy instead of first order as in the global case. Furthermore, the Detection Algorithm works with triangulation methods, and we show their utility for the approximation of the fault lines. The output of our method provides polygonal curves which can be used for the purpose of constrained surface approximation.
TL;DR: An adaptive algorithm where the approximations of the extreme eigenvalues that are needed to obtain upper bounds are computed when running CG leading to an improvement of the upper bounds for the norm of the error.
Abstract: In this paper we consider computing estimates of the norm of the error in the conjugate gradient (CG) algorithm. Formulas were given in a paper by Golub and Meurant (1997). Here, we first prove that these expressions are indeed upper and lower bounds for the A-norm of the error. Moreover, starting from these formulas, we investigate the computation of the l
2-norm of the error. Finally, we define an adaptive algorithm where the approximations of the extreme eigenvalues that are needed to obtain upper bounds are computed when running CG leading to an improvement of the upper bounds for the norm of the error. Numerical experiments show the effectiveness of this algorithm.
TL;DR: The computer-assisted study of two-dimensional (un)stable manifolds of steady states and saddle-type limit cycles for flows in Rn is illustrated and a number of computational issues arising are highlighted.
Abstract: We illustrate and discuss the computer-assisted study (approximation and visualization) of two-dimensional (un)stable manifolds of steady states and saddle-type limit cycles for flows in Rn Our investigation highlights a number of computational issues arising in this task, along with our solutions and “quick-fixes” for some of these problems Two examples illustrative of both successes and shortcomings of our current approach are presented Representative “snapshots” demonstrate the dependence of two-dimensional invariant manifolds on a bifurcation parameter as well as their interactions Such approximation and visualization studies are a necessary component of the computer-assisted study and understanding of global bifurcations
TL;DR: A new representation for diagonally implicit multistage integration methods (DIMSIMs) is derived in which the vector of external stages directly approximates the Nordsieck vector, which is zero-stable for any choice of variable mesh.
Abstract: A new representation for diagonally implicit multistage integration methods (DIMSIMs) is derived in which the vector of external stages directly approximates the Nordsieck vector. The methods in this formulation are zero-stable for any choice of variable mesh. They are also easy to implement since changing step-size corresponds to a simple rescaling of the vector of external approximations. The paper contains an analysis of local truncation error and of error accumulation in a variable step-size situation.
TL;DR: A careful analysis of the algorithm allows a reformulation which overcomes problems without losing efficiency, and a concise vectorized MATLAB-implementation is given.
Abstract: The most elegant way to evaluate box-splines is by using their recursive definition. However, a straightforward implementation reveals numerical difficulties. A careful analysis of the algorithm allows a reformulation which overcomes these problems without losing efficiency. A concise vectorized MATLAB-implementation is given.
TL;DR: A multisplitting two-stage AOR method is presented, which particularly uses the AOR-like iteration as inner iteration and is substantially a relaxed variant of the afore-presented method.
Abstract: The finite difference or the finite element discretizations of many differential or integral equations often result in a class of systems of weakly nonlinear equations. In this paper, by reasonably applying both the multisplitting and the two-stage iteration techniques, and in accordance with the special properties of this system of weakly nonlinear equations, we first propose a general multisplitting two-stage iteration method through the two-stage multiple splittings of the system matrix. Then, by applying the accelerated overrelaxation (AOR) technique of the linear iterative methods, we present a multisplitting two-stage AOR method, which particularly uses the AOR-like iteration as inner iteration and is substantially a relaxed variant of the afore-presented method. These two methods have a forceful parallel computing function and are much more suitable to the high-speed multiprocessor systems. For these two classes of methods, we establish their local convergence theories, and precisely estimate their asymptotic convergence factors under some suitable assumptions when the involved nonlinear mapping is only directionally differentiable. When the system matrix is either an H-matrix or a monotone matrix, and the nonlinear mapping is a P-bounded mapping, we thoroughly set up the global convergence theories of these new methods. Moreover, under the assumptions that the system matrix is monotone and the nonlinear mapping is isotone, we discuss the monotone convergence properties of the new multisplitting two-stage iteration methods, and investigate the influence of the multiple splittings as well as the relaxation parameters upon the convergence behaviours of these methods. Numerical computations show that our new methods are feasible and efficient for parallel solving of the system of weakly nonlinear equations.
TL;DR: In this article, two methods for computing the coefficients of the asymptotic series near the transition point are discussed, and auxiliary functions that can be computed more efficiently than the coefficients in the first method, and do not need the tabulation of many coefficients.
Abstract: Airy-type asymptotic representations of a class of special functions are considered from a numerical point of view. It is well known that the evaluation of the coefficients of the asymptotic series near the transition point is a difficult problem. We discuss two methods for computing the asymptotic series. One method is based on expanding the coefficients of the asymptotic series in Maclaurin series. In the second method we consider auxiliary functions that can be computed more efficiently than the coefficients in the first method, and we do not need the tabulation of many coefficients. The methods are quite general, but the paper concentrates on Bessel functions, in particular on the differential equation of the Bessel functions, which has a turning point character when order and argument of the Bessel functions are equal.
TL;DR: In this article, a successive continuation method for locating connecting orbits in parametrized systems of autonomous ODEs is considered, and a local convergence analysis is presented and several illustrative numerical examples are given.
Abstract: A successive continuation method for locating connecting orbits in parametrized systems of autonomous ODEs is considered. A local convergence analysis is presented and several illustrative numerical examples are given.
TL;DR: A general theorem for establishing the existence of a true periodic orbit near a numerically computed pseudoperiodic orbit of an autonomous system of ordinary differential equations is presented and a numerical method for estimating the Lyapunov exponents of thetrue periodic orbit is given.
Abstract: A general theorem for establishing the existence of a true periodic orbit near a numerically computed pseudoperiodic orbit of an autonomous system of ordinary differential equations is presented. For practical applications, a Newton method is devised to compute appropriate pseudoperiodic orbits. Then numerical considerations for checking the hypotheses of the theorem in terms of quantities which can be computed directly from the pseudoperiodic orbit and the vector field are addressed. Finally, a numerical method for estimating the Lyapunov exponents of the true periodic orbit is given. The theory and computations are designed to be applicable for unstable periodic orbits with long periods. The existence of several such periodic orbits of the Lorenz equations is exhibited.
TL;DR: Fourth order finite-difference algorithms for a semilinear singularly perturbed reaction–diffusion problem are discussed and compared and it is shown that the Bakhvalov mesh produces much better numerical results.
Abstract: Fourth order finite-difference algorithms for a semilinear singularly perturbed reaction–diffusion problem are discussed and compared both theoretically and numerically. One of them is the method of Sun and Stynes (1995) which uses a piecewise equidistant discretization mesh of Shishkin type. Another one is a simplified version of Vulanovic's method (1993), based on a discretization mesh of Bakhvalov type. It is shown that the Bakhvalov mesh produces much better numerical results.
TL;DR: The present algorithm is a generalization of the Fast Poisson Solver developed in the previous paper, based on the Fourier method in combination with a subtraction procedure.
Abstract: We describe high order numerical algorithms for the solution of second order elliptic equations in rectangular domains. These algorithms are based on the Fourier method in combination with a subtraction procedure. The singularities at the corner points, arising due to non-smoothness of the boundaries, are treated explicitly using properly constructed singular corner functions. The present algorithm is a generalization of the Fast Poisson Solver developed in our previous paper.
TL;DR: The algorithms proposed use a method based on selecting well-conditioned pairs of neighbouring polynomials, and the method is equivalent to going round the blocks instead of going across them, as is done in the well-known look-ahead methods.
Abstract: Two algorithms for the solution of a large sparse linear system of equations are proposed. The first is a modification of Lanczos' method and the second is based on one of Brezinski's methods. Both the latter methods are iterative and they can break down. In practical situations, serious numerical error is far more likely to occur because an ill-conditioned pair of polynomials is (implicitly) used in the calculation rather than complete breakdown arising because a large square block of exactly defective polynomials is encountered. The algorithms proposed use a method based on selecting well-conditioned pairs of neighbouring polynomials (in the associated Pade table), and the method is equivalent to going round the blocks instead of going across them, as is done in the well-known look-ahead methods.
TL;DR: For a mathematical model of an error per unit step or error per step adaptive Runge–Kutta algorithm, it may be shown that in a certain Probabilistic sense, with respect to a measure on the space of initial data, the small time-step heuristics are valid with probability one, leading to a probabilistic convergence result for the global error as τ→0.
Abstract: The numerical solution of initial value problems for ordinary differential equations is frequently performed by means of adaptive algorithms with user-input tolerance τ. The time-step is then chosen according to an estimate, based on small time-step heuristics, designed to try and ensure that an approximation to the local error commited is bounded by τ. A question of natural interest is to determine how the global error behaves with respect to the tolerance τ. This has obvious practical interest and also leads to an interesting problem in mathematical analysis. The primary difficulties arising in the analysis are that: (i) the time-step selection mechanisms used in practice are discontinuous as functions of the specified data; (ii) the small time-step heuristics underlying the control of the local error can break down in some cases. In this paper an analysis is presented which incorporates these two difficulties.
For a mathematical model of an error per unit step or error per step adaptive Runge–Kutta algorithm, it may be shown that in a certain probabilistic sense, with respect to a measure on the space of initial data, the small time-step heuristics are valid with probability one, leading to a probabilistic convergence result for the global error as τ → 0. The probabilistic approach is only valid in dimension m > 1; this observation is consistent with recent analysis concerning the existence of spurious steady solutions of software codes which highlights the difference between the cases m = 1 and m > 1. The breakdown of the small time-step heuristics can be circumvented by making minor modifications to the algorithm, leading to a deterministic convergence proof for the global error of such algorithms as τ → 0. An underlying theory is developed and the deterministic and probabilistic convergence results proved as particular applications of this theory.
TL;DR: The divergence-free finite element method (DFFEM) is a method to find an approximate solution of the Navier–Stokes equations in a divergence- free space and significantly reduces the dimension of the system to be solved at each time step.
Abstract: The divergence-free finite element method (DFFEM) is a method to find an approximate solution of the Navier–Stokes equations in a divergence-free space That is, the continuity equation is satisfied a priori DFFEM eliminates the pressure from the calculations and significantly reduces the dimension of the system to be solved at each time step For the standard 9-node velocity and 4-node pressure DFFEM, a basis for the weakly divergence-free subspace is constructed such that each basis function has nonzero support on at most 4 contiguous elements Given this basis, weakly divergence-free macroelements are constructed
TL;DR: Applying the method to spline conversion problems specifies new weights for knot removal and degree reduction, and the scheme is local and the approximation order is optimal, also constrained approximation problems can be solved efficiently.
Abstract: A new method for approximating functions by uniform B-splines is presented. It is based on the orthogonality relations for uniform B-splines in weighted Sobolev spaces, as introduced in (Reif, 1997). The scheme is local and the approximation order is optimal. Moreover, also constrained approximation problems can be solved efficiently; the size of the linear system to be solved is given by the number of constraints. Applying the method to spline conversion problems specifies new weights for knot removal and degree reduction.
TL;DR: A numerical method for evaluating the conjugacy at periodic and homoclinic symbol sequences in a systematic way is developed and it is shown that the resulting nonlinear systems are well conditioned uniformly with respect to the characteristic length of the symbol sequence.
Abstract: Transversal homoclinic orbits of maps are known to generate a Cantor set on which a power of the map conjugates to the Bernoulli shift on two symbols. This conjugacy may be regarded as a coding map, which for example assigns to a homoclinic symbol sequence a point in the Cantor set that lies on a homoclinic orbit of the map with a prescribed number of humps. In this paper we develop a numerical method for evaluating the conjugacy at periodic and homoclinic symbol sequences in a systematic way. The approach combines our previous method for computing the primary homoclinic orbit with the constructive proof of Smale's theorem given by Palmer. It is shown that the resulting nonlinear systems are well conditioned uniformly with respect to the characteristic length of the symbol sequence and that Newton's method converges uniformly too when started at a proper pseudo orbit. For the homoclinic symbol sequences an error analysis is given. The method works in arbitrary dimensions and it is illustrated by examples.
TL;DR: The geometric constraints and mission critical issues are discussed to give the numerical dynamical systems community some insight into the practical considerations and important problems of interest to the space mission designers.
Abstract: Orbits about the Sun–Earth L
1 and L
2 libration points have become very popular for space physics and astrophysics missions. Due to the hyperbolicity of the region of phase space where these orbits are found, they are difficult to compute and analyze. The invariant manifold structure provided by dynamical systems theory has been useful to compute transfer trajectories between orbits. These methods are very promising and require further development. This systematic approach is a great improvement from the difficult and labor-intensive numerical search methods currently popular in the astrodynamics community for studying these orbits. The geometric constraints and mission critical issues are discussed to give the numerical dynamical systems community some insight into the practical considerations and important problems of interest to the space mission designers. It is hoped that this communication will lead to more fruitful exchanges between the two communities.
TL;DR: This work examines the breakdown mechanism for a 2-torus associated to a system of coupled oscillators, and relies on computation of the Lyapunov-type numbers, as defined by Fenichel, for some general results on these numbers.
Abstract: In this work we examine the breakdown mechanism for a 2-torus associated to a system of coupled oscillators. To this end, we rely on computation of the Lyapunov-type numbers, as defined by Fenichel in N. Fenichel, Persistence and smoothness of invariant manifolds for flows, Indiana Univ. Math. J. 21 (1971) 193–226. We give some general results on these numbers, and some specific constructive results for the case in which there is phase-locking on the torus. This is the situation for the system of coupled oscillators we consider, and it leads to great simplifications in the computation of the Lyapunov-type numbers.