A Fast Finite Difference Method for a Continuous Static Linear Bond-Based Peridynamics Model of Mechanics

The peridynamic nonlocal continuum model for solid mechanics is an integro-differential equation that does not involve spatial derivatives of the displacement field. Several numerical methods such as finite element method and collocation method have been developed and analyzed in many articles. However, there is no theory to give a finite difference method because the model does not involve spatial derivatives of the displacement field. Here, we consider a finite difference scheme to solve a continuous static bond-based peridynamics model of mechanics based on its equivalent partial integro-differential equations. Furthermore, we present a fast solution technique to accelerate Toeplitz matrix-vector multiplications arising from finite difference discretization respectively. This fast solution technique is based on a fast Fourier transform and depends on the special structure of coefficient matrices, and it helps to reduce the computational work from \(O(N^{3})\) required by traditional methods to O(Nlog\(^{2}N)\) and the memory requirement from \(O(N^{2})\) to O(N) without using any lossy compression, where N is the number of unknowns. Moreover, the applicability and accuracy of the scheme are demonstrated by numerical experiments to support our theoretical analysis.     

Fast Iterative Method with a Second-Order Implicit Difference Scheme for Time-Space Fractional Convection–Diffusion Equation

In this paper we intend to establish fast numerical approaches to solve a class of initial-boundary problem of time-space fractional convection–diffusion equations. We present a new unconditionally stable implicit difference method, which is derived from the weighted and shifted Grünwald formula, and converges with the second-order accuracy in both time and space variables. Then, we show that the discretizations lead to Toeplitz-like systems of linear equations that can be efficiently solved by Krylov subspace solvers with suitable circulant preconditioners. Each time level of these methods reduces the memory requirement of the proposed implicit difference scheme from \({\mathcal {O}}(N^2)\) to \({\mathcal {O}}(N)\) and the computational complexity from \({\mathcal {O}}(N^3)\) to \({\mathcal {O}}(N\log N)\) in each iterative step, where N is the number of grid nodes. Extensive numerical examples are reported to support our theoretical findings and show the utility of these methods over traditional direct solvers of the implicit difference method, in terms of computational cost and memory requirements.     

A Mixed Discontinuous Galerkin Method Without Interior Penalty for Time-Dependent Fourth Order Problems

A novel discontinuous Galerkin (DG) method is developed to solve time-dependent bi-harmonic type equations involving fourth derivatives in one and multiple space dimensions. We present the spatial DG discretization based on a mixed formulation and central interface numerical fluxes so that the resulting semi-discrete schemes are \(L^2\) stable even without interior penalty. For time discretization, we use Crank–Nicolson so that the resulting scheme is unconditionally stable and second order in time. We present the optimal \(L^2\) error estimate of \(O(h^{k+1})\) for polynomials of degree k for semi-discrete DG schemes, and the \(L^2\) error of \(O(h^{k+1} +(\Delta t)^2)\) for fully discrete DG schemes. Extensions to more general fourth order partial differential equations and cases with non-homogeneous boundary conditions are provided. Numerical results are presented to verify the stability and accuracy of the schemes. Finally, an application to the one-dimensional Swift–Hohenberg equation endowed with a decay free energy is presented.     

Quantum Fourier transform in computational basis

The quantum Fourier transform, with exponential speed-up compared to the classical fast Fourier transform, has played an important role in quantum computation as a vital part of many quantum algorithms (most prominently, Shor’s factoring algorithm). However, situations arise where it is not sufficient to encode the Fourier coefficients within the quantum amplitudes, for example in the implementation of control operations that depend on Fourier coefficients. In this paper, we detail a new quantum scheme to encode Fourier coefficients in the computational basis, with fidelity \(1 - \delta \) and digit accuracy \(\epsilon \) for each Fourier coefficient. Its time complexity depends polynomially on \(\log (N)\), where N is the problem size, and linearly on \(1/\delta \) and \(1/\epsilon \). We also discuss an application of potential practical importance, namely the simulation of circulant Hamiltonians.     

Recursive Computation of Logarithmic Derivatives,Ratios, and Products of Spheroidal Harmonics and Modified Bessel Functions and Applications

Spheroidal harmonics and modified Bessel functions have wide applications in scientific and engineering computing. Recursive methods are developed to compute the logarithmic derivatives, ratios, and products of the prolate spheroidal harmonics (\(P_n^m(x)\), \(Q_n^m(x)\), \(n\ge m\ge 0\), \(x>1\)), the oblate spheroidal harmonics (\(P_n^m(ix)\), \(Q_n^m(ix)\), \(n\ge m\ge 0\), \(x>0\)), and the modified Bessel functions (\(I_n(x)\), \(K_n(x)\), \(n\ge 0\), \(x>0\)) in order to avoid direct evaluation of these functions that may easily cause overflow/underflow for high degree/order and for extreme argument. Stability analysis shows the proposed recursive methods are stable for realistic degree/order and argument values. Physical examples in electrostatics are given to validate the recursive methods.     

A unified framework for two-grid methods for a class of nonlinear problems

In this article, we present a unified error analysis of two-grid methods for a class of nonlinear problems. We first study the two-grid method of Xu by recasting the methodology in the abstract framework of Brezzi, Rappaz, and Raviart (BRR) for approximation of branches of nonsingular solutions and derive a priori error estimates. Our convergence results indicate that the correct scaling between fine and coarse meshes is given by \(h={{\mathcal {O}}}(H^2)\) for all the nonlinear problems which can be written in and applied to the BRR framework. Next, a correction step can be added to the two-grid algorithm, which allows the choice \(h={\mathcal O}(H^3)\). On the other hand, the particular BRR framework with duality pairing, if it is applied to a semilinear problem, allows a higher order relation \(h={{\mathcal {O}}}(H^4)\). Furthermore, even the choice \(h={{\mathcal {O}}}(H^5)\) is possible with the correction step either on fine mesh or coarse mesh. In addition, elliptic problems with gradient nonlinearities and the Naiver–Stokes equations are considered to illustrate our unified theory. Finally, numerical experiments are conducted to confirm our theoretical findings. Numerical results indicate that the correction step used as a simple postprocessing enhances the solution accuracy, particularly for problems with layers.     

Efficient Preconditioning of <Emphasis Type="Italic">hp</Emphasis>-FEM Matrices by Hierarchical Low-Rank Approximations

We introduce a preconditioner based on a hierarchical low-rank compression scheme of Schur complements. The construction is inspired by standard nested dissection, and relies on the assumption that the Schur complements can be approximated, to high precision, by Hierarchically-Semi-Separable matrices. We build the preconditioner as an approximate \(LDM^t\) factorization of a given matrix A, and no knowledge of A in assembled form is required by the construction. The \(LDM^t\) factorization is amenable to fast inversion, and the action of the inverse can be determined fast as well. We investigate the behavior of the preconditioner in the context of DG finite element approximations of elliptic and hyperbolic problems, with respect to both the mesh size and the order of approximation.     

Numerical Methods and Comparison for the Dirac Equation in the Nonrelativistic Limit Regime

We analyze rigorously error estimates and compare numerically spatial/temporal resolution of various numerical methods for the discretization of the Dirac equation in the nonrelativistic limit regime, involving a small dimensionless parameter \(0<\varepsilon \ll 1\) which is inversely proportional to the speed of light. In this limit regime, the solution is highly oscillatory in time, i.e. there are propagating waves with wavelength \(O(\varepsilon ^2)\) and O(1) in time and space, respectively. We begin with several frequently used finite difference time domain (FDTD) methods and obtain rigorously their error estimates in the nonrelativistic limit regime by paying particular attention to how error bounds depend explicitly on mesh size h and time step \(\tau \) as well as the small parameter \(\varepsilon \). Based on the error bounds, in order to obtain ‘correct’ numerical solutions in the nonrelativistic limit regime, i.e. \(0<\varepsilon \ll 1\), the FDTD methods share the same \(\varepsilon \)-scalability on time step and mesh size as: \(\tau =O(\varepsilon ^3)\) and \(h=O(\sqrt{\varepsilon })\). Then we propose and analyze two numerical methods for the discretization of the Dirac equation by using the Fourier spectral discretization for spatial derivatives combined with the symmetric exponential wave integrator and time-splitting technique for temporal derivatives, respectively. Rigorous error bounds for the two numerical methods show that their \(\varepsilon \)-scalability is improved to \(\tau =O(\varepsilon ^2)\) and \(h=O(1)\) when \(0<\varepsilon \ll 1\). Extensive numerical results are reported to support our error estimates.     

The Partial Fast Fourier Transform

An efficient algorithm for computing the one-dimensional partial fast Fourier transform \(f_j=\sum _{k=0}^{c(j)}e^{2\pi ijk/N} F_k\) is presented. Naive computation of the partial fast Fourier transform requires \({\mathcal O}(N^2)\) arithmetic operations for input data of length N. Unlike the standard fast Fourier transform, the partial fast Fourier transform imposes on the frequency variable k a cutoff function c(j) that depends on the space variable j; this prevents one from directly applying standard FFT algorithms. It is shown that the space–frequency domain can be partitioned into rectangular and trapezoidal subdomains over which efficient algorithms can be developed. As in the previous work of Ying and Fomel (Multiscale Model Simul 8(1):110–124, 2009), the contribution from rectangular regions can be reduced to a series of fractional-phase Fourier transforms over squares, each of which can be reduced to a convolution. In this work, we demonstrate that the partial Fourier transform over trapezoidal domains can also be reduced to a convolution. Since the computational complexity of a dealiased convolution of N inputs is \({\mathcal O}(N\log N)\), a fast algorithm for the partial Fourier transform is achieved, with a lower overall coefficient than obtained by Ying and Fomel.     

A greedy approximation algorithm for minimum-gap scheduling

We consider scheduling of unit-length jobs with release times and deadlines, where the objective is to minimize the number of gaps in the schedule. Polynomial-time algorithms for this problem are known, yet they are rather inefficient, with the best algorithm running in time \(O(n^4)\) and requiring \(O(n^3)\) memory. We present a greedy algorithm that approximates the optimum solution within a factor of 2 and show that our analysis is tight. Our algorithm runs in time \(O(n^2 \log n)\) and needs only O(n) memory. In fact, the running time is \(O(n (g^*+1)\log n)\), where \(g^*\) is the minimum number of gaps.     

Hybridized Discontinuous Galerkin Method for Elliptic Interface Problems: Error Estimates Under Low Regularity Assumptions of Solutions

New hybridized discontinuous Galerkin (HDG) methods for the interface problem for elliptic equations are proposed. Unknown functions of our schemes are \(u_h\) in elements and \(\hat{u}_h\) on inter-element edges. That is, we formulate our schemes without introducing the flux variable. We assume that subdomains \(\Omega _1\) and \(\Omega _2\) are polyhedral domains and that the interface \(\Gamma =\partial \Omega _1\cap \partial \Omega _2\) is polyhedral surface or polygon. Moreover, \(\Gamma \) is assumed to be expressed as the union of edges of some elements. We deal with the case where the interface is transversely connected with the boundary of the whole domain \(\overline{\Omega }=\overline{\Omega _1\cap \Omega _2}\). Consequently, the solution u of the interface problem may not have a sufficient regularity, say \(u\in H^2(\Omega )\) or \(u|_{\Omega _1}\in H^2(\Omega _1)\), \(u|_{\Omega _2}\in H^2(\Omega _2)\). We succeed in deriving optimal order error estimates in an HDG norm and the \(L^2\) norm under low regularity assumptions of solutions, say \(u|_{\Omega _1}\in H^{1+s}(\Omega _1)\) and \(u|_{\Omega _2}\in H^{1+s}(\Omega _2)\) for some \(s\in (1/2,1]\), where \(H^{1+s}\) denotes the fractional order Sobolev space. Numerical examples to validate our results are also presented.     

A Two-Grid Block-Centered Finite Difference Method for the Nonlinear Time-Fractional Parabolic Equation

In this article, a two-grid block-centered finite difference scheme is introduced and analyzed to solve the nonlinear time-fractional parabolic equation. This method is considered where the nonlinear problem is solved only on a coarse grid of size H and a linear problem is solved on a fine grid of size h. Stability results are proven rigorously. Error estimates are established on non-uniform rectangular grid which show that the discrete \(L^{\infty }(L^2)\) and \(L^2(H^1)\) errors are \(O(\triangle t^{2-\alpha }+h^2+H^3)\). Finally, some numerical experiments are presented to show the efficiency of the two-grid method and verify that the convergence rates are in agreement with the theoretical analysis.     

A Weak Galerkin Method with an Over-Relaxed Stabilization for Low Regularity Elliptic Problems

A new weak Galerkin (WG) finite element method is developed and analyzed for solving second order elliptic problems with low regularity solutions in the Sobolev space \(W^{2,p}(\Omega )\) with \(p\in (1,2)\). A WG stabilizer was introduced by Wang and Ye (Math Comput 83:2101–2126, 2014) for a simpler variational formulation, and it has been commonly used since then in the WG literature. In this work, for the purpose of dealing with low regularity solutions, we propose to generalize the stabilizer of Wang and Ye by introducing a positive relaxation index to the mesh size h. The relaxed stabilization gives rise to a considerable flexibility in treating weak continuity along the interior element edges. When the norm index \(p\in (1,2]\), we strictly derive that the WG error in energy norm has an optimal convergence order \(O(h^{l+1-\frac{1}{p}-\frac{p}{4}})\) by taking the relaxed factor \(\beta =1+\frac{2}{p}-\frac{p}{2}\), and it also has an optimal convergence order \(O(h^{l+2-\frac{2}{p}})\) in \(L^2\) norm when the solution \(u\in W^{l+1,p}\) with \(p\in [1,1+\frac{2}{p}-\frac{p}{2}]\) and \(l\ge 1\). It is recovered for \(p=2\) that with the choice of \(\beta =1\), error estimates in the energy and \(L^2\) norms are optimal for the source term in the sobolev space \(L^2\). Weak variational forms of the WG method give rise to desirable flexibility in enforcing boundary conditions and can be easily implemented without requiring a sufficiently large penalty factor as in the usual discontinuous Galerkin methods. In addition, numerical results illustrate that the proposed WG method with an over-relaxed factor \(\beta (\ge 1)\) converges at optimal algebraic rates for several low regularity elliptic problems.     

Adjoint-Based,Superconvergent Galerkin Approximations of Linear Functionals

We propose a new technique for computing highly accurate approximations to linear functionals in terms of Galerkin approximations. We illustrate the technique on a simple model problem, namely, that of the approximation of J(u), where \(J(\cdot )\) is a very smooth functional and u is the solution of a Poisson problem; we assume that the solution u and the solution of the adjoint problem are both very smooth. It is known that, if \(u_h\) is the approximation given by the continuous Galerkin method with piecewise polynomials of degree \(k>0\), then, as a direct consequence of its property of Galerkin orthogonality, the functional \(J(u_h)\) converges to J(u) with a rate of order \(h^{2k}\). We show how to define approximations to J(u), with a computational effort about twice of that of computing \(J(u_h)\), which converge with a rate of order \(h^{4k}\). The new technique combines the adjoint-recovery method for providing precise approximate functionals by Pierce and Giles (SIAM Rev 42(2):247–264, 2000), which was devised specifically for numerical approximations without a Galerkin orthogonality property, and the accuracy-enhancing convolution technique of Bramble and Schatz (Math Comput 31(137):94–111, 1977), which was devised specifically for numerical methods satisfying a Galerkin orthogonality property, that is, for finite element methods like, for example, continuous Galerkin, mixed, discontinuous Galerkin and the so-called hybridizable discontinuous Galerkin methods. For the latter methods, we present numerical experiments, for \(k=1,2,3\) in one-space dimension and for \(k=1,2\) in two-space dimensions, which show that \(J(u_h)\) converges to J(u) with order \(h^{2k+1}\) and that the new approximations converges with order \(h^{4k}\). The numerical experiments also indicate, for the p-version of the method, that the rate of exponential convergence of the new approximations is about twice that of \(J(u_h)\).     

A locking-free stabilized mixed finite element method for linear elasticity: the high order case

In this paper, we propose a locking-free stabilized mixed finite element method for the linear elasticity problem, which employs a jump penalty term for the displacement approximation. The continuous piecewise k-order polynomial space is used for the stress and the discontinuous piecewise \((k-1)\)-order polynomial space for the displacement, where we require that \(k\ge 3\) in the two dimensions and \(k\ge 4\) in the three dimensions. The method is proved to be stable and k-order convergent for the stress in \(H(\mathrm {div})\)-norm and for the displacement in \(L^2\)-norm. Further, the convergence does not deteriorate in the nearly incompressible or incompressible case. Finally, the numerical results are presented to illustrate the optimal convergence of the stabilized mixed method.     

Close to linear space routing schemes

Let \(G=(V,E)\) be an unweighted undirected graph with n vertices and m edges, and let \(k>2\) be an integer. We present a routing scheme with a poly-logarithmic header size, that given a source s and a destination t at distance \(\varDelta \) from s, routes a message from s to t on a path whose length is \(O(k\varDelta +m^{1/k})\). The total space used by our routing scheme is \(mn^{O(1/\sqrt{\log n})}\), which is almost linear in the number of edges of the graph. We present also a routing scheme with \(n^{O(1/\sqrt{\log n})}\) header size, and the same stretch (up to constant factors). In this routing scheme, the routing table of every \(v\in V\) is at most \(kn^{O(1/\sqrt{\log n})}deg(v)\), where deg(v) is the degree of v in G. Our results are obtained by combining a general technique of Bernstein (2009), that was presented in the context of dynamic graph algorithms, with several new ideas and observations.     

Cost-efficient design of a quantum multiplier–accumulator unit

This paper proposes a cost-efficient quantum multiplier–accumulator unit. The paper also presents a fast multiplication algorithm and designs a novel quantum multiplier device based on the proposed algorithm with the optimum time complexity as multiplier is the major device of a multiplier–accumulator unit. We show that the proposed multiplication technique has time complexity \(O((3 {\hbox {log}}_{2}n)+1)\), whereas the best known existing technique has \(O(n{\hbox {log}}_{2} n)\), where n is the number of qubits. In addition, our design proposes three new quantum circuits: a circuit representing a quantum full-adder, a circuit known as quantum ANDing circuit, which performs the ANDing operation and a circuit presenting quantum accumulator. Moreover, the proposed quantum multiplier–accumulator unit is the first ever quantum multiplier–accumulator circuit in the literature till now, which has reduced garbage outputs and ancillary inputs to a great extent. The comparative study shows that the proposed quantum multiplier performs better than the existing multipliers in terms of depth, quantum gates, delays, area and power with the increasing number of qubits. Moreover, we design the proposed quantum multiplier–accumulator unit, which performs better than the existing ones in terms of hardware and delay complexities, e.g., the proposed (\(n\times n\))—qubit quantum multiplier–accumulator unit requires \(O(n^{2})\) hardware and \(O({\hbox {log}}_{2}n)\) delay complexities, whereas the best known existing quantum multiplier–accumulator unit requires \(O(n^{3})\) hardware and \(O((n-1)^{2} +1+n)\) delay complexities. In addition, the proposed design achieves an improvement of 13.04, 60.08 and 27.2% for \(4\times 4\), 7.87, 51.8 and 27.1% for \(8\times 8\), 4.24, 52.14 and 27% for \(16\times 16\), 2.19, 52.15 and 27.26% for \(32 \times 32\) and 0.78, 52.18 and 27.28% for \(128 \times 128\)-qubit multiplications over the best known existing approach in terms of number of quantum gates, ancillary inputs and garbage outputs, respectively. Moreover, on average, the proposed design gains an improvement of 5.62% in terms of area and power consumptions over the best known existing approach.     

Numerical Algorithm With High Spatial Accuracy for the Fractional Diffusion-Wave Equation With Neumann Boundary Conditions

A fourth-order compact algorithm is discussed for solving the time fractional diffusion-wave equation with Neumann boundary conditions. The \(L1\) discretization is applied for the time-fractional derivative and the compact difference approach for the spatial discretization. The unconditional stability and the global convergence of the compact difference scheme are proved rigorously, where a new inner product is introduced for the theoretical analysis. The convergence order is \(\mathcal{O }(\tau ^{3-\alpha }+h^4)\) in the maximum norm, where \(\tau \) is the temporal grid size and \(h\) is the spatial grid size, respectively. In addition, a Crank–Nicolson scheme is presented and the corresponding error estimates are also established. Meanwhile, a compact ADI difference scheme for solving two-dimensional case is derived and the global convergence order of \(\mathcal{O }(\tau ^{3-\alpha }+h_1^4+h_2^4)\) is given. Then extension to the case with Robin boundary conditions is also discussed. Finally, several numerical experiments are included to support the theoretical results, and some comparisons with the Crank–Nicolson scheme are presented to show the effectiveness of the compact scheme.     

Unconditionally Optimal Error Analysis of Crank–Nicolson Galerkin FEMs for a Strongly Nonlinear Parabolic System

In this paper, we present unconditionally optimal error estimates of linearized Crank–Nicolson Galerkin finite element methods for a strongly nonlinear parabolic system in \(\mathbb {R}^d\ (d=2,3)\). However, all previous works required certain time-step conditions that were dependent on the spatial mesh size. In order to overcome several entitative difficulties caused by the strong nonlinearity of the system, the proof takes two steps. First, by using a temporal-spatial error splitting argument and a new technique, optimal \(L^2\) error estimates of the numerical schemes can be obtained under the condition \(\tau \ge h\), where \(\tau \) denotes the time-step size and h is the spatial mesh size. Second, we obtain the boundedness of numerical solutions by mathematical induction and inverse inequality when \(\tau \le h\). Then, optimal \(L^2\) and \(H^1\) error estimates are proved in a different way for such case. Numerical results are given to illustrate our theoretical analyses.     

Discontinuous Galerkin Methods for a Stationary Navier–Stokes Problem with a Nonlinear Slip Boundary Condition of Friction Type

In this work, several discontinuous Galerkin (DG) methods are introduced and analyzed to solve a variational inequality from the stationary Navier–Stokes equations with a nonlinear slip boundary condition of friction type. Existence, uniqueness and stability of numerical solutions are shown for the DG methods. Error estimates are derived for the velocity in a broken \(H^1\)-norm and for the pressure in an \(L^2\)-norm, with the optimal convergence order when linear elements for the velocity and piecewise constants for the pressure are used. Numerical results are reported to demonstrate the theoretically predicted convergence orders, as well as the capability in capturing the discontinuity, the ability in handling the shear layers, the capacity in dealing with the advection-dominated problem, and the application to the general polygonal mesh of the DG methods.