Showing papers in "Mathematics of Computation in 2015"

A class of second order difference approximations for solving space fractional diffusion equations

Abstract: A class of second order approximations, called the weighted and shifted Grunwald difference (WSGD) operators, are proposed for Riemann-Liouville fractional derivatives, with their effective applications to numerically solving space fractional diffusion equations in one and two dimensions. The stability and convergence of our difference schemes for space fractional diffusion equations with constant coefficients in one and two dimensions are theoretically established. Several numerical examples are implemented to test the efficiency of the numerical schemes and confirm the convergence order, and the numerical results for variable coefficients problem are also presented.

Generalized Jacobi functions and their applications to fractional differential equations

Abstract: In this paper, we consider spectral approximation of fractional dif- ferential equations (FDEs). A main ingredient of our approach is to define a new class of generalized Jacobi functions (GJFs), which is intrinsically related to fractional calculus and can serve as natural basis functions for properly de- signed spectral methods for FDEs. We establish spectral approximation results for these GJFs in weighted Sobolev spaces involving fractional derivatives. We construct efficient GJF-Petrov-Galerkin methods for a class of prototypical fractional initial value problems (FIVPs) and fractional boundary value prob- lems (FBVPs) of general order, and we show that with an appropriate choice of the parameters in GJFs, the resulting linear systems are sparse and well- conditioned. Moreover, we derive error estimates with convergence rates only depending on the smoothness of data, so true spectral accuracy can be attained if the data are smooth enough. The ideas and results presented in this paper will be useful in dealing with more general FDEs involving Riemann-Liouville or Caputo fractional derivatives.

Numerical approximation of fractional powers of elliptic operators

Abstract: We present and study a novel numerical algorithm to approximate the action of $T^\beta:=L^{-\beta}$ where $L$ is a symmetric and positive definite unbounded operator on a Hilbert space $H_0$. The numerical method is based on a representation formula for $T^{-\beta}$ in terms of Bochner integrals involving $(I+t^2L)^{-1}$ for $t\in(0,\infty)$. To develop an approximation to $T^\beta$, we introduce a finite element approximation $L_h$ to $L$ and base our approximation to $T^\beta$ on $T_h^\beta:= L_h^{-\beta}$. The direct evaluation of $T_h^{\beta}$ is extremely expensive as it involves expansion in the basis of eigenfunctions for $L_h$. The above mentioned representation formula holds for $T_h^{-\beta}$ and we propose three quadrature approximations denoted generically by $Q_h^\beta$. The two results of this paper bound the errors in the $H_0$ inner product of $T^\beta-T_h^\beta\pi_h$ and $T_h^\beta-Q_h^\beta$ where $\pi_h$ is the $H_0$ orthogonal projection into the finite element space. We note that the evaluation of $Q_h^\beta$ involves application of $(I+(t_i)^2L_h)^{-1}$ with $t_i$ being either a quadrature point or its inverse. Efficient solution algorithms for these problems are available and the problems at different quadrature points can be straightforwardly solved in parallel. Numerical experiments illustrating the theoretical estimates are provided for both the quadrature error $T_h^\beta-Q_h^\beta$ and the finite element error $T^\beta-T_h^\beta\pi_h$.

Variational formulation of problems involving fractional order differential operators

Abstract: In this work, we consider boundary value problems involving either Caputo or Riemann-Liouville fractional derivatives of order α ∈ (1, 2) on the unit interval (0, 1). These fractional derivatives lead to nonsymmetric boundary value problems, which are investigated from a variational point of view. The variational problem for the Riemann-Liouville case is coercive on the space Hα/2 0 (0, 1) but the solutions are less regular, whereas that for the Caputo case involves different test and trial spaces. The numerical analysis of these problems requires the so-called shift theorems which show that the solutions of the variational problem are more regular. The regularity pickup enables one to establish convergence rates of the finite element approximations. The analytical theory is then applied to the Sturm-Liouville problem involving a fractional derivative in the leading term. Finally, extensive numerical results are presented to illustrate the error estimates for the source problem and eigenvalue problem.

An accelerated randomized Kaczmarz algorithm

Abstract: The randomized Kaczmarz ($\RK$) algorithm is a simple but powerful approach for solving consistent linear systems $Ax=b$. This paper proposes an accelerated randomized Kaczmarz ($\ARK$) algorithm with better convergence than the standard $\RK$ algorithm on ill conditioned problems. The per-iteration cost of $\RK$ and $\ARK$ are similar if $A$ is dense, but $\RK$ is much more able to exploit sparsity in $A$ than is $\ARK$. To deal with the sparse case, an efficient implementation for $\ARK$, called $\SARK$, is proposed. A comparison of convergence rates and average per-iteration complexities among $\RK$, $\ARK$, and $\SARK$ is given, taking into account different levels of sparseness and conditioning. Comparisons with the leading deterministic algorithm --- conjugate gradient applied to the normal equations --- are also given. Finally, the analysis is validated via computational testing.

Linear regression MDP scheme for discrete backward stochastic differential equations under general conditions

Abstract: We design a numerical scheme for solving the Multi step-forward Dynamic Programming (MDP) equation arising from the time-discretization of backward stochastic differential equations. The generator is assumed to be locally Lipschitz, which includes some cases of quadratic drivers. When the large sequence of conditional expectations is computed using empirical least-squares regressions, under general conditions we establish an upper bound error as the average, rather than the sum, of local regression errors only, suggesting that our error estimation is tight. Despite the nested regression problems, the interdependency errors are justified to be at most of the order of the statistical regression errors (up to logarithmic factor). Finally, we optimize the algorithm parameters, depending on the dimension and on the smoothness of value functions, in the limit as the time mesh size goes to zero and compute the complexity needed to achieve a given accuracy. Numerical experiments are presented illustrating theoretical convergence estimates.

Convergence analysis of a fully discrete finite difference scheme for the Cahn-Hilliard-Hele-Shaw equation

Abstract: We present an error analysis for an unconditionally energy stable, fully discrete finite difference scheme for the Cahn-Hilliard-Hele-Shaw equa- tion, a modified Cahn-Hilliard equation coupled with the Darcy flow law. The scheme, proposed by S. M. Wise, is based on the idea of convex splitting. In this paper, we rigorously prove first order convergence in time and second or- der convergence in space. Instead of the (discrete)L ∞ (0,T;L 2 )∩L 2 (0,T;H 2 h ) error estimate, which would represent the typical approach, we provide a dis- creteL ∞ (0,T;H 1 )∩L 2 (0,T;H 3 h )e rror estimate for the phase variable, which allows us to treat the nonlinear convection term in a straightforward way. Our convergence is unconditionalin the sense that the time stepsis in no way constrained by the mesh spacingh .T his is accomplished with the help of anL 2 (0,T;H 3 h ) bound of the numerical approximation of the phase variable. To facilitate both the stability and convergence analyses, we establish a finite difference analog of a Gagliardo-Nirenberg type inequality.

Discrete and conforming smooth de Rham complexes in three dimensions

Abstract: Conforming discrete de Rham complexes consisting of finite element spaces with extra smoothness are constructed. In particular, we develop H2, H1(curl), H1 and L2 conforming finite element spaces and show that an exactness property is satisfied. These results naturally lead to discretizations for Stokes and Brinkman type problems as well as conforming approximations to fourth order curl problems. In addition, we reduce the question of stability of the three-dimensional Scott-Vogelius finite element to a simply stated conjecture.

Weak convergence for a spatial approximation of the nonlinear stochastic heat equation

Abstract: We find the weak rate of convergence of approximate solutions of the nonlinear stochastic heat equation, when discretized in space by a standard finite element method. Both multiplicative and additive noise is considered under different assumptions. This extends an earlier result of Debussche in which time discretization is considered for the stochastic heat equation perturbed by white noise. It is known that this equation only has a solution in one space dimension. In order to get results for higher dimensions, colored noise is considered here, besides the white noise case where considerably weaker assumptions on the noise term is needed. Integration by parts in the Malliavin sense is used in the proof. The rate of weak convergence is, as expected, essentially twice the rate of strong convergence.

A blob method for the aggregation equation

Abstract: Author(s): Craig, Katy; Bertozzi, Andrea L | Abstract: Motivated by classical vortex blob methods for the Euler equations, we develop a numerical blob method for the aggregation equation. This provides a counterpoint to existing literature on particle methods. By regularizing the velocity field with a mollifier or "blob function", the blob method has a faster rate of convergence and allows a wider range of admissible kernels. In fact, we prove arbitrarily high polynomial rates of convergence to classical solutions, depending on the choice of mollifier. The blob method conserves mass and the corresponding particle system is both energy decreasing for a regularized free energy functional and preserves the Wasserstein gradient flow structure. We consider numerical examples that validate our predicted rate of convergence and illustrate qualitative properties of the method.

Strong rate of convergence for the Euler-Maruyama approximation of stochastic differential equations with irregular coefficients

Abstract: We consider the Euler-Maruyama approximation for multi-dimensional stochastic differential equations with irregular coefficients. We provide the rate of strong convergence where the possibly discontinuous drift coefficient satisfies a one-sided Lipschitz condition and the diffusion coefficient is Holder continuous.

Rigorous numerics for analytic solutions of differential equations: the radii polynomial approach

Abstract: Judicious use of interval arithmetic, combined with careful pen and paper estimates, leads to effective strategies for computer assisted analysis of nonlinear operator equations. The method of radii polynomials is an efficient tool for bounding the smallest and largest neighborhoods on which a Newton-like operator associated with a nonlinear equation is a contraction mapping. The method has been used to study solutions of ordinary, partial, and delay differential equations such as equilibria, periodic orbits, solutions of initial value problems, heteroclinic and homoclinic connecting orbits in the C category of functions. In the present work we adapt the method of radii polynomials to the analytic category. For ease of exposition we focus on studying periodic solutions in Cartesian products of infinite sequence spaces. We derive the radii polynomials for some specific application problems, and give a number of computer assisted proofs in the analytic framework.

Computing periods of rational integrals

Abstract: A period of a rational integral is the result of integrating, with respect to one or several variables, a rational function over a closed path. This work focuses particularly on periods depending on a parameter: in this case the period under consideration satisfies a linear differential equation, the Picard-Fuchs equation. I give a reduction algorithm that extends the Griffiths-Dwork reduction and apply it to the computation of Picard-Fuchs equations. The resulting algorithm is elementary and has been successfully applied to problems that were previously out of reach.

Convergence of a nonlinear entropy diminishing Control Volume Finite Element scheme for solving anisotropic degenerate parabolic equations

Abstract: In this paper, we propose and analyze a Control Volume Finite Elements (CVFE) scheme for solving possibly degenerated parabolic equations. This scheme does not require the introduction of the so-called Kirchhoff transform in its definition. We prove that the discrete solution obtained \emph{via} the scheme remains in the physical range, and that the natural entropy of the problem decreases with time. The convergence of the method is proved as the discretization steps tend to $0$. Finally, numerical examples illustrate the efficiency of the method.

Conditional bounds for the least quadratic non-residue and related problems

Abstract: This paper studies explicit and theoretical bounds for several inter- esting quantities in number theory, conditionally on the Generalized Riemann Hypothesis. Specically, we improve the existing explicit bounds for the least qua- dratic non-residue and the least prime in an arithmetic progression. We also rene the classical conditional bounds of Littlewood for L-functions at s = 1. In par- ticular, we derive explicit upper and lower bounds for L(1; ) and (1 +it), and deduce explicit bounds for the class number of imaginary quadratic elds. Finally, we improve the best known theoretical bounds for the least quadratic non-residue, and more generally, the least k-th power non-residue.

Optimal error estimates of the direct discontinuous Galerkin method for convection-diffusion equations

Abstract: . In this paper, we present the optimal L2-error estimate ofO(hk+1) for polynomial elements of degree k of the semidiscrete direct discontinuous Galerkin method for convection-diffusion equations. The main technical difficulty lies in the control of the inter-element jump terms which arise because of the convection and the discontinuous nature of numerical solutions. The main idea is to use some global projections satisfying interface conditions dictated by the choice of numerical fluxes so that trouble terms at the cell interfaces are eliminated or controlled. In multi-dimensional case, the orders of k + 1 hinge on a superconvergence estimate when tensor product polynomials of degree k are used on Cartesian grids. A collection of projection errors in both oneand multi-dimensional cases is established.

Exploring the tree of numerical semigroups

Abstract: In this paper we describe an algorithm visiting all numerical semigroups up to a given genus using a well suited representation. The interest of this algorithm is that it fits particularly well the architecture of modern computers allowing very large optimizations: we obtain the number of numerical semigroups of genus g 67 and we confirm the Wilf conjecture for g 60.

A fully well-balanced, positive and entropy-satisfying Godunov-type method for the shallow-water equations

Abstract: This work is devoted to the derivation of a fully well-balanced numerical scheme for the well-known shallow-water model. During the last two decades, several well-balanced strategies have been introduced with a special attention to the exact capture of the stationary states associated with the so-called lake at rest. By fully well-balanced, we mean here that the proposed Godunov-type method is also able to preserve stationary states with non zero velocity. The numerical procedure is shown to preserve the positiveness of the water height and satisfies a discrete entropy inequality.

On high order finite element spaces of differential forms

Abstract: We show how the high order finite element spaces of differential forms due to Raviart-Thomas-Nedelec-Hiptmair fit into the framework of finite element systems, in an elaboration of the finite element exterior calculus of Arnold-Falk-Winther. Based on observations by Bossavit, we provide new low order degrees of freedom. As an alternative to existing choices of bases, we provide canonical resolutions in terms of scalar polynomials and Whitney forms.

Estimating $\pi(x)$ and related functions under partial RH assumptions

Abstract: The aim of this paper is to give a direct interpretation of the validity of the Riemann hypothesis up to a certain height $T$ in terms of the prime-counting function $\pi(x)$. This is done by proving the well-known explicit Schoenfeld bound on the RH to hold as long as $4.92 \sqrt{x/\log(x)} \leq T$. Similar statements are proven for the Riemann prime-counting function and the Chebyshov functions $\psi(x)$ and $\vartheta(x)$. Apart from that, we also improve some of the existing bounds of Chebyshov type for the function $\psi(x)$.

Spectral analysis and spectral symbol of matrices in isogeometric collocation methods

Abstract: A linear full elliptic second order Partial Dierential Equation (PDE), defined on a d-dimensional domain , is approximated by the isogeometric Galerkin method based on uniform tensor-product Bsplines of degrees (p1;:::; pd). The considered approximation process leads to a d-level stiness matrix, banded in a multilevel sense. This matrix is close to a d-level Toeplitz structure when the PDE coecients are constant and the physical domain is just the hypercube (0; 1) d without using any geometry map. In such a simplified case, a detailed spectral analysis of the stiness matrices has been carried out in a previous work. In this paper, we complete the picture by considering non-constant PDE coecients and an arbitrary domain , parameterized with a non-trivial geometry map. We compute and study the spectral symbol of the related stiness matrices. This symbol describes the asymptotic eigenvalue distribution when the fineness parameters tend to zero (so that the matrix-size tends to infinity). The mathematical technique used for computing the symbol is based on the theory of Generalized Locally Toeplitz (GLT) sequences.

Analogues of Vélu’s formulas for isogenies on alternate models of elliptic curves

Abstract: Isogenies of elliptic curves have been well-studied, in part because there are several cryptographic applications. Using Velu’s formula, isogenies can be evaluated explicitly given their kernel. However, Velu’s formula applies to elliptic curves given by a Weierstrass equation. In this paper we show how to similarly evaluate isogenies on Edwards curves and Huff curves. Edwards and Huff curves are new normal forms for elliptic curves, different than the traditional Weierstrass form.

Eventual linear convergence of the Douglas-Rachford iteration for basis pursuit

Abstract: We provide a simple analysis of the Douglas-Rachford splitting algorithm in the context of ‘ 1 minimization with linear constraints, and quantify the asymptotic linear convergence rate in terms of principal angles between relevant vector spaces. In the compressed sensing setting, we show how to bound this rate in terms of the restricted isometry constant. More general iterative schemes obtained by ‘ 2 -regularization and over-relaxation including the dual split Bregman method [27] are also treated, which answers the question how to choose the relaxation and soft-thresholding parameters to accelerate the asymptotic convergence rate. We make no attempt at characterizing the transient regime preceding the onset of linear convergence.

Enumeration of MOLS of small order

Abstract: We report the results of a computer investigation of sets of mutually orthogonal latin squares (MOLS) of small order. For $n\le9$ we 1. Determine the number of orthogonal mates for each species of latin square of order $n$. 2. Calculate the proportion of latin squares of order $n$ that have an orthogonal mate, and the expected number of mates when a square is chosen uniformly at random. 3. Classify all sets of MOLS of order $n$ up to various different notions of equivalence. We also provide a triple of latin squares of order 10 that is the closest to being a set of MOLS so far found.

An explicit density estimate for Dirichlet $L$-series

Abstract: Dirichlet L-series L(s, χ) = ∑ n≥1 χ(n)n −s associated to primitive Dirichlet characters χ are one of the keys to the distribution of primes. Even the simple case χ = 1 which corresponds to the Riemann zeta-function contains many informations on primes and on the Farey dissection. There have been many generalizations of these notions, and they all have arithmetical properties and/or applications, see [45, 29, 33] for instance. Investigations concerning these functions range over many directions, see [14] or [43]. We note furthermore that Dirichlet characters have been the subject of numerous studies, see [2, 50, 4]; Dirichlet series in themselves are still mysterious, see [3] and [6]. One of the main problem concerns the location of the zeroes of these functions in the strip 0 <

Superconvergence of Local Discontinuous Galerkin methods for one-dimensional linear parabolic equations

Abstract: In this paper, we study the superconvergence properties of the LDG method for the one-dimensional linear Schrodinger equation. We build a special interpolation function by constructing a correction function, and prove the numerical solution is superclose to the interpolation function in the \(L^{2}\) norm. The order of superconvergence is \(2k+1\), when the polynomials of degree at most k are used. Even though the linear Schrodinger equation involves only second order spatial derivative, it is actually a wave equation because of the coefficient i. It is not coercive and there is no control on the derivative for later time based on the initial condition of the solution itself, as for the parabolic case. In our analysis, the special correction functions and special initial conditions are required, which are the main differences from the linear parabolic equations. We also rigorously prove a \((2k+1)\)-th order superconvergence rate for the domain, cell averages, and the numerical fluxes at the nodes in the maximal and average norm. Furthermore, we prove the function value and the derivative approximation are superconvergent with a rate of \((k+2)\)-th order at the Radau points. All theoretical findings are confirmed by numerical experiments.

On the accuracy of finite element approximations to a class of interface problems

Abstract: The jump is defined as [∇u · n] = ∇u · n + ∇u ·n where u = u|Ω± and n is the unit outward pointing normal to Ω (see figure 1). Also, we denote [u] = u − u. Many numerical methods have been developed for problem (1.1). Perhaps the most notable ones are the finite difference method of Peskin [18] (i.e., immersed boundary method) and the method of LeVeque and Li [11] (i.e., the immersed interface method ; see also the method of Mayo [14, 15, 16]) .The method of LeVeque and Li [11] was developed for the more general problem with discontinuous diffusion coefficients, while the method of Peskin [18] was developed for fluid flow problems with an immersed boundary. Although the method of Peskin [18] is formulated with a force function F that incorporates the elastic force of the immersed boundary Γ, it was shown in [19] that it can be re-formulated as an interface problem (with α = 0) where β encodes the elastic force. Since the two important papers [18, 11] there have been many articles extending or improving these methods. In particular, finite element versions of these methods have appeared; see for example [3, 9, 6, 2]. For the above problem (α = 0), it is well known that the method of Peskin [18] is only first-order accurate whereas the method of LeVeque and Li [11] is second-order

The use of bad primes in rational reconstruction

Abstract: A standard method for computing a rational number from its values modulo a collection of primes is to determine its value modulo the product of the primes via Chinese Remaindering, and then use Farey sequences for rational reconstruction. This method is guaranteed to work if we restrict ourselves to ”good” primes. Depending on the particular application, however, there is often no efficient way of finding good primes. This note shows that in most situations, we can simply ignore this problem. With regard to applications, we are particularly interested in the design of modular and, thus, parallel versions of algorithms in commutative algebra and algebraic geometry. Here, typically, the final result consists of one or several a priori unknown ideals which are found via constructions yielding the (reduced) Grobner bases of the ideals.