Showing papers in "Mathematics of Computation in 2008"

The completion of locally refined simplicial partitions created by bisection

Abstract: Recently, in [Found. Comput. Math., 7(2) (2007), 245-269], we proved that an adaptive finite element method based on newest vertex bisection in two space dimensions for solving elliptic equations, which is essentially the method from [SINUM, 38 (2000), 466-488] by Morin, Nochetto, and Siebert, converges with the optimal rate. The number of triangles N in the output partition of such a method is generally larger than the number M of triangles that in all intermediate partitions have been marked for bisection, because additional bisections are needed to retain conforming meshes. A key ingredient to our proof was a result from [Numer. Math., 97(2004), 219-268] by Binev, Dahmen and DeVore saying that N - N-0 <= CM for some absolute constant C, where N-0 is the number of triangles from the initial partition that have never been bisected. In this paper, we extend this result to bisection algorithms of n-simplices, with that generalizing the result concerning optimality of the adaptive finite element method to general space dimensions.

On splitting methods for Schrödinger-Poisson and cubic nonlinear Schrödinger equations

Abstract: We give an error analysis of Strang-type splitting integrators for nonlinear Schrodinger equations. For Schrodinger-Poisson equations with an H 4 -regular solution, a first-order error bound in the H 1 norm is shown and used to derive a second-order error bound in the L 2 norm. For the cubic Schrodinger equation with an H 4 -regular solution, first-order convergence in the H 2 norm is used to obtain second-order convergence in the L 2 norm. Basic tools in the error analysis are Lie-commutator bounds for estimating the local error and H m -conditional stability for error propagation, where m = 1 for the Schrodinger-Poisson system and m = 2 for the cubic Schrodinger equation.

A superconvergent ldg-hybridizable galerkin method for second-order elliptic problems

Abstract: We identify and study an LDG-hybridizable Galerkin method, which is not an LDG method, for second-order elliptic problems in several space dimensions with remarkable convergence properties. Unlike all other known discontinuous Galerkin methods using polynomials of degree k ≥ 0 for both the potential as well as the flux, the order of convergence in L 2 of both unknowns is k + 1. Moreover, both the approximate potential as well as its numerical trace superconverge in L 2 -like norms, to suitably chosen projections of the potential, with order k + 2. This allows the application of element-by-element postprocessing of the approximate solution which provides an approximation of the potential converging with order k+2 in L 2 . The method can be thought to be in between the hybridized version of the Raviart-Thomas and that of the Brezzi-Douglas-Marini mixed methods.

Finite elements for symmetric tensors in three dimensions

Abstract: We construct finite element subspaces of the space of symme- tric tensors with square-integrable divergence on a three-dimensional domain. These spaces can be used to approximate the stress field in the classical Hellinger-Reissner mixed formulation of the elasticty equations, when stan- dard discontinous finite element spaces are used to approximate the displace- ment field. These finite element spaces are defined with respect to an arbitrary simplicial triangulation of the domain, and there is one for each positive value of the polynomial degree used for the displacements. For each degree, these provide a stable finite element discretization. The construction of the spaces is closely tied to discretizations of the elasticity complex, and can be viewed as the three-dimensional analogue of the triangular element family for plane elasticity previously proposed by Arnold and Winther.

Weak order for the discretization of the stochastic heat equation

Abstract: In this paper we study the approximation of the distribution of $X_t$ Hilbert--valued stochastic process solution of a linear parabolic stochastic partial differential equation written in an abstract form as $$ dX_t+AX_t \, dt = Q^{1/2} d W_t, \quad X_0=x \in H, \quad t\in[0,T], $$ driven by a Gaussian space time noise whose covariance operator $Q$ is given. We assume that $A^{-\alpha}$ is a finite trace operator for some $\alpha>0$ and that $Q$ is bounded from $H$ into $D(A^\beta)$ for some $\beta\geq 0$. It is not required to be nuclear or to commute with $A$. The discretization is achieved thanks to finite element methods in space (parameter $h>0$) and implicit Euler schemes in time (parameter $\Delta t=T/N$). We define a discrete solution $X^n_h$ and for suitable functions $\varphi$ defined on $H$, we show that $$ |\E \, \varphi(X^N_h) - \E \, \varphi(X_T) | = O(h^{2\gamma} + \Delta t^\gamma) $$ oindent where $\gamma<1- \alpha + \beta$. Let us note that as in the finite dimensional case the rate of convergence is twice the one for pathwise approximations.

The Euler implicit/explicit scheme for the 2D time-dependent Navier-Stokes equations with smooth or non-smooth initial data

Abstract: This paper considers the stability and convergence results for the Euler implicit/explicit scheme applied to the spatially discretized two-dimensional (2D) time-dependent Navier-Stokes equations. A Galerkin finite element spatial discretization is assumed, and the temporal treatment is implicit/explict scheme, which is implicit for the linear terms and explicit for the nonlinear term. Here the stability condition depends on the smoothness of the initial data u 0 ∈ H α , i.e., the time step condition is τ < C 0 in the case of a = 2, τ| log h| ≤ C 0 in the case of a = 1 and Th -2 ≤ C 0 in the case of a = 0 for mesh size h and some positive constant C 0 . We provide the H 2 -stability of the scheme under the stability condition with a = 0,1, 2 and obtain the optimal H 1 - L 2 error estimate of the numerical velocity and the optimal L 2 error estimate of the numerical pressure under the stability condition with a = 1, 2.

Superconvergence of mixed finite element methods for optimal control problems

Abstract: T. In this paper, we investigate the superconvergence property of the numerical solution of a quadratic convex optimal control problem by using rectangular mixed finite element methods. The state and co-state variables are approximated by the lowest order Raviart-Thomas mixed finite element spaces and the control variable is approximated by piecewise constant functions. Some realistic regularity assumptions are presented and applied to error estimation by using an operator interpolation technique. We derive L 2 superconvergence properties for the flux functions along the Gauss lines and for the scalar functions at the Gauss points via mixed projections. Moreover, global L 2 superconvergence results are obtained by virtue of an interpolation postprocessing technique. Thus, based on these superconvergence estimates, some asymptotic exactness a posteriori error estimators are presented for the mixed finite element methods. Finally, some numerical examples are given to demonstrate the practical side of the theoretical results about superconvergence.

The complexity of class polynomial computation via floating point approximations

Abstract: We analyse the complexity of computing class polynomials, that are an important ingredient for CM constructions of elliptic curves, via complex floating point approximations of their roots. The heart of the algorithm is the evaluation of modular functions in several arguments. The fastest one of the presented approaches uses a technique devised by Dupont to evaluate modular functions by Newton iterations on an expression involving the arithmetic-geometric mean. It runs in time $O (|D| \log^5 |D| \log \log |D|) = O (|D|^{1 + \epsilon}) = O ( h^{2 + \epsilon})$ for any $\epsilon > 0$, where $D$ is the CM discriminant and $h$ is the degree of the class polynomial. Another fast algorithm uses multipoint evaluation techniques known from symbolic computation; its asymptotic complexity is worse by a factor of $\log |D|$. Up to logarithmic factors, this running time matches the size of the constructed polynomials. The estimate also relies on a new result concerning the complexity of enumerating the class group of an imaginary-quadratic order and on a rigorously proven upper bound for the height of class polynomials.

Filtering in Legendre spectral methods

Abstract: We discuss the impact of modal filtering in Legendre spectral methods, both on accuracy and stability. For the former, we derive sufficient conditions on the filter to recover high order accuracy away from points of discontinuity. Computational results confirm that less strict necessary conditions appear to be adequate. We proceed to discuss a instability mechanism in polynomial spectral methods and prove that filtering suffices to ensure stability. The results are illustrated by computational experiments.

A globally convergent BFGS method for nonlinear monotone equations without any merit functions

Abstract: r. Since 1965, there has been significant progress in the theoretical study on quasi-Newton methods for solving nonlinear equations, especially in the local convergence analysis. However, the study on global convergence of quasi-Newton methods is relatively fewer, especially for the BFGS method. To ensure global convergence, some merit function such as the squared norm merit function is typically used. In this paper, we propose an algorithm for solving nonlinear monotone equations, which combines the BFGS method and the hyperplane projection method. We also prove that the proposed BFGS method converges globally if the equation is monotone and Lipschitz continuous without differentiability requirement on the equation, which makes it possible to solve some nonsmooth equations. An attractive property of the proposed method is that its global convergence is independent of any merit function.We also report some numerical results to show efficiency of the proposed method.

Fast algorithms for computing isogenies between elliptic curves

Abstract: We survey algorithms for computing isogenies between elliptic curves defined over a field of characteristic either 0 or a large prime. We introduce a new algorithm that computes an isogeny of degree $\ell$ ($\ell$ different from the characteristic) in time quasi-linear with respect to $\ell$. This is based in particular on fast algorithms for power series expansion of the Weierstrass $\wp$-function and related functions.

The hyperdeterminant and triangulations of the 4-cube

Abstract: PETER HUGGINS, BERND STURMFELS,JOSEPHINE YU AND DEBBIE S. YUSTERAbstract. The hyperdeterminant of format 2×2×2×2 is a polynomialof degree 24 in 16 unknowns which has 2894276 terms. We compute theNewton polytope of this polynomial and the secondary polytope of the 4-cube. The 87959448 regular triangulations of the 4-cube are classified into25448 D-equivalence classes, one for each vertex of the Newton polytope.The 4-cube has 80876 coarsest regular subdivisions, one for each facet of thesecondary polytope, but only 268 of them come from the hyperdeterminant.

Local projection stabilization of equal order interpolation applied to the Stokes problem

Abstract: The local projection stabilization allows us to circumvent the Babuska-Brezzi condition and to use equal order interpolation for discretizing the Stokes problem. The projection is usually done in a two-level approach by projecting the pressure gradient onto a discontinuous finite element space living on a patch of elements. We propose a new local projection stabilization method based on (possibly) enriched finite element spaces and discontinuous projection spaces defined on the same mesh. Optimal order of convergence is shown for pairs of approximation and projection spaces satisfying a certain inf-sup condition. Examples are enriched simplicial finite elements and standard quadrilateral/hexahedral elements. The new approach overcomes the problem of an increasing discretization stencil and, thus, is simple to implement in existing computer codes. Numerical tests confirm the theoretical convergence results which are robust with respect to the user-chosen stabilization parameter.

A multiscale method for highly oscillatory ordinary differential equations with resonance

Abstract: A multiscale method for computing the effective behavior of a class of stiff and highly oscillatory ordinary differential equations (ODEs) is presented. The oscillations may be in resonance with on ...

Superconvergence of quadratic finite elements on mildly structured grids

Abstract: Superconvergence estimates are studied in this paper on quadratic finite element discretizations for second order elliptic boundary value problems on mildly structured triangular meshes. For a large class of practically useful grids, the finite element solution u h is proven to be superclose to the inter-polant u I and as a result a postprocessing gradient recovery scheme for u h can be devised. The analysis is based on a number of carefully derived identities. In addition to its own theoretical interests, the result in this paper can be used for deriving asymptotically exact a posteriori error estimators for quadratic finite element methods.

The construction of good extensible rank-1 lattices

Abstract: It has been shown by Hickernell and Niederreiter that there exist generating vectors for integration lattices which yield small integration errors for n = p, p 2 ,... for all integers p ≥ 2. This paper provides algorithms for the construction of generating vectors which are finitely extensible for n = p, p 2 ,... for all integers p ≥ 2. The proofs which show that our algorithms yield good extensible rank-1 lattices are based on a sieve principle. Particularly fast algorithms are obtained by using the fast component-by-component construction of Nuyens and Cools. Analogous results are presented for generating vectors with small weighted star discrepancy.

Robust a-posteriori estimator for advection-diffusion-reaction problems

Abstract: We propose an almost-robust residual-based a-posteriori estimator for the advection-diffusion-reaction model problem. The theory is developed in the one-dimensional setting. The numerical error is measured with respect to a norm which was introduced by the author in 2005 and somehow plays the role that the energy norm has with respect to symmetric and coercive differential operators. In particular, the mentioned norm possesses features that allow us to obtain a meaningful a-posteriori estimator, robust up to a √log(Pe) factor, where Pe is the global Peclet number of the problem. Various numerical tests are performed in one dimension, to confirm the theoretical results and show that the proposed estimator performs better than the usual one known in literature. We also consider a possible two-dimensional extension of our result and only present a few basic numerical tests, indicating that the estimator seems to preserve the good features of the one-dimensional setting.

Convergence of classes of high-order semi-Lagrangian schemes for the Vlasov–Poisson system

Abstract: In this paper we present some classes of high-order semi-Lagran- gian schemes for solving the periodic one-dimensional Vlasov-Poisson system in phase-space on uniform grids. We prove that the distribution function $ f(t,x,v)$ and the electric field $ E(t,x)$ converge in the $ L^2$ norm with a rate of $\displaystyle \mathcal{O}\left(\Delta t^2 +h^{m+1}+ \frac{h^{m+1}}{\Delta t}\right),$ where $ m$ is the degree of the polynomial reconstruction, and $ \Delta t$ and $ h$ are respectively the time and the phase-space discretization parameters

On Schönhage’s algorithm and subquadratic integer gcd computation

Abstract: We describe a new subquadratic left-to-right GCD algorithm, inspired by Schonhage's algorithm for reduction of binary quadratic forms, and compare it to the first subquadratic GCD algorithm discovered by Knuth and Schonhage, and to the binary recursive GCD algorithm of Stehle and Zimmer-mann. The new GCD algorithm runs slightly faster than earlier algorithms, and it is much simpler to implement. The key idea is to use a stop condition for HGCD that is based not on the size of the remainders, but on the size of the next difference. This subtle change is sufficient to eliminate the back-up steps that are necessary in all previous subquadratic left-to-right GCD algorithms. The subquadratic GCD algorithms all have the same asymptotic running time, O(n(log n)(2) log log n).

Inverse-Type Estimates on hp-Finite Element Spaces and Applications.

Abstract: This work is concerned with the development of inverse-type inequalities for piecewise polynomial functions and, in particular, functions belonging to hp-finite element spaces. The cases of positive and negative Sobolev norms are considered for both continuous and discontinuous finite element functions.The inequalities are explicit both in the local polynomial degree and the local mesh size. The assumptions on the hp-finite element spaces are very weak, allowing anisotropic (shape-irregular) elements and varying polynomial degree across elements. Finally, the new inverse-type inequalities are used to derive bounds for the condition number of symmetric stiffness matrices of hp-boundary element method discretisations of integral equations, with element-wise discontinuous basis functions constructed via scaled tensor products of Legendre polynomials.

Optimal logarithmic energy points on the unit sphere

Abstract: We study minimum energy point charges on the unit sphere S d in R d+1 , d ≥ 3, that interact according to the logarithmic potential log (1/r), where r is the Euclidean distance between points. Such optimal N-point configurations are uniformly distributed as N →∞. We quantify this result by estimating the spherical cap discrepancy of optimal energy configurations. The estimate is of order O (N- 1/(d+2) ). Essential is an improvement of the lower bound of the optimal logarithmic energy which yields the second term (1/d)(log N)/N in the asymptotical expansion of the optimal energy. Previously, this was known for the unit sphere in R 3 only. Furthermore, we present an upper bound for the error of integration for an equally-weighted numerical integration rule Q N with the N nodes forming an optimal logarithmic energy configuration. For polynomials p of degree at most n this bound is Cd(N 1 /d/n) -d/2 ||p||∞ as n/N 1/d → 0. For continuous functions f of S d satisfying a Lipschitz condition with constant C f the bound is (12dC f + C' d ||f||∞)O(N -1/(d+2) ) as N → ∞.

On Estimators for Eigenvalue/Eigenvector Approximations

Abstract: We consider a large class of residuum based a posteriori eigen-value/eigenvector estimates and present an abstract framework for proving their asymptotic exactness. Equivalence of the estimator and the error is also established. To demonstrate the strength of our abstract approach we present a detailed study of hierarchical error estimators for Laplace eigenvalue problems in planar polygonal regions. To this end we develop new error analysis for the Galerkin approximation which avoids the use of the strengthened Cauchy-Schwarz inequality and the saturation assumption, and gives reasonable and explicitly computable upper bounds on the discretization error. A brief discussion is also given concerning the design of estimators which are in the same spirit, but are based on different a posteriori techniques-notably, those of gradient recovery type.

An interpolated stochastic algorithm for quasi-linear PDEs

Abstract: In this paper, we improve the forward-backward algorithm for quasi-linear PDEs introduced in a previous work. The new discretization scheme takes advantage of the standing regularity properties of the true solution through an interpolation procedure. For the convergence analysis, we also exploit the optimality of the square Gaussian quantization used to approximate the conditional expectations involved. The resulting bound for the error is closely related to the Holder exponent of the second order spatial derivatives of the true solution and turns out to be more satisfactory than the one previously established.

Sign changes in sums of the Liouville function

Abstract: The Liouville function A(n) is the completely multiplicative function whose value is -1 at each prime. We develop some algorithms for computing the sum T(n) = Σ n k=1 λ(k)/k, and use these methods to determine the smallest positive integer n where T(n) < 0. This answers a question originating in some work of Turan, who linked the behavior of T(n) to questions about the Riemann zeta function. We also study the problem of evaluating Polya's sum L(n) = Σ n k=1 λ(k), and we determine some new local extrema for this function, including some new positive values.

Analysis of a finite element PML approximation for the three dimensional time-harmonic Maxwell problem

Abstract: In our paper [Math. Comp. 76, 2007, 597-614] we considered the acoustic and electromagnetic scattering problems in three spatial dimensions. In particular, we studied a perfectly matched layer (PML) approximation to an electromagnetic scattering problem. We demonstrated both the solvability of the continuous PML approximations and the exponential convergence of the resulting solution to the solution of the original acoustic or electromagnetic problem as the layer increased. In this paper, we consider finite element approximation of the truncated PML electromagnetic scattering problem. Specifically, we consider approximations which result from the use of Nedelec (edge) finite elements. We show that the resulting finite element problem is stable and gives rise to quasi-optimal convergence when the mesh size is sufficiently small.

Global convergence of the Polak-Ribière-Polyak conjugate gradient method with an Armijo-type inexact line search for nonconvex unconstrained optimization problems

Abstract: We propose two algorithms for nonconvex unconstrained optimization problems that employ Polak-Ribiere-Polyak conjugate gradient formula and new inexact line search techniques. We show that the new algorithms converge globally if the function to be minimized has Lipschitz continuous gradients. Preliminary numerical results show that the proposed methods for particularly chosen line search conditions are very promising.

Adaptive wavelet algorithms for elliptic PDE's on product domains

Abstract: With standard isotropic approximation by (piecewise) polynomials of fixed order in a domain D subset of R-d, the convergence rate in terms of the number N of degrees of freedom is inversely proportional to the space dimension d. This so-called curse of dimensionality can be circumvented by applying sparse tensor product approximation, when certain high order mixed derivatives of the approximated function happen to be bounded in L-2. It was shown by Nitsche (2006) that this regularity constraint can be dramatically reduced by considering best N-term approximation from tensor product wavelet bases. When the function is the solution of some well-posed operator equation, dimension independent approximation rates can be practically realized in linear complexity by adaptive wavelet algorithms, assuming that the infinite stiffness matrix of the operator with respect to such a basis is highly compressible. Applying piecewise smooth wavelets, we verify this compressibility for general, non-separable elliptic PDEs in tensor domains. Applications of the general theory developed include adaptive Galerkin discretizations of multiple scale homogenization problems and of anisotropic equations which are robust, i.e., independent of the scale parameters, resp. of the size of the anisotropy.

Error analysis for the elastic flow of parametrized curves

Abstract: We analyze a semidiscrete numerical scheme for approximating the evolution of parametric curves by elastic flow in R n . The fourth order equation is split into two coupled second order problems, which are approximated by linear finite elements. We prove error bounds for the resulting scheme and present numerical test calculations that confirm our analysis.

New expansions of numerical eigenvalues for -Δ= by nonconforming elements

Abstract: The paper explores new expansions of the eigenvalues for -Δu = Λpu in S with Dirichlet boundary conditions by the bilinear element (denoted Q 1 ) and three nonconforming elements, the rotated bilinear element (denoted Q 1 rot ), the extension of Q 1 rot (denoted EQ 1 rot ) and Wilson's elements. The expansions indicate that Q 1 and Q 1 rot provide upper bounds of the eigenvalues, and that EQ 1 rot and Wilson's elements provide lower bounds of the eigenvalues. By extrapolation, the O(h 4 ) convergence rate can be obtained, where h is the maximal boundary length of uniform rectangles. Numerical experiments are carried out to verify the theoretical analysis made.