Search

Differential-algebraic equations (DAEs) are widely used for modeling of dynamical systems. In numerical analysis of DAEs, consistent initialization and index reduction are important preprocessing prior to numerical integration. Existing DAE solvers commonly adopt structural preprocessing methods based on combinatorial optimization. Unfortunately, the structural methods fail if the DAE has numerical or symbolic cancellations. For such DAEs, methods have been proposed to modify them to other DAEs to which the structural methods are applicable, based on the combinatorial relaxation technique. Existing modification methods, however, work only for a class of DAEs that are linear or close to linear. This paper presents two new modification methods for nonlinear DAEs: the substitution method and the augmentation method. Both methods are based on the combinatorial relaxation approach and are applicable to a large class of nonlinear DAEs. The substitution method symbolically solves equations for some derivatives based on the implicit function theorem and substitutes the solution back into the system. Instead of solving equations, the augmentation method modifies DAEs by appending new variables and equations. The augmentation method has advantages that the equation solving is not needed and the sparsity of DAEs is retained. It is shown in numerical experiments that both methods, especially the augmentation method, successfully modify high-index DAEs that the DAE solver in MATLAB cannot handle., Comment: A preliminary version of this paper is to appear in Proceedings of the 44th International Symposium on Symbolic and Algebraic Computation (ISSAC 2019), Beijing, China, July 2019

The paper focuses on the analysis of the Euler projection Galerkin finite element method (FEM) for the dynamics of magnetization in ferromagnetic materials, described by the Landau–Lifshitz equation with the point-wise constraint $|{\textbf{m}}|=1$. The method is based on a simple sphere projection that projects the numerical solution onto a unit sphere at each time step, and the method has been used in many areas in the past several decades. However, error analysis for the commonly used method has not been done since the classical energy approach cannot be applied directly. In this paper we present an optimal $\textbf{L}^2$ error analysis of the backward Euler sphere projection method by using quadratic or higher order finite elements under a time step condition $\tau =O(\epsilon _0 h)$ with some small $\epsilon _0>0$. The analysis is based on more precise estimates of the extra error caused by the sphere projection in both $\textbf{L}^2$ and $\textbf{H}^1$ norms, and the classical estimate of dual norm. Numerical experiment is provided to confirm our theoretical analysis.

This paper is concerned with optimal error estimates and recovery technique of a classical mixed finite element method for the thermistor problem, which is governed by a parabolic/elliptic system with strong nonlinearity and coupling. The method is based on a popular combination of the lowest-order Raviart–Thomas mixed approximation for the electric potential/field $(\phi , \boldsymbol{\theta })$ and the linear Lagrange approximation for the temperature $u$. A common question is how the first-order approximation influences the accuracy of the second-order approximation to the temperature in such a strongly coupled system, while previous work only showed the first-order accuracy $O(h)$ for all three components in a traditional way. In this paper, we prove that the method produces the optimal second-order accuracy $O(h^2)$ for $u$ in the spatial direction, although the accuracy for the potential/field is in the order of $O(h)$. And more importantly, we propose a simple one-step recovery technique to obtain a new numerical electric potential/field of second-order accuracy. The analysis presented in this paper relies on an $H^{-1}$-norm estimate of the mixed finite element methods and analysis on a nonclassical elliptic map. We provide numerical experiments in both two- and three-dimensional spaces to confirm our theoretical analyses.

In this paper we revisit the problem of performing a $QZ$ step with a so-called ‘perfect shift’, which is an ‘exact’ eigenvalue of a given regular pencil $\lambda B-A$ in unreduced Hessenberg triangular form. In exact arithmetic, the $QZ$ step moves that eigenvalue to the bottom of the pencil, while the rest of the pencil is maintained in Hessenberg triangular form, which then yields a deflation of the given eigenvalue. But in finite precision the $QZ$ step gets ‘blurred’ and precludes the deflation of the given eigenvalue. In this paper we show that when we first compute the corresponding eigenvector to sufficient accuracy, then the $QZ$ step can be constructed using this eigenvector, so that the deflation is also obtained in finite precision. An important application of this technique is the computation of the index of a system of differential algebraic equations, since an exact deflation of the infinite eigenvalues is needed to impose correctly the algebraic constraints of such differential equations.

We propose and analyse numerical algorithms based on weighted least squares for the approximation of a real-valued function on a general bounded domain $\Omega \subset \mathbb{R}^d$. Given any $n$-dimensional approximation space $V_n \subset L^2(\Omega)$, the analysis in [6] shows the existence of stable and optimally converging weighted least-squares estimators, using a number of function evaluations $m$ of the order $n \log n$. When an $L^2(\Omega)$-orthonormal basis of $V_n$ is available in analytic form, such estimators can be constructed using the algorithms described in [6,Section 5]. If the basis also has product form, then these algorithms have computational complexity linear in $d$ and $m$. In this paper we show that, when $\Omega$ is an irregular domain such that the analytic form of an $L^2(\Omega)$-orthonormal basis is not available, stable and quasi-optimally weighted least-squares estimators can still be constructed from $V_n$, again with $m$ of the order $n \log n$, but using a suitable surrogate basis of $V_n$ orthonormal in a discrete sense. The computational cost for the calculation of the surrogate basis depends on the Christoffel function of $\Omega$ and $V_n$. Numerical results validating our analysis are presented., Comment: Version of the paper accepted for publication

In a series of papers published since 2017 the author introduced a simple alternative of the $n$-simplex type, to enhance the accuracy of approximations of second-order boundary value problems subject to Dirichlet boundary conditions, posed on smooth curved domains. This technique is based upon trial functions consisting of piecewise polynomials defined on straight-edged triangular or tetrahedral meshes, interpolating the Dirichlet boundary conditions at points of the true boundary. In contrast, the test functions are defined by the standard degrees of freedom associated with the underlying method for polytopic domains. While the mathematical analysis of the method for Lagrange and Hermite methods for two-dimensional second- and fourth-order problems was carried out in earlier paper by the author this paper is devoted to the study of the three-dimensional case. Well-posedness, uniform stability and optimal a priori error estimates in the energy norm are proved for a tetrahedron-based Lagrange family of finite elements. Novel error estimates in the $L^2$-norm, for the class of problems considered in this work, are also proved. A series of numerical examples illustrates the potential of the new technique. In particular, its superior accuracy at equivalent cost, as compared to the isoparametric technique, is highlighted.

Nonlinear elastic models are widely used to describe the elastic response of crystalline solids, for example, the well-known Cauchy–Born model. While the Cauchy–Born model only depends on the strain, effects of higher-order strain gradients are significant and higher-order continuum models are preferred in various applications such as defect dynamics and modeling of carbon nanotubes. In this paper we rigorously derive a higher-order nonlinear elasticity model for crystals from its atomistic description in one dimension. We show that, compared to the second-order accuracy of the Cauchy–Born model, the higher-order continuum model in this paper is of fourth-order accuracy with respect to the interatomic spacing in the thermal dynamic limit. In addition we discuss the key issues for the derivation of higher-order continuum models in more general cases. The theoretical convergence results are demonstrated by numerical experiments.

In this paper we analyze a class of trace finite element methods for the discretization of vector-Laplace equations. A key issue in the finite element discretization of such problems is the treatment of the constraint that the unknown vector field must be tangential to the surface (‘tangent condition’). We study three different natural techniques for treating the tangent condition, namely a consistent penalty method, a simpler inconsistent penalty method and a Lagrange multiplier method. The main goal of the paper is to present an analysis that reveals important properties of these three different techniques for treating the tangent constraint. A detailed error analysis is presented that takes the approximation of both the geometry of the surface and the solution of the partial differential equation into account. Error bounds in the energy norm are derived that show how the discretization error depends on relevant parameters such as the degree of the polynomials used for the approximation of the solution, the degree of the polynomials used for the approximation of the level set function that characterizes the surface, the penalty parameter and the degree of the polynomials used for the approximation of the Lagrange multiplier.

First, a new convergence analysis is given for the semidiscrete (finite elements in space) numerical method that is used in Le et al. (2016, Numerical solution of the time-fractional Fokker–Planck equation with general forcing. SIAM J. Numer. Anal.,54 1763–1784) to solve the time-fractional Fokker–Planck equation on a domain $\varOmega \times [0,T]$ with general forcing, i.e., where the forcing term is a function of both space and time. Stability and convergence are proved in a fractional norm that is stronger than the $L^2(\varOmega )$ norm used in the above paper. Furthermore, unlike the bounds proved in Le et al., the constant multipliers in our analysis do not blow up as the order of the fractional derivative $\alpha $ approaches the classical value of $1$. Secondly, for the semidiscrete (L1 scheme in time) method for the same Fokker–Planck problem, we present a new $L^2(\varOmega )$ convergence proof that avoids a flaw in the analysis of Le et al.’s paper for the semidiscrete (backward Euler scheme in time) method.

The augmented Lagrangian method (ALM) is fundamental in solving convex programming problems with linear constraints. The proximal version of ALM, which regularizes ALM’s subproblem over the primal variable at each iteration by an additional positive-definite quadratic proximal term, has been well studied in the literature. In this paper we show that it is not necessary to employ a positive-definite quadratic proximal term for the proximal ALM and the convergence can be still ensured if the positive definiteness is relaxed to indefiniteness by reducing the proximal parameter. An indefinite proximal version of the ALM is thus proposed for the generic setting of convex programming problems with linear constraints. We show that our relaxation is optimal in the sense that the proximal parameter cannot be further reduced. The consideration of indefinite proximal regularization is particularly meaningful for generating larger step sizes in solving ALM’s primal subproblems. When the model under discussion is separable in the sense that its objective function consists of finitely many additive function components without coupled variables, it is desired to decompose each ALM’s subproblem over the primal variable in Jacobian manner, replacing the original one by a sequence of easier and smaller decomposed subproblems, so that parallel computation can be applied. This full Jacobian splitting version of the ALM is known to be not necessarily convergent, and it has been studied in the literature that its convergence can be ensured if all the decomposed subproblems are further regularized by sufficiently large proximal terms. But how small the proximal parameter could be is still open. The other purpose of this paper is to show the smallest proximal parameter for the full Jacobian splitting version of ALM for solving multi-block separable convex minimization models.

In this paper we present a derivation of the surface Helmholtz decomposition, discuss its relation to the surface Hodge decomposition and derive a well-posed stream function formulation of a class of surface Stokes problems. We consider a $C^2$ connected (not necessarily simply connected) oriented hypersurface $\varGamma \subset \mathbb{R}^3$ without boundary. The surface gradient, divergence, curl and Laplace operators are defined in terms of the standard differential operators of the ambient Euclidean space $\mathbb{R}^3$. These representations are very convenient for the implementation of numerical methods for surface partial differential equations. We introduce surface $\mathbf H({\mathop{\rm div}}_{\varGamma})$ and $\mathbf H({\mathop{\rm curl}}_{\varGamma})$ spaces and derive useful properties of these spaces. A main result of the paper is the derivation of the Helmholtz decomposition, in terms of these surface differential operators, based on elementary differential calculus. As a corollary of this decomposition we obtain that for a simply connected surface to every tangential divergence-free velocity field there corresponds a unique scalar stream function. Using this result the variational form of the surface Stokes equation can be reformulated as a well-posed variational formulation of a fourth-order equation for the stream function. The latter can be rewritten as two coupled second-order equations, which form the basis for a finite element discretization. A particular finite element method is explained and the results of a numerical experiment with this method are presented.

This paper concerns the adaptive finite element method for elliptic Dirichlet boundary control problems in the energy space. The contribution of this paper is twofold. First, we rigorously derive efficient and reliable a posteriori error estimates for finite element approximations of Dirichlet boundary control problems. As a by-product, a priori error estimates are derived in a simple way by introducing appropriate auxiliary problems and establishing certain norm equivalence. Secondly, for the coupled elliptic partial differential system that resulted from the first-order optimality system, we prove that the sequence of adaptively generated discrete solutions including the control, the state and the adjoint state, guided by our newly derived a posteriori error indicators, converges to the true solution along with the convergence of the error estimators. We give some numerical results to confirm our theoretical findings.

Variants of the coordinate descent approach for minimizing a nonlinear function are distinguished in part by the order in which coordinates are considered for relaxation. Three common orderings are cyclic (CCD), in which we cycle through the components of $x$ in order; randomized (RCD), in which the component to update is selected randomly and independently at each iteration; and random-permutations cyclic (RPCD), which differs from CCD only in that a random permutation is applied to the variables at the start of each cycle. Known convergence guarantees are weaker for CCD and RPCD than for RCD, though in most practical cases, computational performance is similar among all these variants. There is a certain type of quadratic function for which CCD is significantly slower than for RCD; a recent paper by Sun & Ye (2016, Worst-case complexity of cyclic coordinate descent: $O(n^2)$ gap with randomized version. Technical Report. Stanford, CA: Department of Management Science and Engineering, Stanford University. arXiv:1604.07130) has explored the poor behavior of CCD on functions of this type. The RPCD approach performs well on these functions, even better than RCD in a certain regime. This paper explains the good behavior of RPCD with a tight analysis.

This paper provides a unified error analysis for nonconforming space discretizations of linear wave equations in the time domain. We propose a framework that studies wave equations as first-order evolution equations in Hilbert spaces and their space discretizations as differential equations in finite-dimensional Hilbert spaces. A lift operator maps the semidiscrete solution from the approximation space to the continuous space. Our main results are a priori error bounds in terms of interpolation, data and conformity errors of the method. Such error bounds are the key to the systematic derivation of convergence rates for a large class of problems. To show that this approach significantly eases the proof of new convergence rates, we apply it to an isoparametric finite element discretization of the wave equation with acoustic boundary conditions in a smooth domain. Moreover, our results reproduce known convergence rates for already investigated conforming and nonconforming space discretizations in a concise and unified way. The examples discussed in this paper comprise discontinuous Galerkin discretizations of Maxwell’s equations and finite elements with mass lumping for the acoustic wave equation.

In this paper we focus on the analysis of energy-preserving schemes for solving high-dimensional nonlinear Klein–Gordon equations. A novel energy-preserving scheme is developed based on the discrete gradient method and the Duhamel principle. The local error, global convergence and nonlinear stability of the new scheme are analysed in detail. Numerical experiments are implemented to compare with existing numerical methods in the literature, and the numerical results show the remarkable efficiency of the new energy-preserving scheme presented in this paper.

In this paper, we propose two low-order nonconforming finite element methods (FEMs) for the three-dimensional Stokes flow that generalize the nonconforming FEM of Kouhia & Stenberg (1995, A linear nonconforming finite element method for nearly incompressible elasticity and Stokes flow. Comput. Methods Appl. Mech. Eng, 124, 195–212). The finite element spaces proposed in this paper consist of two globally continuous components (one piecewise affine and one enriched component) and one component that is continuous at the midpoints of interior faces. We prove that the discrete Korn inequality and a discrete inf–sup condition hold uniformly in the mesh size and also for a nonempty Neumann boundary. Based on these two results, we show the well-posedness of the discrete problem. Two counterexamples prove that there is no direct generalization of the Kouhia–Stenberg FEM to three space dimensions: the finite element space with one nonconforming and two conforming piecewise affine components does not satisfy a discrete inf–sup condition with piecewise constant pressure approximations, while finite element functions with two nonconforming and one conforming component do not satisfy a discrete Korn inequality.

In this paper we suggest a method for transforming a vector subdivision scheme generating $C^{\ell}$ limits to another such scheme of the same dimension, generating $C^{\ell+1}$ limits. In scalar subdivision, it is well known that a scheme generating $C^{\ell}$ limit curves can be transformed to a new scheme producing $C^{\ell+1}$ limit curves by multiplying the scheme's symbol with the smoothing factor $\tfrac{z+1}{2}$. We extend this approach to vector and Hermite subdivision schemes, by manipulating symbols. The algorithms presented in this paper allow to construct vector (Hermite) subdivision schemes of arbitrarily high regularity from a convergent vector scheme (from a Hermite scheme whose Taylor scheme is convergent with limit functions of vanishing first component)., 28 pages, 4 figures. Corrected typos, updated contact information

n this paper, optimal control problems subject to a nonlinear scalar conservation law are studied. Such optimal control problems are challenging both at the continuous and at the discrete level since the control-to-state operator poses difficulties as it is, e.g., not differentiable. Therefore discretization of the underlying optimal control problem should be designed with care. Here the discretize-then-optimize approach is employed where first the full discretization of the objective function as well as the underlying PDE is considered. Then, the derivative of the reduced objective is obtained by using an adjoint calculus. In this paper total variation diminishing Runge-Kutta (TVD-RK) methods for the time discretization of such problems are studied. TVD-RK methods, also called strong stability preserving (SSP), are originally designed to preserve total variation of the discrete solution. It is proven in this paper that providing an SSP state scheme, is enough to ensure stability of the discrete adjoint. However requiring SSP for both discrete state and adjoint is too strong. Also approximation properties that the discrete adjoint inherits from the discretization of the state equation are studied. Moreover order conditions are derived. In addition, optimal choices with respect to CFL constant are discussed and numerical experiments are presented.

This paper analyses discontinuous Galerkin finite element methods (DGFEM) to approximate a regular solution to the von Karman equations defined on a polygonal domain. A discrete inf–sup condition sufficient for the stability of the discontinuous Galerkin discretization of a well-posed linear problem is established, and this allows the proof of local existence and uniqueness of a discrete solution to the nonlinear problem with a Banach fixed point theorem. The Newton scheme is locally second-order convergent and appears to be a robust solution strategy up to machine precision. A comprehensive a priori and a posteriori energy-norm error analysis relies on one sufficiently large stabilization parameter and sufficiently fine triangulations. In case the other stabilization parameter degenerates towards infinity, the DGFEM reduces to a novel C0-interior penalty method (IPDG). In contrast to the known C0-IPDG dueto Brenner et al., (2016, A C0 interior penalty method for a von Karman plate. Numer. Math., 1–30), the overall discrete formulation maintains symmetry of the trilinear form in the first two components—despite the general nonsymmetry of the discrete nonlinear problems. Moreover, a reliable and efficient a posteriori error analysis immediately follows for the DGFEM of this paper, while the different norms in the known C0-IPDG lead to complications with some nonresidual-type remaining terms. Numerical experiments confirm the best-approximation results and the equivalence of the error and the error estimators. A related adaptive mesh-refining algorithm leads to optimal empirical convergence rates for a nonconvex domain.

In this paper, we establish the optimal convergence for a second-order exponential-type integrator from Hofmanová & Schratz (2017, An exponential-type integrator for the KdV equation. Numer. Math., 136, 1117–1137) for solving the Korteweg–de Vries equation with rough initial data. The scheme is explicit and efficient to implement. By rigorous error analysis, we show that the scheme provides second-order accuracy in $H^\gamma $ for initial data in $H^{\gamma +4}$ for any $\gamma \geq 0$, where the regularity requirement is lower than for classical methods. The result is confirmed by numerical experiments, and comparisons are made with the Strang splitting scheme.

We consider the diffraction of an electromagnetic plane wave by a biperiodic structure. This paper is concerned with a numerical solution of the diffraction grating problem for three-dimensional Maxwell’s equations. Based on the Dirichlet-to-Neumann (DtN) operator, an equivalent boundary value problem is formulated in a bounded domain by using a transparent boundary condition. An a posteriori error estimate-based adaptive edge finite element method is developed for the variational problem with the truncated DtN operator. The estimate takes account of both the finite element approximation error and the truncation error of the DtN operator, where the former is used for local mesh refinements and the latter is shown to decay exponentially with respect to the truncation parameter. Numerical experiments are presented to demonstrate the competitive behaviour of the proposed method.

In recent work it has been established that deep neural networks (DNNs) are capable of approximating solutions to a large class of parabolic partial differential equations without incurring the curse of dimension. However, all this work has been restricted to problems formulated on the whole Euclidean domain. On the other hand, most problems in engineering and in the sciences are formulated on finite domains and subjected to boundary conditions. The present paper considers an important such model problem, namely the Poisson equation on a domain $D\subset \mathbb {R}^d$ subject to Dirichlet boundary conditions. It is shown that DNNs are capable of representing solutions of that problem without incurring the curse of dimension. The proofs are based on a probabilistic representation of the solution to the Poisson equation as well as a suitable sampling method.

This paper focuses on the stability analysis of a general class of linear multistep methods for decoupled forward–backward stochastic differential equations (FBSDEs). The general linear multistep methods we consider contain many well-known linear multistep methods from the ordinary differential equation framework, such as Adams, Nyström, Milne--Simpson and backward differentiation formula methods. Under the classical root condition, we prove that general linear multistep methods are mean-square (zero) stable for decoupled FBSDEs with generator function related to both y and z. Based on the stability result, we further establish a fundamental convergence theorem.

The Gaver–Stehfest algorithm is widely used for numerical inversion of the Laplace transform. In this paper we provide the first rigorous study of the rate of convergence of the Gaver–Stehfest algorithm. We prove that Gaver–Stehfest approximations converge exponentially fast if the target function is analytic in a neighbourhood of a point and they converge at a rate $o(n^{-k})$ if the target function is $(2k+3)$-times differentiable at a point.

In this paper we discuss spectral properties of operators associated with the least-squares finite-element approximation of elliptic partial differential equations. The convergence of the discrete eigenvalues and eigenfunctions towards the corresponding continuous eigenmodes is studied and analyzed with the help of appropriate $L^2$ error estimates. A priori and a posteriori estimates are proved.

This paper is concerned with the numerical approximation of stochastic mechanical systems with nonlinear holonomic constraints. Such systems are described by second order stochastic differential-algebraic equations involving an implicitly given Lagrange multiplier process. The explicit representation of the Lagrange multiplier leads to an underlying stochastic ordinary differential equation, the drift coefficient of which is typically not globally one-sided Lipschitz continuous. We investigate a half-explicit drift-truncated Euler scheme which fulfills the constraint exactly. Pathwise uniform $L_p$-convergence is established. The proof is based on a suitable decomposition of the discrete Lagrange multipliers and on norm estimates for the single components, enabling the verification of consistency, semi-stability and moment growth properties of the scheme. To the best of our knowledge, the presented result is the first strong convergence result for a constraint-preserving scheme in the considered setting., Comment: 39 pages, 1 figure

In this paper, we present a fully discrete analysis of an H(div)-conforming DG method with semi-implicit time-marching for the evolutionary incompressible Navier–Stokes equations. We use a semi-implicit time-discrete scheme in which the convection velocity is treated explicitly for the convection term. A stability analysis and a priori error estimates are given, in which the constants are independent of the negative powers of the viscosity. For inf-sup stable H(div)-conforming finite element pairs $BDM_k/P_{k-1}$ and $RT_k/P_k$, the rate of convergence $k+1/2$ is proved for the $L^2$ error of the velocity in the case of $\nu \leq C h$, where $k$ is the degree of the polynomials in the velocity approximation. In particular, for the inf-sup stable finite element pair $RT_k/P_k$, the convergence rate of the pressure is also $k+1/2$ when $\nu \leq C h$. The numerical experiments verify the analytical results.

Non divergence-free discretisations for the incompressible Stokes problem may suffer from a lack of pressure-robustness characterised by large discretisations errors due to irrotational forces in the momentum balance. This paper argues that also divergence-free virtual element methods (VEM) on polygonal meshes are not really pressure-robust as long as the right-hand side is not discretised in a careful manner. To be able to evaluate the right-hand side for the test functions, some explicit interpolation of the virtual test functions is needed that can be evaluated pointwise everywhere. The standard discretisation via an $L^2$-best approximation does not preserve the divergence and so destroys the orthogonality between divergence-free test functions and possibly eminent gradient forces in the right-hand side. To repair this orthogonality and restore pressure-robustness another divergence-preserving reconstruction is suggested based on Raviart-Thomas approximations on local subtriangulations of the polygons. All findings are proven theoretically and are demonstrated numerically in two dimensions. The construction is also interesting for hybrid high-order methods on polygonal or polyhedral meshes., 18 pages, 6 figures, 1 table, submitted to SINUM

The goal of the paper is to study the convergence of finite volume approximations of the Navier–Stokes–Fourier system describing the motion of compressible, viscous and heat-conducting fluids. The numerical flux uses upwinding with an additional numerical diffusion of order $\mathcal O(h^{ \varepsilon +1})$, $0

In this paper we study the numerical approximation of the stochastic Cahn–Hilliard–Navier–Stokes system on a bounded polygonal domain of $\mathbb{R}^{d}$, $d=2,3$. We propose and analyze an algorithm based on the finite element method and a semiimplicit Euler scheme in time for a fully discretization. We prove that the proposed numerical scheme satisfies the discrete mass conservative law, has finite energies and constructs a weak martingale solution of the stochastic Cahn–Hilliard–Navier–Stokes system when the discretization step (both in time and in space) tends to zero.

Observing that the recent developments of spatial discretization schemes based on recursive (product) quantization can be applied to a wide family of discrete time Markov chains, including all standard time discretization schemes of diffusion processes, we establish in this paper a generic strong error bound for such quantized schemes under a Lipschitz propagation assumption. We also establish a marginal weak error estimate that is entirely new to our best knowledge. As an illustration of their generality, we show how to recursively quantize the Euler scheme of a jump diffusion process, including details on the algorithmic aspects grid computation, transition weight computation, etc. Finally, we test the performances of the recursive quantization algorithm by pricing a European put option in a jump Merton model.

In this paper, we propose a two-level additive Schwarz domain decomposition preconditioner for the symmetric interior penalty Galerkin method for a second-order elliptic boundary value problem with highly heterogeneous coefficients. A specific feature of this preconditioner is that it is based on adaptively enriching its coarse space with functions created by solving generalized eigenvalue problems on thin patches covering the subdomain interfaces. It is shown that the condition number of the underlined preconditioned system is independent of the contrast if an adequate number of functions are used to enrich the coarse space. Numerical results are provided to confirm this claim.

This paper is concerned with a fast finite element method for the three-dimensional Poisson equation in infinite domains. Both the exterior problem and the strip-tail problem are considered. Exact Dirichlet-to-Neumann (DtN)-type artificial boundary conditions (ABCs) are derived to reduce the original infinite-domain problems to suitable truncated-domain problems. Based on the best relative Chebyshev approximation for the square-root function, a fast algorithm is developed to approximate exact ABCs. One remarkable advantage is that one need not compute the full eigensystem associated with the surface Laplacian operator on artificial boundaries. In addition, compared with the modal expansion method and the method based on Pad$\acute{\textrm{e}}$ approximation for the square-root function, the computational cost of the DtN mapping is further reduced. An error analysis is performed and numerical examples are presented to demonstrate the efficiency of the proposed method.

This paper is concerned with new discretization methods for the numerical treatment of elliptic partial differential equations. We derive an adaptive approximation scheme that is based on frames of quarkonial type, which can be interpreted as a wavelet version of $hp$ finite element dictionaries. These new frames are constructed from a finite set of functions via translation, dilation and multiplication by monomials. By using nonoverlapping domain decomposition ideas, we establish quarkonial frames on domains that can be decomposed into the union of parametric images of unit cubes. We also show that these new representation systems are stable in a certain range of Sobolev spaces. The construction is performed in such a way that, similar to the wavelet setting, the frame elements, the so-called quarklets, possess a certain number of vanishing moments. This enables us to generalize the basic building blocks of adaptive wavelet algorithms to the quarklet case. The applicability of the new approach is demonstrated by numerical experiments for the Poisson equation on $L$-shaped domains.

This paper investigates numerical schemes for stochastic differential equations driven by multi-dimensional fractional Brownian motions (fBms) with Hurst parameter $H\in (\frac 12,1)$. Based on the continuous dependence of numerical solutions on the driving noises, we propose the order conditions of Runge–Kutta methods for the strong convergence rate $2H-\frac 12$, which is the optimal strong convergence rate for approximating the Lévy area of fBms. We provide an alternative way to analyse the convergence rate of explicit schemes by adding ‘stage values’ such that the schemes are interpreted as Runge–Kutta methods. Taking advantage of this technique the strong convergence rate of simplified step-$N$ Euler schemes is obtained, which gives an answer to a conjecture in Deya et al. (2012) when $H\in (\frac 12,1)$. Numerical experiments verify the theoretical convergence rate.

Non-inf-sup-stable finite element approximations to the incompressible Navier–Stokes equations based on equal-order spaces for velocity and pressure are studied in this paper. To account for the violation of the discrete inf-sup condition, different types of symmetric pressure stabilization terms are considered. It is shown in the numerical analysis that these terms also improve stabilization of dominating convection in the following sense: error bounds with constants independent of inverse powers of the viscosity are derived. For proving the bound for the $L^2$ error of the pressure the choice of a suitable initial approximation for the velocity is essential.

How to choose the step size of gradient descent method has been a popular subject of research. In this paper we propose a modified limited memory steepest descent method (MLMSD). In each iteration we propose a selection rule to pick a unique step size from a candidate set, which is calculated by Fletcher’s limited memory steepest descent method (LMSD), instead of going through all the step sizes in a sweep, as in Fletcher’s original LMSD algorithm. MLMSD is motivated by an inexact super-linear convergence rate analysis. The R-linear convergence of MLMSD is proved for a strictly convex quadratic minimization problem. Numerical tests are presented to show that our algorithm is efficient and robust.

In this paper we propose a locally conservative, lowest-order staggered discontinuous Galerkin method for the coupled Stokes–Darcy problem on general quadrilateral and polygonal meshes. This model is composed of Stokes flow in the fluid region and Darcy flow in the porous media region, coupling together through mass conservation, balance of normal forces and the Beavers–Joseph–Saffman condition. Stability of the proposed method is proved. A new regularization operator is constructed to show the discrete trace inequality. Optimal convergence estimates for all the approximations covering low regularity are achieved. Numerical experiments are given to illustrate the performances of the proposed method. The numerical results indicate that the proposed method can be flexibly applied to rough grids such as the trapezoidal grid and $h$-perturbation grid.

Hybrid high-order (HHO) methods for elliptic diffusion problems have been originally formulated for loads in the Lebesgue space $L^2(\varOmega )$. In this paper we devise and analyse a variant thereof, which is defined for any load in the dual Sobolev space $H^{-1}(\varOmega )$. The main feature of the present variant is that its $H^1$-norm error can be bounded only in terms of the $H^1$-norm best error in a space of broken polynomials. We establish this estimate with the help of recent results on the quasi-optimality of nonconforming methods. We prove also an improved error bound in the $L^2$-norm by duality. Compared to previous works on quasi-optimal nonconforming methods the main novelties are that HHO methods handle pairs of unknowns and not a single function and, more crucially, that these methods employ a reconstruction that is one polynomial degree higher than the discrete unknowns. The proposed modification affects only the formulation of the discrete right-hand side. This is obtained by properly mapping discrete test functions into $H^1_0(\varOmega )$.

In this paper, we develop and study algorithms for approximately solving linear algebraic systems: ${{\mathcal{A}}}_h^\alpha u_h = f_h$, $ 0< \alpha

In this paper we introduce and analyze a fully mixed formulation for the nonlinear problem given by the coupling of the Stokes and Darcy–Forchheimer equations with the Beavers–Joseph–Saffman condition on the interface. This new approach yields non-Hilbert normed spaces and a twofold saddle point structure for the corresponding operator equation, whose continuous and discrete solvabilities are analyzed by means of a suitable abstract theory developed for this purpose. In particular, feasible choices of finite element subspaces include PEERS of the lowest order for the stress of the fluid, Raviart–Thomas of the lowest order for the Darcy velocity, piecewise constants for the pressures and continuous piecewise linear elements for the vorticity. An a priori error estimates and associated rates of convergence are derived, and several numerical results illustrating the good performance of the method are reported.

This paper deals with the numerical analysis of a system of second order in time partial differential equations modeling the vibrations of a coupled system that consists of an elastic solid in contact with an inviscid compressible fluid. We analyze a weak formulation with the unknowns in both media being the respective displacement fields. For its numerical approximation, we propose first a semidiscrete in space discretization based on standard Lagrangian elements in the solid and Raviart–Thomas elements in the fluid. We establish its well-posedness and derive error estimates in appropriate norms for the proposed scheme. In particular, we obtain an $\mathrm L^{\infty }(\mathrm L^2)$ optimal rate of convergence under minimal regularity assumptions of the solution, which are proved to hold for appropriate data of the problem. Then we consider a fully discrete approximation based on a family of implicit finite difference schemes in time, from which we obtain optimal error estimates for sufficiently smooth solutions. Finally, we report some numerical results, which allow us to assess the performance of the method. These results also show that the numerical solution is not polluted by spurious modes as is the case with other alternative approaches.

This paper deals with a priori pointwise error estimates for the finite element solution of boundary value problems with Neumann boundary conditions in polygonal domains. Due to the corners of the domain, the convergence rate of the numerical solutions can be lower than in the case of smooth domains. As a remedy, the use of local mesh refinement near the corners is considered. In order to prove quasi-optimal a priori error estimates, regularity results in weighted Sobolev spaces are exploited. This is the first work on the Neumann boundary value problem where both the regularity of the data is exactly specified and the sharp convergence order $h^{2} \lvert \ln h \rvert $ in the case of piecewise linear finite element approximations is obtained. As an extension we show the same rate for the approximate solution of a semilinear boundary value problem. The proof relies in this case on the supercloseness between the Ritz projection to the continuous solution and the finite element solution.

We provide the first significant extension of the unified transform (also known as the Fokas method) applied to elliptic boundary value problems, namely, we extend the method to curvilinear polygons and partial differential equations (PDEs) with variable coefficients. This is used to solve the generalized Dirichlet-to-Neumann map. The central component of the unified transform is the coupling of certain integral transforms of the given boundary data and of the unknown boundary values. This has become known as the global relation and, in the case of constant coefficient PDEs, simply links the Fourier transforms of the Dirichlet and Neumann boundary values. We extend the global relation to PDEs with variable coefficients and to domains with curved boundaries. Furthermore, we provide a natural choice of global relations for separable PDEs. These generalizations are numerically implemented using a method based on Chebyshev interpolation for efficient and accurate computation of the integral transforms that appear in the global relation. Extensive numerical examples are provided, demonstrating that the method presented in this paper is both accurate and fast, yielding exponential convergence for sufficiently smooth solutions. Furthermore, the method is straightforward to use, involving just the construction of a simple linear system from the integral transforms, and does not require knowledge of Green’s functions of the PDE. Details on the implementation are discussed at length.

In this paper, we present a multi-level mixed element scheme for the Helmholtz transmission eigenvalue problem on polygonal domains that are not necessarily able to be covered by rectangle grids. We first construct an equivalent linear mixed formulation of the transmission eigenvalue problem and then discretize it with Lagrangian finite elements of low regularities. The proposed scheme admits a natural nested discretization, based on which we construct a multi-level scheme. Optimal convergence rate and optimal com- putational cost can be obtained with the scheme., 20 pages, 13 figures

This paper considers a matrix optimization problem where the objective function is continuously differentiable and the constraints involve a semidefinite-box constraint and a rank constraint. We first replace the rank constraint by adding a non-Lipschitz penalty function in the objective and prove that this penalty problem is exact with respect to the original problem. Next, for the penalty problem we present a nonmonotone proximal gradient (NPG) algorithm whose subproblem can be solved by Newton’s method with globally quadratic convergence. We also prove the convergence of the NPG algorithm to a first-order stationary point of the penalty problem. Furthermore, based on the NPG algorithm, we propose an adaptive penalty method (APM) for solving the original problem. Finally, the efficiency of an APM is shown via numerical experiments for the sensor network localization problem and the nearest low-rank correlation matrix problem.

We study virtual element methods (VEMs) for solving the obstacle problem, which is a representative elliptic variational inequality of the first kind. VEMs can be regarded as a generalization of standard finite element methods with the addition of some suitable nonpolynomial functions, and the degrees of freedom are carefully chosen so that the stiffness matrix can be computed without actually computing the nonpolynomial functions. With this special design, VEMS can easily deal with complicated element geometries. In this paper we establish a priori error estimates of VEMs for the obstacle problem. We prove that the lowest-order ($k=1$) VEM achieves the optimal convergence order, and suboptimal order is obtained for the VEM with $k=2$. Two numerical examples are reported to show that VEM can work on very general polygonal elements, and the convergence orders in the $H^1$ norm agree well with the theoretical prediction.

In this paper, a partitioned time-stepping scheme for transient electromagnetically and thermally driven flow is analysed. The flow is modeled by coupled evolutionary magneto-hydrodynamic (MHD) equations with temperature-dependent coefficients. The partitioned scheme requires solving only one uncoupled MHD equation and one heat transfer equation per time step. It is based on Crank–Nicolson discretization in time and extrapolated treatment of the coupling and nonlinear terms such that skew-symmetry properties of the nonlinear terms are preserved. We prove that the proposed time-stepping scheme is unconditionally stable and derive error estimates for the fully discretized scheme using finite element spatial discretization in suitable norms. Numerical examples are presented that illustrate the accuracy and efficiency of the scheme.

Many applications of computational fluid dynamics require multiple simulations of a flow under different input conditions. In this paper, a numerical algorithm is developed to efficiently determine a set of such simulations in which the individually independent members of the set are subject to different viscosity coefficients, initial conditions, and/or body forces. The proposed scheme applied to the flow ensemble leads to need to solve a single linear system with multiple right-hand sides, and thus is computationally more efficient than solving for all the simulations separately. We show that the scheme is nonlinearly and long-term stable under certain conditions on the time-step size and a parameter deviation ratio. Rigorous numerical error estimate shows the scheme is of first-order accuracy in time and optimally accurate in space. Several numerical experiments are presented to illustrate the theoretical results., Comment: 20 pages, 3 figures