Showing papers in "Ima Journal of Numerical Analysis in 2003"
TL;DR: A new method is introduced for large-scale convex constrained optimization and is a generalization of the Spectral Projected Gradient method (SPG), but can be used when projections are difficult to compute.
Abstract: A new method is introduced for large-scale convex constrained optimization. The general model algorithm involves, at each iteration, the approximate minimization of a convex quadratic on the feasible set of the original problem and global convergence is obtained by means of nonmonotone line searches. A specific algorithm, the Inexact Spectral Projected Gradient method (ISPG), is implemented using inexact projections computed by Dykstra's alternating projection method and generates interior iterates. The ISPG method is a generalization of the Spectral Projected Gradient method (SPG), but can be used when projections are difficult to compute. Numerical results for constrained least-squares rectangular matrix problems are presented.
169 citations
TL;DR: In this paper, the authors prove the convergence of a particular iterative scheme for one factor uncertain volatility models and demonstrate how non-monotone discretization schemes (such as standard Crank-Nicolson timestepping) can converge to incorrect solutions.
Abstract: The pricing equations derived from uncertain volatility models in finance are often cast in the form of nonlinear partial differential equations. Implicit timestepping leads to a set of nonlinear algebraic equations which must be solved at each timestep. To solve these equations, an iterative approach is employed. In this paper, we prove the convergence of a particular iterative scheme for one factor uncertain volatility models. We also demonstrate how non-monotone discretization schemes (such as standard Crank-Nicolson timestepping) can converge to incorrect solutions, or lead to instability. Numerical examples are provided.
157 citations
TL;DR: In this paper, a coupled system of two singularly perturbed linear reaction-diffusion two-point boundary value problems is examined, where the leading term is multiplied by a small positive parameter, but these parameters may have different magnitudes.
Abstract: A coupled system of two singularly perturbed linear reaction-diffusion two-point boundary value problems is examined. The leading term of each equation is multiplied by a small positive parameter, but these parameters may have different magnitudes. The solutions to the system have boundary layers that overlap and interact. The structure of these layers is analysed, and this leads to the construction of a piecewise-uniform mesh that is a variant of the usual Shishkin mesh. On this mesh central differencing is proved to be almost first-order accurate, uniformly in both small parameters. Supporting numerical results are presented for a test problem.
130 citations
TL;DR: In this paper, the authors considered discretization in time of an inhomogeneous parabolic equation in a Banach space setting, using a representation of the solution as an integral along a smooth curve in the complex left half-plane which, after transformation to a finite interval, is then evaluated to high accuracy by a quadrature rule.
Abstract: We consider the discretization in time of an inhomogeneous parabolic equation in a Banach space setting, using a representation of the solution as an integral along a smooth curve in the complex left half-plane which, after transformation to a finite interval, is then evaluated to high accuracy by a quadrature rule. This reduces the problem to a finite set of elliptic equations with complex coefficients, which may be solved in parallel. The paper is a further development of earlier work by the authors, where we treated the homogeneous equation in a Hilbert space framework. Special attention is given here to the treatment of the forcing term. The method is combined with finite-element discretization in spatial variables.
126 citations
TL;DR: In this article, the authors investigated the order of convergence of a general class of two-step hybrid methods for systems of differential equations of the special form y = f(x, y).
Abstract: The theory of B-series is used to investigate the order of convergence of a general class of two-step hybrid methods for systems of differential equations of the special form y = f(x, y). The main result is a set of order conditions, analogous to those for Runge-Kutta methods, offering an alternative to the customary ad hoc Taylor expansions. A byproduct is a remarkably simple formula from which the order of dispersion of such methods is easily determined. Conditions under which the two-step methods are symmetric are established, and particular examples are considered.
117 citations
TL;DR: In this paper, the authors generalized Smale's α theory to the context of intrinsic Newton iteration on geodesically complete analytic Riemannian and Hermitian manifolds.
Abstract: In this paper, Smale's α theory is generalized to the context of intrinsic Newton iteration on geodesically complete analytic Riemannian and Hermitian manifolds. Results are valid for analytic mappings from a manifold to a linear space of the same dimension, or for analytic vector fields on the manifold. The invariant γ is defined by means of high-order covariant derivatives. Bounds on the size of the basin of quadratic convergence are given. If the ambient manifold has negative sectional curvature, those bounds depend on the curvature. A criterion of quadratic convergence for Newton iteration from the information available at a point is also given.
109 citations
TL;DR: In this article, a modified steepest descent method, whose stepsizes alternately minimize the function value and the gradient norm along the line of gradient descent was proposed, which is known as the alternate minimization (AM) gradient method.
Abstract: It is well known that the minimization of a smooth function f(x) is equivalent to minimizing its gradient norm ‖g(x)‖
2 in some sense. In this paper, we propose a modified steepest descent method, whose stepsizes alternately minimize the function value and the gradient norm along the line of steepest descent. Hence the name ‘alternate minimization (AM) gradient method’. For strictly convex quadratics, the AM method is proved to be Q‐superlinearly convergent in two dimensions, and Q‐linearly convergent in any dimension. Numerical experiments are presented for symmetric and positive definite linear systems. They suggest that the AM method is much better than the classical steepest descent (SD) method and comparable with some existing gradient methods. They also show that the AM method is an efficient alternative if a solution with a low precision is required. Two variants of the AM method, named shortened SD step gradient methods, are also presented and analysed in this paper. By designing a new kind of line search, the two variants are extended to the field of unconstrained optimization.
107 citations
TL;DR: A restrictively preconditioned conjugate gradient method is presented for solving a large sparse system of linear equations and is applied to solve the linear systems of block two-by-two structures.
Abstract: A restrictively preconditioned conjugate gradient method is presented for solving a large sparse system of linear equations. This new method originates from the classical conjugate gradient method and its restrictively preconditioned variant, and covers many standard Krylov subspace iteration methods such as the conjugate gradient, conjugate residual, CGNR, CGNE and the corresponding preconditioned variants. In particular, the new method is applied to solve the linear systems of block two-by-two structures, and consequently, special but quite efficient preconditioned conjugate gradient-like methods are obtained for these problems. Numerical computations show that the new methods are more effective than some standard preconditioned Krylov subspace methods like the preconditioned conjugate gradient and the preconditioned MINRES.
100 citations
TL;DR: An algorithm for computing the ‘pseudospectral abscissa’, which is the largest real part of such an eigenvalue, measures the robust stability of A, and proves global and local quadratic convergence.
Abstract: x = Ax is robustly stable when all eigenvalues of complex matrices within a given distance of the square matrix A lie in the left half-plane. The ‘pseudospectral abscissa’, which is the largest real part of such an eigenvalue, measures the robust stability of A .W epresent an algorithm for computing the pseudospectral abscissa, prove global and local quadratic convergence, and discuss numerical implementation. As with analogous methods for calculating H∞ norms, our algorithm depends on computing the eigenvalues of associated Hamiltonian matrices.
81 citations
TL;DR: It is shown that, in fact, there exists an elegant proof of this feature independent of the space dimension, and superconvergence for dimensions four and up is proved simultaneously.
Abstract: Superconvergence of the gradient for the linear simplicial finite-element method applied to elliptic equations is a well known feature in one, two, and three space dimensions. In this paper we show that, in fact, there exists an elegant proof of this feature independent of the space dimension. As a result, superconvergence for dimensions four and up is proved simultaneously. The key ingredient will be that we embed the gradients of the continuous piecewise linear functions into a larger space for which we describe an orthonormal basis having some useful symmetry properties. Since gradients and rotations of standard finite-element functions are in fact the rotation-free and divergence-free elements of Raviart-Thomas and Nedelec spaces in three dimensions, we expect our results to have applications also in those contexts.
74 citations
TL;DR: It is shown that the proposed fully discrete stabilized finite-element method results in the optimal order error bounds for the velocity and the pressure.
Abstract: A fully discrete stabilized finite-element method is presented for the two-dimensional time-dependent Navier-Stokes problem. The spatial discretization is based on a finite-element space pair (X h , M h ) for the approximation of the velocity and the pressure, constructed by using the Q 1 - P 0 quadrilateral element or the P 1 - P 0 triangular element; the time discretization is based on the Euler semi-implicit scheme. It is shown that the proposed fully discrete stabilized finite-element method results in the optimal order error bounds for the velocity and the pressure.
TL;DR: In this paper, mean-square and weak quasi-symplectic methods for Langevin type equations with separable Hamiltonians and additive noise were constructed. And the methods derived are based on symplectic schemes for stochastic Hamiltonian systems.
Abstract: Langevin type equations are an important and fairly large class of systems close to Hamiltonian ones. The constructed mean-square and weak quasi-symplectic methods for such systems degenerate to symplectic methods when a system degenerates to a stochastic Hamiltonian one. In addition, quasi-symplectic methods' law of phase volume contractivity is close to the exact law. The methods derived are based on symplectic schemes for stochastic Hamiltonian systems. Mean-square symplectic methods were obtained in Milstein et al. (2002, SIAM J. Numer. Anal., 39, 2066-2088; 2003, SIAM J. Numer. Anal., 40, 1583-1604) while symplectic methods in the weak sense are constructed in this paper. Special attention is paid to Hamiltonian systems with separable Hamiltonians and with additive noise. Some numerical tests of both symplectic and quasi-symplectic methods are presented. They demonstrate superiority of the proposed methods in comparison with standard ones.
TL;DR: A Strang error estimate containing the consistency terms arising from the approximation of the continuous operators involved is deduced and applied to analyse a fully discrete Galerkin scheme for a twofold saddle point formulation of a nonlinear elliptic boundary value problem in divergence form.
Abstract: We provide a general abstract theory for the solvability and Galerkin approximation of nonlinear twofold saddle point problems. In particular, a Strang error estimate containing the consistency terms arising from the approximation of the continuous operators involved is deduced. Then we apply these results to analyse a fully discrete Galerkin scheme for a twofold saddle point formulation of a nonlinear elliptic boundary value problem in divergence form. Some numerical results are also presented.
TL;DR: The paper deals with a Nitsche-type finite-element method for treating non-matching meshes at the interface of some domain decomposition for transmission problems of the plane with Dirichlet boundary conditions entailing singularities at the corners or endpoints of the polygonal interface.
Abstract: The paper deals with a Nitsche-type finite-element method for treating non-matching meshes at the interface of some domain decomposition. The method is extended to transmission (or interface) problems of the plane with Dirichlet boundary conditions entailing singularities at the corners or endpoints of the polygonal interface. In a natural way, the interface of the transmission problem is taken as the interface of the domain decomposition and of the non-matching meshes. Properties of the finite-element scheme and error estimates are proved. For appropriate local mesh refinement, optimal convergence rates as known for the classical finite-element method and regular solution are derived. Some numerical tests illustrate the approach and confirm the theoretical results.
TL;DR: In this paper, the existence and uniqueness of M-matrix solutions of the quadratic matrix equation X 2 - EX - F = 0 were discussed and the Schur method was used to find the unique stabilizing or almost stabilizing solution.
Abstract: We study the quadratic matrix equation X 2 - EX - F = 0, where E is diagonal and F is an M-matrix. Quadratic matrix equations of this type arise in noisy Wiener-Hopf problems for Markov chains. The solution of practical interest is a particular M-matrix solution. The existence and uniqueness of M-matrix solutions and numerical methods for finding the desired M-matrix solution are discussed by transforming the equation into an equation that belongs to a special class of non-symmetric algebraic Riccati equations (AREs). We also discuss the general non-symmetric ARE and describe how we can use the Schur method to find the stabilizing or almost stabilizing solution if it exists. The desired M-matrix solution of the quadratic matrix equation (a special non-symmetric ARE by itself) turns out to be the unique stabilizing or almost stabilizing solution.
TL;DR: In this article, the stability properties of numerical methods for this kind of equation have been studied by numerous authors who have mainly considered the constant coefficient case and integration over meshes with constant stepsize.
Abstract: where a(t) and b(t) are continuous complex-valued functions and q ∈ (0, 1) is a fixed constant. We also extend our analysis to the neutral equation y � (t) = a(t) y(t) + b(t) y( qt ) + c(t) y � ( qt ), where c(t) is also a continuous complex-valued function. In recent years, stability properties of numerical methods for this kind of equation have been studied by numerous authors who have mainly considered the constant coefficient case and integration over meshes with constant stepsize. In general, the developed techniques give rise to non-standard recurrence relations. Instead, in the present paper we study constrained variable stepsize schemes. In Bellen et al. (1997, Appl. Numer. Math.,24, 279–293) the behaviour of the class of Θ-methods was analysed with relevance for the constant coefficient version of the above equation. As with the constant coefficient case, we show here that the methods in this class are stable if and only if 1 Θ > 1/2. The basic tool we use in our analysis is the spectral radius of a family of matrices, computed by means of polytope extremal norms. An extension of the stability results to a more general class of equations is also discussed in the last section.
TL;DR: In this paper, the shape-preservation property is secured by adjusting 'tension' parameters that arise upon relaxing parametric continuity to geometric continuity, and a simpler and cheaper alternative is also introduced.
Abstract: The interpolation of a planar sequence of points p 0 ,..., p N by shape-preserving G 1 or G 2 PH quintic splines with specified end conditions is considered. The shape-preservation property is secured by adjusting 'tension' parameters that arise upon relaxing parametric continuity to geometric continuity. In the G 2 case, the PH spline construction is based on applying Newton-Raphson iterations to a global system of equations, commencing with a suitable initialization strategy-this generalizes the construction described previously in Numerical Algorithms 27, 35-60 (2001). As a simpler and cheaper alternative, a shape-preserving G 1 PH quintic spline scheme is also introduced. Although the order of continuity is lower, this has the advantage of allowing construction through purely local equations.
TL;DR: In this article, lower and upper bounds for the optimal p-norm condition number achievable by two-sided diagonal scalings are given, and a class of new lower bounds for componentwise distance to the nearest singular matrix is given.
Abstract: In this note we give lower and upper bounds for the optimal p-norm condition number achievable by two-sided diagonal scalings. There are no assumptions on the irreducibility of certain matrices. The bounds are shown to be optimal for the 2-norm. For the 1-norm and inf-norm the (known) exact value of the optimal condition number is conflrmed. We give means how to calculate the minimizing diagonal matrices. Furthermore, a class of new lower bounds for the componentwise distance to the nearest singular matrix is given. They are shown to be superior to existing ones.
TL;DR: In this article, the authors constructed an approximate solution to one-dimensional nonlinear drift-diffusion equations by a finite-volume scheme, which yields the global existence of solutions to the equations.
Abstract: We construct an approximate solution to one-dimensional nonlinear drift-diffusion equations by a finite-volume scheme. The convergence of the scheme is proved, which yields the global existence of solutions to the equations. This result is valid even in the presence of vacuum state: the densities of free carriers of charge vanish for some times and positions, which implies that the continuity equations are degenerate parabolic equations. Finally, numerical simulations are performed by the finite-volume scheme.
TL;DR: In this paper, the authors consider the problem of constructing spatial finite-difference approximations on an arbitrary fixed grid which preserve any number of integrals of the partial differential equation and preserve some of its symmetries.
Abstract: We consider the problem of constructing spatial finite-difference approximations on an arbitrary fixed grid which preserve any number of integrals of the partial differential equation and preserve some of its symmetries. A basis for the space of such finite-difference operators is constructed; most cases of interest involve a single such basis element. (The 'Arakawa' Jacobian is such an element, as are discretizations satisfying 'summation by parts' identities.) We show how the grid, its symmetries, and the differential operator interact to affect the complexity of the finite difference.
TL;DR: Numerical results indicate that the AL e gendre Galerkin–Chebyshev collocation method for Burgers-like equations is as stable and accurate as the standard Legendre collocations method, and as efficient and easy to implement as thestandard ChebysheV collocation methods.
Abstract: AL e gendre Galerkin–Chebyshev collocation method for Burgers-like equations is developed. This method is based on the Legendre–Galerkin variational form, but the nonlinear term and the right-hand term are treated by Chebyshev–Gauss interpolation. Error estimates of the semi-discrete scheme and the fully discrete scheme are given in the L 2 -norm. Numerical results indicate that our method is as stable and accurate as the standard Legendre collocation method, and as efficient and easy to implement as the standard Chebyshev collocation method.
TL;DR: In this article, the authors present adaptive procedures for the numerical study of positive solutions of the following problem: u t = u xx (x,t) ∈ (0,1) x [0,T], u x (t) = 0 t ∈ [ 0,T], u x 1,t = up(1,t), u x 0,0) = u 0 (x) x ∈ 0, 1, with p > 1.
Abstract: In this paper we present adaptive procedures for the numerical study of positive solutions of the following problem: u t = u xx (x,t) ∈ (0,1) x [0,T), u x (0,t) = 0 t ∈ [0,T], u x (1,t) = up(1,t) t ∈ [0, T), u(x,0) = u 0 (x) x ∈ (0, 1), with p > 1. We describe two methods. The first one refines the mesh in the region where the solution becomes bigger in a precise way that allows us to recover the blow-up rate and the blow-up set of the continuous problem. The second one combines the ideas used in the first one with moving mesh methods and moves the last points when necessary. This scheme also recovers the blow-up rate and set. Finally, we present numerical experiments to illustrate the behaviour of both methods.
TL;DR: It is shown that the PCG method is at least as accurate (fast) as the SOR-like method.
Abstract: Golub et al. (2001, BIT, 41, 71n85) gave a generalized successive over-relaxation method for the augmented systems. In this paper, the connection between the SOR-like method and the preconditioned conjugate gradient (PCG) method for the augmented systems is investigated. It is shown that the PCG method is at least as accurate (fast) as the SOR-like method. Numerical examples demonstrate that the PCG method is much faster than the SOR-like method.
TL;DR: In this paper, locally supported wavelets on manifolds Γ are given as the closure of a disjoint union of general smooth parametric images of an n-simplex, and the wavelets are proven to generate Riesz bases for Sobolev spaces H s (Γ) when s ∈ (-1, 3/2), if not limited by the global smoothness of F.
Abstract: We construct locally supported, continuous wavelets on manifolds Γ that are given as the closure of a disjoint union of general smooth parametric images of an n-simplex. The wavelets are proven to generate Riesz bases for Sobolev spaces H s (Γ) when s ∈ (-1, 3/2), if not limited by the global smoothness of F. These results generalize the findings from Dahmen & Stevenson (1999) SIAM J. Numer. Anal., 37, 319-352, where it was assumed that each parametrization has a constant Jacobian determinant. The wavelets can be arranged to satisfy the cancellation property of, in principle, any order, except for wavelets with supports that extend to different patches, which generally satisfy the cancellation property of only order 1.
TL;DR: In this paper, smoothing properties and approximation of time derivatives for time discretization schemes with constant time steps for a homogeneous parabolic problem formulated as an abstract initial value problem in a Banach space were studied.
Abstract: We study smoothing properties and approximation of time derivatives for time discretization schemes with constant time steps for a homogeneous parabolic problem formulated as an abstract initial-value problem in a Banach space. The time stepping schemes are based on using rational functions r(z) e–z which are A()-stable for suitable [0, /2] and satisfy |r()| < 1, and the approximations of time derivatives are based on using difference quotients in time. Both smooth and non-smooth data error estimates of optimal order for the approximation of time derivatives are proved. Further, we apply the results to obtain error estimates of time derivatives in the supremum norm for fully discrete methods based on discretizing the spatial variable by a finite-element method.
TL;DR: In this article, a new parameter-free alternating direction implicit Laplace-modified method was proposed, in which the standard central difference formula was used for the time approximation and orthogonal spline collocation is used for spatial discretization.
Abstract: On a rectangular region, we consider a linear second-order hyperbolic initial-boundary value problem involving a mixed derivative term, continuous variable coefficients and nonhomogeneous Dirichlet boundary conditions. In comparison to the alternating direction implicit Laplace-modified method of Fernandes (1997), we formulate and analyse a new parameter-free alternating direction implicit scheme in which the standard central difference formula is used for the time approximation and orthogonal spline collocation is used for the spatial discretization. We establish unconditional stability of the scheme, and its optimal order in the discrete maximum norm in time and the H 1 norm in space. Numerical experiments indicate that the new scheme, which has the same order as the method of Fernandes (1997, Numer. Math., 77, 223-241), is more accurate. We also show that the new scheme is easily generalized to the second-order hyperbolic problems on rectangular polygons. Extensions of the scheme to problems with discontinuous coefficients, nonlinear problems, and problems with other boundary conditions are also discussed.
TL;DR: In this article, the authors give a proof of convergence for the approximate solution of an elliptic-hyperbolic system, describing the conservation of two immiscible incompressible phases flowing in a porous medium.
Abstract: This paper gives a proof of convergence for the approximate solution of an elliptic-hyperbolic system, describing the conservation of two immiscible incompressible phases flowing in a porous medium. The approximate solution is obtained by a mixed finite-element method on a large class of meshes for the elliptic equation and a finite-volume method for the hyperbolic equation. Since the considered meshes are not necessarily structured, the proof uses a weak total variation inequality, which cannot yield a BV-estimate. We thus prove, under an L∞ estimate, the weak convergence of the finite-volume approximation. The strong convergence proof is then sketched under regularity assumptions which ensure that the flux is Lipschitz continuous.
TL;DR: In this article, a numerical solution of a λ-rational Sturm-Liouville problem using the shooting technique is presented, in which some classical procedures (boundary value methods and shooting with the Magnus method) are applied.
Abstract: This paper deals with the numerical solution of a λ-rational Sturm-Liouville problem using the shooting technique. Some classical procedures (boundary value methods and shooting with the Magnus method) are applied in order to examine their effectiveness in the presence of an eigenvalue embedded in the essential spectrum. It is shown that all these methods experience a degeneration of their order of accuracy: none, regardless of the order it may have on a classical Sturm-Liouville problem, has order greater than 2. Moreover, this situation cannot be rectified by applying the regularizing transformation of Niessen and Zettl. Nevertheless, it is possible to generalize a very efficient correction technique, first suggested by the authors for classical Sturm-Liouville problems, to the present case.
TL;DR: This paper designs high-order (non)local artificial boundary conditions (ABCs) which are different from those proposed by Han, H. & Bao, W. for incompressible materials, and presents error bounds for the finite-element approximation of the exterior Stokes equations in two dimensions.
Abstract: In this paper we design high-order (non)local artificial boundary conditions (ABCs) which are different from those proposed by Han, H. & Bao, W. (1997 Numer. Math., 77, 347-363) for incompressible materials, and present error bounds for the finite-element approximation of the exterior Stokes equations in two dimensions. The finite-element approximation (especially its corresponding stiff matrix) becomes much simpler (sparser) when it is formulated in a bounded computational domain using the new (non)local approximate ABCs. Our error bounds indicate how the errors of the finite-element approximations depend on the mesh size, terms used in the approximate ABCs and the location of the artificial boundary. Numerical examples of the exterior Stokes equations outside a circle in the plane are presented. Numerical results demonstrate the performance of our error bounds.
TL;DR: An a posteriori error estimate is presented that allows the estimation of the difference between continuous and semidiscrete solutions by quantities that can be calculated from the approximation and given data.
Abstract: A time discrete scheme is used to approximate the solution to a phase field system of Penrose-Fife type with a non-conserved order parameter. An a posteriori error estimate is presented that allows the estimation of the difference between continuous and semidiscrete solutions by quantities that can be calculated from the approximation and given data.