Search

Probabilistic domain decomposition (PDD) is an alternative paradigm for solving boundary value problems (BVPs) in parallel with excellent scalability properties, thanks to its reliance on stochastic representations of the BVP. However, there are cases when the latter is less numerically convenient, or unknown. Semilinear elliptic BVPs and the Helmholtz equation are prominent examples of either class. In this paper, we overcome this issue by designing suitable iterative schemes for either problem. These schemes not only retain the desirable properties of PDD but also are optimally suited for pathwise variance reduction, resulting in a systematic, nearly cost-free reduction of the statistical error through the iterations. Numerical tests carried out on the supercomputer Marconi100 are presented. [ABSTRACT FROM AUTHOR]

In this paper new innovative fourth order compact schemes for Robin and Neumann boundary conditions have been developed for boundary value problems of elliptic PDEs in two and three dimensions. Different from the traditional finite difference operator approach, which may not work for the flux type of boundary conditions, carefully designed undetermined coefficient methods are utilized in developing high order compact (HOC) schemes. The new methods not only can be utilized to design HOC schemes for flux type boundary conditions but also can be applied to general elliptic PDEs including Poisson, Helmholtz, diffusion-advection, and anisotropic equations with linear boundary conditions. In the newly developed HOC methods, the coefficient matrices are generally M-matrices, which guarantee the discrete maximum principle for well-posed problems, so the convergence of the HOC methods. The developed HOC methods are versatile and can cover most HOC schemes in the literature. The HOC methods for Robin boundary conditions and for anisotropic diffusion and advection equations with Robin or even Dirichlet boundary conditions are likely the first ones that have ever been developed. With the help of pseudoinverse or singular value decomposition, we have also observed that the developed HOC methods usually have smaller error constants compared with traditional HOC methods when applicable. Nontrivial examples with large wave numbers and oscillatory solutions are presented to confirm the performance of the new HOC methods. [ABSTRACT FROM AUTHOR]

The a posteriori error estimate and greedy algorithms play central roles in the reduced basis method (RBM). In [M. Yano, Comput. Methods Appl. Mech. Engrg., 287 (2015), pp. 290--309; ESAIM Math. Model. Numer. Anal., 50 (2016), pp. 163--185; SIAM J. Sci. Comput., 40 (2018), pp. A388--A420], several versions of RBMs based on exact error certifications and greedy algorithms with spatio-parameter adaptivities are developed. In this paper, with the parametric symmetric coercive elliptic boundary value problem as an example of the primal-dual variational problems satisfying the strong duality, we develop primal-dual RBMs (PD-RBM) with robust true error certifications and discuss three versions of greedy algorithms to balance the finite element error, the exact RB error, and the adaptive mesh refinements. For a class of convex minimization variational problems which has corresponding dual problems satisfying the strong duality, the primal-dual gap between the primal and dual functionals can be used as an a posteriori error estimator. This primal-dual gap error estimator is robust with respect to the parameters of the problem, and it can be used for both mesh refinements of finite element methods and the true RB error certification. With the help of an integration by parts formula, the primal-dual variational theory is developed for the symmetric coercive elliptic boundary value problems with nonhomogeneous boundary conditions by both the conjugate function and Lagrangian theories. A generalized Prager--Synge identity, which is the primal-dual gap error representation for this specific problem, is developed. RBMs for both the primal and dual problems with robust error estimates are developed. The dual variational problem often can be viewed as a constraint optimization problem. In the paper, different from the standard saddlepoint finite element approximation, the dual RBM is treated as a Galerkin projection by constructing RB spaces satisfying the homogeneous constraint. Inspired by the greedy algorithm with spatioparameter adaptivity of [M. Yano, SIAM J. Sci. Comput., 40 (2018), pp. A388--A420], adaptive balanced greedy algorithms with primal-dual finite element and RB error estimators are discussed. Numerical tests are presented to test the PD-RBM with adaptive balanced greedy algorithms. [ABSTRACT FROM AUTHOR]

In this paper, a parallel probabilistic method using the Feynman--Kac formula of killed Brownian motions is proposed to solve the mixed boundary value problems (BVPs) of 3D Laplace equations. To avoid using reflecting Brownian motions and the calculation of their local time L(t) in the Feynman--Kac representation of solutions for Neumann and Robin BVPs, the proposed method uses an iterative approach to approximate the solutions where each iteration will use the Feynman--Kac formula to solve a pure Dirichlet problem, thus only involving killed Brownian motions. First, the boundary of the domain is decomposed with overlapping local patches formed by the intersection of hemispheres superimposed on the domain boundary. The iteration starts with an arbitrary initial guess for the Dirichlet data on the Neumann and Robin boundaries; then, using the Feynman--Kac formula for a pure Dirichlet problem with the current available Dirichlet data on the whole boundary, the solution over the hemispheres can be obtained by the Feynman--Kac formula for the killed Brownian motions, sampled by a walk-on-spheres (WOS) algorithm. Second, by solving a local boundary integral equation (BIE) over each hemisphere and a local patch on the domain boundary, the Dirichlet data on the Neumann and Robin boundaries can be updated. By continuing this process, the proposed hybrid probabilistic and deterministic BIE-WOS method gives a highly parallel algorithm for the global solution of any mixed-type BVPs of the Laplace equations. Numerical results of various mixed interior and exterior BVPs demonstrate the parallel efficiency and accuracy of the proposed method. [ABSTRACT FROM AUTHOR]

In this paper, we propose a numerical method to solve the classic L²-optimal transport problem. Our algorithm is based on the use of multiple shooting, in combination with a continuation procedure, to solve the boundary value problem associated to the transport problem. Based on the viewpoint of Wasserstein Hamiltonian flow with initial and target densities, our algorithm reflects the Hamiltonian structure of the underlying problem and exploits it in the numerical discretization. Several numerical examples are presented to illustrate the performance of the method. [ABSTRACT FROM AUTHOR]

In anisotropic media, the standard perfectly matched layer (PML) technique suffers irrevocable instability in terminating the unbounded problem domains. It remains an open question whether a stable PML-like absorbing boundary condition exists. For wave scattering in a layered orthotropic medium, this question is affirmatively answered for the first time in this paper. In each orthotropic medium, the permittivity tensor uniquely determines a change of coordinates that transforms the governing anisotropic Helmholtz equation into an isotropic Helmholtz equation in the new coordinate system. This leads us to propose a novel Sommerfeld radiation condition (SRC) to rigorously characterize outgoing waves in the layered orthotropic medium. Naturally, the SRC motivates a regionalized PML (RPML) to truncate the scattering problem, in the sense that a standard PML is set up in the new coordinate system in each orthotropic region. It is revealed that the RPML is unconditionally stable compared with the unstable uniaxial PML. A high-accuracy boundary integral equation (BIE) method is developed to solve the resulting boundary value problem. Numerical experiments are carried out to validate the stability of the RPML and the accuracy of the BIE method, showing exponentially decaying truncation errors as the RPML parameters increase. [ABSTRACT FROM AUTHOR]

In this paper, we address the temporal energy growth associated with numerical approximations of the perfectly matched layer (PML) for Maxwell's equations in first order form. In the literature, several studies have shown that a numerical method which is stable in the absence of the PML can become unstable when the PML is introduced. We demonstrate in this paper that this instability can be directly related to numerical treatment of boundary conditions in the PML. First, at the continuous level, we establish the stability of the constant coefficient initial boundary value problem for the PML. To enable the construction of stable numerical boundary procedures, we derive energy estimates for the PML in the Laplace space. Second, we develop a high order accurate and stable numerical approximation for the PML using summation-by-parts finite difference operators to approximate spatial derivatives and weak enforcement of boundary conditions using penalties. By constructing analogous discrete energy estimates we show discrete stability and convergence of the numerical method. Numerical experiments verify the theoretical results. [ABSTRACT FROM AUTHOR]

In a series of recent papers [B. Adcock, A. C. Hansen, and C. Poon, Appl. Comput. Harmon. Anal., 36 (2014), pp. 387-415; B. Adcock, M. Gataric, and A. C. Hansen, SIAM J. Imaging Sci., 7 (2014), pp. 1690-1723; Adcock et al., SIAM J. Math. Anal., 47 (2015), pp. 1196-1233], it was shown that one can optimally recover the wavelet coefficients of an unknown compactly supported function from pointwise evaluations of its Fourier transform via the method of generalized sampling. While these papers focused on the optimality of generalized sampling in terms of its stability and error bounds, the current paper explains how this optimal method can be implemented to yield a computationally efficient algorithm. In particular, we show that generalized sampling has a computational complexity of O (M(N) logN) when recovering the first N boundary-corrected wavelet coefficients of an unknown compactly supported function from M(N) Fourier samples. Therefore, due to the linear correspondences between the number of samples M and the number of coefficients N shown previously, generalized sampling offers a computationally optimal way of recovering wavelet coefficients from Fourier data. [ABSTRACT FROM AUTHOR]

This paper is concerned with the boundary integral equation method for solving the exterior Neumann boundary value problem of dynamic poroelasticity in two dimensions. The main contribution of this work consists of two aspescts: the proposal of a novel regularized boundary integral equation, and the presentation of new regularized formulations of the strongly singular and hyper-singular boundary integral operators. First, turning to the spectral properties of the doublelayer operator and the corresponding Calder\'on relation of the poroelasticity, we propose the novel low-GMRES-iteration integral equation whose eigenvalues are bounded away from zero and infinity. Second, with the help of the classical idea of using G\"unter derivatives and integration-by-parts, we reformulate the strongly singular and hyper-singular integral operators into combinations of the weakly singular operators and the tangential derivatives. The accuracy and efficiency of the proposed methodology are demonstrated through several numerical examples. [ABSTRACT FROM AUTHOR]

In many applications of practical interest, solutions of partial differential equation models arise as critical points of an underlying (energy) functional. If such solutions are saddle points, rather than being maxima or minima, then the theoretical framework is nonstandard, and the development of suitable numerical approximation procedures turns out to be highly challenging. In this paper, our aim is to present an iterative discretization methodology for the numerical solution of nonlinear variational problems with multiple (saddle point) solutions. In contrast to traditional numerical approximation schemes, which typically fail in such situations, the key idea of the current work is to employ a simultaneous interplay of a previously dev eloped local minimax approach and adaptive Galerkin discretizations. We thereby derive an adaptive local minimax Galerkin (LMMG) method, which combines the search for saddle point solutions and their approximation in finitedimensional spaces in a highly effective way. Under certain assumptions, we will prove that the generated sequence of approximate solutions converges to the solution set of the variational problem. This general framework will be applied to the specific context of finite element discretizations of (singularly perturbed) semilinear elliptic boundary value problems, and a series of numerical experiments will be presented. [ABSTRACT FROM AUTHOR]

In this paper, a new localized radial basis function (RBF) method based on partition of unity (PU) is proposed for solving boundary and initial-boundary value problems. The new method benefits from a direct discretization approach and is called the "direct RBF partition of unity (D-RBF-PU)" method. Thanks to avoiding all derivatives of PU weight functions as well as all lower derivatives of local approximants, the new method is faster and simpler than the standard RBF-PU method. Besides, the discontinuous PU weight functions can now be utilized to develop the method in a more efficient and less expensive way. Alternatively, the new method is an RBF-generated finite difference (RBF-FD) method in a PU setting which is much faster and in some situations more accurate than the original RBF-FD. The polyharmonic splines are used for local approximations, and the error and stability issues are considered. Some numerical experiments on irregular two- and three-dimensional domains, as well as cost comparison tests, are performed to support the theoretical analysis and to show the efficiency of the new method. [ABSTRACT FROM AUTHOR]

In this paper, two finite difference/element approaches for the time-fractional subdiffusion equation with Dirichlet boundary conditions are developed, in which the time direction is approximated by the fractional linear multistep method and the space direction is approximated by the finite element method. The two methods are unconditionally stable and convergent of order O(τq +hr+1) in the L² norm, where q = 2-β or 2 when the analytical solution to the subdiffusion equation is sufficiently smooth, β (0 < β < 1) is the order of the fractional derivative, τ and h are the step sizes in time and space, respectively, and r is the degree of the polynomial space. The corresponding schemes for the subdiffusion equation with Neumann boundary conditions are presented as well, where the stability and convergence are shown. Numerical examples are provided to verify the theoretical analysis. Comparisons between the algorithms derived in this paper and the existing algorithms are given, which show that our numerical schemes exhibit better performances than the existing ones. [ABSTRACT FROM AUTHOR]

We introduce a new higher order scheme for computing a tangentially stabilized curve shortening flow with a driving force represented by an intrinsic partial differential equation for an evolving curve position vector. Our new scheme is a combination of the explicit forward Euler and the fully implicit backward Euler schemes. At any discrete time step, the solution is found efficiently using a few semi-implicit iterations. Basic properties of the new scheme axe proved in the paper, and its precision is tested by comparing the results with known analytical solutions. For any choice of the time step, the new higher order scheme gives exact radius of evolving uniformly discretized circles in case of flow by curvature and in case of rotation by a constant tangential velocity. Such properties do not hold for other schemes solving flow by mean curvature like the classical explicit, semi-implicit, or fully implicit schemes. In general, the scheme is second order accurate, which is shown by comparing a numerically evolving encompassed area with known analytical expression. The behavior of the scheme is discussed on representative examples, and its advantages with respect to the balance between CPU time and precision are shown. [ABSTRACT FROM AUTHOR]

This paper introduces a fast and numerically stable algorithm for the solution of fourth-order linear boundary value problems on an interval. This type of equation arises in a variety of settings in physics and signal processing. Our method reformulates the equation as a collection of second-kind integral equations defined on local subdomains. Each such equation can be stably discretized and solved. The boundary values of these local solutions are matched by solving a banded linear system. The method of deferred corrections is then used to increase the accuracy of the scheme. Deferred corrections require applying the integral operator to a function on the entire domain, for which we provide an algorithm with linear cost. We illustrate the performance of our method on several numerical examples. [ABSTRACT FROM AUTHOR]

In the present paper we describe a class of algorithms for the solution of Laplace's equation on polygonal domains with Dirichlet and Neumann boundary conditions. It is well known that in such cases the solutions may have singularities near the corners, which poses a challenge for many existing methods. We present a high-order solver for Laplace's equation on polygonal domains requiring relatively few degrees of freedom to resolve the behavior near corners accurately. Our approach is based on the observation that if the boundary data are smooth on each edge of the polygon then in the vicinity of each corner the solution to the corresponding boundary integral equation has an expansion in terms of certain (analytically available) singular powers. We construct a set of discretization nodes and weights which accurately integrate and interpolate these singular powers, enabling the construction of high-order Nyström schemes requiring relatively few discretization nodes. Our results are illustrated with several numerical examples. [ABSTRACT FROM AUTHOR]

Recent work has given rise to a novel and simple numerical technique for solving elliptic boundary value problems formulated in convex polygons in two dimensions. The method, based on the unified transform, involves expanding the unknown boundary values in a Legendre basis and determining the expansion coeffcients by evaluating the so-called global relation at appropriate points in the complex Fourier plane (spectral collocation). In this paper we provide a significant advancement of this numerical technique by providing a fast and effcient method to evaluate the solution in the domain interior. The use of a Legendre basis allows the relevant integrals to be computed effciently and accurately using Chebyshev interpolation, even for large degree. For the particular case of the Laplace equation this allows an explicit expansion in the domain interior in terms of hypergeometric functions. Evaluation in the interior is found to converge more rapidly than the approximation of the unknown boundary values, allowing accurate approximation of solutions with weak corner singularities. For stronger singularities, the method can be combined with global singular functions for rapid convergence. Numerical examples are provided, showing that the method compares well against standard spectral methods and opens up the possibility of applying the method to general curvilinear domains. [ABSTRACT FROM AUTHOR]

In this paper, we study the asymptotic analysis and numerical method for singularly perturbed eigenvalue problems (SPEPs). We first make a close asymptotic analysis of the SPEP and prove that, for piecewise smooth potential functions V (x), the eigenvalues converge to the minimum value of V (x), and the eigenfunctions are concentrated in the immediate vicinity of the minimal points of V (x) as ∞ 0+. Then we propose some new schemes based on the tailored finite point method (TFPM) for numerical solutions of SPEPs with higher accuracy. Our numerical examples verify our theory and show the feasibility and efficiency of our TFPM. [ABSTRACT FROM AUTHOR]

We are concerned with the numerical solution of a nonlocal wave equation in an infinite two-dimensional space. The contribution of this paper is the derivation of an absorbing boundary condition which allows the wave field defined on the finite computational domain to retain the same feature as that defined on the original infinite domain. We resort to the idea of a first-kind integral equation method and develop a solution formulation in terms of a potential summation on a surrounding ghost region. This new formulation can be taken as an absorbing boundary condition of generalized Dirichlet-to-Dirichlet type. The accuracy and effectiveness of our approach are illustrated by some numerical examples. [ABSTRACT FROM AUTHOR]

This work is on the numerical approximation of incoming solutions to Maxwell's equations with dissipative boundary conditions, whose energy decays exponentially with time. Such solutions are called asymptotically disappearing (ADS) and they play an important role in inverse back-scattering problems. The existence of ADS is a difficult mathematical problem. For the exterior of a sphere, such solutions have been constructed analytically by Colombini, Petkov, and Rauch [Proc. Amer. Math. Soc., 139 (2011), pp. 2163-2173] by specifying appropriate initial conditions. However, for general domains of practical interest (such as Lipschitz polyhedra), the existence of such solutions is not evident. This paper considers a finite-element approximation of Maxwell's equations in the exterior of a polyhedron, whose boundary approximates the sphere. Standard Nédélec--Raviart--Thomas elements are used with a Crank--Nicolson scheme to approximate the electric and magnetic fields. Discrete initial conditions interpolating the ones chosen by Colombini, Petkov, and Rauch are modified so that they are (weakly) divergence-free. We prove that with such initial conditions, the approximation to the electric field is weakly divergence-free for all time. Finally, we show numerically that the finite-element approximations of the ADS also approximates this exponential decay (quadratically) with time when the mesh size and the time step become small. [ABSTRACT FROM AUTHOR]

In this paper, we present an accurate, efficient and stable numerical method for the incompressible Navier-Stokes equations (NSEs). The method is based on (1) an equivalent pressure Poisson equation formulation of the NSE with proper pressure boundary conditions, which facilitates the design of high-order and stable numerical methods, and (2) the Krylov deferred correction (KDC) accelerated method of lines transpose (MOLT), which is very stable, efficient, and of arbitrary order in time. Numerical tests with known exact solutions in three dimensions show that the new method is spectrally accurate in time, and a numerical order of convergence 9 was observed. Two-dimensional computational results of flow past a cylinder and flow in a bifurcated tube are also reported. [ABSTRACT FROM AUTHOR]

This paper addresses the construction of different families of absorbing boundary conditions for the one- and two-dimensional Schrödinger equations with a general variable nonlinear potential. Various semidiscrete time schemes are built for the associated initial boundary value problems. Finally, some numerical simulations give a comparison of the various absorbing boundary conditions and associated schemes to analyze their accuracy and efficiency. [ABSTRACT FROM AUTHOR]

In this paper we address the error growth in time for hyperbolic problems in first order form in bounded domains. The energy method is used to study when an error growth or a fixed error bound is obtained. It is shown that the choice of boundary procedure is a crucial point. Numerical experiments corroborate the theoretical findings. [ABSTRACT FROM AUTHOR]

Numerical simulation of the solution to the sine-Gordon equation on the whole real axis is considered in this paper. Based on nonlinear spectral analysis, exact nonreflecting boundary conditions are derived at two artificially introduced boundary points. Then a numerical scheme of second order is proposed to approximate the solution. In the end, some numerical examples are provided to demonstrate the effectiveness of the proposed scheme. [ABSTRACT FROM AUTHOR]

This paper is concerned with the development of robust and efficient iterative techniques for solving the coupled systems of linear equations which arise from a mixed finite element discretization of dissipative acoustic equations. The system of linearized, time-harmonic acoustic equations we consider is based on the Navier-Stokes, continuity, and energy equations accompanied by suitable equations of state. In this paper, we focus on the design of preconditioners that can be used in combination with the Krylov subspace methods. Our starting point is the Schur complement formulation of the mixed problem. One of the difficulties encountered with the approach we consider is generally associated with the issue of defining boundary conditions for the associated Schur complement problem. To derive suitable boundary conditions for this problem, we first consider an iterative solution method which is based on the direct use of the Schur complement formulation. This development naturally leads us to address the key question of what boundary conditions should be adjoined to the Schur complement preconditioner. We give numerical examples which demonstrate the effectiveness and robustness of the proposed preconditioner. The examples include both twoand three-dimensional problems with boundary conditions of various types. [ABSTRACT FROM AUTHOR]

This paper is concerned with the numerical approximation of the solutions of a macroscopic model for the description of the flow of pedestrians. Solutions of the associated Riemann problem are known to be possibly nonclassical in the sense that the underlying discontinuities may well violate the Lax inequalities, which makes their numerical approximation very sensitive. This study proposes to extend the transport-equilibrium strategy proposed in [C. Chalons, C. R. Acad. Sci. Paris Ser. I, 342 (2006), pp. 623-626] and [C. Chalons, Transport-equilibrium schemes for computing nonclassical shocks. I. Scalar conservation laws, submitted] for computing the nonclassical solutions of scalar conservation laws with either a concave-convex or a convex-concave flux function and supplemented with an invertible kinetic function. These strong properties are not fulfilled in the present setting, since the flow function admits several inflection points and the kinetic function is not invertible. We nevertheless succeed in obtaining an efficient numerical scheme. Numerical evidence is provided. [ABSTRACT FROM AUTHOR]

Creating in parallel guaranteed quality large unstructured meshes is a challenging problem. Parallel mesh generation procedures decompose the original mesh generation problem into smaller subproblems that can be solved in parallel. The subproblems can be treated as either completely or partially coupled, or they can be treated as completely decoupled. In this paper we present a parallel guaranteed quality Delaunay method for 2-dimensional domains which is based on the complete decoupling of the subproblems. As a result the method eliminates the communication and the synchronization during the meshing of the subproblems. Moreover, it achieves 100% code reuse of existing, fine-tuned, and well-tested sequential mesh generators. The approach we describe in this paper presents for the first time an effective way to create in parallel guaranteed quality meshes with billions of elements in a few hundreds of seconds, and at the same time demonstrates that these meshes can be generated in an efficient and scalable way. Our performance data indicate superlinear speedups. [ABSTRACT FROM AUTHOR]

Optimization problems constrained by nonlinear partial differential equations have been the focus of intense research in scientific computing lately. Current methods for the parallel numerical solution of such problems involve sequential quadratic programming (SQP), with either reduced or full space approaches. In this paper we propose and investigate a class of parallel full space SQP Lagrange-Newton-Krylov-Schwarz (LNKSz) algorithms. In LNKSz, a Lagrangian functional is formed and differentiated to obtain a Karush-Kuhn-Tucker (KKT) system of nonlinear equations. An inexact Newton method with line search is then applied. At each Newton iteration the linearized KKT system is solved with a Schwarz preconditioned Krylov subspace method. We apply LNKSz to the parallel numerical solution of some boundary control problems of two-dimensional incompressible Navier-Stokes equations. Numerical results are reported for different combinations of Reynolds number, mesh size, and number of parallel processors. We also compare the application of the LNKSz method to flow control problems against the application of the Newton-Krylov-Schwarz (NKSz) method to flow simulation problems. [ABSTRACT FROM AUTHOR]

Iterated deferred correction is a widely used approach to the numerical solution of first order systems of nonlinear two-point boundary value problems. Normally the orders of accuracy of the various methods used in a deferred correction scheme differ by 2, and, as a direct result, each time a deferred correction is applied the order of the overall scheme is increased by a maximum of 2. In this paper we consider the construction of mono-implicit Runge­Kutta (MIRK) methods where an increase of four orders of accuracy is obtained for each deferred correction. We develop a very powerful yet rather straightforward theory which allows us to identify the appropriate Runge­Kutta formulae for inclusion in such schemes. In particular, we will focus on the construction of pairs of MIRK formulae of order 4 and 8 which will allow this superconvergence to be realized. We will further show that it is possible to derive formulae of this type for which high order interpolants and accurate error estimates are readily available. [ABSTRACT FROM AUTHOR]

Two first-order system least squares (FOSLS) methods based on L2 norms are applied to various boundary value problems of planar linear elasticity. Both use finite element discretization and multigrid solution methods. They are two-stage algorithms that solve first for the displacement flux variable (the gradient of displacement, which easily yields the deformation and stress variables), then for the displacement variable itself. As a complement to a companion theoretical paper, this paper focuses on numerical results, including finite element accuracy and multigrid convergence estimates that confirm uniform optimal performance-even as the material tends to the incompressible limit. [ABSTRACT FROM AUTHOR]

A long-standing open question associated with the use of collocation methods for boundary value ordinary differential equations is concerned with the development of a high order continuous solution approximation to augment the high order discrete solution approximation, obtained at the mesh points which subdivide the problem interval. It is well known that the use of collocation at Gauss points leads to solution approximations at the mesh points for which the global error is O(h2k), where k is the number of collocation points used per subinterval and h is the subinterval size. This discrete solution is said to be superconvergent. The collocation solution also yields a C0 continuous solution approximation that has a global error of O(hk+1). In this paper, we show how to efficiently augment the superconvergent discrete collocation solution to obtain C¹ continuous "superconvergent" interpolants whose global errors are O(h2k). The key ideas are to use the theoretical framework of continuous Runge{Kutta schemes and to augment the collocation solution with inexpensive monoimplicit Runge-Kutta stages. Specific schemes are derived for k = 1, 2, 3, and 4. Numerical results are provided to support the theoretical analysis. [ABSTRACT FROM AUTHOR]

We consider an unbounded steady-state flow of viscous fluid over a three-dimensional finite body or configuration of bodies. For the purpose of solving this flow numerically, we discretize the governing equations (Navier{Stokes) on a finite-difference grid. Prior to the discretization, we obviously need to truncate the original unbounded domain by introducing an artificial computational boundary at a finite distance from the body; otherwise, the number of discrete variables will not be finite. This artificial boundary is typically the external boundary of the domain covered by the grid. The flow problem (both continuous and discretized) formulated on the finite computational domain is clearly subdefinite unless supplemented by some artificial boundary conditions (ABCs) at the external computational boundary. In this paper, we present an innovative approach to constructing highly accurate ABCs for three-dimensional flow computations. The approach extends our previous technique developed for the two-dimensional case; it employs the finite-difference counterparts to Calderón's pseudodifferential boundary projections calculated in the framework of the difference potentials method (DPM) of Ryaben'kii. The resulting ABCs appear spatially nonlocal but are particularly easy to implement along with the existing flow solvers. The new boundary conditions have been successfully combined with the NASA-developed production code TLNS3D and used for the analysis of wing-shaped configurations in subsonic and transonic flow regimes. As demonstrated by the computational experiments and comparison with the standard local methods, the DPM-based ABCs allow one to greatly reduce the size of the computational domain while still maintaining high accuracy of the numerical solution. Moreover, they may provide for a noticeable speedup of multigrid convergence. [ABSTRACT FROM AUTHOR]

In this paper, we design robust and efficient linear solvers for the numerical approximation of solutions to Maxwell's equations with dissipative boundary conditions. We consider a structure-preserving finite-element approximation with standard Nédélec-Raviart-Thomas elements in space and a Crank-Nicolson scheme in time to approximate the electric and magnetic fields. We focus on two types of block preconditioners. The first type is based on the well-posedness results of the discrete problem. The second uses an exact block factorization of the linear system, for which the structure-preserving discretization yields sparse Schur complements. We prove robustness and optimality of these block preconditioners and provide supporting numerical tests. [ABSTRACT FROM AUTHOR]

We consider the identification of spatially distributed parameters under H¹ regularization. Solving the associated minimization problem by Gauss-Newton iteration results in linearized problems to be solved in each step that can be cast as boundary value problems involving a low-rank modification of the Laplacian. Using an algebraic multigrid as a fast Laplace solver, the Sherman-Morrison-Woodbury formula can be employed to construct a preconditioner for these linear problems which exhibits excellent scaling w.r.t. the relevant problem parameters. We first develop this approach in the functional setting, thus obtaining a consistent methodology for selecting boundary conditions that arise from the H¹ regularization. We then construct a method for solving the discrete linear systems based on combining any fast Poisson solver with the Woodbury formula. The efficacy of this method is then demonstrated with scaling experiments. These are carried out for a common nonlinear parameter identification problem arising in electrical resistivity tomography. [ABSTRACT FROM AUTHOR]

In this work, we show that solvers of elliptic boundary value problems in d dimensions can be approximated to accuracy ε from only O(log(N) logd (N/ε)) matrix-vector products with carefully chosen vectors (right-hand sides). The solver is only accessed as a black box, and the under-lying operator may be unknown and of an arbitrarily high order. Our algorithm (1) has complexity O(N log²(N) log2d (N/ε)) and represents the solution operator as a sparse Cholesky factorization with O(N log (N) logd (N/ε)) nonzero entries, (2) allows for embarrassingly parallel evaluation of the solution operator and the computation of its log-determinant, (3) allows for O(log (N) logd (N/ε)) complexity computation of individual entries of the matrix representation of the solver that, in turn, enables its recompression to an O(N logd (N/ε)) complexity representation. As a byproduct, our compression scheme produces a homogenized solution operator with near-optimal approximation accuracy. By polynomial approximation, we can also approximate the continuous Green's function (in operator and Hilbert-Schmidt norm) to accuracy ε from O(log1+d(ε-1)) solutions of the PDE. We include rigorous proofs of these results. To the best of our knowledge, our algorithm achieves the best known trade-off between accuracy ε and the number of required matrix-vector products. [ABSTRACT FROM AUTHOR]

This paper describes a exible C++ software framework, called GRINS, for simulating complex multiphysics systems of partial differential equations using the finite element method. GRINS is designed to facilitate the rapid development and reuse of finite element formulation code as well as other features of complex simulations, including quantity-of-interest functionals, solvers, boundary conditions, initial conditions, and postprocessing. GRINS is built on the FEMSystem framework that is part of the libMesh finite element library. GRINS, through FEMSystem and libMesh, supports a number of contemporary algorithms and computing features including hybrid distributed and shared memory parallelism, parallel adaptive mesh refinement on unstructured grids, goal-oriented refinement driven by dual-weighted residual error estimates in quantities of interest, and adjoint-based parameter sensitivities. GRINS is deployed both as an application binary, allowing immediate use of existing features, and as a library to which the user may build and link a stand-alone application, extending the framework to their needs. Several examples demonstrate the effectiveness of the design, including uncertainty quantification applications. Planned developments are also discussed. [ABSTRACT FROM AUTHOR]

Solutions for many problems of interest exhibit singular behaviors at domain corners or points where the boundary condition changes type. For these types of problems, direct spectral methods with the usual polynomial basis functions do not lead to a satisfactory convergence rate. We develop in this paper a Müntz-Galerkin method which is based on specially tuned Müntz polynomials to deal with the singular behaviors of the underlying problems. By exploring the relations between Jacobi polynomials and Müntz polynomials, we develop efficient implementation procedures for the Müntz-Galerkin method, and provide optimal error estimates. As an example of applications, we consider the Poisson equation with mixed Dirichlet-Neumann boundary conditions, whose solution behaves like O(r1/2) near the singular point, and demonstrate that the Müntz-Galerkin method greatly improves the rates of convergence of the usual spectral method. [ABSTRACT FROM AUTHOR]

Partial differential equations (PDEs) are commonly used to model a wide variety of physical phenomena. A PDE model of a physical problem is typically described by conservation laws, constitutive laws, material properties, boundary conditions, boundary data, and geometry. In most practical applications, however, the PDE model is only an approximation to the real physical problem due to both (i) the deliberate mathematical simplification of the model to keep it tractable and (ii) the inherent uncertainty of the physical parameters. In such cases, the PDE model may not produce a good prediction of the true state of the underlying physical problem. In this paper, we introduce a functional regression method that incorporates observations into a deterministic linear PDE model to improve its prediction of the true state. Our method is devised as follows. First, we augment the PDE model with a random Gaussian functional which serves to represent various sources of uncertainty in the model. We next derive a linear regression model for the Gaussian functional by utilizing observations and adjoint states. This allows us to determine the posterior distribution of the Gaussian functional and the posterior distribution for our estimate of the true state. Furthermore, we consider the problem of experimental design in this setting, wherein we develop an algorithm for designing experiments to efficiently reduce the variance of our state estimate. We provide several examples from the heat conduction, the convection-diffusion equation, and the reduced wave equation, all of which demonstrate the performance of the proposed methodology. [ABSTRACT FROM AUTHOR]

We present boundary integral representations of several initial boundary value problems related to the heat equation. A Galerkin discretization with piecewise constant functions in time and piecewise linear functions in space leads to optimal a priori error estimates, provided that the meshwidths in space and time satisfy ht = O(hx²). Each time step involves the solution of a linear system, whose spectral condition number is independent of the refinement under the same assumption on the mesh. We show that if the parabolic multipole method is used to apply parabolic boundary integral operators, the overall complexity of the scheme is log-linear while preserving the convergence of the Galerkin discretization method. The theoretical estimates are confirmed numerically at the end of the paper. [ABSTRACT FROM AUTHOR]

For diffusion-reaction equations employing a splitting procedure is attractive as it reduces the computational demand and facilitates a parallel implementation. Moreover, it opens up the possibility to construct second-order integrators that preserve positivity independent of the time step size used. However, for boundary conditions that are neither periodic nor of homogeneous Dirichlet type, order reduction limits its usefulness. In the situation described the Strang splitting procedure is not more accurate than Lie splitting. In this paper, we propose a modified Lie/Strang splitting procedure that, while retaining all the favorable properties of the original method, does not suffer from order reduction. That is, the modified Strang splitting is second-order accurate in time. We demonstrate our results by conducting numerical simulations in one and two space dimensions with inhomogeneous and time-dependent Dirichlet boundary conditions. In addition, a mathematically rigorous convergence analysis is conducted that confirms the results observed in the numerical simulations. [ABSTRACT FROM AUTHOR]

A type of parallel augmented subspace scheme for eigenvalue problems is proposed by using coarse space in the multigrid method. With the help of coarse space in the multigrid method, solving the eigenvalue problem in the finest space is decomposed into solving the standard linear boundary value problems and very-low-dimensional eigenvalue problems. The computational efficiency can be improved since there is no direct eigenvalue solving in the finest space and the multigrid method can act as the solver for the deduced linear boundary value problems. Furthermore, for different eigenvalues, the corresponding boundary value problem and low-dimensional eigenvalue problem can be solved in the parallel way since they are independent of each other and there exists no data exchanging. This property means that we do not need to do the orthogonalization in the highest-dimensional spaces. This is the main aim of this paper since avoiding orthogonalization can improve the scalability of the proposed numerical method. Some numerical examples are provided to validate the proposed parallel augmented subspace method. [ABSTRACT FROM AUTHOR]

In this paper we propose a local orthogonal decomposition method (LOD) for elliptic partial differential equations with inhomogeneous Dirichlet and Neumann boundary conditions. For this purpose, we present new boundary correctors which preserve the common convergence rates of the LOD, even if the boundary condition has a rapidly oscillating fine scale structure. We prove a corresponding a priori error estimate and present numerical experiments. We also demonstrate numerically that the method is reliable with respect to thin conductivity channels in the diffusion matrix. Accurate results are obtained without resolving these channels by the coarse grid and without using patches that contain the channels. [ABSTRACT FROM AUTHOR]

Large-scale networks arise in many applications. It is often of interest to be able to identify the most important nodes of a network or to ascertain the ease of traveling between nodes. These and related quantities can be determined by evaluating expressions of the form uT f(A)w, where A is the adjacency matrix that represents the graph of the network, f is a nonlinear function, such as the exponential function, and u and w are vectors, for instance, axis vectors. This paper describes a novel technique for determining upper and lower bounds for expressions uT f(A) w when A is symmetric and bounds for many vectors u and w are desired. The bounds are computed by first evaluating a low-rank approximation of A, which is used to determine rough bounds for the desired quantities for all nodes. These rough bounds indicate for which vectors and more accurate bounds should be computed with the aid of Gauss-type quadrature rules. This hybrid approach is cheaper than only using Gauss-type rules to determine accurate upper and lower bounds in the common situation when it is not known a priori for which vectors u and w accurate bounds for uT f(A)w should be computed. Several computed examples, including an application to software engineering, illustrate the performance of the hybrid method. [ABSTRACT FROM AUTHOR]

We describe a high-order accurate methodology for the numerical simulation of time-harmonic waves governed by the Helmholtz equation. Our approach combines compact finite difference schemes that provide an inexpensive venue toward high-order accuracy with the method of difference potentials developed by Ryaben'kii. The latter can be interpreted as a generalized discrete version of the method of Calderon's operators in the theory of partial differential equations. The method of difference potentials can accommodate nonconforming boundaries on regular structured grids with no loss of accuracy due to staircasing. It introduces a universal framework for treating boundary conditions of any type. A significant advantage of this method is that changing the boundary condition within a fairly broad variety does not require any major changes to the algorithm and is computationally inexpensive. In this paper, we address various types of boundary conditions using the method of difference potentials. We demonstrate the resulting numerical capabilities by solving a range of nonstandard boundary value problems for the Helmholtz equation. These include problems with variable coefficient Robin boundary conditions (including discontinuous coefficients) and problems with mixed (Dirichlet/Neumann) boundary conditions. In all our simulations, we use a Cartesian grid and a circular boundary curve. For those test cases where the overall solution was smooth, our methodology has consistently demonstrated the design fourth-order rate of grid convergence, whereas when the regularity of the solution was not sufficient, the convergence slowed down, as expected. We also show that every additional boundary condition requires only an incremental additional expense. [ABSTRACT FROM AUTHOR]

In this paper we develop a new limiter for linear reconstruction on non-coordinatealigned meshes in two space dimensions, with focus on Cartesian embedded boundary grids. Our limiter is inherently two dimensional and linearity preserving. It separately limits the x and y components of the gradient, as opposed to a scalar limiter which limits all components simultaneously with one scalar. The limiter is based on solving a tiny linear program (LP) on each cell, using a very efficient version of the simplex method. A variety of computational results on triangular and embedded boundary meshes are presented. They demonstrate that the LP limiter successfully removes oscillations and significantly increases solution accuracy compared to a scalar limiter. [ABSTRACT FROM AUTHOR]

In this paper, a hybrid approach for solving the Laplace equation in general three-dimensional (3-D) domains is presented. The approach is based on a local method for the Dirichletto-Neumann (DtN) mapping of a Laplace equation by combining a deterministic (local) boundary integral equation (BIE) method and the probabilistic Feynman-Kac formula for solutions of elliptic partial differential equations. This hybridization produces a parallel algorithm where the bulk of the computation has no need for data communication between processors. Given Dirichlet data of the solution on a domain boundary, a local BIE is established over the boundary of a local region formed by a hemisphere superimposed on the domain boundary. By using a homogeneous Dirichlet Green's function for the whole sphere, the resulting BIE involves only the Dirichlet data (the solution value) over the surface of the hemisphere, while over the patch of the domain boundary intersected by the hemisphere, both Dirichlet and Neumann data are used. Then, first, the solution value on the surface of the hemisphere is computed by the Feynman-Kac formula, which is implemented by a Monte Carlo walk-on-spheres algorithm. Second, a boundary collocation method is employed to solve the integral equation on the aforementioned local patch of the domain boundary to yield the required Neumann data there. As a result, a local method of finding the DtN mapping is obtained, which can be used to find all Neumann data on the whole domain boundary in a parallel manner. Finally, the potential solution in the whole space can be computed by an integral representation using both the Dirichlet and Neumann data over the domain boundary. [ABSTRACT FROM AUTHOR]

Evolution equations of the joint response-excitation probability density function (REPDF) generalize the existing PDF evolution equations and enable us to compute the PDF of the solution of stochastic systems driven by colored random noise. This paper aims at developing an efficient numerical method for this evolution equation of REPDF by considering the response and excitation spaces separately. For the response space, a nonconforming adaptive discontinuous Galerkin method is used to resolve both local and discontinuous dynamics while a probabilistic collocation method is used for the excitation space. We propose two fundamentally different adaptive schemes for the response space using either the local variance combined with the boundary flux difference or using particle trajectories. The effectiveness of the proposed new algorithm is demonstrated in two prototype applications dealing with randomly forced nonlinear oscillators. We first study the stochastic pendulum problem and compare the resulting PDF against the one obtained from Monte Carlo simulation. Subsequently, we study the Duffing oscillator for two different types of stochastic forcing and random initial conditions. We observe both oscillatory and chaotic dynamics and compare the results against the solution of the effective Fokker--Planck equation. The framework we develop here is general and can be readily extended to stochastic PDEs subject to random boundary conditions, random initial conditions, or random forcing terms. [ABSTRACT FROM AUTHOR]

We present a model and the related discretization for phase change problems. In particular, we are focused on the evaporation of water. The governing equations inside the do-mains and conditions on the moving interface are derived. Afterward the numerical methods for the discretization are presented. We use a level set method to capture the interface motion and use composite finite elements (CFEs) to solve the required equations in the whole domain. CFEs are a special kind of finite elements, allowing a fast calculation, because they use structured grids and respect the geometry by adapting the basis functions in the neighborhood of an interface or at the domain boundary. For the special construction of CFEs used in this paper, we present a method to take into account Dirichlet boundary conditions on the complicated domain boundary. Also, the method used for the calculation of the interface conditions within the CFE-grid is presented. We tested the Dirichlet boundary condition method for the CFEs by solving an elliptic test problem with a precomputed right-hand side. The phase change discretization is tested by verifying the d2 law, a widely used test in the phase change and moving boundary context. With the CFE discretization, we are able to solve the full three-dimensional (3D) problem. Afterward we give a short introduction to radio-frequency ablation and present a possibility for integrating the presented phase change model into the simulation of the ablation. [ABSTRACT FROM AUTHOR]

The problem of optimal mass transport arises in numerous applications, including image registration, mesh generation, reflector design, and astrophysics. One approach to solving this problem is via the Monge--Ampère equation. While recent years have seen much work in the development of numerical methods for solving this equation, very little has been done on the implementation of the transport boundary condition. In this paper, we propose a method for solving the transport problem by iteratively solving a Monge--Ampère equation with Neumann boundary conditions. To enable mappings between variable densities, we extend an earlier discretization of the equation to allow for right-hand sides that depend on gradients of the solution [B. D. Froese and A. M. Oberman, SIAM J. Numer. Anal., 49 (2011), pp. 1692-1714]. This discretization provably converges to the viscosity solution. The resulting system is solved efficiently with Newton's method. We provide several challenging computational examples that demonstrate the effectiveness 1.3 and effciency (O (M)--O (M1.3) time) of the proposed method. [ABSTRACT FROM AUTHOR]

In this paper, in order to numerically solve for multiple positive solutions to a singularly perturbed Neumann boundary value problem in mathematical biology and other applications, a local minimax method is modified with new local mesh refinement and other strategies. Algorithm convergence and other related properties are verified. Motivated by the numerical algorithm and convinced by the numerical results, a Morse index approach is used to identify the Morse index of the root solution u¹ε = 1 at any perturbation value, its bifurcation points, and then the critical perturbation value. Many interesting numerical solutions are computed for the first time and displayed with their contours and mesh profiles to illustrate the theory and method. [ABSTRACT FROM AUTHOR]

In this paper, a time-space fractional diffusion equation in two dimensions (TSFDE-2D) with homogeneous Dirichlet boundary conditions is considered. The TSFDE-2D is obtained from the standard diffusion equation by replacing the first-order time derivative with the Caputo fractional derivative Due to image rights restrictions, multiple line equation(s) cannot be graphically displayed., γ ∈ (0,1), and the second-order space derivatives with the fractional Laplacian -(-Δ)α/2, α ∈ (1,2]. Using the matrix transfer technique proposed by Ilić et al. [M. Ilić, F. Liu, I. Turner, and V. Anh, Fract. Calc. Appl. Anal, 9 (2006), pp. 333-349], the TSFDE-2D is transformed into a time fractional differential system as Due to image rights restrictions, multiple line equation(s) cannot be graphically displayed., where A is the approximate matrix representation of (-Δ). Traditional approximation of Aα/2 requires diagonalization of A, which is very time-consuming for large sparse matrices. The novelty of our proposed numerical schemes is that, using either the finite difference method or the Laplace transform to handle the Caputo time fractional derivative, the solution of the TSFDE-2D is written in terms of a matrix function vector product f(A)b at each time step, where b is a suitably defined vector. Depending on the method used to generate the matrix A, the product f(A)b can be approximated using either the preconditioned Lanczos method when A is symmetric or the M-Lanzcos method when A is nonsymmetric, which are powerful techniques for solving large linear systems. We give error bounds for the new methods and illustrate their roles in solving the TSFDE-2D. We also derive the analytical solution of the TSFDE-2D in terms of the Mittag-Leffler function. Finally, numerical results are presented to verify the proposed numerical solution strategies. [ABSTRACT FROM AUTHOR]