Search

This paper establishes a posteriori error analysis for the Stokes equations discretized by an interior penalty type method using H(div) finite elements. The a posteriori error estimator is then employed for designing two grid refinement strategies; one is locally based and the other is globally based. The locally based refinement technique is believed to be able to capture local singularities in the numerical solution. The numerical formulations for the Stokes problem make use of H(div) conforming elements of Raviart-Thomas type. Therefore, the finite element solution features a full satisfaction of the continuity equation (mass conservation). The result of this paper provides a rigorous analysis for the method's reliability and efficiency. In particular, an H¹ norm a posteriori error estimator is obtained, together with upper and lower bound estimates. Numerical results are presented to verify the new theory of a posteriori error estimators. [ABSTRACT FROM AUTHOR]

This paper presents applications of an optimization-based multigrid method to state constrained aerodynamic shape optimization problems. The one-shot simultaneous pseudo-timestepping method has been applied successfully to such problems, reducing the cost of computation up to 70% from that of a traditional gradient method. However, the number of optimization iterations is comparatively large in this method since we update the design parameters in each time-step of state and costate runs. In this paper we use an optimization-based multigrid strategy to reduce the total number of optimization iterations while keeping the idea of the one-shot pseudo-time-stepping method. Design examples of drag reduction with constant lift, together with geometrical constraints, for an RAE2822 airfoil are included. The total effort required for convergence of the optimization problem is between four and six times that of forward simulation runs depending on computational grid size. [ABSTRACT FROM AUTHOR]

As a solution to the frequently reported problem of the high frequencies reconstructed sometimes many iterations after the low frequencies in tomography, we proposed in a recent paper [V. Israel-Jost, Ph. Choquet, and A. Constantinesco, Int. J. Biomed. Imaging, 2006, paper 34043.] a general adaptation that can apply to many algorithms: the incomplete backprojection. We stressed the fact that this type of frequency adapted (FA) algorithm works only when the linear problem of tomography to be solved involves deblurring, i.e., when the point spread functions (PSFs) associated with each voxel are large. Thus, depending on the desired level of detail, a parameter is set to achieve the reconstruction within a small number of iterations. The aim of this paper is to give a mathematical analysis of both the acceleration obtained and the convergence of an FA algorithm derived from the simultaneous algebraic reconstruction technique (SART) algorithm of Andersen and Kak, the so-called FA-SART algorithm. [ABSTRACT FROM AUTHOR]

The paper contains a noniterative solver for the Helmholtz and the modified Helmholtz equations in a hexahedron. The solver is based on domain decomposition. The solution domain is divided into mostly parallelepiped subdomains. In each subdomain a particular solution of the nonhomogeneous Helmholtz equation is first computed by a fast spectral 3-D method which was developed in our earlier papers (see, for example, SIAM J. Sci. Comput., 20 (1999), pp. 2237--2260). This method is based on the application of the discrete Fourier transform accompanied by a subtraction technique. For high accuracy the subdomain boundary conditions must be compatible with the specified inhomogeneous right-hand side at the edges of all the interfaces. In the following steps the partial solutions are hierarchically matched. At each step pairs of adjacent subdomains are merged into larger units. The paper describes in detail the matching algorithm for two boxes which is a basis for the domain decomposition scheme. The hierarchical approach is convenient for parallelization and can minimize the global communication. The algorithm requires O (N3, log, N) operations, where N is the number of grid points in each direction. [ABSTRACT FROM AUTHOR]

We develop a novel method for measuring the impact of random jumps in the expectation of functionals of jump-diffusion processes, without simulating sample paths of jumpdiffusion processes for Monte Carlo estimation. The proposed method constructs deterministic upper and lower bounding functions for the expectation of making full use of a similar structure in the infinitesimal generator of diffusion and jump-diffusion processes. Our hard bounding functions have the form of Markov-type inequalities, parametrized by the jump intensity parameter, so that the upper and lower bounding functions converge to each other when the jump intensity tends to zero. Hard bounding functions are derived also for the sensitivities of the expectation with respect to the jump intensity. The proposed method requires only well-developed numerical methods for PDEs (not partial integro-differential equations) and some elementary numerical integration of smooth functions. Numerical results are presented throughout the paper to support our theoretical developments and to illustrate the effectiveness of our methodology. [ABSTRACT FROM AUTHOR]

In this paper we present an improved DQDS algorithm for computing all the singular values of a bidiagonal matrix to high relative accuracy. There are two key contributions: a novel deflation strategy that improves the convergence for badly scaled matrices, and some modifications to certain shift strategies that accelerate the convergence for most bidiagonal matrices. These techniques together ensure linear worst case complexity of the improved algorithm (denoted by V5). Our extensive numerical experiments indicate that V5 is typically 1.2x-4x faster than DLASQ (the LAPACK-3.4.0 implementation of DQDS) without any degradation in accuracy. On matrices for which DLASQ shows very slow convergence, V5 can be 3x-10x faster. We develop a hybrid algorithm (HDLASQ) by combining our improvements with the aggressive early deflation strategy (AggDef2 in [Y. Nakatsukasa, K. Aishima, and I. Yamazaki, SIAM J. Matrix Anal. Appl., 33 (2012), pp. 22--51]). Numerical results show that HDLASQ is the fastest among these different versions. [ABSTRACT FROM AUTHOR]

In this work, we discuss the approximation of elliptic problems with log-normal random coefficients using the Wick product and the Mikulevicius-Rozovskii formula. The main idea is that the multiplication between the log-normal coefficient and the gradient of the solution can be regarded as a Taylor-like expansion in terms of the Wick product and the Malliavin derivative. For the classical model, the coefficients of Wiener chaos expansion are fully coupled together in the uncertainty propagator, while the Wick-Malliavin model yields an uncertain propagator with weak coupling in the upper-triangular part for a relatively small truncation order in the Mikulevicius-Rozovskii formula. In this paper we focus on the difference between the classical model and the Wick-Malliavin model with respect to the standard deviation of the underlying Gaussian random process. Both theoretical and numerical discussions are presented. [ABSTRACT FROM AUTHOR]

Many complex real-world systems can be accurately modeled by simulations. However, high-fidelity simulations may take hours or even days to compute. Because this can be impractical, a surrogate model is often used to approximate the dynamic behavior of the original simulator. This model can then be used as a cheap, drop-in replacement for the simulator. Because simulations can be very expensive, the data points, which are required to build the model, must be chosen as optimally as possible. Sequential design strategies offer a huge advantage over one-shot experimental designs because they can use information gathered from previous data points in order to determine I lie location of new data points. Each sequential design strategy must perform a trade-off between exploration and exploitation, where the former involves selecting data points in unexplored regions of the design space, while the latter suggests adding data points in regions which were previously identified to be interesting (for example, highly nonlinear regions). In this paper, a novel hybrid sequential design strategy is proposed which uses a Monte Carlo-based approximation of a Voronoi tessellation for exploration and local linear approximations of the simulator for exploitation. The advantage of this method over other sequential design methods is that it is independent of the model type, and can therefore be used in heterogeneous modeling environments, where multiple model types are used at the same time. The new method is demonstrated on a number of test problems, showing that it is a robust, competitive, and efficient sequential design strategy. [ABSTRACT FROM AUTHOR]

This paper is concerned with the numerical simulations for the dynamics of the molecular beam epitaxy (MBE) model. The numerical simulations of the MBE models require long time computations, and therefore large time-stepping methods become necessary. In this work, we consider some unconditionally energy stable finite difference schemes, which will be used in the time adaptivity strategies. It is found that the use of the time adaptivity cannot only resolve the steady-state solutions but also the dynamical changes of the solution accurately and efficiently. The adaptive time step is selected based on the energy variation or the change of the roughness of the solution. The numerical experiments demonstrated that the CPU time is significantly saved for long time simulations. [ABSTRACT FROM AUTHOR]

Estimation of the stiffness of a biological soft tissue is useful for the detection of pathologies such as tumors or atherosclerotic plaques. Elastography is a method based on the comparison between two images before and after a forced deformation of the tissue of interest. An inverse elasticity problem is then solved for Young's modulus estimation. In the case of arteries. no forced deformation is required, since vessels naturally move under the action of blood. Young's modulus can therefore be estimated by solving a coupled inverse fluid-structure interaction problem. In this paper we locus on the mathematical properties of this problem and its numerical solution. We give some well posedness analysis and some preliminary results based on a synthetic data set, i.e., test cases where the exact Young's modulus is known and the displacement dataset is numerically generated by solving a forward fluid-structure interaction problem. We address the problem of the presence of the noise in the measured displacement and of the proper sampling frequency for obtaining reliable estimates. [ABSTRACT FROM AUTHOR]

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

Position-dependent smoothness-increasing accuracy-conserving (SIAC) filtering is a promising technique not only ill improving the order of the numerical solution obtained by a discontinuous Galerkin (DG) method but also in increasing the smoothness of the field and improving the magnitude of the errors. This was initially established as an accuracy enhancement technique by Cockburn et al. for linear hyperbolic equations to handle smooth solutions [Math. Comp., 72 (2003), pp. 577 606]. By implementing this technique, the quality of the solution can be improved from order k + 1 to order 2k + 1 in the L2-norm. Ryan and Shu used these ideas to extend this technique to be able to handle postprocessing near boundaries ms well as discontinuities [Methods Appl. Anal., 10 (2003), pp. 295-307]. However, this presented difficulties as the resulting error had a stair-stepping effect and the errors themselves were not improved over those of the DG solution unless the mesh was suitably refined. In this paper, we discuss an improved filter for enhancing DG solutions that easily switches between one-sided postprocessing to handle boundaries or discontinuities and symmetric postprocessing for smooth regions. We numerically demonstrate that the magnitude of the errors using the modified postprocessor is roughly the same ms that of the errors for the symmetric postprocessor itself, regardless of the boundary conditions. [ABSTRACT FROM AUTHOR]

We develop in this paper some efficient algorithms which are essential to implementations of spectral methods on the sparse grid by Smolyak's construction based on a nested quadrature. More precisely, we develop a fast algorithm for the discrete transform between the values at the sparse grid and the coefficients of expansion in a hierarchical basis; and by using the aforementioned fast transform, we construct two very efficient sparse spectral-Galerkin methods for a model elliptic equation. In particular, the Chebyshev Legendre-Galerkin method leads to a sparse matrix with a much lower number of nonzero elements than that of low-order sparse grid methods based on finite elements or wavelets, and can be efficiently solved by a suitable sparse solver. Ample numerical results are presented to demonstrate the efficiency and accuracy of our algorithms. [ABSTRACT FROM AUTHOR]

In this paper a local Fourier analysis technique for multigrid methods on triangular grids is presented. The analysis is based on an expression of the Fourier transform in new coordinate systems, both in space variables and in frequency variables, associated with reciprocal bases. This tool makes it possible to study different components of the multigrid method in a very similar way to that of rectangular grids. Different smoothers for the discrete Laplace operator obtained with linear finite elements are analyzed. A new three-color smoother has been studied and has proven to be the best choice for "near" equilateral triangles. It is also shown that the block-line smoothers are more appropriate for irregular triangles. Numerical test calculations validate the theoretical predictions.. [ABSTRACT FROM AUTHOR]

In systems biology, the stochastic description of biochemical reaction kinetics is increasingly being employed to model gene regulatory networks and signaling pathways. Mathematically speaking, such models require the numerical solution of the underlying evolution equation, known as the chemical master equation (CME). Until now, the CME has primarily been treated by Monte Carlo techniques, the most prominent of which is the stochastic simulation algorithm [D. T. Gillespie, J. Comput. Phys., 22 (1976), pp. 403-434]. The paper presents an alternative, which focuses on the discrete partial differential equation (PDE) structure of the CME. This allows us to adopt ideas from adaptive discrete Galerkin methods as first suggested by Deuflhard and Wulkow [IMPACT Comput. Sci. Engrg., 1 (1989), pp. 269-301] for polyreaction kinetics and independently developed by Engblom. From the two different options for discretizing the CME as a discrete PDE, Engblom chose the method of lines approach (first space, then time), whereas we strongly advocate use of the Rothe method (first time, then space) for clear theoretical and algorithmic reasons. Numerical findings at two rather challenging problems illustrate the promising features of the proposed method and, at the same time, indicate lines of necessary further improvement of the method worked out here. [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 introduces a strategy for automatically generating a block preconditioner for solving the incompressible Navier–Stokes equations. We consider the ''pressure convection–diffusion preconditioners" proposed by Kay, Loghin, and Wathen [SIAM J. Sci. Comput., 24 (2002), pp. 237-256] and Silvester, Elman, Kay, and Wathen [J. Comput. Appl. Math., 128 (2001), pp. 261- 279]. Numerous theoretical and numerical studies have demonstrated mesh independent convergence on several problems and the overall efficacy of this methodology. A drawback, however, is that it requires the construction of a convection–diffusion operator (denoted Fp) projected onto the discrete pressure space. This means that integration of this idea into a code that models incompressible flow requires a sophisticated understanding of the discretization and other implementation issues, something often held only by the developers of the model. As an alternative, we consider automatic ways of computing Fp based on purely algebraic considerations. The new methods are closely related to the ''BFBt preconditioner'' of Elman [SIAM J. Sci. Comput., 20 (1999), pp. 1299-1316]. We use the fact that the preconditioner is derived from considerations of commutativity between the gradient and convection–diffusion operators, together with methods for computing sparse approximate inverses, to generate the required matrix Fp automatically. We demonstrate that with this strategy the favorable convergence properties of the preconditioning methodology are retained. [ABSTRACT FROM AUTHOR]

This paper presents the results of a comparison of three parallel algebraic multigrid (AMG) preconditioners for structural mechanics applications. In particular, we are interested in investigating both the scalability and robustness of the preconditioners. Numerical results are given for a range of structural mechanics problems with various degrees of difficulty, and they show that algebraic multigrid methods can be applied successfully even to difficult problems in structural mechanics. [ABSTRACT FROM AUTHOR]

As the development of the PSOR method (a point parallel SOR method by mesh domain partitioning proposed in [SIAM J. Sci. Comput., 20 (1999), pp. 2261-2281], this paper introduces a new mesh domain partition and ordering (the multitype partition and ordering), and proposes a new block parallel SOR (BPSOR) method for numerically solving 2-dimensional (2D) or three-dimensional (3D) elliptic boundary problems. A general mathematical analysis shows that the BPSOR method can have the same asymptotic convergence rate as the corresponding sequential block SOR method if the coefficient matrix of the block linear system is ''consistently ordered.'' It also shows that the original sequential ordering can be maintained in the parallel implementation of the BPSOR method so that the BPSOR method can be effectively applied to solve complex elliptic boundary problems. Furthermore, three particular multitype orderings are proposed based on strip and block mesh partitions, which lead to three effective BPSOR methods for solving the five-point like linear systems (in two dimensions) and the seven-point like linear systems (in three dimensions). In addition, it is shown that the PSOR method can be generated from BPSOR if each block equation is solved approximately by only one point SOR iteration. Thus, the PSOR method is well defined without involving any interprocessor data communication operations, and extended to solving three-dimensional (3D) problems. Finally, numerical results are presented which confirm the theoretical results and show that BPSOR has a good parallel performance on a parallel MIMD computer. [ABSTRACT FROM AUTHOR]

In this paper we present optimization algorithms for image restoration based on the total variation (TV) minimization framework of Rudin, Osher, and Fatemi (ROF). Our approach formulates TV minimization as a second-order cone program which is then solved by interior-point algorithms that are efficient both in practice (using nested dissection and domain decomposition) and in theory (i.e., they obtain solutions in polynomial time). In addition to the original ROF minimization model, we show how to apply our approach to other TV models, including ones that are not solvable by PDE-based methods. Numerical results on a varied set of images are presented to illustrate the effectiveness of our approach. [ABSTRACT FROM AUTHOR]

The mesh assessment problem is investigated in this paper by taking into account the shape and size of elements and the solution behavior. Three elementwise mesh quality measures characterizing the shape, alignment, and adaptation features of elements are introduced according to the estimates of interpolation error developed on a general mesh. An adaptive mesh is assessed by an overall quality measure defined as a weighted Lebesgue norm of a product of the three elementwise quality measures. It is shown that the overall quality of a mesh is good if the overall mesh quality measure is small or significantly smaller than the so-called roughness measure of the solution, defined as the ratio of two Lebesgue norms of a derivative of the solution. The definition of the overall mesh quality measure comes in such a way that the measure appears in the underlying error bound as the only factor depending substantially on the mesh. As an immediate result, the task of mesh adaptation becomes to control the overall mesh quality. This idea is applied to variational mesh adaptation to develop two functionals, one new and the other related to an existing functional recently developed using the regularity and equidistribution arguments. Numerical experiments are given to demonstrate the ability of the functionals to generate adaptive meshes of good quality. [ABSTRACT FROM AUTHOR]

In this paper, we are interested in the iterative solution of ill-conditioned Toeplitz systems generated by continuous nonnegative real-valued functions f with a finite number of zeros. We construct new w-circulant preconditioners without explicit knowledge of the generating function f by approximating f by its convolution f * KN with a suitable positive reproducing kernel KN. By the restriction to positive kernels we obtain positive definite preconditioners. Moreover, if f has only zeros of even order ≤ 2s, then we can prove that the property [This symbol cannot be presented in ASCII format]π-π t2k KN(t)dt ≤ CN-2k (k = 0, … ,s) of the kernel is necessary and sufficient to ensure the convergence of the PCG method in a number of iteration steps independent of the dimension N of the system. Our theoretical results were confirmed by numerical tests. [ABSTRACT FROM AUTHOR]

In this paper, a strategy is proposed for alternative computations of the residual vectors in Krylov subspace methods, which improves the agreement of the computed residuals and the true residuals to the level of O(u). A. . x. . Building on earlier ideas on residual replacement and on insights in the finite precision behavior of the Krylov subspace methods, computable error bounds are derived for iterations that involve occasionally replacing the computed residuals by the true residuals, and they are used to monitor the deviation of the two residuals and hence to select residual replacement steps, so that the recurrence relations for the computed residuals, which control the convergence of the method, are perturbed within safe bounds. Numerical examples are presented to demonstrate the effectiveness of this new residual replacement scheme. [ABSTRACT FROM AUTHOR]

In this paper, we present a meshless discretization technique for instationary convec- tion-diffusion problems. It is based on operator splitting, the method of characteristics, and a generalized partition of unity method. We focus on the discretization process and its quality. The method may be used as an h-version or a p-version. Even for general particle distributions, the convergence behavior of the different versions corresponds to that of the respective version of the finite element method on a uniform grid. We discuss the implementational aspects of the proposed method. Furthermore, we present the results of numerical examples, where we considered instationary convection-diffusion, instationary diffusion, linear advection, and elliptic problems. [ABSTRACT FROM AUTHOR]

This paper formulates several algorithms for the direct computation of periodic orbits as solutions of boundary value problems. The algorithms emphasize the use of coarse meshes and high orders of accuracy. Convergence theorems are given in the limit of increasing order with a fixed mesh. The algorithms are implemented with the use of MATLAB and ADOL-C, a software package for automatic differentiation. Automatic differentiation enables accurate computation of high-order derivatives of functions without the truncation errors inherent in finite difference calculations. We embed the algorithms in a continuation framework and extend them to compute saddle-node bifur- cations of periodic orbits directly. We present data from numerical studies of four test problems, making some comparisons with other methods for computing periodic orbits. These results demonstrate that high-order methods based upon automatic differentiation are capable of high precision with small meshes. [ABSTRACT FROM AUTHOR]

In this paper we formulate a time dependent model to approximate the solution to the nonlinear total variation optimization problem for deblurring and noise removal introduced by Rudin and Osher [Total variation based image restoration with free local constraints, in Proceedings IEEE Internat. Conf. Imag. Proc., IEEE Press, Piscataway, NJ, (1994), pp. 31–35] and Rudin, Osher, and Fatemi [Phys. D, 60 (1992), pp. 259–268], respectively. Our model is based on level set motion whose steady state is quickly reached by means of an explicit procedure based on Roe's scheme [J. Comput. Phys., 43 (1981), pp. 357–372], used in fluid dynamics. We show numerical evidence of the speed of resolution and stability of this simple explicit procedure in some representative 1D and 2D numerical examples. [ABSTRACT FROM AUTHOR]