Fast direct solvers for Poisson equation on 2D polar and spherical geometries

A simple and efficient class of FFT‐based fast direct solvers for Poisson equation on 2D polar and spherical geometries is presented. These solvers rely on the truncated Fourier series expansion, where the differential equations of the Fourier coefficients are solved by the second‐ and fourth‐order finite difference discretizations. Using a grid by shifting half mesh away from the origin/poles, and incorporating with the symmetry constraint of Fourier coefficients, the coordinate singularities can be easily handled without pole condition. By manipulating the radial mesh width, three different boundary conditions for polar geometry including Dirichlet, Neumann, and Robin conditions can be treated equally well. The new method only needs O(MN log₂ N) arithmetic operations for M × N grid points. © 2002 John Wiley & Sons, Inc. Numer Methods Partial Differential Eq 18: 56–68, 2002     

Preordering saddle‐point systems for sparse LDLT factorization without pivoting

This paper focuses on efficiently solving large sparse symmetric indefinite systems of linear equations in saddle‐point form using a fill‐reducing ordering technique with a direct solver. Row and column permutations partition the saddle‐point matrix into a block structure constituting a priori pivots of order 1 and 2. The partitioned matrix is compressed by treating each nonzero block as a single entry, and a fill‐reducing ordering is applied to the corresponding compressed graph. It is shown that, provided the saddle‐point matrix satisfies certain criteria, a block LDL^T factorization can be computed using the resulting pivot sequence without modification. Numerical results for a range of problems from practical applications using a modern sparse direct solver are presented to illustrate the effectiveness of the approach.     

Fast Iterative Solution of Saddle Point Problems in Optimal Control Based on Wavelets

In this paper, wavelet techniques are employed for the fast numerical solution of a control problem governed by an elliptic boundary value problem with boundary control. A quadratic cost functional involving natural norms of the state and the control is to be minimized. Firstly the constraint, the elliptic boundary value problem, is formulated in an appropriate weak form that allows to handle varying boundary conditions explicitly: the boundary conditions are treated by Lagrange multipliers, leading to a saddle point problem. This is combined with a fictitious domain approach in order to cover also more complicated boundaries.Deviating from standard approaches, we then use (biorthogonal) wavelets to derive an equivalent infinite discretized control problem which involves only ₂-norms and -operators. Classical methods from optimization yield the corresponding optimality conditions in terms of two weakly coupled (still infinite) saddle point problems for which a unique solution exists. For deriving finite-dimensional systems which are uniformly invertible, stability of the discretizations has to be ensured. This together with the ₂-setting circumvents the problem of preconditioning: all operators have uniformly bounded condition numbers independent of the discretization.In order to numerically solve the resulting (finite-dimensional) linear system of the weakly coupled saddle point problems, a fully iterative method is proposed which can be viewed as an inexact gradient scheme. It consists of a gradient algorithm as an outer iteration which alternatingly picks the two saddle point problems, and an inner iteration to solve each of the saddle point problems, exemplified in terms of the Uzawa algorithm. It is proved here that this strategy converges, provided that the inner systems are solved sufficiently well. Moreover, since the system matrix is well-conditioned, it is shown that in combination with a nested iteration strategy this iteration is asymptotically optimal in the sense that it provides the solution on discretization level J with an overall amount of arithmetic operations that is proportional to the number of unknows N _J on that level.Finally, numerical results are provided.     

Matrix decomposition algorithms for elliptic boundary value problems: a survey

We provide an overview of matrix decomposition algorithms (MDAs) for the solution of systems of linear equations arising when various discretization techniques are applied in the numerical solution of certain separable elliptic boundary value problems in the unit square. An MDA is a direct method which reduces the algebraic problem to one of solving a set of independent one-dimensional problems which are generally banded, block tridiagonal, or almost block diagonal. Often, fast Fourier transforms (FFTs) can be employed in an MDA with a resulting computational cost of O(N ² logN) on an N × N uniform partition of the unit square. To formulate MDAs, we require knowledge of the eigenvalues and eigenvectors of matrices arising in corresponding two–point boundary value problems in one space dimension. In many important cases, these eigensystems are known explicitly, while in others, they must be computed. The first MDAs were formulated almost fifty years ago, for finite difference methods. Herein, we discuss more recent developments in the formulation and application of MDAs in spline collocation, finite element Galerkin and spectral methods, and the method of fundamental solutions. For ease of exposition, we focus primarily on the Dirichlet problem for Poisson’s equation in the unit square, sketch extensions to other boundary conditions and to more involved elliptic problems, including the biharmonic Dirichlet problem, and report extensions to three dimensional problems in a cube. MDAs have also been used extensively as preconditioners in iterative methods for solving linear systems arising from discretizations of non-separable boundary value problems.     

Projected Schur complement method for solving non‐symmetric systems arising from a smooth fictitious domain approach

This paper deals with a fast method for solving large‐scale algebraic saddle‐point systems arising from fictitious domain formulations of elliptic boundary value problems. A new variant of the fictitious domain approach is analyzed. Boundary conditions are enforced by control variables introduced on an auxiliary boundary located outside the original domain. This approach has a significantly higher convergence rate; however, the algebraic systems resulting from finite element discretizations are typically non‐symmetric. The presented method is based on the Schur complement reduction. If the stiffness matrix is singular, the reduced system can be formulated again as another saddle‐point problem. Its modification by orthogonal projectors leads to an equation that can be efficiently solved using a projected Krylov subspace method for non‐symmetric operators. For this purpose, the projected variant of the BiCGSTAB algorithm is derived from the non‐projected one. The behavior of the method is illustrated by examples, in which the BiCGSTAB iterations are accelerated by a multigrid strategy. Copyright © 2007 John Wiley & Sons, Ltd.     

Direct solution method for the equilibrium problem for elastic stents

We are interested in numerical methods for approximating vector‐valued functions on a metric graph. As a model problem, we formulate and analyze a numerical method for the solution of the stationary problem for the one‐dimensional elastic stent model. The approximation is built using the mixed finite element method. The discretization matrix is a symmetric saddle‐point matrix, and we discuss sparse direct methods for the fast and robust solution of the associated equilibrium system. The convergence of the numerical method is proven and the error estimate is obtained. Numerical examples confirm the theoretical estimates.     

Matrix decomposition algorithms for the <Emphasis Type="Italic">C</Emphasis><Superscript>0</Superscript>-quadratic finite element Galerkin method

Explicit expressions for the eigensystems of one-dimensional finite element Galerkin (FEG) matrices based on C ⁰ piecewise quadratic polynomials are determined. These eigensystems are then used in the formulation of fast direct methods, matrix decomposition algorithms (MDAs), for the solution of the FEG equations arising from the discretization of Poisson’s equation on the unit square subject to several standard boundary conditions. The MDAs employ fast Fourier transforms and require O(N ²log N) operations on an N×N uniform partition. Numerical results are presented to demonstrate the efficacy of these algorithms.     

Interplay between the Inverse Scattering Method and Fokas's Unified Transform with an Application

It is known that the initial‐boundary value problem for certain integrable Partial Differential Equations (PDEs) on the half‐line with integrable boundary conditions can be mapped to a special case of the inverse scattering method (ISM) on the full‐line. This can also be established within the so‐called unified transform (UT) of Fokas for initial‐boundary value problems with linearizable boundary conditions. In this paper, we show a converse to this statement within the Ablowitz‐Kaup‐Newell‐Segur (AKNS) scheme: the ISM on the full‐line can be mapped to an initial‐boundary value problem with linearizable boundary conditions. To achieve this, we need a matrix version of the UT that was introduced by the author to study integrable PDEs on star‐graphs. As an application of the result, we show that the new, nonlocal reduction of the AKNS scheme introduced by Ablowitz and Musslimani to obtain the nonlocal nonlinear Schrödinger (NLS) equation can be recast as an old, local reduction, thus putting the nonlocal NLS and the NLS equations on equal footing from the point of view of the reduction group theory of Mikhailov.     

A fast Fourier transform based direct solver for the Helmholtz problem

This article is devoted to the efficient numerical solution of the Helmholtz equation in a two‐ or three‐dimensional (2D or 3D) rectangular domain with an absorbing boundary condition (ABC). The Helmholtz problem is discretized by standard bilinear and trilinear finite elements on an orthogonal mesh yielding a separable system of linear equations. The main key to high performance is to employ the fast Fourier transform (FFT) within a fast direct solver to solve the large separable systems. The computational complexity of the proposed FFT‐based direct solver is $\mathrm{??}(N\log N)$ operations. Numerical results for both 2D and 3D problems are presented confirming the efficiency of the method discussed.     

An effective direct solution method for certain boundary element equations in 3D

The boundary element method for the Dirichlet problem in a three-dimensional rotational domain leads to a system of linear equations with a full dense matrix having a special block structure. A direct solution method for such systems is presented, which requires O(N^3/2 ln N) arithmetical operations only, using a Fast Fourier Transformation (FFT), where N denotes the number of unknowns on the boundary surface.     

Fast Wavelet Transform for Toeplitz Matrices and Property Analysis

Fast wavelet transform algorithms for Toeplitz matrices are proposed in this paper. Distinctive from the well known discrete trigonometric transforms, such as the discrete cosine transform （DCT） and the discrete Fourier transform （DFT） for Toeplitz matrices, the new algorithms are achieved by compactly supported wavelet that preserve the character of a Toeplitz matrix after transform, which is quite useful in many applications involving a Toeplitz matrix. Results of numerical experiments show that the proposed method has good compression performance similar to using wavelet in the digital image coding. Since the proposed algorithms turn a dense Toeplitz matrix into a band-limited form, the arithmetic operations required by the new algorithms are O（N） that are reduced greatly compared with O（N log N） by the classical trigonometric transforms.     

Fast Fourier Transform Solvers and Preconditioners for Quadratic Spline Collocation

Quadratic Spline Collocation (QSC) methods of optimal order of convergence have been recently developed for the solution of elliptic Partial Differential Equations (PDEs). In this paper, linear solvers based on Fast Fourier Transforms (FFT)are developed for the solution of the QSC equations. The complexity of the FFT solvers is O(N ² log N), where N is the gridsize in one dimension. These direct solvers can handle PDEs with coefficients in one variable or constant, and Dirichlet, Neumann, alternating Dirichlet-Neumann or periodic boundary conditions, along at least one direction of a rectangular domain. General variable coefficient PDEs are handled by preconditioned iterative solvers. The preconditioner is the QSC matrix arising from a constant coefficient PDE. The convergence analysis of the preconditioner is presented. It is shown that, under certain conditions, the convergence rate is independent of the gridsize. The preconditioner is solved by FFT techniques, and integrated with one-step or acceleration methods, giving rise to asymptotically almost optimal linear solvers, with complexity O(N ² log N). Numerical experiments verify the effectiveness of the solvers and preconditioners, even on problems more general than the analysis assumes. The development and analysis of FFT solvers and preconditioners is extended to QSC equations corresponding to systems of elliptic PDEs.     

Adaptive wavelet schemes for an elliptic control problem with Dirichlet boundary control

An adaptive algorithm based on wavelets is proposed for the fast numerical solution of control problems governed by elliptic boundary value problems with Dirichlet boundary control. A quadratic cost functional representing Sobolev norms of the state and a regularization in terms of the control is to be minimized subject to linear constraints in weak form. In particular, the constraints are formulated as a saddle point problem that allows to handle the varying boundary conditions explicitly. In the framework of (biorthogonal) wavelets, a representer for the functional is derived in terms of ₂-norms of wavelet expansion coefficients and the constraints are written in form of an ₂ automorphism. Standard techniques from optimization are then used to deduce the resulting first order necessary conditions as a (still infinite) system in ₂. Applying the machinery developed in [8,9] which has been extended to control problems in [14], an adaptive method is proposed which can be interpreted as an inexact gradient method for the control. In each iteration step, in turn the primal and the adjoint saddle point system are solved up to a prescribed accuracy by an adaptive iterative Uzawa algorithm for saddle point problems which has been proposed in [10]. Under these premises, it can be shown that the adaptive algorithm containing now three layers of iterations is asymptotically optimal. This means that the convergence rate achieved for computing the solution up to a desired target tolerance is asymptotically the same as the wavelet-best N-term approximation of the solution, and the total computational work is proportional to the number of computational unknowns. AMS subject classification 65K10, 65N99, 93B40Angela Kunoth: This work has been supported partly by the Deutsche Forschungsgemeinschaft (SFB 611) at the Universität Bonn and by the European Communitys Human Potential Programme under contract HPRN-CT-2002-00286 Breaking Complexity.     

On realization of the Kreĭn–Langer class Nκ of matrix‐valued functions in Pontryagin spaces

In this paper the realization problems for the Kre?n–Langer class N_κ of matrix‐valued functions are being considered. We found the criterion when a given matrix‐valued function from the class N_κ can be realized as linear‐fractional transformation of the transfer function of canonical conservative system of the M. Livsic type (Brodskii–Livsic rigged operator colligation) with the main operator acting on a rigged Pontryagin space Π_κ with indefinite metric. We specify three subclasses of the class N_κ (R) of all realizable matrix‐valued functions that correspond to different properties of a realizing system, in particular, when the domains of the main operator of a system and its conjugate coincide, when the domain of the hermitian part of a main operator is dense in Π_κ . Alternatively we show that the class N_κ (R) can be realized as transfer matrix‐functions of some canonical impedance systems with self‐adjoint main operators in rigged spaces Π_κ . The case of scalar functions of the class N_κ (R) is considered in details and some examples are presented. (© 2008 WILEY‐VCH Verlag GmbH & Co. KGaA, Weinheim)     

Real valued iterative methods for solving complex symmetric linear systems

Complex valued systems of equations with a matrix R + 1S where R and S are real valued arise in many applications. A preconditioned iterative solution method is presented when R and S are symmetric positive semi‐definite and at least one of R, S is positive definite. The condition number of the preconditioned matrix is bounded above by 2, so only very few iterations are required. Applications when solving matrix polynomial equation systems, linear systems of ordinary differential equations, and using time‐stepping integration schemes based on Padé approximation for parabolic and hyperbolic problems are also discussed. Numerical comparisons show that the proposed real valued method is much faster than the iterative complex symmetric QMR method. Copyright © 2000 John Wiley & Sons, Ltd.     

Wavelet Techniques for the Fictitious-Domain-Lagrange-Multiplier-Approach

We consider second order elliptic boundary value problems when essential boundary conditions are enforced with the aid of Lagrange multipliers. This is combined with a fictitious domain approach into which the physical domain is embedded. The resulting saddle point problem will be discretized in terms of wavelets, resulting in an operator equation in ₂. Stability of the discretization and consequently the uniform boundedness of the condition number of the finite-dimensional operator independent of the discretization is guaranteed by an appropriate LBB condition. For the iterative solution of the saddle point system, an incomplete Uzawa algorithm is employed. It can be shown that the iterative scheme combined with a nested iteration strategy is asymptotically optimal in the sense that it provides the solution up to discretization error on discretization level J in an overall amount of iterations of order O(N _J ), where N _J is the number of unknowns on level J. Finally, numerical results are provided.     

Inverse spectral problem of fourth‐order boundary value problems with finite spectrum

In the present paper, we consider a class of inverse spectral problem of fourth‐order boundary value problems. Under the so‐called “Atkinson type” conditions, the problem has finite spectrum and corresponding matrix representations. By using the method of inverse matrix eigenvalue problems of two‐banded matrix, the leading coefficient and potential functions are reconstructed from the given three sets of interlacing real numbers satisfying certain conditions.     

Multilevel preconditioning on the refined interface and optimal boundary solvers for the Laplace equation

In this paper we propose and analyze some strategies to construct asymptotically optimal algorithms for solving boundary reductions of the Laplace equation in the interior and exterior of a polygon. The interior Dirichlet or Neumann problems are, in fact, equivalent to a direct treatment of the Dirichlet-Neumann mapping or its inverse, i.e., the Poincaré-Steklov (PS) operator. To construct a fast algorithm for the treatment of the discrete PS operator in the case of polygons composed of rectangles and regular right triangles, we apply the Bramble-Pasciak-Xu (BPX) multilevel preconditioner to the equivalent interface problem in theH ^1/2-setting. Furthermore, a fast matrix-vector multiplication algorithm is based on the frequency cutting techniques applied to the local Schur complements associated with the rectangular substructures specifying the nonmatching decomposition of a given polygon. The proposed compression scheme to compute the action of the discrete interior PS operator is shown to have a complexity of the orderO(N log ^q N),q [2, 3], with memory needsO(N log² N), whereN is the number of degrees of freedom on the polygonal boundary under consideration. In the case of exterior problems we propose a modification of the standard direct BEM whose implementation is reduced to the wavelet approximation applied to either single layer or hypersingular harmonic potentials and, in addition, to the matrix-vector multiplication for the discrete interior PS operator.     

Partial Fourier approximation of the Lamé equations in axisymmetric domains

In this paper, we study the partial Fourier method for treating the Lamé equations in three‐dimensional axisymmetric domains subjected to non‐axisymmetric loads. We consider the mixed boundary value problem of the linear theory of elasticity with the displacement û , the body force f̂ ϵ (L₂)³ and homogeneous Dirichlet and Neumann boundary conditions. The partial Fourier decomposition reduces, without any error, the three‐dimensional boundary value problem to an infinite sequence of two‐dimensional boundary value problems, whose solutions û _n (n = 0, 1, 2,…) are the Fourier coefficients of û . This process of dimension reduction is described, and appropriate function spaces are given to characterize the reduced problems in two dimensions. The trace properties of these spaces on the rotational axis and some properties of the Fourier coefficients û _n are proved, which are important for further numerical treatment, e.g. by the finite‐element method. Moreover, generalized completeness relations are described for the variational equation, the stresses and the strains. The properties of the resulting system of two‐dimensional problems are characterized. Particularly, a priori estimates of the Fourier coefficients û _n and of the error of the partial Fourier approximation are given. Copyright © 1999 John Wiley & Sons, Ltd.