Search

In this paper, we develop, analyze and numerically test an invariant-preserving three-level linearized implicit difference scheme for a rotation-two-component Camassa--Holm system [L. Fan, H. Gao, and Y. Liu, Adv. Math. 291 (2016), pp. 59-89], which contains strongly nonlinear terms and high-order derivative terms. We prove that the numerical scheme is uniquely solvable and second-order convergent for both the spatial and temporal discretizations. Optimal error estimates for the velocity in the L∞-norm and for the surface elevation in the L2-norm are obtained under a suitable step-size ratio restriction based on energy analysis. In particular, the analysis for the strongly nonlinear term is carried out by splitting it into several recursive expressions. Moreover, the numerical scheme preserves at least two conservation invariants: mass and energy. Extensive numerical experiments including long time simulations for zero/nonzero rotation parameters demonstrate solution behavior and verify the convergence results as well as mass/energy conservation. [ABSTRACT FROM AUTHOR]

Image registration is one of the most challenging problems in image processing, where ill-posedness arises due to noisy data as well as nonuniqueness, and hence the choice of regularization is crucial. This paper presents hyperelasticity as a regularizer and introduces a new and stable numerical implementation. On one hand, hyperelastic registration is an appropriate model for large and highly nonlinear deformations, for which a linear elastic model needs to fail. On the other hand, the hyperelastic regularizer yields very regular and diffeomorphic transformations. While hyperelasticity might be considered as just an additional outstanding regularization option for some applications, it becomes inevitable for applications involving higher order distance measures like mass-preserving registration. The paper gives a short introduction to image registration and hyperelasticity. The hyperelastic image registration problem is phrased in a variational setting, and an existence proof is provided. The focus of the paper, however, is on a robust numerical scheme. A key challenge is an unbiased discretization of hyperelasticity, which enables the numerical monitoring of variations of length, surface, and volume of infinitesimal reference elements. We resolve this issue by using a nodal-based discretization with a special tetrahedral partitioning. The potential of the hyperelastic registration is demonstrated in a direct comparison with a linear elastic registration in an academic example. The paper also presents a real life application from three-dimensional positron emission tomography of the human heart which requires mass preservation, and thus hyperelastic registration is the only option. [ABSTRACT FROM AUTHOR]

In this paper we propose a novel computational technique to solve the Eikonal equation efficiently on parallel architectures. The proposed method manages the list of active nodes and iteratively updates the solutions on those nodes until they converge. Nodes are added to or removed from the list based on a convergence measure, but the management of this list does not entail an extra burden of expensive ordered data structures or special updating sequences. The proposed method has suboptimal worst-case performance but, in practice, on real and synthetic datasets, runs faster than guaranteed-optimal alternatives. Furthermore, the proposed method uses only local, synchronous updates and therefore has better cache coherency, is simple to implement, and scales efficiently on parallel architectures. This paper describes the method, proves its consistency, gives a performance analysis that compares the proposed method against the state-of-the-art Eikonal solvers, and describes the implementation on a single instruction multiple datastream (SIMD) parallel architecture. [ABSTRACT FROM AUTHOR]

We study high order numerical quadratures to one dimensional delta function integrals in this paper. This is motivated by the fact that traditional numerical quadratures give only first order accuracy in general. We provide criteria on discrete delta functions and support size formulas which ensure any desired accuracy of the numerical quadratures. Discrete delta functions and support size formulas satisfying these criteria are designed. Numerical examples are presented to verify the performed analysis and the high order accuracy of the proposed numerical quadratures. [ABSTRACT FROM AUTHOR]

This paper analyzes the process of merging groups of many Gaussian basis functions into a single basis function in vortex simulations. Analysis of the equations governing this process yields fundamental parameters and uniform estimates of the merger-induced error. In a merging event, the uniform error bound depends only on relationships between nearby computational elements permitting fast and effective merging schemes. In this paper, one such algorithm is proposed and demonstrated. [ABSTRACT FROM AUTHOR]

This paper presents some guidance for choosing multigrid preconditioners for systems of partial differential equations (PDEs). Rather than developing a black-box algorithm for arbitrary systems of PDEs, the strategy in this paper is to classify systems with multigrid preconditioners. From this classification, a new preconditioner is developed, and it is shown to be suitable for a particular class of strongly cross-coupled systems, solving the whole cross-coupled system almost as effectively as solving just the diagonal system. Also, this paper gives an explanation of why block diagonal preconditioners are somewhat effective for certain strongly cross-coupled systems that at first glance appear to be intractable for multigrid preconditioning. [ABSTRACT FROM AUTHOR]

In a recent paper, the authors introduced the recovery method (local energy matching principle) for solving large systems of lattice equations. The idea is to construct a partial differential equation along with a finite element discretization such that the arising system of linear equations has equivalent energy as the original system of lattice equations. Since a vast variety of efficient solvers is available for solving large systems of finite element discretizations of elliptic PDEs, these solvers may serve as preconditioners for the system of lattice equations. In this paper, we will focus on both the theoretical and the numerical dependence of the method on various mesh-dependent parameters, which can be easily computed and monitored during the solution process. Systematic parameter tests have been performed which underline (a) the robustness and the efficiency of the recovery method and (b) the reliability of the control parameters, which are computed in a preprocessing step to predict the performance of the preconditioner based on the recovery method. [ABSTRACT FROM AUTHOR]

In this paper, we propose iterative algorithms for solving image restoration problems. The iterative algorithms are based on decoupling of deblurring and denoising steps in the restoration process. In the deblurring step, an efficient deblurring method using fast transforms can be employed. In the denoising step, effective methods such as the wavelet shrinkage denoising method or the total variation denoising method can be used. The main advantage of this proposal is that the resulting algorithms can be very efficient and can produce better restored images in visual quality and signalto-noise ratio than those by the restoration methods using the combination of a data-fitting term and a regularization term. The convergence of the proposed algorithms is shown in the paper. Numerical examples are also given to demonstrate the effectiveness of these algorithms. [ABSTRACT FROM AUTHOR]

In this paper, we present a family of domain decomposition based on an Aitken-like acceleration of the Schwarz method seen as an iterative procedure with a linear rate of convergence. This paper is a generalization of the method first introduced at the 12th International Conference on Domain Decomposition that was restricted to regular Cartesian grids. The potential of this method to provide scalable parallel computing on a geographically broad grid of parallel computers was demonstrated for some linear and nonlinear elliptic problems discretized by finite differences on a Cartesian mesh. The main purpose of this paper is to present a generalization of the method to nonuniform Cartesian meshes. The salient feature of the method consists of accelerating the sequence of traces on the artificial interfaces generated by the Schwarz procedure using a good approximation of the main eigenvectors of the trace transfer operator. For linear separable elliptic operators, our solver is a direct solver. For nonlinear operators, we use an approximation of the eigenvectors of the Jacobian of the trace transfer operator. The acceleration is then applied to the sequence generated by the Schwarz algorithm applied directly to the nonlinear operator. [ABSTRACT FROM AUTHOR]

In this paper, we propose a stochastic mixed multiscale finite element method. The proposed method solves the stochastic porous media flow equation on the coarse grid using a set of precomputed basis functions. The precomputed basis functions are constructed based on selected realizations of the stochastic permeability field, and furthermore the solution is projected onto the finite-dimensional space spanned by these basis functions. We employ multiscale methods using limited global information since the permeability fields do not have apparent scale separation. The proposed approach does not require any interpolation in stochastic space and can easily be coupled with interpolation-based approaches to predict the solution on the coarse grid. Numerical results are presented for permeability fields with normal and exponential variograms. [ABSTRACT FROM AUTHOR]

We present a local Fourier analysis of multigrid methods for the two-dimensional curl-curl formulation of Maxwell's equations. Both the hybrid smoother proposed by Hiptmair and the overlapping block smoother proposed by Arnold, Falk, and Winther are considered. The key to our approach is the identification of two-dimensional eigenspaces of the discrete curl-curl problem by decoupling the Fourier modes for edges with different orientations. This procedure is used to quantify the smoothing properties of the considered smoothers and the convergence behavior of the multigrid methods. Additionally, we identify the Helmholtz splitting in Fourier space. This allows several well known properties to be recovered in Fourier space, such as the commutation properties of the classical Nédélec prolongator and the equivalence of the curl-curl operator and the vector Laplacian for divergence-free vectors. We show how the approach used in this paper can be generalized to twoand three-dimensional problems in H(curl) and H(div) and to other types of regular meshes. [ABSTRACT FROM AUTHOR]

In this paper we derive explicit formulas for the normal form coefficients to verify the nondegeneracy of eight codimension two bifurcations of fixed points with one or two critical eigenvalues. These include all strong resonances, as well as the degenerate flip and Neimark--Sacker bifurcations. Applying our results to n-dimensional maps, one avoids the computation of the center manifold, but one can directly evaluate the critical normal form coefficients in the original basis. The formulas remain valid also for n=2 and allow one to avoid the transformation of the linear part of the map into Jordan form. The developed techniques are tested on two examples: (1) a three-dimensional map appearing in adaptive control; (2) a periodically forced epidemic model. We compute numerically the critical normal form coefficients for several codim 2 bifurcations occurring in these models. To compute the required derivatives of the Poincaré map for the epidemic model, the automatic differentiation package ADOL-C is used. [ABSTRACT FROM AUTHOR]

In this paper we apply the generalized Taylor­Galerkin finite element model to simulate bore wave propagation in a domain of two dimensions. For stability and accuracy reasons, we generalize the model through the introduction of four free parameters. One set of parameters is rig- orously determined to obtain the high-order finite element solution. The other set of free parameters is determined from the underlying discrete maximum principle to obtain the monotonic solutions. The resulting two models are used in combination through the flux correct transport technique of Zalesak, thereby constructing a finite element model which has the ability to capture hydraulic dis- continuities. In addition, this paper highlights the implementation of two Krylov subspace iterative solvers, namely, the bi-conjugate gradient stabilized (Bi-CGSTAB) and the generalized minimum residual (GMRES) methods. For the sake of comparison, the multifrontal direct solver is also con- sidered. The performance characteristics of the investigated solvers are assessed using results of a standard test widely used as a benchmark in hydraulic modeling. Based on numerical results, it is shown that the present finite element method can render the technique suitable for solving shallow water equations with sharply varying solution profiles. Also, the GMRES solver is shown to have a much better convergence rate than the Bi-CGSTAB solver, thereby saving much computing time compared to the multifrontal solver. [ABSTRACT FROM AUTHOR]

The Schwarz alternating method can be used to solve linear elliptic boundary value problems on domains which consist of two or more overlapping subdomains. The solution is approx- imated by an infinite sequence of functions which result from solving a sequence of elliptic boundary value problems in each of the subdomains. This paper considers four Schwarz alternating methods for the N-dimensional, steady, viscous, incompressible Navier­Stokes equations, N ≤ 4. It is shown that the Schwarz sequences converge to the true solution provided that the Reynolds number is sufficiently small. [ABSTRACT FROM AUTHOR]

In this paper, we examine different methods using techniques of blocking, buffering, and padding for efficient implementations of bit-reversals. We evaluate the merits and limits of each technique and its application and architecture-dependent conditions for developing cache-optimal methods. Besides testing the methods on different uniprocessors, we conducted both simulation and measurements on two commercial symmetric multiprocessors (SMP) to provide architectural insights into the methods and their implementations. We present two contributions in this paper: (1) Our integrated blocking methods, which match cache associativity and translation-lookaside buffer (TLB) cache size and which fully use the available registers, are cache-optimal and fast. (2) We show that our padding methods outperform other software-oriented methods, and we believe they are the fastest in terms of minimizing both CPU and memory access cycles. Since the padding methods are almost independent of hardware, they could be widely used on many uniprocessor workstations and multiprocessors. [ABSTRACT FROM AUTHOR]

This paper studies the solution of quadratic eigenvalue problems by the quadratic residual iteration method. The focus is on applications arising from finite-element simulations in acoustics. One approach is the shift-and-invert Arnoldi method applied to the linearized problem. When more than one eigenvalue is wanted, it is advisable to use locking or deflation of converged eigenvectors (or Schur vectors). In order to avoid unlimited growth of the subspace dimension, one can restart the method by purging unwanted eigenvectors (or Schur vectors). Both locking and restarting use the partial Schur form. The disadvantage of this approach is that the dimension of the linearized problem is twice that of the quadratic problem. The quadratic residual iteration and Jacobi­Davidson methods directly solve the quadratic problem. Unfortunately, the Schur form is not defined, nor are locking and restarting. This paper shows a link between methods for solving quadratic eigenvalue problems and the linearized problem. It aims to combine the benefits of the quadratic and the linearized approaches by employing a locking and restarting scheme based on the Schur form of the linearized problem in quadratic residual iteration and Jacobi­Davidson. Numerical experiments illustrate quadratic residual iteration and Jacobi­Davidson for computing the linear Schur form. It also makes a comparison with the shift-and-invert Arnoldi method. [ABSTRACT FROM AUTHOR]

In this paper we consider domain decomposition methods for solving the radial basis function interpolation equations. There are three interwoven threads to the paper. The first thread provides good ways of setting up and solving small- to medium-sized radial basis function interpolation problems. These may occur as subproblems in a domain decomposition solution of a larger interpolation problem. The usual formulation of such a problem can suffer from an unfortunate scale dependence not intrinsic in the problem itself. This scale dependence occurs, for instance, when fitting polyharmonic splines in even dimensions. We present and analyze an alternative formulation, available for all strictly conditionally positive definite basic functions, which does not suffer from this drawback, at least for the very important example previously mentioned. This formulation changes the problem into one involving a strictly positive definite symmetric system, which can be easily and efficiently solved by Cholesky factorization. The second section considers a natural domain decomposition method for the interpolation equations and views it as an instance of von Neumann's alternating projection algorithm. Here the underlying Hilbert space is the reproducing kernel Hilbert space induced by the strictly conditionally positive definite basic function. We show that the domain decomposition method presented converges linearly under very weak nondegeneracy conditions on the possibly overlapping subdomains. The last section presents some algorithmic details and numerical results of a domain decomposition interpolatory code for polyharmonic splines in 2 and 3 dimensions. This code has solved problems with 5 million centers and can fit splines with 10,000 centers in approximately 7 seconds on very modest hardware. [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]

Dynamical systems can often be decomposed into loosely coupled subsystems. The system of ordinary differential equations (ODEs) modelling such a problem can then be partitioned corresponding to the subsystems, and the loose couplings can be exploited by special integration methods to solve the problem using a parallel computer or just solve the problem more efficiently than by standard methods. This paper presents accuracy analysis of methods for the numerical integration of stiff partitioned systems of ODEs. The discretization formulas are based on the implicit Euler formula and the second order implicit backward differentiation formula (BDF2). Each subsystem of the partitioned problem is discretized independently, and the couplings to the other subsystems are based on solution values from previous time steps. Applied this way, the discretization formulas are called decoupled. The stability properties of the decoupled implicit Euler formula are well understood. This paper presents error bounds and asymptotic error expansions to be used in controlling step size, relaxation between subsystems and the validity of the partitioning. The decoupled BDF2 formula is analyzed within the same framework. Finally, the analysis is used in the design of a decoupled numerical integration algorithm with variable step size to control the local error and adaptive selection of partitionings. Two versions of the algorithm with decoupled implicit Euler and BDF2, respectively, are used in examples where a realistic problem is solved. The examples compare the results from the decoupled implicit Euler and BDF2 formulas, and compare them with results from the corresponding classical formulas. [ABSTRACT FROM AUTHOR]

In this paper, we consider a discontinuous Galerkin finite element method for first order initial-final value problems arising from optimal control of launch vehicles. In particular, we derive an a posteriori error estimate for this method that is subsequently implemented to provide a computational error estimate used for adaptive error control. We discuss several practical issues in the implementation of the computational error estimate as well. We test the theory on a real-life trajectory problem and find that the estimates are reliable and accurate while the adaptive error control leads to a significant gain in efficiency. [ABSTRACT FROM AUTHOR]

This paper describes a new approach for efficient Jacobian calculation using automatic differentiation. This approach obviates the need for vector accumulation by reducing certain parts of the computational graph representing the system of equations, allowing the Jacobian to be accumulated through a series of scalar operations. The advantages of this new approach include low spatial complexity, the ability to adapt to available memory if resources are limited, the ability to efficiently handle linear equations in a convenient manner, and low computational complexity. In addition, this approach is particularly well suited for evaluations within an interpretive environment. The approach can be applied to both the forward and reverses modes of automatic differentiation.

Considered are numerical integration schemes for nondissipative dynamical systems in which multiple time scales are present. It is assumed that one can do an explicit separation of the RHS ‘forces’ into fast forces and slow forces such that (i) the fast forces contain the high frequency part of the solution, (ii) the fast forces are conservative, and (iii) the reduced problem consisting only of the fast forces can be integrated much more cheaply than the full problem. The fast forces are allowed to have low frequency components. Particular applications for which the schemes are intended include N-body problems (for which most of the forces are slow) and nonlinear wave phenomena (for which the fast forces can be propagated by spectral methods). The assumption of cheap integration of fast forces implies that the overall cost of integration is primarily determined by the step size used to sample the slow forces. A long-time-step method is one in which this step size exceeds half the period of the fastest normal mode present in the full system. An existing method that comes close to qualifying is the ‘impulse’ method, also known as Verlet-I and r-RESPA. It is shown that it might fail, though, for a couple of reasons. First, it suffers a serious loss of accuracy if the step size is near a multiple of the period of a normal mode, and, second, it is unstable if the step size is near a multiple of half the period of a normal mode. Proposed in this paper is a ‘mollified’ impulse method having an error bound that is independent of the frequency of the fast forces. It is also shown to possess superior stability properties. Theoretical results are supplemented by numerical experiments. The method is efficient and reasonably easy to implement. [ABSTRACT FROM AUTHOR]

In this paper we consider how to compute the singular value decomposition (SVD) A = UΣVτ of A = [a1, a2] ∈ Rm×2 accurately in floating point arithmetic. It is shown how to compute the Jacobi rotation V (the right singular vector matrix) and how to compute AV = UΣ even if the floating point representation of V is the identity matrix. In the case . a1. 2 >> . a2. 2, underflow can produce the identity matrix as the floating point value of V, even for a1, a2 that are far from being mutually orthogonal. This can cause loss of accuracy and failure of convergence of the floating point implementation of the Jacobi method for computing the SVD. The modified Jacobi method recommended in this paper can be implemented as a reliable and highly accurate procedure for computing the SVD of general real matrices whenever the exact singular values do not exceed the underflow or overflow limits. [ABSTRACT FROM AUTHOR]

An approximate method for solving the Riemann problem is needed to construct Godunov schemes for relativistic hydrodynamical equations. Such an approximate Riemann solver is presented in this paper which treats all waves emanating from an initial discontinuity as themselves discontinuous. Therefore, jump conditions for shocks are approximately used for rarefaction waves. The solver is easy to implement in a Godunov scheme and converges rapidly for relativistic hydrodynamics. The fast convergence of the solver indicates the potential of a higher performance of a Godunov scheme in which the solver is used. [ABSTRACT FROM AUTHOR]

The computation of derivatives for optimizing time-dependent flow problems is often based on the integration of the adjoint differential equation. For this purpose, the knowledge of the complete forward solution is required. In the area of flow control, especially for three-dimensional problems, it may be impossible to keep track of the full forward solution due to the lack of storage capacities. One usual method to overcome this problem is checkpointing. Using a checkpointing strategy, only parts of the forward solution are kept in the main memory and additional time steps have to be performed. If one extends this approach such that some checkpoints can be stored on disc too, one can reduce the number of additional time steps. On the other hand, one has to take the access cost to one checkpoint into account. On parallel machines, one may use parallel I/O facilities to store and retrieve checkpoints. In these cases, the access cost of the checkpoints is also no longer negligible. Therefore, in this paper, the write and read counts for each checkpoint in a binomial checkpointing approach are examined. This way, the overall access cost to checkpoints is minimized. Numerical results illustrate the derived checkpointing approaches. They also show that checkpointing techniques may reduce the overall computing time despite the required recalculations. [ABSTRACT FROM AUTHOR]

This paper proposes new finite element spaces that can be constructed for agglomerates of standard elements that have certain regular structure. The main requirement is that the agglomerates share faces that have closed boundaries composed of 1-d edges. The spaces resulting from the agglomerated elements are subspaces of the original de Rham sequence of H1-conforming, H(curl)-conforming, H(div)-conforming, and piecewise constant spaces associated with an unstructured "fine" mesh. The procedure can be recursively applied so that a sequence of nested de Rham complexes can be constructed. As an illustration we generate coarser spaces from the sequence corresponding to the lowest-order Nédélec spaces, lowest-order Raviart-Thomas spaces, and for piecewise linear H1-conforming spaces, all in three dimensions. The resulting V-cycle multigrid methods used in preconditioned conjugate gradient iterations appear to perform similar to those of the geometrically refined case. [ABSTRACT FROM AUTHOR]

The linear sampling method is a qualitative procedure for the visualization of both impenetrable and inhomogeneous scatterers, which requires the regularized solution of a linear illposed integral equation of the first kind. An open issue in this technique is the one of determining the optimal scatterer profile from the visualization maps in an automatic manner. In the present paper this problem is addressed in two steps. First, linear sampling is optimized by using a new regularization algorithm for the solution of the integral equation, which provides more accurate maps for different levels of the noise affecting the data. Then an edge detection technique based on active contours is applied to the optimized maps. Our computation exploits a recently introduced implementation of the linear sampling method, which enhances both the accuracy and the numerical effectiveness of the approach. [ABSTRACT FROM AUTHOR]

In this paper, we consider the unilateral frictional contact problem of a hyperelastic body in the case of large displacements and small strains. In order to retain the linear elasticity framework, we decompose the deformation into a large global rotation and a small elastic displacement. This corotational approach is combined with a primal-dual active set strategy to tackle the contact problem. The resulting algorithm preserves both energy and angular momentum. [ABSTRACT FROM AUTHOR]

BPCONT, a driver for the software package AUTO for the numerical continuation of simple branch points of algebraic and boundary-value problems, is described in detail. Simple branch points are points in the continuation space where two solution branches intersect transversally. Detection and accurate location of branch points along the continuation of a solution branch and switching to the continuation of the other branch are implemented in AUTO and other software packages, but branch point continuation has not been fully supported. BPCONT fills this gap. It handles both generic problems, where branch point continuation is performed in two extra parameters, and nongeneric cases, where one extra parameter is typically involved due to problem-specific symmetries. This paper presents the extended algebraic and boundary-value problems which define branch point continuation, discusses their initialization at branch point detection and the symmetry-breaking technique used to automatically handle nongeneric cases, and describes the BPCONT implementation. Several examples, including generic and nongeneric algebraic problems, generic and nongeneric periodic boundary-value problems, and a nongeneric nonperiodic boundary-value problem, are also presented. [ABSTRACT FROM AUTHOR]

Spectral element methods on tetrahedra with symmetric collocation points can be accelerated by factorizing the discrete operators according to Hesthaven and Teng [SIAM J. Sci. Comput., 21 (2000), pp. 2352-2380]. While these authors focused on first-order conservation laws, the present paper provides an extension to second-order problems. Though factorization is easily accomplished for planar elements, difficulties arise from the presence of variable metric coefficients in curved tetrahedra. Two approaches are considered to cope with this peculiarity: (i) approximation of the metric terms by collocation projection, (ii) Gauss quadrature based on axisymmetric point sets. The first method achieves a separation of the metric terms such that the discrete operators can be reduced to factorizable standard matrices. As a consequence, the performance is comparable to the planar case, whereas the accuracy is limited by the projection step. The second approach maintains accuracy since all terms are evaluated individually in the quadrature points. Nonetheless, complete factorization is achieved by exploiting the symmetry in the quadrature points. Performance analysis shows that the curved element operator is less than three times as costly as the planar element counterpart. [ABSTRACT FROM AUTHOR]

In this paper we present a new incomplete factorization of a square matrix into triangular factors in which we get standard LU or LDLT factors (direct factors) and their inverses (inverse factors) at the same time. Algorithmically, we derive this method from the approach based on the Sherman-Morrison formula [R. Bru, J. Cerdán, J. Marín, and J. Mas, SIAM J. Sci. Comput., 25 (2003), pp. 701-715]. In contrast to the robust incomplete decomposition (RIF) algorithm [M. Benzi and M. Tůma, Numer. Linear Algebra Appl., 10 (2003), pp. 385-400] the direct and inverse factors here directly influence each other throughout the computation. Consequently, the algorithm to compute the approximate factors may mutually balance dropping in the factors and control their conditioning in this way. For the symmetric positive definite case, we derive the theory and present an algorithm for computing the incomplete LDLT factorization, and we discuss experimental results. We call this new approximate LDLT factorization the balanced incomplete factorization (BIF). Our experimental results confirm that this factorization is very robust and may be useful in solving difficult ill conditioned problems by preconditioned iterative methods. Moreover, the internal coupling of the computation of direct and inverse factors results in much shorter setup times (times to compute approximate decomposition) than RIF, a method of a similar and very high level of robustness. We also derive and present the theory for the general nonsymmetric case, but do not discuss its implementation. [ABSTRACT FROM AUTHOR]

In this paper, we propose an efficient and accurate numerical method for computing the ground state of spin-1 Bose-Einstein condensates (BECs) by using the normalized gradient flow or imaginary time method. The key idea is to find a third projection or normalization condition based on the relation between the chemical potentials so that the three projection parameters used in the projection step of the normalized gradient flow are uniquely determined by this condition as well as the other two physical conditions given by the conservation of total mass and total magnetization. This allows us to successfully extend the most popular and powerful normalized gradient flow or imaginary time method for computing the ground state of a single-component BEC to compute the ground state of spin-1 BECs. An efficient and accurate discretization scheme, the backward-forward Euler sine-pseudospectral method, is proposed to discretize the normalized gradient flow. Extensive numerical results on ground states of spin-1 BECs with ferromagnetic/antiferromagnetic interaction and harmonic/optical lattice potential in one/three dimensions are reported to demonstrate the efficiency of our new numerical method. [ABSTRACT FROM AUTHOR]

In this paper the numerical solutions of nonlinear neutral differential equations with constant delay and with proportional delay, two typical examples of neutral functional differential equations, are investigated. The focus is on the stability of numerical solutions produced by θ- methods. It is proved that if ½ ≤ θ? ≤ 1 the numerical solution to neutral differential equations with constant delay is globally stable and asymptotically stable under some conditions which are sufficient for stability and asymptotic stability of the exact solution. For neutral differential equations with proportional delay, two discretization approaches, the approach to transforming it into neutral equations with constant delay and the approach to directly discretizing it based on full-geometric mesh, are discussed. The integration on the first mesh interval [0, T0], which needs to be considered for both approaches, is also investigated. By taking into account the order of the methods, we construct two algorithms for the integration on [0, T0]. It is also proved that under some assumptions the numerical solutions to the neutral equation with proportional delay produced by the two approaches are globally stable and asymptotically stable if ½ ≤ θ ≤ 1. [ABSTRACT FROM AUTHOR]

We describe an asynchronous parallel derivative-free algorithm for linearly constrained optimization. Generating set search (GSS) is the basis of our method. At each iteration, a GSS algorithm computes a set of search directions and corresponding trial points and then evaluates the objective function value at each trial point. Asynchronous versions of the algorithm have been developed in the unconstrained and bound-constrained cases which allow the iterations to continue (and new trial points to be generated and evaluated) as soon as any other trial point completes. This enables better utilization of parallel resources and a reduction in overall run time, especially for problems where the objective function takes minutes or hours to compute. For linearly constrained GSS, the convergence theory requires that the set of search directions conforms to the nearby boundary. This creates an immediate obstacle for asynchronous methods where the definition of nearby is not well defined. In this paper, we develop an asynchronous linearly constrained GSS method that overcomes this difficulty and maintains the original convergence theory. We describe our implementation in detail, including how to avoid function evaluations by caching function values and using approximate lookups. We test our implementation on every CUTEr test problem with general linear constraints and up to 1000 variables. Without tuning to individual problems, our implementation was able to solve 95% of the test problems with 10 or fewer variables, 73% of the problems with 11-100 variables, and nearly half of the problems with 100-1000 variables. To the best of our knowledge, these are the best results that have ever been achieved with a derivative-free method for linearly constrained optimization. Our asynchronous parallel implementation is freely available as part of the APPSPACK software. [ABSTRACT FROM AUTHOR]

In this paper, we develop a numerical scheme for the interface problem in the planar symmetric radiative transfer equation with isotropic scattering. Such problems arise in the modeling of the propagation of energy density for waves in heterogeneous media with weak random fluctuation in the high frequency regime. The idea, following the earlier work of Jin and Wen for regular transmission and reflection, is to build the interface condition, which characterizes the reflection and transmission, into the numerical flux. The new contribution of this article is to deal with the diffuse transmission at the interface. The positivity of the scheme is proven. Moreover, we show, via numerical examples, that this new scheme is able to capture the correct (regular or diffuse) transmission and reflection through the interface, and, as the mean free path goes to zero, capture the diffusion limit with the correct interface conditions. [ABSTRACT FROM AUTHOR]

We discuss in this paper the numerical approximation of fluid-structure interaction (FSI) problems dealing with strong added-mass effect. We propose new semi-implicit algorithms based on inexact block-LU factorization of the linear system obtained after the space-time discretization and linearization of the FSI problem. As a result, the fluid velocity is computed separately from the coupled pressure-structure velocity system at each iteration, reducing the computational cost. We investigate explicit-implicit decomposition through algebraic splitting techniques originally designed for the FSI problem. This approach leads to two different families of methods which extend to FSI the algebraic pressure correction method and the Yosida method, two schemes that were previously adopted for pure fluid problems. Furthermore, we have considered the inexact factorization of the fluid-structure system as a preconditioner. The numerical properties of these methods have been tested on a model problem representing a blood-vessel system. [ABSTRACT FROM AUTHOR]

This paper presents a preconditioning method based on combining two-sided permutations with a multilevel approach. The nonsymmetric permutation exploits a greedy strategy to put large entries of the matrix in the diagonal of the upper leading submatrix. The method can be regarded as a complete pivoting version of the incomplete LU factorization. This leads to an effective incomplete factorization preconditioner for general nonsymmetric, irregularly structured, sparse linear systems. [ABSTRACT FROM AUTHOR]

This paper gives an overview of the use of polynomial chaos (PC) expansions to represent stochastic processes in numerical simulations. Several methods are presented for performing arithmetic on, as well as for evaluating polynomial and nonpolynomial functions of variables represented by PC expansions. These methods include Taylor series, a newly developed integration method, as well as a sampling-based spectral projection method for nonpolynomial function evaluations. A detailed analysis of the accuracy of the PC representations, and of the different methods for nonpolynomial function evaluations, is performed. It is found that the integration method offers a robust and accurate approach for evaluating nonpolynomial functions, even when very high-order information is present in the PC expansions. [ABSTRACT FROM AUTHOR]

In this paper we present an efficient branch-following procedure that can be used not only to compute branches of periodic solutions of periodically forced dynamical systems but also to determine the stability of the periodic solutions. The procedure combines Broyden's method with a subspace iteration method to determine the dominant eigenvalues. The method has connections with the hybrid Newton-Picard methods developed by Lust et al. in [SIAM 1. Sci. Comput., 19 (1998) pp. 1188-12091. A convergence analysis of the procedure is presented. The method is applied to the computation of periodic states of a reverse flow reactor, and its performance is compared with two variants of the hybrid Newton-Picard method. [ABSTRACT FROM AUTHOR]

In this paper, we develop a numerical scheme for the Smoluchowski coagulation equation, which relies on a conservative formulation and a finite volume approach. Several numerical simulations are performed to test the validity of the scheme and the expected behavior of the model. In particular the gelation phenomenon and the long time behavior of the solution are numerically studied. [ABSTRACT FROM AUTHOR]

Substantial effort has been focused over the last two decades on developing multi-level iterative methods capable of solving the large linear systems encountered in engineering practice. These systems often arise from discretizing partial differential equations over unstructured meshes, and the particular parameters or geometry of the physical problem being discretized may be unavailable to the solver. Algebraic multigrid (AMG) and multilevel domain decomposition methods of algebraic type have been of particular interest in this context because of their promises of optimal performance without the need for explicit knowledge of the problem geometry. These methods construct a hierarchy of coarse problems based on the linear system itself and on certain assumptions about the smooth components of the error. For smoothed aggregation (SA) methods applied to discretizations of elliptic problems, these assumptions typically consist of knowledge of the near-nullspace of the weak form. This paper introduces an extension of the SA method in which good convergence properties are achieved in situations where explicit knowledge of the near-nullspace components is unavailable. This extension is accomplished by using the method itself to determine near-nullspace components and adjusting the coarsening processes accordingly. [ABSTRACT FROM AUTHOR]

This paper is about the construction of a BGK Navier-Stokes (BGK-NS) solver in the discontinuous Galerkin (DG) framework. Since in the DG formulation the conservative variables and their slopes can be updated simultaneously, the flow evolution in each element involves only the flow variables in the nearest neighboring cells. Instead of using the semidiscrete approach in the Runge-Kutta discontinuous Galerkin (RKDG) method, the current DG-BGK method integrates the governing equations in time as well. Due to the coupling of advection and dissipative terms in the gas-kinetic formulation, the DG-BGK method solves the viscous governing equations directly. Numerical examples for the one-dimensional compressible Navier-Stokes solutions will be presented. [ABSTRACT FROM AUTHOR]

This paper addresses the problem of one-dimensional partitioning of structurally unsymmetric square and rectangular sparse matrices for parallel matrix-vector and matrix-transpose-vector multiplies. The objective is to minimize the communication cost while maintaining the balance on computational loads of processors. Most of the existing partitioning models consider only the total message volume hoping that minimizing this communication-cost metric is likely to reduce other metrics. However, the total message latency (start-up time) may be more important than the total message volume. Furthermore, the maximum message volume and latency handled by a single processor are also important metrics. We propose a two-phase approach that encapsulates all these four communication-cost metrics. The objective in the first phase is to minimize the total message volume while maintaining the computational-load balance. The objective in the second phase is to encapsulate the remaining three communication-cost metrics. We propose communication-hypergraph and partitioning models for the second phase. We then present several methods for partitioning communication hypergraphs. Experiments on a wide range of test matrices show that the proposed approach yields very effective partitioning results. A parallel implementation on a PC cluster verifies that the theoretical improvements shown by partitioning results hold in practice. [ABSTRACT FROM AUTHOR]

This paper presents a solver based on the HLLC (Harten-Lax-van Leer contact wave) approximate nonlinear Riemann solver for gas dynamics for the ideal magnetohydrodynamics (MHD) equations written in conservation form. It is shown how this solver also can be considered a modification of Linde's "adequate" solver. This approximation method is intended to resolve slow, Alfvén, and contact waves better than the original HLL (Harten-Lax-van Leer) solver. Compared to exact nonlinear solvers and Roe's solver, this new solver is computationally inexpensive. In addition, the method will exactly resolve isolated contacts and fast shocks. The method also preserves positive density and pressure with two caveats: first, the numerical signal velocities (the eigenvalues of the Roe average matrix) do not underestimate the physical signal velocities, and second, in a very few cases it may be required to change the wavespeeds of the Riemann fan for the underlying HLL method to guarantee positive pressures. These conditions are less restrictive on the definitions of the wavespeeds than the conditions needed to make the HLLC method positively conservative for gas dynamics. While the method is intended for a three-dimensional MHD problem, the simulation results concentrate on one-dimensional test cases. [ABSTRACT FROM AUTHOR]

This paper presents an approach to obtain a deformation which matches two images acquired from different medical imaging modalities. This problem arises in the investigation of human brains. Two distance functionals for the images are proposed with different pros and cons. These functionals are to be minimized. We add a smoothing term to the minimization problem which retains certain desired elastic features in the solution. At each minimization step an approximate solution for the linearized problem is computed with a multigrid method as an inner iteration. Furthermore, we use a multiresolution minimization approach to obtain a suitable initial guess. Finally, we present some experimental results for registration problems of synthetic images and for a real computer tomography (CT)­magnetic resonance imaging (MRI) registration. [ABSTRACT FROM AUTHOR]

This paper studies the numerical solution of two-dimensional nonlinear degenerate reaction-diffusion differential equations with singular forcing terms over rectangular domains. The equations considered may generate strong quenching singularities. This investigation focuses on a variable time step Peaceman­Rachford splitting method for the aforementioned problem. The time adaptation is implemented based on arc-length estimations of the first time derivative of the solution. The two-dimensional problem is split into several one-dimensional problems so that the computational cost is significantly reduced. The monotonicity and localized linear stability of the variable step scheme are investigated. We give some numerical examples to illustrate our results as well as to demonstrate the viability and efficiency of the method over existing methods for the quenching problem. It is also shown that the numerical solution obtained preserves important properties of the physical solution of the given problem. [ABSTRACT FROM AUTHOR]

In the Davidson method, any preconditioner can be exploited for the iterative computation of eigenpairs. However, the convergence of the eigenproblem solver may be poor for a high quality preconditioner. Theoretically, this counter-intuitive phenomenon with the Davidson method is remedied by the Jacobi­Davidson approach, where the preconditioned system is restricted to appropriate subspaces of codimension one. However, it is not clear how the restricted system can be solved accurately and efficiently in the case of a good preconditioner. The obvious approach introduces instabilities that hamper convergence. In this paper, we show how incomplete decomposition based on multilevel approaches can be used in a stable way. We also show how these preconditioners can be efficiently improved when better approximations for the eigenvalue of interest become available. The additional costs for updating the preconditioners are negligible. Furthermore, our approach leads to a good initial guess for the wanted eigenpair and to high quality preconditioners for nearby eigenvalues. We illustrate our ideas for the MRILU preconditioner. [ABSTRACT FROM AUTHOR]

In this paper, we present an efficient transformation technique for accurate numerical evaluation of weakly singular integrals with interior singularities. The present transformation technique does not require any division of the integration interval. It is composed of two parts of a sigmoidal transformation whose tails coincide with a singular point smoothly up to the order of the sigmoidal transformation employed. Therefore, in using the standard Gauss quadrature rule, the present transformation has a feature in which the integration points are clustered into the interior singular point very closely, as in the case of endpoint singular integrals. The results of some numerical examples show the superiority of the present method over the existing methods. [ABSTRACT FROM AUTHOR]

We introduce spectral element-based algebraic multigrid (ρAMGe), a new algebraic multigrid method for solving systems of algebraic equations that arise in Ritz-type finite element discretizations of partial differential equations. The method requires access to the element stiffness matrices, which enables accurate approximation of algebraically “smooth” vectors (i.e., error components that relaxation cannot effectively eliminate). Most other algebraic multigrid methods are based in some manner on predefined concepts of smoothness. Coarse-grid selection and prolongation, for example, are often defined assuming that smooth errors vary slowly in the direction of “strong” connections (relatively large coefficients in the operator matrix). One aim of ρAMGe is to broaden the range of problems to which the method can be successfully applied by avoiding any implicit premise about the nature of the smooth error. ρAMGe uses the spectral decomposition of small collections of element stiffness matrices to determine local representations of algebraically smooth error components. This provides a foundation for generating the coarse level and for defining effective interpolation. This paper presents a theoretical foundation for ρAMGe along with numerical experiments demonstrating its robustness. [ABSTRACT FROM AUTHOR]