Numerical Taylor expansions of invariant manifolds in large dynamical systems

Summary. In this paper we develop a numerical method for computing higher order local approximations of invariant manifolds, such as stable, unstable or center manifolds near steady states of a dynamical system. The underlying system is assumed to be large in the sense that a large sparse Jacobian at the equilibrium occurs, for which only a linear (black box) solver and a low dimensional invariant subspace is available, but for which methods like the QR–Algorithm are considered to be too expensive. Our method is based on an analysis of the multilinear Sylvester equations for the higher derivatives which can be solved under certain nonresonance conditions. These conditions are weaker than the standard gap conditions on the spectrum which guarantee the existence of the invariant manifold. The final algorithm requires the solution of several large linear systems with a bordered Jacobian. To these systems we apply a block elimination method recently developed by Govaerts and Pryce [12, 14]. Received March 12, 1996 / Revised version reveiced August 8, 1997     

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     

Global conservative solutions for a model equation for shallow water waves of moderate amplitude

In this paper, we study the continuation of solutions to an equation for surface water waves of moderate amplitude in the shallow water regime beyond wave breaking (in [11], Constantin and Lannes proved that this equation accommodates wave breaking phenomena). Our approach is based on a method proposed by Bressan and Constantin [2]. By introducing a new set of independent and dependent variables, which resolve all singularities due to possible wave breaking, the evolution problem is rewritten as a semilinear system. Local existence of the semilinear system is obtained as fixed points of a contractive transformation. Moreover, this formulation allows one to continue the solution after collision time, giving a global conservative solution where the energy is conserved for almost all times. Finally, returning to the original variables, we obtain a semigroup of global conservative solutions, which depend continuously on the initial data.     

The strength of nonstationary iteration

This is the second of three papers in which we study global convergence of iterations using linear information for the solution of nonlinear equations. In Wasilkowski [6] we proved that for the class of all analytic scalar complex functions having only simple zeros there exists no globally convergentstationary iteration using linear information. Here we exhibit anonstationary iteration using linear information which is globally convergent even for the multivariate and abstract cases. This demonstrates the strength of nonstationary iteration. In Wasilkowski [7] we shall prove that any globally convergent iteration using linear information hasinfinite complexity even for the class of scalar complex polynomials having only simple zeros.     

On the convergence of splittings for semidefinite linear systems

Recently, Lee et al. [Young-ju Lee, Jinbiao Wu, Jinchao Xu, Ludmil Zikatanov, On the convergence of iterative methods for semidefinite linear systems, SIAM J. Matrix Anal. Appl. 28 (2006) 634-641] introduce new criteria for the semi-convergence of general iterative methods for semidefinite linear systems based on matrix splitting. The new conditions generalize the classical notion of P-regularity introduced by Keller [H.B. Keller, On the solution of singular and semidefinite linear systems by iterations, SIAM J. Numer. Anal. 2 (1965) 281-290]. In view of their results, we consider here stipulations on a splitting A=M-N, which lead to fixed point systems such that, the iterative scheme converges to a weighted Moore-Penrose solution to the system Ax=b. Our results extend the result of Lee et al. to a more general case and we also show that it requires less restrictions on the splittings than Keller’s P-regularity condition to ensure the convergence of iterative scheme.     

A multi-grid method with a priori and a posteriori level choice for the regularization of nonlinear ill-posed problems

Summary In this paper we study a multi-grid method for the numerical solution of nonlinear systems of equations arising from the discretization of ill-posed problems, where the special eigensystem structure of the underlying operator equation makes it necessary to use special smoothers. We provide uniform contraction factor estimates and show that a nested multigrid iteration together with an a priori or a posteriori chosen stopping index defines a regularization method for the ill-posed problem, i.e., a stable solution method, that converges to an exact solution of the underlying infinite-dimensional problem as the data noise level goes to zero, with optimal rates under additional regularity conditions. Supported by the Fonds zur F?rderung der wissenschaftlichen Forschung under grant T 7-TEC and project F1308 within Spezialforschungsbereich 13     

Collocation Methods and Their Modifications for Cauchy Singular Integral Equations on the Interval

In this paper we present polynomial collocation methods and their modi.cations for the numerical solution of Cauchy singular integral equations over the interval [-1, 1]. More precisely, the operators of the integral equations have the form with piecewise continuous coefficients a and b, and with a Jacobi weight . Using the splitting property of the singular values of the collocation methods, we obtain enough stable approximate methods to .nd the least square solution of our integral equation. Moreover, the modifications of the collocation methods enable us to compute kernel and cokernel dimensions of operators from a C*-algebra, which is generated by operators of the Cauchy singular integral equations.     

Computation and analysis for a constrained entropy optimization problem in finance

In [T. Coleman, C. He, Y. Li, Calibrating volatility function bounds for an uncertain volatility model, Journal of Computational Finance (2006) (submitted for publication)], an entropy minimization formulation has been proposed to calibrate an uncertain volatility option pricing model (UVM) from market bid and ask prices. To avoid potential infeasibility due to numerical error, a quadratic penalty function approach is applied. In this paper, we show that the solution to the quadratic penalty problem can be obtained by minimizing an objective function which can be evaluated via solving a Hamilton–Jacobian–Bellman (HJB) equation. We prove that the implicit finite difference solution of this HJB equation converges to its viscosity solution. In addition, we provide computational examples illustrating accuracy of calibration.     

Stability of backward Euler multirate methods and convergence of waveform relaxation

For a large class of traditional backward Euler multirate methods we show that stability is preserved when the methods are applied to certain stable (but not necessarily monotonic) non-linear systems. Methods which utilize waveform relaxation sweeps are shown to be stable and converge for certain monotonic systems.     

Kernel B-splines on general lattices

The aim of the paper is to provide a computationally effective way to construct stable bases on general non-degenerate lattices. In particular, we define new stable bases on hexagonal lattices and we give some numerical examples which show their usefulness in applications.     

A posteriori error estimators for nonconforming finite element methods of the linear elasticity problem

In this work we derive and analyze a posteriori error estimators for low-order nonconforming finite element methods of the linear elasticity problem on both triangular and quadrilateral meshes, with hanging nodes allowed for local mesh refinement. First, it is shown that equilibrated Neumann data on interelement boundaries are simply given by the local weak residuals of the numerical solution. The first error estimator is then obtained by applying the equilibrated residual method with this set of Neumann data. From this implicit estimator we also derive two explicit error estimators, one of which is similar to the one proposed by Dörfler and Ainsworth (2005) [24] for the Stokes problem. It is established that all these error estimators are reliable and efficient in a robust way with respect to the Lamé constants. The main advantage of our error estimators is that they yield guaranteed, i.e., constant-free upper bounds for the energy-like error (up to higher order terms due to data oscillation) when a good estimate for the inf-sup constant is available, which is confirmed by some numerical results.     

Semi-analytic time-marching algorithms for semi-linear parabolic equations

New time marching algorithms for numerical solution of semi-linear parabolic equations are described. They are based on the approximation method proposed by the first author. An important feature of the algorithms is that they are both explicit and stable under mild restrictions to the time step, which come from the non-linear part of the equation.     

Iterative refinement enhances the stability ofQR factorization methods for solving linear equations

Iterative refinement is a well-known technique for improving the quality of an approximate solution to a linear system. In the traditional usage residuals are computed in extended precision, but more recent work has shown that fixed precision is sufficient to yield benefits for stability. We extend existing results to show that fixed precision iterative refinement renders anarbitrary linear equations solver backward stable in a strong, componentwise sense, under suitable assumptions. Two particular applications involving theQR factorization are discussed in detail: solution of square linear systems and solution of least squares problems. In the former case we show that one step of iterative refinement suffices to produce a small componentwise relative backward error. Our results are weaker for the least squares problem, but again we find that iterative refinement improves a componentwise measure of backward stability. In particular, iterative refinement mitigates the effect of poor row scaling of the coefficient matrix, and so provides an alternative to the use of row interchanges in the HouseholderQR factorization. A further application of the results is described to fast methods for solving Vandermonde-like systems.     

Multilevel methods for mixed finite elements in three dimensions

In this paper we consider second order scalar elliptic boundary value problems posed over three–dimensional domains and their discretization by means of mixed Raviart–Thomas finite elements [18]. This leads to saddle point problems featuring a discrete flux vector field as additional unknown. Following Ewing and Wang [26], the proposed solution procedure is based on splitting the flux into divergence free components and a remainder. It leads to a variational problem involving solenoidal Raviart–Thomas vector fields. A fast iterative solution method for this problem is presented. It exploits the representation of divergence free vector fields as s of the –conforming finite element functions introduced by Nédélec [43]. We show that a nodal multilevel splitting of these finite element spaces gives rise to an optimal preconditioner for the solenoidal variational problem: Duality techniques in quotient spaces and modern algebraic multigrid theory [50, 10, 31] are the main tools for the proof. Received November 4, 1996 / Revised version received February 2, 1998     

Spline collocation methods for linear multi-term fractional differential equations

Some regularity properties of the solution of linear multi-term fractional differential equations are derived. Based on these properties, the numerical solution of such equations by piecewise polynomial collocation methods is discussed. The results obtained in this paper extend the results of Pedas and Tamme (2011) [15] where we have assumed that in the fractional differential equation the order of the highest derivative of the unknown function is an integer. In the present paper, we study the attainable order of convergence of spline collocation methods for solving general linear fractional differential equations using Caputo form of the fractional derivatives and show how the convergence rate depends on the choice of the grid and collocation points. Theoretical results are verified by some numerical examples.     

Boundary element preconditioners for a hypersingular integral equation on an interval

We propose an almost optimal preconditioner for the iterative solution of the Galerkin equations arising from a hypersingular integral equation on an interval. This preconditioning technique, which is based on the single layer potential, was already studied for closed curves [11,14]. For a boundary element trial space, we show that the condition number is of order (1 + | log h _min|)², where h _min is the length of the smallest element. The proof requires only a mild assumption on the mesh, easily satisfied by adaptive refinement algorithms.     

An Improved Theta-method for Systems of Ordinary Differential Equations

The / -method of order 1 or 2 (if / =1/2) is often used for the numerical solution of systems of ordinary differential equations. In the particular case of linear constant coefficient stiff systems the constraint 1/2 h / h 1, which excludes the explicit forward Euler method, is essential for the method to be A -stable. Moreover, unless / =1/2, this method is not elementary stable in the sense that its fixed-points do not display the linear stability properties of the fixed-points of the involved differential equation. We design a non-standard version of the / -method of the same order. We prove a result on the elementary stability of the new method, irrespective of the value of the parameter / ] [0,1]. Some absolute elementary stability properties pertinent to stiffness are discussed.     

Algorithms for the computation of functionals defined on the solution of a discrete ill-posed problem

We study a linear, discrete ill-posed problem, by which we mean a very ill-conditioned linear least squares problem. In particular we consider the case when one is primarily interested in computing a functional defined on the solution rather than the solution itself. In order to alleviate the ill-conditioning we require the norm of the solution to be smaller than a given constant. Thus we are lead to minimizing a linear functional subject to two quadratic constraints. We study existence and uniqueness for this problem and show that it is essentially equivalent to a least squares problem with a linear and a quadratic constraint, which is easier to handle computationally. Efficient algorithms are suggested for this problem.     

The life-span of backward error analysis for numerical integrators

Summary. Backward error analysis is a useful tool for the study of numerical approximations to ordinary differential equations. The numerical solution is formally interpreted as the exact solution of a perturbed differential equation, given as a formal and usually divergent series in powers of the step size. For a rigorous analysis, this series has to be truncated. In this article we study the influence of this truncation to the difference between the numerical solution and the exact solution of the perturbed differential equation. Results on the long-time behaviour of numerical solutions are obtained in this way. We present applications to the numerical phase portrait near hyperbolic equilibrium points, to asymptotically stable periodic orbits and Hopf bifurcation, and to energy conservation and approximation of invariant tori in Hamiltonian systems. Received October 18, 1995 / Revised version received February 28, 1996     

Discontinuous Least-Squares finite element method for the Div-Curl problem

In this paper, we consider the div-curl problem posed on nonconvex polyhedral domains. We propose a least-squares method based on discontinuous elements with normal and tangential continuity across interior faces, as well as boundary conditions, weakly enforced through a properly designed least-squares functional. Discontinuous elements make it possible to take advantage of regularity of given data (divergence and curl of the solution) and obtain convergence also on nonconvex domains. In general, this is not possible in the least-squares method with standard continuous elements. We show that our method is stable, derive a priori error estimates, and present numerical examples illustrating the method.