Search

In this paper, we discuss the numerical solution of certain matrix-valued PDEs. Such PDEs arise, for example, when constructing a Riemannian contraction metric for a dynamical system given by an autonomous ODE. We develop and analyze a new meshfree discretization scheme using kernel-based approximation spaces. However, since these approximation spaces have now to be matrix-valued, the kernels we need to use are fourth-order tensors. We will review and extend recent results on even more general reproducing kernel Hilbert spaces. We will then apply this general theory to solve a matrix-valued PDE and derive error estimates for the approximate solution. The paper ends with applications to typical examples from dynamical systems. [ABSTRACT FROM AUTHOR]

This paper presents a bottom-up approach to derivation of orientation tensor closures for fiber suspension flow models. To begin with, we consider polynomial approximations based on the two-dimensional (2D) versions of the linear, quadratic, natural, and orthotropic smooth closures for reconstruction of the fourth-order orientation tensor. A numerical study is performed for simple flows. The investigation of planar closures provides new insights and boundary conditions for the design of orthotropic closures in three dimensions. The proposed extensions use finite element shape functions to interpolate the data at principal orientation states and additional points. The results for 3D simple flows indicate that natural closures based on (extended) quadratic and piecewise-linear interpolation provide a far better description of the 3D orientation dynamics than any other orthotropic closure considered in this study. [ABSTRACT FROM AUTHOR]

In this paper, we develop a numerical multiscale method to solve the fractional Laplacian with a heterogeneous diffusion coefficient. When the coefficient is heterogeneous, this adds to the computational costs. Moreover, the fractional Laplacian is a nonlocal operator in its standard form; however, the Caffarelli{Silvestre extension allows for a localization of the equations. This adds a complexity of an extra spacial dimension and a singular/degenerate coefficient depending on the fractional order. Using a subgrid correction method, we correct the basis functions in a natural weighted Sobolev space and show that these corrections are able to be truncated to design a computationally efficient scheme with optimal convergence rates. A key ingredient of this method is the use of quasi-interpolation operators to construct the fine scale spaces. Since the solution of the extended problem on the critical boundary is our main interest, we construct a projective quasi-interpolation that has both d and d + 1 dimensional averages over subsets in the spirit of the Scott{Zhang operator. We show that this operator satisfies local stability and local approximation properties in weighted Sobolev spaces. We further show that we can obtain a greater rate of convergence for sufficient smooth forces, utilizing a global L² projection on the critical boundary. We present some numerical examples, utilizing our projective quasi-interpolation in dimension 2 + 1 for analytic and heterogeneous cases to demonstrate the rates and effectiveness of the method. [ABSTRACT FROM AUTHOR]

In this paper, we study the problem of computing the effective diffusivity for a particle moving in chaotic and stochastic flows. In addition, we numerically investigate the residual diffusion phenomenon in chaotic advection. The residual diffusion refers to the nonzero effective (homogenized) diffusion in the limit of zero molecular diffusion as a result of chaotic mixing of the streamlines. In this limit, traditional numerical methods typically fail since the solutions of the advection-diffusion equations develop sharp gradients. Instead of solving the Fokker-Planck equation in the Eulerian formulation, we compute the motion of particles in the Lagrangian formulation, which is modeled by stochastic differential equations (SDEs). We propose an effective numerical integrator based on a splitting method to solve the corresponding SDEs in which the deterministic subproblem is symplectic preserving while the random subproblem can be viewed as a perturbation. We provide rigorous error analysis for the new numerical integrator using the backward error analysis technique and show that our method outperforms standard Euler-based integrators. Numerical results are presented to demonstrate the accuracy and efficiency of the proposed method for several typical chaotic and stochastic flow problems of physical interest. The existence of residual diffusivity for these flow problems is also investigated. [ABSTRACT FROM AUTHOR]

In this paper, we investigate the damped stochastic nonlinear Schrödinger (NLS) equation with multiplicative noise and its splitting-based approximation. When the damped effect is large enough, we prove that the solutions of both the damped stochastic NLS equation and the splitting scheme are exponentially stable and possess some exponential integrability. These properties show that the strong order of the scheme is 1 and independent of time. Additionally, we analyze the regularity of the Kolmogorov equation with respect to the stochastic NLS equation. As a consequence, the weak order of the scheme is shown to be 1 and independent of time. [ABSTRACT FROM AUTHOR]

We propose a new discretization method for the Stokes equations. The method is an improved version of the method recently presented in [C. Lehrenfeld and J. Schöberl, Comp. Meth. Appl. Mech. Eng., 361 (2016)] which is based on an H(div)-conforming finite element space and a hybrid discontinuous Galerkin (HDG) formulation of the viscous forces. H(div)-conformity results in favorable properties such as pointwise divergence-free solutions and pressure robustness. However, for the approximation of the velocity with a polynomial degree k, it requires unknowns of degree k on every facet of the mesh. In view of the superconvergence property of other HDG methods, where only unknowns of polynomial degree k-1 on the facets are required to obtain an accurate polynomial approximation of order k (possibly after a local postprocessing), this is suboptimal. The key idea in this paper is to slightly relax the H(div)-conformity so that only unknowns of polynomial degree k-1 are involved for normal continuity. This allows for optimality of the method also in the sense of superconvergent HDG methods. In order not to lose the benefits of H(div)-conformity, we introduce a cheap reconstruction operator which restores pressure robustness and pointwise divergence-free solutions and suits well to the finite element space with relaxed H(div)-conformity. We present this new method, carry out a thorough h-version error analysis, and demonstrate the performance of the method on numerical examples. [ABSTRACT FROM AUTHOR]

The improved Boussinesq equation is studied in this paper. Control properties for this equation posed on a bounded interval are first considered. When the control acts through the Dirichlet boundary condition the linearized system is proved to be approximately but not spectrally controllable. In a second part, the equation is posed on the one-dimensional torus and distributed moving controls are considered. Under some condition on the velocity at which the control moves, exact controllability results for both linear and nonlinear improved Boussinesq equations are obtained applying the moment method and a fixed point argument. [ABSTRACT FROM AUTHOR]

Let F be a compact set of a Banach space X. This paper analyzes the "generalized empirical interpolation method," which, given a function f ∊ F, builds an interpolant Jn[f] in an n dimensional subspace Xn ⊂ X with the knowledge of n outputs (σi(f))ni=1, where σi ∊ X' and X' is the dual space of X. The space Xn is built with a greedy algorithm that is adapted to F in the sense that it is generated by elements of F itself. The algorithm also selects the linear functionals (σi)ni=1 from a dictionary Σ ⊂ X'. In this paper, we study the interpolation error maxf∊F ‖f - Jn[f]‖X by comparing it with the best possible performance on an n dimensional space, i.e., the Kolmogorov n-width of F in X, dn(F,X). For polynomial or exponential decay rates of dn(F,X), we prove that the interpolation error has the same behavior modulo the norm of the interpolation operator. Sharper results are obtained in the case where X is a Hilbert space. [ABSTRACT FROM AUTHOR]

In this paper, a numerical solution of partial differential equations on the unit sphere is given by using a kernel trial approximation in combination with a special Petrov-Galerkin test discretization. The solvability of the scheme is proved, and the error bounds are obtained for functions in appropriate Sobolev spaces. The condition number of the final system is estimated in terms of discretization parameters. The method is meshless because in the trial side the numerical solution parameterizes entirely in terms of scattered points and in the test side everything breaks down to simple numerical integrations over independent spherical caps. This means that no connected background mesh is required for either approximation or integration. [ABSTRACT FROM AUTHOR]

As is well known, the smallest possible ratio between the spectral norm and the Frobenius norm of an m x n matrix with m x n is 1= p m and is (up to scalar scaling) attained only by matrices having pairwise orthonormal rows. In the present paper, the smallest possible ratio between spectral and Frobenius norms of n1 find tensors of order d, also called the best rank-one approximation ratio in the literature, is investigated. The exact value is not known for most configurations of n1 nd. Using a natural definition of orthogonal tensors over the real field (resp., unitary tensors over the complex field), it is shown that the obvious lower bound 1= p n1 nd1 is attained if and only if a tensor is orthogonal (resp., unitary) up to scaling. Whether or not orthogonal or unitary tensors exist depends on the dimensions n1; : : : ; nd and the field. A connection between the (non)existence of real orthogonal tensors of order three and the classical Hurwitz problem on composition algebras can be established: existence of orthogonal tensors of size ` x m x n is equivalent to the admissibility of the triple [`; m; n] to the Hurwitz problem. Some implications for higher-order tensors are then given. For instance, real orthogonal n n tensors of order d x 3 do exist, but only when n = 1; 2; 4; 8. In the complex case, the situation is more drastic: unitary tensors of size ` x m x n with ` x m x n exist only when `m x n. Finally, some numerical illustrations for spectral norm computation are presented. [ABSTRACT FROM AUTHOR]

This paper describes a method that can be used either to compute exactly or toapproximate, in an appropriate norm, the distance of a matrix pair (A;B) from a pair ( ~ A; ~B) for which the product ~ A~B T has rank strictly less than the maximum of the ranks of A or B. If A and B both have less than full column rank, then the exact distance is given by the smallest singular value of an augmented matrix that is parameterized by a single real parameter. The particular value of the parameter that gives the distance can be determined by a line search. Even when A or B has full column rank, it is still always possible to at least approximate the distance. The method and its properties are demonstrated with several small numerical examples. [ABSTRACT FROM AUTHOR]

This paper focuses on the numerical analysis for three-dimensional Bean's critical-state model in type-II superconductivity. We derive hyperbolic mixed variational inequalities of the second kind for the evolution Maxwell equations with Bean's constitutive law between the electric Field and the current density. On the basis of the variational inequality in the magnetic induction formulation, a semidiscrete Ritz-Galerkin approximation problem is rigorously analyzed, and a strong convergence result is proven. Thereafter, we propose a concrete realization of the Ritz-Galerkin approximation through a mixed Finite element method based on edge elements of Nédélec's First family, Raviart-Thomas face elements, divergence-free Raviart-Thomas face elements, and piecewise constant elements. As a Final result, we prove error estimates for the proposed mixed Finite element method. [ABSTRACT FROM AUTHOR]

This paper derives a posteriori error estimates for conforming numerical approximations of the Laplace eigenvalue problem with a homogeneous Dirichlet boundary condition. In particular, upper and lower bounds for an arbitrary simple eigenvalue are given. These bounds are guaranteed, fully computable, and converge with optimal speed to the given exact eigenvalue. They are valid without restrictions on the computational mesh or on the approximate eigenvector; we only need to assume that the approximate eigenvalue is separated from the surrounding smaller and larger exact ones, which can be checked in practice. Guaranteed, fully computable, optimally convergent, and polynomial-degree robust bounds on the energy error in the approximation of the associated eigenvector are derived as well, under the same hyp otheses. Remarkably, there appears no unknown (solution-, regularity-, or polynomial-degree-dependent) constant in our theory, and no convexity/regularity assumption on the computational domain/exact eigenvector(s) is needed. The multiplicative constant appearing in our estimates depends on (computable estimates of) the gaps to the surrounding exact eigenvalues. Its two improvements are presented. First, it is reduced by a fixed factor under an explicit, a posteriori calculable condition on the mesh and on the approximate eigenvector--eigenvalue pair. Second, when an elliptic regularity assumption on the corresponding source problem is satisfied with known constants, this multiplicative constant can be brought to the optimal value of one. Inexact algebraic solvers are taken into account; the estimates are valid on each iteration and can serve for the design of adaptive stopping criter ia. The application of our framework to conforming finite element approximations of arbitrary polynomial degree is provided, along with a numerical illustration on a set of test problems. [ABSTRACT FROM AUTHOR]

Multisecant quasi-Newton methods have been shown to be particularly suited to solve nonlinear fixed-point equations that arise from partitioned multiphysics simulations where the exact Jacobian is inaccessible. In all these methods, the underdetermined multisecant equation for the approximate (inverse) Jacobian is enhanced by a norm minimization condition. The standard choice is the minimization of the Frobenius norm of the approximate inverse Jacobian. In this setting, it is well known that transient fluid-structure simulations typically require the use of secant information also from previous time steps to achieve a small enough number of iterations per implicit time step. The number of these time steps highly depends on the application, the physical parameters, the used solvers, and the mesh resolution. Using too few leads to a relatively high number of iterations, while using too many leads not only to a computational overhead but also to an increase in the number of iterations as well. Determining the optimal number requires a costly trial-and-error process. In this paper, we present results for two different approaches to overcome this issue: The first approach is to use a modified method (presented in [F. Lindner, M. Mehl, K. Scheufele, and B. Uekermann, "A Comparison of Various Quasi-Newton Schemes for Partitioned Fluid-Structure Interaction," in Proceedings of ECCOMAS Coupled Problems, Venice, 2015, pp. 1-12]) that minimizes the Frobenius norm of the difference between the current approximate (inverse) Jacobian and that of the previous time step. Thus, previous time step information is taken into account in an implicit and automatized way without magic parameters. The second approach is to use an unrestricted number of previous time steps in combination with a suitable filtering algorithm automatically removing secant information that is outdated (thus, slowing down convergence) or contradicting newer information (deteriorating the condition of the multisecant equation system). We present a novel algorithm to realize the first idea with linear complexity in the number of coupling surface unknowns (note that already storing the approximate inverse Jacobian would induce quadratic complexity) and the efficient parallel implementation for both approaches. This results in highly efficient, parallelizable, and robust iterative solvers applicable for surface coupling in many types of multiphysics simulations using black-box solver software. In addition, our numerical results for the fluid-structure benchmark (FSI3) from Turek and Hron ["Proposal for Numerical Benchmarking of Fluid-Structure Interaction between an Elastic Object and Laminar Incompressible Flow," in Fluid-Structure Interaction, Springer, Berlin, 2006, pp. 371-385] and for a flow through a flexible tube with a large variety of parameter settings prove the robustness and numerical efficiency of the first approach in particular. The second approach can be shown to be highly sensitive to the choice of the filtering method for the secant information and not always robust for the filters we currently use. [ABSTRACT FROM AUTHOR]

We seek discrete approximations to solutions u : Ω → ℝ of semilinear elliptic PDE of the form Δu + fs (u) = 0, where fs is a one-parameter family of nonlinear functions and Ω is a domain in ℝd. The main achievement of this paper is the approximation of solutions to the PDE on the cube Ω = (0,π)³ ⊆ ℝ³. There are 323 possible isotropy subgroups of functions on the cube, which fall into 99 conjugacy classes. The bifurcations with symmetry in this problem are quite interesting, including many with 3-dimensional critical eigenspaces. Our automated symmetry analysis is necessary with so many isotropy subgroups and bifurcations among them, and it allows our code to follow one branch in each equivalence class that is created at a bifurcation point. Our most complicated result is the complete analysis of a degenerate bifurcation with a 6-dimensional critical eigenspace. This article extends our work in [Internat. J. Bifur. Chaos Appl. Sci. Engrg., 19 (2009), pp. 2531-2556], wherein we combined symmetry analysis with modified implementations of the gradient Newton-Galerkin algorithm (GNGA [J. M. Neuberger and J. W. Swift, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 11 (2001), pp. 801-820]) to automatically generate bifurcation diagrams and solution graphics for small, discrete problems with large symmetry groups. The code described in the current paper is efficiently implemented in parallel, allowing us to investigate a relatively fine-mesh discretization of the cube. We use the methodology and corresponding library presented in our previous paper [Int. J. Parallel Program., 40 (2012), pp. 443-464]. [ABSTRACT FROM AUTHOR]

This paper sets out to provide a general framework for the pricing of average-type options via lower and upper bounds. This class of options includes Asian, basket, and options on the volume-weighted average price. The use of lower and upper bounds is proposed in response to the inherent difficulty in finding analytical representations for the true price of these options and the requirement for numerical procedures to be fast and efficient. We demonstrate that in some cases lower bounds allow for the dimensionality of the problem to be reduced and that these methods provide reasonable approximations to the price of the option. [ABSTRACT FROM AUTHOR]

In this paper we derive a framework for multiscale approximation of elliptic problems on standard and mixed form. The method presented is based on a splitting into coarse and fine scales together with a systematic technique for approximation of the fine scale part, based on the solution of decoupled localized subgrid problems. The fine scale approximation is then used to modify the coarse scale equations. A key feature of the method is that symmetry in the bilinear form is preserved in the discrete system. Other key features are a posteriori error bounds and adaptive algorithms based on these bounds. The adaptive algorithms are used for automatic tuning of the method parameters. In the last part of the paper we present numerical examples where we apply the framework to a problem in oil reservoir simulation. [ABSTRACT FROM AUTHOR]

Let A ∈ ℂN x N be a nonsingular complex matrix and b and c be complex vectors of length N. The goal of this paper is to investigate approaches for efficient approximations of the bilinear form c*A-1b. Equivalently, we wish to approximate the scalar value c*x, where x solves the linear system Ax = b. Here the matrix A can be very large or its elements can be too costly to compute so that A is not explicitly available and it is used only in the form of the matrix-vector product. Therefore a direct method is not an option. For A Hermitian positive definite, b*A-1b can be efficiently approximated as a by-product of the conjugate-gradient iterations, which is mathematically equivalent to the matching moment approximations computed via the Gauss-Christoffel quadrature. In this paper we propose a new method using the biconjugate gradient iterations which is applicable to the general complex case. The proposed approach will be compared with existing ones using analytic arguments and numerical experiments. [ABSTRACT FROM AUTHOR]

This paper continues the convergence analysis in [M. Giles and S. Ulbrich, SIAM J. Numer. Anal., 48 (2010), pp. 882-904] of discrete approximations to the linearized and adjoint equations arising from an unsteady one-dimensional hyperbolic equation with a convex flux function. We consider a simple modified Lax-Friedrichs discretization on a uniform grid, and a key point is that the numerical smoothing increases the number of points across the nonlinear discontinuity as the grid is refined. It is proved that there is convergence in the discrete approximation of linearized output functionals even for Dirac initial perturbations and pointwise convergence almost everywhere for the solution of the adjoint discrete equations. In particular, the adjoint approximation converges to the correct uniform value in the region in which characteristics propagate into the discontinuity. Moreover, it is shown that the results of [M. Giles and S. Ulbrich, SIAM J. Numer. Anal., 48 (2010), pp. 882-904] and the present paper hold also for quite general nonlinear initial data which contain multiple shocks and for which shocks form at a later time and/or merge. [ABSTRACT FROM AUTHOR]

This paper studies a new recovery-based a posteriori error estimator for the conforming linear finite element approximation to elliptic interface problems. Instead of recovering the gradient in the continuous finite element space, the flux is recovered through a weighted L2 projection onto H(div) conforming finite element spaces. The resulting error estimator is analyzed by establishing the reliability and efficiency bounds and is supported by numerical results. This paper also proposes an adaptive finite element method based on either the recovery-based estimators or the edge estimator through local mesh refinement and establishes its convergence. In particular, it is shown that the reliability and efficiency constants as well as the convergence rate of the adaptive method are independent of the size of jumps. [ABSTRACT FROM AUTHOR]

This paper studies a new recovery-based a posteriori error estimator for the conforming linear finite element approximation to elliptic interface problems. Instead of recovering the gradient in the continuous finite element space, the flux is recovered through a weighted L2 projection onto H(div) conforming finite element spaces. The resulting error estimator is analyzed by establishing the reliability and efficiency bounds and is supported by numerical results. This paper also proposes an adaptive finite element method based on either the recovery-based estimators or the edge estimator through local mesh refinement and establishes its convergence. In particular, it is shown that the reliability and efficiency constants as well as the convergence rate of the adaptive method are independent of the size of jumps. [ABSTRACT FROM AUTHOR]

In (SIAM J. Matrix Anal. Appl., 23 (2002), pp. 706–727, we developed numerical algorithms for computing sparse low-rank approximations of matrices, and we also provided a detailed error analysis of the proposed algorithms together with some numerical experiments. The low-rank approximations are constructed in a certain factored form with the degree of sparsity of the factors controlled by some user-specified parameters. In this paper, we cast the sparse low-rank approximation problem in the framework of penalized optimization problems. We discuss various approximation schemes for the penalized optimization problem which are more amenable to numerical computations. We also include some analysis to show the relations between the original optimization problem and the reduced one. We then develop a globally convergent discrete Newton-like iterative method for solving the approximate penalized optimization problems. We also compare the reconstruction errors of the sparse low-rank approximations computed by our new methods with those obtained using the methods in the earlier paper and several other existing methods for computing sparse low-rank approximations. Numerical examples show that the penalized methods are more robust and produce approximations with factors which have fewer columns and are sparser. [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]

A collection of algorithms is described for numerically computing with smooth functions defined on the unit sphere. Functions are approximated to essentially machine precision by using a structure-preserving iterative variant of Gaussian elimination together with the double Fourier sphere method. We show that this procedure allows for stable differentiation, reduces the oversampling of functions near the poles, and converges for certain analytic functions. Operations such as function evaluation, differentiation, and integration are particularly efficient and can be computed by essentially one-dimensional algorithms. A highlight is an optimal complexity direct solver for Poisson's equation on the sphere using a spectral method. Without parallelization, we solve Poisson's equation with 100 million degrees of freedom in 1 minute on a standard laptop. Numerical results are presented throughout. In a companion paper (part II) we extend the ideas presented here to computing with functions on the disk. [ABSTRACT FROM AUTHOR]

Max-plus methods have previously been used to solve deterministic control problems. The methods are based on max-plus (or min-plus) expansions and can yield curse-of-dimensionalityfree numerical methods. In this paper, we explore min-plus methods for continuous-time stochastic control on a finite-time horizon. We first approximate the original value function via timediscretization. By generalizing the min-plus distributive property to continuum spaces, we obtain an algorithm for recursive computation of the time-discretized values, which we refer to as the idempotent distributed dynamic programming principle (IDDPP). Under the IDDPP, the value function at each step can be represented as an infimum of functions in a certain class. This is a min-plus expansion for the value function. For the specific class of problems considered here, we see that the class can be taken as that consisting of the quadratic functions. A means for reducing the numbers of constituent quadratic functions is discussed. [ABSTRACT FROM AUTHOR]

In the field of uncertainty quantification (UQ), propagation of uncertainty is one of the most important tasks, as it is essential to almost all UQ analysis. While many numerical methods have been developed for uncertainty propagation in probabilistic framework, much less has been discussed in the nonprobabilistic framework, which is important whenever one does not possess sufficient data or knowledge of the underlying systems. In this paper, we focus on the use of fuzzy sets to model uncertainty and their propagation through physical systems. In particular, we develop a numerical strategy that can efficiently propagate fuzzy sets in complex systems. The method utilizes an accurate approximation model for the solution over the support of the input fuzzy sets and then retrieves the output fuzzy set information via the approximation model. By doing so the method becomes highly efficient, as the only simulation cost is in the construction of the approximation model. In particular, we discuss the use of orthogonal polynomials for the approximation model. This is very similar to the probabilistic method of generalized polynomial chaos. However, no probability is assumed and the method remains nonprobabilistic. We prove that the fuzzy sets propagated by the approximation model can be accurate and establish the convergence via a new distance measure between fuzzy sets. Several numerical examples are provided to demonstrate the effectiveness of the method and verify the convergence analysis. [ABSTRACT FROM AUTHOR]

In this paper we are concerned with the solution of compressed sensing (CS) problems where the signals to be recovered are sparse in coherent and redundant dictionaries. We extend the primal-dual Newton Conjugate Gradient method (pdNCG) in [T. F. Chan, G. H. Golub, and P. Mulet, SIAM J. Sci. Comput., 20 (1999), pp. 1964-1977] to CS problems. We provide an inexpensive and provably effective preconditioning technique for linear systems using pdNCG. Numerical results are presented on CS problems which demonstrate the performance of pdNCG with the proposed preconditioner compared to state-of-the-art existing solvers. [ABSTRACT FROM AUTHOR]

This paper concerns the accurate evaluation of the principal value integral governing axisymmetric vortex sheet motion. Previous quadrature rules for this integral lose accuracy near the axis of symmetry. An approximation by de Bernadinis and Moore (dBM) that converges pointwise at the rate of O(h3) has maximal errors near the axis that are O(h). As a result, the discretization error is not smooth. It contains high wavenumber frequencies that make it difficult to resolve the vortex sheet motion. This paper explains the reason for the degeneracy near the axis and proposes a modified quadrature rule that is uniformly O(h3). The results are based on an analytic approximation of the integrand, whose integral can be precomputed. The modification is implemented at negligible additional cost per timestep. As an example, it is applied to compute the evolution of an initially spherical vortex sheet. [ABSTRACT FROM AUTHOR]

Generalized Gaussian quadratures appear to have been introduced by Markov late in the last century and have been studied in great detail as a part of modern analysis. They have not been widely used as a computational tool, in part due to an absence of effective numerical schemes for their construction. Recently, a numerical scheme for the design of such quadratures was introduced by Ma et al.; numerical results presented in their paper indicate that such quadratures dramatically reduce the computational cost of the evaluation of integrals under certain conditions. In this paper, we modify their approach, improving the stability of the scheme and extending its range of applicability. The performance of the method is illustrated with several numerical examples. [ABSTRACT FROM AUTHOR]

In recent years the growing popularity of supercomputers has fostered the development of algorithms able to take advantage of the massive parallelism offered by multiple processors. Direct methods, though robust and computationally efficient, hardly exploit high degrees of parallelism. By contrast, Krylov methods preconditioned by Factored Sparse Approximate Inverses (FSAI) provide, at least in principle, a perfectly parallel approach but are often thwarted by an excessive set-up cost. In this paper we extend the concept of supernode from sparse LU factorizations to approximate inverses, and use it to accelerate the computation of an FSAI-type preconditioner. The numerical experiments on real-world problems show that the overall FSAI efficiency can be significantly increased while preserving its intrinsic parallelism. [ABSTRACT FROM AUTHOR]

We consider the one-dimensional one phase Stefan problems with temperature or thermal flux boundary conditions. Exact and closed-form solutions for these two free boundary problems are very limited and highly restricted to particular initial and boundary conditions, and so are the similarity solutions. In general, only numerical methods allow one to completely solve such nonlinear problems. In this paper, the quasi-similarity solutions are defined and, for each problem, an analytical approximation of the interfacial position is given as function of the time-dependent boundary datum. [ABSTRACT FROM AUTHOR]

All solutions of a matrix interpolation problem with purely imaginary interpolation points are given, by a descriptor system realization, if necessary and sufficient conditions are satisfied. The first step is to reduce the interpolation problem to a minimal realization. No assumptions are needed in the paper. [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]

In this paper quasi-Monte Carlo (QMC) methods are applied to a class of elliptic partial differential equations (PDEs) with random coefficients, where the random coefficient is parametrized by a countably infinite number of terms in a Karhunen--Loeve expansion. Models of this kind appear frequently in numerical models of physical systems, and in uncertainty quantification. The method uses a QMC method to estimate expected values of linear functionals of the exact or approximate solution of the PDE, with the expected value considered as an infinite dimensional integral in the parameter space corresponding to the randomness induced by the random coefficient. The error analysis, arguably the first rigorous application of the modern theory of QMC in weighted spaces, exploits the regularity with respect to both the physical variables (the variables in the physical domain) and the parametric variables (the parameters corresponding to randomness). In the weighted-space theory of QMC methods, "weights," describing the varying difficulty of different subsets of the variables, are introduced in order to make sure that the high-dimensional integration problem is tractable. It turns out that the weights arising from the present analysis are of a nonstandard kind, being of neither product nor order dependent form, but instead a hybrid of the two--we refer to these as "product and order dependent weights," or "POD weights" for short. These POD weights are of a simple enough form to permit a component-by-component construction of a randomly shifted lattice rule that has optimal convergence properties for the given weighted space setting. If the terms in the expansion for the random coefficient have an appropriate decay property, and if we choose POD weights that minimize a certain upper bound on the error, then the solution of the PDE belongs to the joint function space needed for the analysis, and the QMC error (in the sense of a root-mean-square error averaged over shifts) is of order O(N-1-δ) for arbitrary δ > 0, where N denotes the number of sampling points in the parameter space. Moreover, this convergence rate and slower convergence rates under more relaxed conditions are achieved under conditions similar to those found recently by Cohen, De Vore, and Schwab [Found. Comput. Math., 10 (2010), pp. 615-646] for the same model by best N-term approximation. We analyze the impact of a finite element (FE) discretization on the overall efficiency of the scheme, in terms of accuracy versus overall cost, with results that are comparable to those of the best N-term approximation. [ABSTRACT FROM AUTHOR]

The present paper is concerned with the unilateral contact model in linear elastostatics, the so-called Signorini problem. (Our results can also be applied to the scalar Signorini problem.) A standard continuous linear finite element approximation is first chosen to approach the two-dimensional problem. We develop a new error analysis in the H¹-norm using estimates on Poincaré constants with respect to the size of the areas of the noncontact sets. In particular we do not assume any additional hypothesis on the finiteness of the set of transition points between contact and noncontact. This approach allows us to establish better error bounds under sole Hτ assumptions on the solution: if 3/2 < τ < 2 we improve the existing rate by a factor h(τ-3/2)² and if τ = 2 the existing rate (h¾) is improved by a new rate of h √"ln(h)|. Using the same finite element spaces as previously we then consider another discrete approximation of the (nonlinear) contact condition in which the same kind of analysis leads to the same convergence rates as for the first approximation. [ABSTRACT FROM AUTHOR]

In this paper we construct an "interior penalty" discontinuous Galerkin method to approximate the minimizer of a variational problem related to the p(x)-Laplacian. The function p : Ω → [p1,p2] is log-Hölder continuous and 1 < p1 ≤ p2 < ∞. We prove that the minimizers of the discrete functional converge to the solution. We also make some numerical experiments in dimension one to compare this method with the conforming Galerkin method, in the case where p1 is close to one. This example is motivated by its applications to image processing. [ABSTRACT FROM AUTHOR]

In this paper we present several numerical results performed with a fully deterministic scheme to discretize the Boltzmann equation of rarefied gas dynamics in a bounded domain for multiscale problems. Periodic, specular reflection and diffusive boundary conditions are discussed and investigated numerically. The collision operator is treated by a Fourier approximation of the collision integral, which guarantees spectral accuracy in velocity with a computational cost of MN log(N), where N is the number of degrees of freedom in velocity space andM represents the number of discrete angles of the collision kernel. This algorithm is coupled with a second order finite volume scheme in space and a time discretization allowing us to deal for rarefied regimes as well as their hydrodynamic limit. Our numerical results show that the proposed approach significantly improves the near-wall nonstationary flow accuracy of standard numerical methods over a wide range of Knudsen numbers, in particular when the solution to the Boltzmann equation is close to the local equilibrium and for slow motion flows. [ABSTRACT FROM AUTHOR]

In this paper we present PDE and finite element analyses for a system of PDEs consisting of the Darcy equation and the Cahn--Hilliard equation, which arises as a diffuse interface model for the two-phase Hele-Shaw flow. In the model the two sets of equations are coupled through an extra phase induced force term in the Darcy equations and a fluid induced transport term in the Cahn--Hilliard equation. We propose a fully discrete implicit finite element method for approximating the PDE system, which consists of the implicit Euler method combined with a convex splitting energy strategy for the temporal discretization, the standard finite element discretization for the pressure, and a split (or mixed) finite element discretization for the fourth-order Cahn--Hilliard equation. It is shown that the proposed numerical method satisfies a mass conservation law in addition to a discrete energy law that mimics the basic energy law for the Darcy--Cahn--Hilliard phase field model and holds uniformly in the phase field parameter ε. With the help of the discrete energy law, we first prove that the fully discrete finite method is unconditionally energy stable and uniquely solvable at each time step. We then show that, using the compactness method, the finite element solution has an accumulation point that is a weak solution of the PDE system. As a result, the convergence result also provides a constructive proof of the existence of global-in-time weak solutions to the Darcy--Cahn--Hilliard phase field model in both two and three dimensions. Numerical experiments based on the overall solution method of combining the proposed finite element discretization and a nonlinear multigrid solver are presented to validate the theoretical results and to show the effectiveness of the proposed fully discrete finite element method for approximating the Darcy--Cahn--Hilliard phase field model. [ABSTRACT FROM AUTHOR]

This paper deals with the approximation of the eigenvalues of evolution operators for linear retarded functional differential equations through the reduction to finite dimensional operators by a pseudospectral collocation. Fundamental applications such as determination of asymptotic stability of equilibria and periodic solutions of nonlinear autonomous retarded functional differential equations follow at once. Numerical tests are provided. [ABSTRACT FROM AUTHOR]

This paper is on the construction of energy-minimizing coarse spaces that obey certain functional constraints and can thus be used, for example, to build robust coarse spaces for elliptic problems with large variations in the coefficients. In practice they are built by patching together solutions to appropriate local saddle point or eigenvalue problems. We develop an abstract framework for such constructions, akin to an abstract Bramble--Hilbert-type lemma, and then apply it in the design of coarse spaces for discretizations of PDEs with highly varying coefficients. The stability and approximation bounds of the constructed interpolant are in the weighted L2 norm and are independent of the variations in the coefficients. Such spaces can be used, for example, in two-level overlapping Schwarz algorithms for elliptic PDEs with large coefficient jumps generally not resolved by a standard coarse grid or for numerical upscaling purposes. Some numerical illustration is provided. [ABSTRACT FROM AUTHOR]

This paper presents two methods based on domain decomposition concepts for determining the diagonal of the inverse of specific matrices. The first uses a divide-and-conquer principle and the Sherman-Morrison-Woodbury formula and assumes that the matrix can be decomposed into a 2 x 2 block-diagonal matrix and a low-rank matrix. The second method is a standard domain decomposition approach in which local solves are combined with a global correction. Both methods can be successfully combined with iterative solvers and sparse approximation techniques. The efficiency of the methods usually depends on the specific implementation, which should be fine-tuned for different test problems. Preliminary results for some two-dimensional (2D) problems are reported to illustrate the proposed methods. [ABSTRACT FROM AUTHOR]

The paper deals with the effective conductivity tensor Kef of anisotropic random media subject to mean uniform flux. The hydraulic conductivity K field is modeled as a collection of spheroidal disjoint inclusions of different, isotropic and independent Y = ln K; the latter is a random variable with given distribution of variance Due to image rights restrictions, multiple line equation(s) cannot be graphically displayed. Inclusions are embedded in homogeneous background of anisotropic conductivity K0. The K field is anisotropic, characterized by the anisotropy ratio , ratio of the vertical and horizontal integral scales of K. We derive Kef by accurate numerical simulations; the numerical model for anisotropic media is presented here for the first time, and it generalizes a previously developed model for isotropic formations [I. Jankovic, A. Fiori, and G. Dagan, Multiscale Model. Simul.,1(2003), pp. 40-56]. The numerical model is capable of solving complex three-dimensional flow problems with high accuracy for any configuration of the spheroidal inclusions and arbitrary K distribution. The numerically derived Kef for a normal is compared with its prediction by (i) the self-consistent solution Ksc, (ii) the first-order approximation in s , and (iii) the exponential conjecture [L. J. Gelhar and C. L. Axness Water. Resour. Res., 19 (1983), pp. 161-180]. It is found that the self-consistent solution Ksc is very accurate for a broad range of the values of the parameters s ,f and for the densest inclusions packing. In contrast, the first-order solution strongly deviates from Kef for large Due to image rights restrictions, multiple line equation(s) cannot be graphically displayed., as expected, and the exponential conjecture is generally unable to correctly represent the effective conductivity. The numerical solution for the potential is expressed as an infinite series of spheroidal harmonics, attached to the interior and exterior of each inclusion, which accounts for the nonlinear interaction between neighboring inclusions. [ABSTRACT FROM AUTHOR]

In this paper we study numerical approximations of continuous solutions to the nonlocal p-Laplacian type diffusion equation, ut(t, x) = ∫Ω J(x - y)|u(t, y) - u(t, x)|p-2(u(t, y) - u(t,x)) dy. First, we find that a semidiscretization ill space of this problem gives rise to an ODE system whose solutions converge uniformly to the continuous one as the mesh size goes to zero. Moreover, the semidiscrete approximation shares some properties of the continuous problem: it preserves the total mass and the solution converges to the mean value of the initial condition as t goes to infinity. Next, we also discretize the time variable and present a totally discrete method which also enjoys the above mentioned properties. In addition, we investigate the limit as p goes to infinity in these approximations and obtain a discrete model for the evolution of a sandpile. Finally, we present some numerical experiments that illustrate our results. [ABSTRACT FROM AUTHOR]

This paper is devoted to the numerical approximation of Lyapunov and Sacker-Sell Spectral intervals for linear differential-algebraic equations (DAEs). The spectral analysis for DAEs is improved and the concepts of leading directions and solution subspaces associated with spectral intervals are extended to DAEs. Numerical methods based on smooth singular value decompositions are introduced for Computing all or only some spectral intervals and their associated leading directions. The numerical algorithms as well as implementation issues are discussed in detail and Numerical examples are presented to illustrate the theoretical results. [ABSTRACT FROM AUTHOR]

This paper investigates backscatter data for the inverse obstacle problem in impedance tomography, when the obstacles are small. It is shown that under this circumstance the backscatter data are close approximations of a rational function that has second or fourth order poles at the locations of the obstacles. Furthermore, a numerical method is presented to locate the obstacles via the poles of certain Laurent-Padé approximations. Numerical experiments explore the potential of this algorithm also for extended obstacles, taking into consideration that the problem is severely ill-posed. [ABSTRACT FROM AUTHOR]

In this paper we develop a novel rule for choosing regularization parameters in non-smooth Tikhonov functional. It is solely based on the value function and applicable to a broad range of nonsmooth models, and it extends one known criterion. A posteriori error estimates of the approximations are derived. An efficient numerical algorithm for computing the minimizer is developed, and its convergence properties are discussed. Numerical results for several common nonsmooth models are presented, including deblurring natural images. The numerical results indicate the rule can yield results comparable with those achieved with the discrepancy principle and the optimal choice, and the algorithm merits a fast and steady convergence. [ABSTRACT FROM AUTHOR]

In this paper we propose a rigorous numerical technique for the computation of symmetric connecting orbits for ordinary differential equations. The idea is to solve a projected boundary value problem (BVP) in a function space via a fixed point argument. The formulation of the projected BVP involves a high order parameterization of the invariant manifolds at the steady states. Using this parameterization, one can obtain explicit exponential asymptotic bounds for the coefficients of the expansion of the manifolds. Combining these bounds with piecewise linear approximations, one can construct a contraction in a function space whose unique fixed point corresponds to the wanted connecting orbit. We have implemented the method to demonstrate its effectiveness, and we have used it to prove the existence of a family of even homoclinic orbits for the Gray-Scott equation. [ABSTRACT FROM AUTHOR]

In this paper, an Ulm-like method is proposed for solving inverse singular value problems. This method has an advantage over Newton's methods since it avoids solving approximate Jacobian equations. Under some mild assumptions, we show that the proposed method converges at least quadratically in the root sense. Our numerical tests, based on comparison with the inexact Newton method given by Bai and Xu [Linear Algebra Appl., 429 (2008), pp. 527-547], demonstrate the effectiveness of the new method. [ABSTRACT FROM AUTHOR]

This paper studies existence and uniqueness results and interlacing properties of nonlinear modifications of small rank of symmetric eigenvalue problems. Approximation properties of the Rayleigh functional are used to design numerical methods, the local convergence of which is quadratic or even cubic. Numerical examples demonstrate their efficiency. [ABSTRACT FROM AUTHOR]

A new algorithm is developed for computing etAB, where A is an n x n matrix and B is n x n0 with n0 « n. The algorithm works for any A, its computational cost is dominated by the formation of products of A with n x n0 matrices, and the only input parameter is a backward error tolerance. The algorithm can return a single matrix etA B or a sequence Due to image rights restrictions, multiple line equation(s) cannot be graphically displayed. on an equally spaced grid of points tk. It uses the scaling part of the scaling and squaring method together with a truncated Taylor series approximation to the exponential. It determines the amount of scaling and the Taylor degree, using the recent analysis of Al-Mohy and Higham [SIAM J. Matrix Anal. Appl., 31 (2009), pp. 970 989], which provides sharp truncation error bounds expressed in terms of the quantities ||Ak