Showing papers in "SIAM Journal on Scientific Computing in 2007"
TL;DR: This paper considers how specially structured tensors allow for efficient storage and computation, and proposes storing sparse tensors using coordinate format and describes the computational efficiency of this scheme for various mathematical operations, including those typical to tensor decomposition algorithms.
Abstract: In this paper, the term tensor refers simply to a multidimensional or $N$-way array, and we consider how specially structured tensors allow for efficient storage and computation. First, we study sparse tensors, which have the property that the vast majority of the elements are zero. We propose storing sparse tensors using coordinate format and describe the computational efficiency of this scheme for various mathematical operations, including those typical to tensor decomposition algorithms. Second, we study factored tensors, which have the property that they can be assembled from more basic components. We consider two specific types: A Tucker tensor can be expressed as the product of a core tensor (which itself may be dense, sparse, or factored) and a matrix along each mode, and a Kruskal tensor can be expressed as the sum of rank-1 tensors. We are interested in the case where the storage of the components is less than the storage of the full tensor, and we demonstrate that many elementary operations can be computed using only the components. All of the efficiencies described in this paper are implemented in the Tensor Toolbox for MATLAB.
505 citations
TL;DR: New convergence results that show superlinear convergence of the parareal algorithm when used on bounded time intervals, and linear convergence for unbounded intervals are shown.
Abstract: The parareal algorithm is a method to solve time-dependent problems parallel in time: it approximates parts of the solution later in time simultaneously to parts of the solution earlier in time. In this paper the relation of the parareal algorithm to space-time multigrid and multiple shooting methods is first briefly discussed. The focus of the paper is on new convergence results that show superlinear convergence of the algorithm when used on bounded time intervals, and linear convergence for unbounded intervals.
413 citations
TL;DR: A new Lagrangian cell-centered scheme for two-dimensional compressible flows with main new feature of the introduction of four pressures on each edge, two for each node on each side of the edge, and a semidiscrete entropy inequality is provided.
Abstract: We present a new Lagrangian cell-centered scheme for two-dimensional compressible flows. The primary variables in this new scheme are cell-centered, i.e., density, momentum, and total energy are defined by their mean values in the cells. The vertex velocities and the numerical fluxes through the cell interfaces are not computed independently, contrary to standard approaches, but are evaluated in a consistent manner due to an original solver located at the nodes. The main new feature of the algorithm is the introduction of four pressures on each edge, two for each node on each side of the edge. This extra degree of freedom allows us to construct a nodal solver which fulfills two properties. First, the conservation of momentum and total energy is ensured. Second, a semidiscrete entropy inequality is provided. In the case of a one-dimensional flow, the solver reduces to the classical Godunov acoustic solver: it can be considered as its two-dimensional generalization. Many numerical tests are presented. They are representative test cases for compressible flows and demonstrate the robustness and the accuracy of this new solver.
362 citations
TL;DR: A new projection method to solve large-scale continuous-time Lyapunov matrix equations based on matrix factorizations, generated as a combination of Krylov subspaces in A and A^{-1}$.
Abstract: In this paper we propose a new projection method to solve large-scale continuous-time Lyapunov matrix equations. The new approach projects the problem onto a much smaller approximation space, generated as a combination of Krylov subspaces in $A$ and $A^{-1}$. The reduced problem is then solved by means of a direct Lyapunov scheme based on matrix factorizations. The reported numerical results show the competitiveness of the new method, compared to a state-of-the-art approach based on the factorized alternating direction implicit iteration.
349 citations
TL;DR: A new method is introduced, RBF-QR, which entirely eliminates such ill-conditioning in the special case when the data points are distributed over the surface of a sphere, and it allows the RBF shape parameter to be optimized without the limitations imposed by stability concerns.
Abstract: When radial basis functions (RBFs) are made increasingly flat, the interpolation error typically decreases steadily until some point when Runge-type oscillations either halt or reverse this trend. Because the most obvious method to calculate an RBF interpolant becomes a numerically unstable algorithm for a stable problem in the case of near-flat basis functions, there will typically also be a separate point at which disastrous ill-conditioning enters. We introduce here a new method, RBF-QR, which entirely eliminates such ill-conditioning, and we apply it in the special case when the data points are distributed over the surface of a sphere. This algorithm works even for thousands of node points, and it allows the RBF shape parameter to be optimized without the limitations imposed by stability concerns. Since interpolation in the flat RBF limit on a sphere is found to coincide with spherical harmonics interpolation, new insights are gained as to why the RBF approach (with nonflat basis functions) often is the more accurate of the two methods.
245 citations
TL;DR: This paper considers the localization of a collection of small, three-dimensional bounded homogeneous inclusions via time-harmonic electromagnetic means, typically using arrays of electric or magnetic dipole transmitters and receivers with given polarization(s) at some distance from the collection, possibly also lying in the far field.
Abstract: In this paper we consider the localization of a collection of small, three-dimensional bounded homogeneous inclusions via time-harmonic electromagnetic means, typically using arrays of electric or magnetic dipole transmitters and receivers with given polarization(s) at some distance from the collection, possibly also lying in the far field. The inclusions, somewhat apart or closely spaced, are buried within a homogeneous medium, and are of arbitrary contrast of permittivity, conductivity, and permeability vis-ag-vis this embedding medium. The problem is formulated as an inverse scattering problem for the full Maxwell equations and it involves a robust asymptotic modeling of the multistatic response matrix. No specific application is studied at this stage, but characterization of obstacles in subsoils, nondestructive evaluation of man-made structures, and medical imaging are primary fields of application envisaged. The proposed approach uses a MUSIC (multiple signal classification)-type algorithm, and it yields fast numbering, accurate localization, and estimates of the electromagnetic and geometric parameters (polarization tensors) of the inclusions. The mathematical machinery is detailed first, some specific attention being given to triaxial ellipsoidal inclusions and degenerate spherical shapes (for the latter known results are retrieved). Then, the viability of this algorithm—which would be easily extended to planarly layered environments by introduction of their Green’s functions—is documented by a variety of numerical results from synthetic, noiseless, and severely noisy field data.
219 citations
TL;DR: A semi-Lagrangian method for solving the HJB equation for a typical gas storage valuation problem is presented and is able to handle a wide class of spot price models that exhibit mean-reverting seasonality dynamics and price jumps.
Abstract: The valuation of a gas storage facility is characterized as a stochastic control problem, resulting in a Hamilton-Jacobi-Bellman (HJB) equation. In this paper, we present a semi-Lagrangian method for solving the HJB equation for a typical gas storage valuation problem. The method is able to handle a wide class of spot price models that exhibit mean-reverting seasonality dynamics and price jumps. We develop fully implicit and Crank-Nicolson timestepping schemes based on a semi-Lagrangian approach and prove the convergence of fully implicit timestepping to the viscosity solution of a modified HJB equation posed on a bounded domain, provided that a strong comparison result holds. The semi-Lagrangian approach avoids Policy-type iterations required by an implicit finite difference method without requiring additional cost. Numerical experiments are presented for several variants of the basic scheme.
161 citations
TL;DR: An error estimator and an adaptive algorithm for efficient solution of parabolic partial differential equations are developed and the space and time discretization errors are equilibrated, leading to an efficient method.
Abstract: In this paper, we develop an error estimator and an adaptive algorithm for efficient solution of parabolic partial differential equations. The error estimator assesses the discretization error with respect to a given quantity of physical interest and separates the influence of the time and space discretizations. This allows us to set up an efficient adaptive strategy producing economical (locally) refined meshes for each time step and an adapted time discretization. The space and time discretization errors are equilibrated, leading to an efficient method.
152 citations
TL;DR: A block-iterative projection method for solving linear equations and/or inequalities that allows diagonal componentwise relaxation in conjunction with orthogonal projections onto the individual hyperplanes of the system, and is thus called diagonally relaxed orthOGonal projections (DROP).
Abstract: We propose and study a block-iterative projection method for solving linear equations and/or inequalities. The method allows diagonal componentwise relaxation in conjunction with orthogonal projections onto the individual hyperplanes of the system, and is thus called diagonally relaxed orthogonal projections (DROP). Diagonal relaxation has proven useful in accelerating the initial convergence of simultaneous and block-iterative projection algorithms, but until now it was available only in conjunction with generalized oblique projections in which there is a special relation between the weighting and the oblique projections. DROP has been used by practitioners, and in this paper a contribution to its convergence theory is provided. The mathematical analysis is complemented by some experiments in image reconstruction from projections which illustrate the performance of DROP.
149 citations
TL;DR: A new directional multilevel algorithm for solving N-body or N-point problems with highly oscillatory kernels which is proved to have $O(N\log N)$ computational complexity for any given accuracy when the points are sampled from a two dimensional surface.
Abstract: This paper introduces a new directional multilevel algorithm for solving $N$-body or $N$-point problems with highly oscillatory kernels. These systems often result from the boundary integral formulations of scattering problems and are difficult due to the oscillatory nature of the kernel and the non-uniformity of the particle distribution. We address the problem by first proving that the interaction between a ball of radius $r$ and a well-separated region has an approximate low rank representation, as long as the well-separated region belongs to a cone with a spanning angle of $O(1/r)$ and is at a distance which is at least $O(r^2)$ away from from the ball. We then propose an efficient and accurate procedure which utilizes random sampling to generate such a separated, low rank representation. Based on the resulting representations, our new algorithm organizes the high frequency far field computation by a multidirectional and multiscale strategy to achieve maximum efficiency. The algorithm performs well on a large group of highly oscillatory kernels. Our algorithm is proved to have $O(N\log N)$ computational complexity for any given accuracy when the points are sampled from a two dimensional surface. We also provide numerical results to demonstrate these properties.
144 citations
TL;DR: Numerical experiments, comparing the implicit midpoint rule with Heun and leapfrog methods on nonlinear equations with additive or multiplicative noise, produce behavior similar to the linear case.
Abstract: We seek numerical methods for second-order stochastic differential equations that reproduce the stationary density accurately for all values of damping. A complete analysis is possible for scalar linear second-order equations (damped harmonic oscillators with additive noise), where the statistics are Gaussian and can be calculated exactly in the continuous-time and discrete-time cases. A matrix equation is given for the stationary variances and correlation for methods using one Gaussian random variable per timestep. The only Runge-Kutta method with a nonsingular tableau matrix that gives the exact steady state density for all values of damping is the implicit midpoint rule. Numerical experiments, comparing the implicit midpoint rule with Heun and leapfrog methods on nonlinear equations with additive or multiplicative noise, produce behavior similar to the linear case.
TL;DR: In this paper, a comprehensive spectral analysis of the Helmholtz operator preconditioned with a shifted Laplacian is presented, and the results of the spectral analysis are combined with an upper bound on the GMRES-residual norm.
Abstract: Shifted Laplace preconditioners have attracted considerable attention as a technique to speed up convergence of iterative solution methods for the Helmholtz equation. In this paper we present a comprehensive spectral analysis of the Helmholtz operator preconditioned with a shifted Laplacian. Our analysis is valid under general conditions. The propagating medium can be heterogeneous, and the analysis also holds for different types of damping, including a radiation condition for the boundary of the computational domain. By combining the results of the spectral analysis of the preconditioned Helmholtz operator with an upper bound on the GMRES-residual norm, we are able to provide an optimal value for the shift and to explain the mesh-dependency of the convergence of GMRES preconditioned with a shifted Laplacian. We illustrate our results with a seismic test problem.
TL;DR: A boundary element method with basis functions that incorporate the asymptotic behavior of the solution at high frequencies and combines this hybrid method with very effective quadrature rules for oscillatory integrals to obtain a sparse discretization matrix for the oscillatory problem.
Abstract: We consider two-dimensional scattering problems, formulated as an integral equation defined on the boundary of the scattering obstacle. The oscillatory nature of high-frequency scattering problems necessitates a large number of unknowns in classical boundary element methods. In addition, the corresponding discretization matrix of the integral equation is dense. We formulate a boundary element method with basis functions that incorporate the asymptotic behavior of the solution at high frequencies. The method exhibits the effectiveness of asymptotic methods at high frequencies with only few unknowns, but retains accuracy for lower frequencies. New in our approach is that we combine this hybrid method with very effective quadrature rules for oscillatory integrals. As a result, we obtain a sparse discretization matrix for the oscillatory problem. Moreover, numerical experiments indicate that the accuracy of the solution actually increases with increasing frequency. The sparse discretization applies to problems where the phase of the solution can be predicted a priori, for example in the case of smooth and convex scatterers.
TL;DR: A simple and computationally inexpensive adaptive strategy that allows us to simultaneously capture the unique entropy solution and to achieve a high resolution of the computed solution is proposed.
Abstract: We discover that the choice of a piecewise polynomial reconstruction is crucial in computing solutions of nonconvex hyperbolic (systems of) conservation laws. Using semidiscrete central-upwind schemes, we illustrate that the obtained numerical approximations may fail to converge to the unique entropy solution or the convergence may be so slow that achieving a proper resolution would require the use of (almost) impractically fine meshes. For example, in the scalar case, all computed solutions seem to converge to solutions that are entropy solutions for some entropy pairs. However, in most applications, one is interested in capturing the unique (Kruzhkov) solution that satisfies the entropy condition for all convex entropies. We present a number of numerical examples that demonstrate the convergence of the solutions, computed with the dissipative second-order minmod reconstruction, to the unique entropy solution. At the same time, more compressive and/or higher-order reconstructions may fail to resolve composite waves, typically present in solutions of nonconvex conservation laws, and thus may fail to recover the Kruzhkov solution. In this paper, we propose a simple and computationally inexpensive adaptive strategy that allows us to simultaneously capture the unique entropy solution and to achieve a high resolution of the computed solution. We use the dissipative minmod reconstruction near the points where convexity changes and utilize a fifth-order weighted essentially nonoscillatory (WENO5) reconstruction in the rest of the computational domain. Our numerical examples (for one- and two-dimensional scalar and systems of conservation laws) demonstrate the robustness, reliability, and nonoscillatory nature of the proposed adaptive method.
TL;DR: It is discussed here how state-of-the-art computational geometry methods make it tractable to solve the problem of finding the extreme rays of a cone with a degenerate vertex at the origin, a difficult problem.
Abstract: We discuss an implementation of a derivative-free generating set search method for linearly constrained minimization with no assumption of nondegeneracy placed on the constraints. The convergence guarantees for generating set search methods require that the set of search directions possesses certain geometrical properties that allow it to approximate the feasible region near the current iterate. In the hard case, the calculation of the search directions corresponds to finding the extreme rays of a cone with a degenerate vertex at the origin, a difficult problem. We discuss here how state-of-the-art computational geometry methods make it tractable to solve this problem in connection with generating set search. We also discuss a number of other practical issues of implementation, such as the careful treatment of equality constraints and the desirability of augmenting the set of search directions beyond the theoretically minimal set. We illustrate the behavior of the implementation on several problems from the CUTEr test suite. We have found it to be successful on problems with several hundred variables and linear constraints.
TL;DR: A projection-based iterative algorithm based on a joint bidiagonalization algorithm appropriate for large-scale problems when it is computationally infeasible to transform the regularized problem to standard form and how estimates of the corresponding optimal regularization parameter can be efficiently obtained is presented.
Abstract: We present a projection-based iterative algorithm for computing general-form Tikhonov regularized solutions to the problem $\min_x\{\| Ax-b \|_2^2+\lambda^2\| Lx \|_2^2\}$, where the regularization matrix $L$ is not the identity. Our algorithm is designed for the common case where $\lambda$ is not known a priori. It is based on a joint bidiagonalization algorithm and is appropriate for large-scale problems when it is computationally infeasible to transform the regularized problem to standard form. By considering the projected problem, we show how estimates of the corresponding optimal regularization parameter can be efficiently obtained. Numerical results illustrate the promise of our projection-based approach.
TL;DR: The objective of this article is to show how an efficient discretization can be achieved by hierarchical approximation as well as asymptotic expansions of the underlying continuous problem.
Abstract: A major challenge in computational finance is the pricing of options that depend on a large number of risk factors. Prominent examples are basket or index options where dozens or even hundreds of stocks constitute the underlying asset and determine the dimensionality of the corresponding degenerate parabolic equation. The objective of this article is to show how an efficient discretization can be achieved by hierarchical approximation as well as asymptotic expansions of the underlying continuous problem. The relation to a number of state-of-the-art methods is highlighted.
TL;DR: A proof for existence and uniqueness of the weak solution in the function space of zero-mean potential functions, using a subtraction approach, and state convergence properties of the finite element (FE) method for the numerical solution to the correction potential.
Abstract: In electroencephalography (EEG) source analysis, a dipole is widely used as the model of the current source. The dipole introduces a singularity on the right-hand side of the governing Poisson-type differential equation that has to be treated specifically when solving the equation toward the electric potential. In this paper, we give a proof for existence and uniqueness of the weak solution in the function space of zero-mean potential functions, using a subtraction approach. The method divides the total potential into a singularity and a correction potential. The singularity potential is due to a dipole in an infinite region of homogeneous conductivity. We then state convergence properties of the finite element (FE) method for the numerical solution to the correction potential. We validate our approach using tetrahedra and regular and geometry-conforming node-shifted hexahedra elements in an isotropic three-layer sphere model and a model with anisotropic middle compartment. Validation is carried out using sophisticated visualization techniques, correlation coefficient (CC), and magnification factor (MAG) for a comparison of the numerical results with analytical series expansion formulas at the surface and within the volume conductor. For the subtraction approach, with regard to the accuracy in the anisotropic three-layer sphere model (CC of 0.998 or better and MAG of 4.3% or better over the whole range of realistic eccentricities) and to the computational complexity, 2mm node-shifted hexahedra achieve the best results. A relative FE solver accuracy of $10^{-4}$ is sufficient for the used algebraic multigrid preconditioned conjugate gradient approach. Finally, we visualize the computed potentials of the subtraction method in realistically shaped FE head volume conductor models with anisotropic skull compartments.
TL;DR: This paper considers thermoacoustic tomography as the inverse problem of determining from lateral Cauchy data the unknown initial conditions in a wave equation with spatially varying coefficients and derives a numerical scheme for the solution of the quasi-reversibility problem by a $B$-spline Galerkin method.
Abstract: In this paper we consider thermoacoustic tomography as the inverse problem of determining from lateral Cauchy data the unknown initial conditions in a wave equation with spatially varying coefficients. This problem also occurs in several applications in the area of medical imaging and nondestructive testing. Using the method of quasi-reversibility, the original ill-posed problem is replaced with a boundary value problem for a fourth order partial differential equation. We find a weak $H^2$ solution of this problem and show that it is a well-posed elliptic problem. Error estimates and convergence of the approximation follow from observability estimates for the wave equation, which are proved using a Carleman estimate. We derive a numerical scheme for the solution of the quasi-reversibility problem by a $B$-spline Galerkin method, for which we give error estimates. Finally, we present numerical results supporting the robustness of this method for the reconstruction of initial conditions from full and limited boundary data.
TL;DR: Two three-level BDDC methods are introduced for solving the coarse problem approximately for problems in three dimensions, an extension of previous work for the two-dimensional case.
Abstract: Balancing domain decomposition by constraints (BDDC) methods are nonoverlapping iterative substructuring domain decomposition methods for the solution of large sparse linear algebraic systems arising from the discretization of elliptic boundary value problems. Their coarse problems are given in terms of a small number of continuity constraints for each subdomain, which are enforced across the interface. The coarse problem matrix is generated and factored by a direct solver at the beginning of the computation and it can ultimately become a bottleneck if the number of subdomains is very large. In this paper, two three-level BDDC methods are introduced for solving the coarse problem approximately for problems in three dimensions. This is an extension of previous work for the two-dimensional case. Edge constraints are considered in this work since vertex constraints alone, which work well in two dimensions, result in a noncompetitive algorithm in three dimensions. Some new technical tools are then needed in the analysis and this makes the three-dimensional case more complicated. Estimates of the condition numbers are provided for two three-level BDDC methods, and numerical experiments are also discussed.
TL;DR: It is demonstrated that a small modification of the multiplicative, additive, and restricted additive Schwarz preconditioner at the algebraic level, motivated by optimized Schwarz methods defined at the continuous level, leads to a significant reduction in the iteration count of the iterative solver.
Abstract: We demonstrate that a small modification of the multiplicative, additive, and restricted additive Schwarz preconditioner at the algebraic level, motivated by optimized Schwarz methods defined at the continuous level, leads to a significant reduction in the iteration count of the iterative solver. Numerical experiments using finite difference and spectral element discretizations of the positive definite Helmholtz problem and an idealized climate simulation illustrate the effectiveness of the new approach.
TL;DR: A variational formulation of combined motion by minus the Laplacian of curvature and mean curvature flow is presented, which has very good properties with respect to the equidistribution of mesh points and, if applicable, area conservation.
Abstract: We present a variational formulation of combined motion by minus the Laplacian of curvature and mean curvature flow, as well as related flows. The proposed scheme covers both the closed curve case and the case of curves that are connected via triple or quadruple junction points or intersect the external boundary. On introducing a parametric finite element approximation, we prove stability bounds and compare our scheme with existing approaches. The presented scheme has very good properties with respect to the equidistribution of mesh points and, if applicable, area conservation.
TL;DR: In this paper, a general purpose algorithm for rapidly computing certain types of oscillatory integrals which frequently arise in problems connected to wave propagation, general hyperbolic equations, and curvilinear tomography is introduced.
Abstract: We introduce a general purpose algorithm for rapidly computing certain types of oscillatory integrals which frequently arise in problems connected to wave propagation, general hyperbolic equations, and curvilinear tomography. The problem is to numerically evaluate a so-called Fourier integral operator (FIO) of the form $\int e^{2\pi i \Phi(x,\xi)} a(x,\xi) \hat{f}(\xi) \d\xi$ at points given on a Cartesian grid. Here, $\xi$ is a frequency variable, $\hat f(\xi)$ is the Fourier transform of the input $f$, $a(x,\xi)$ is an amplitude, and $\Phi(x,\xi)$ is a phase function, which is typically as large as $|\xi|$; hence the integral is highly oscillatory. Because a FIO is a dense matrix, a naive matrix vector product with an input given on a Cartesian grid of size $N$ by $N$ would require $O(N^4)$ operations. This paper develops a new numerical algorithm which requires $O(N^{2.5} \log N)$ operations and as low as $O(\sqrt{N})$ in storage space (the constants in front of these estimates are small). It operates by localizing the integral over polar wedges with small angular aperture in the frequency plane. On each wedge, the algorithm factorizes the kernel $e^{2 \pi i \Phi(x,\xi)} a(x,\xi)$ into two components: (1) a diffeomorphism which is handled by means of a nonuniform FFT and (2) a residual factor which is handled by numerical separation of the spatial and frequency variables. The key to the complexity and accuracy estimates is the fact that the separation rank of the residual kernel is provably independent of the problem size. Several numerical examples demonstrate the numerical accuracy and low computational complexity of the proposed methodology. We also discuss the potential of our ideas for various applications such as reflection seismology.
TL;DR: A multi-resolution scheme is shown to greatly improve the robustness of the Galerkin procedure in presence of steep dependences, but this improvement comes with a higher computational cost which drastically increases with the number of uncertain reaction rates.
Abstract: This paper presents a multi-resolution approach for the propagation of parametric uncertainty in chemical systems. It is motivated by previous studies where Galerkin formulations of Wiener-Hermite expansions were found to fail in the presence of steep dependences of the species concentrations with regard to the reaction rates. The multi-resolution scheme is based on representation of the uncertain concentration in terms of compact polynomial multi-wavelets, allowing for the control of the convergence in terms of polynomial order and resolution level. The resulting representation is shown to greatly improve the robustness of the Galerkin procedure in presence of steep dependences. However, this improvement comes with a higher computational cost which drastically increases with the number of uncertain reaction rates. To overcome this drawback an adaptive strategy is proposed to control locally (in the parameter space) and in time the resolution level. The efficiency of the method is demonstrated for an uncertain chemical system having eight random parameters.
TL;DR: This study considers various ILU-class preconditioners and investigates the parameters that render them safely applicable to common surface integral formulations without increasing the complexity of MLFMA, and concludes that the no-fill ILU(0) preconditionser is an optimal choice for the combined-field integral equation (CFIE).
Abstract: Iterative solution of large-scale scattering problems in computational electromagnetics with the multilevel fast multipole algorithm (MLFMA) requires strong preconditioners, especially for the electric-field integral equation (EFIE) formulation. Incomplete LU (ILU) preconditioners are widely used and available in several solver packages. However, they lack robustness due to potential instability problems. In this study, we consider various ILU-class preconditioners and investigate the parameters that render them safely applicable to common surface integral formulations without increasing the ${\cal O}(n\log n)$ complexity of MLFMA. We conclude that the no-fill ILU(0) preconditioner is an optimal choice for the combined-field integral equation (CFIE). For EFIE, we establish the need to resort to methods depending on drop tolerance and apply pivoting for problems with high condition estimate. We propose a strategy for the selection of the parameters so that the preconditioner can be used as a black-box method. Robustness and efficiency of the employed preconditioners are demonstrated over several test problems.
TL;DR: In the present paper a simple Newton-type procedure for certain piecewise linear systems is derived and shown to have a finite termination property, i.e., it converges to the exact solution in a finite number of steps.
Abstract: The correct formulation of numerical models for free-surface hydrodynamics often requires the solution of special linear systems whose coefficient matrix is a piecewise constant function of the solution itself. In so doing one may prevent the development of unrealistic negative water depths. The resulting piecewise linear systems are equivalent to particular linear complementarity problems whose solutions could be obtained by using, for example, interior point methods. These methods may have a favorable convergence property, but they are purely iterative and convergence to the exact solution is proven only in the limit of an infinite number of iterations. In the present paper a simple Newton-type procedure for certain piecewise linear systems is derived and discussed. This procedure is shown to have a finite termination property, i.e., it converges to the exact solution in a finite number of steps, and, actually, it converges very quickly, as confirmed by a few numerical tests.
TL;DR: Pointwise weighted essentially nonoscillatory (PWENO) interpolation is applied to obtain semi-Lagrangian and flux balance methods that together with splitting techniques form the building blocks of the numerical approach.
Abstract: We demonstrate the ability of nonoscillatory interpolation strategies for solving efficiently the transport phase in kinetic systems with applications to charged particle transport in plasmas and semiconductors. Pointwise weighted essentially nonoscillatory (PWENO) interpolation is applied to obtain semi-Lagrangian and flux balance methods that together with splitting techniques form the building blocks of our numerical approach. These methods do not present the restrictive CFL condition typical of finite-difference methods with explicit time-solvers, and, moreover, they provide reliable results controlling parasite oscillations from classical polynomial interpolation while giving highly accurate approximations of smooth parts of the solutions. We perform and compare these methods in different benchmark problems for Vlasov or collisional models for charged particle transport.
TL;DR: The main result is that under mild assumptions the total complexity for solving the interpolation problem at $M$ arbitrary nodes is of order ${\cal O}(M\log M)$.
Abstract: A fast and reliable algorithm for the optimal interpolation of scattered data on the torus $\mathbb{T}^d$ by multivariate trigonometric polynomials is presented. The algorithm is based on a variant of the conjugate gradient method in combination with the fast Fourier transforms for nonequispaced nodes. The main result is that under mild assumptions the total complexity for solving the interpolation problem at $M$ arbitrary nodes is of order ${\cal O}(M\log M)$. This result is obtained by the use of localized trigonometric kernels where the localization is chosen in accordance with the spatial dimension $d$. Numerical examples show the efficiency of the new algorithm.
TL;DR: This research approaches the eigenproblem from the nonlinear perspective, which helps to develop two nearly optimal methods, one of which extends the recent Jacobi-Davidson conjugate gradient method to JDQMR, improving robustness and efficiency.
Abstract: Large, sparse, Hermitian eigenvalue problems are still some of the most computationally challenging tasks. Despite the need for a robust, nearly optimal preconditioned iterative method that can operate under severe memory limitations, no such method has surfaced as a clear winner. In this research we approach the eigenproblem from the nonlinear perspective, which helps us develop two nearly optimal methods. The first extends the recent Jacobi-Davidson conjugate gradient (JDCG) method to JDQMR, improving robustness and efficiency. The second method, generalized-Davidson+1 (GD+1), utilizes the locally optimal conjugate gradient recurrence as a restarting technique to achieve almost optimal convergence. We describe both methods within a unifying framework and provide theoretical justification for their near optimality. A choice between the most efficient of the two can be made at runtime. Our extensive experiments confirm the robustness, the near optimality, and the efficiency of our multimethod over other state-of-the-art methods.
TL;DR: A numerical method for simulations of nonlinear surface water waves over variable bathymetry based on the reduction of the problem to a lower-dimensional Hamiltonian system involving boundary quantities alone, using the Dirichlet-Neumann operator.
Abstract: We present a numerical method for simulations of nonlinear surface water waves over variable bathymetry. It is applicable to either two- or three-dimensional flows, as well as to either static or moving bottom topography. The method is based on the reduction of the problem to a lower-dimensional Hamiltonian system involving boundary quantities alone. A key component of this formulation is the Dirichlet-Neumann operator which, in light of its joint analyticity properties with respect to surface and bottom deformations, is computed using its Taylor series representation. We present new, stabilized forms for the Taylor terms, each of which is efficiently computed by a pseudospectral method using the fast Fourier transform. Physically relevant applications are displayed to illustrate the performance of the method; comparisons with analytical solutions and laboratory experiments are provided.