A Fast Poisson Solver of Arbitrary Order Accuracy in Rectangular Regions

Abstract

Abstract. In this paper we propose a direct method for the solution of the Poisson equation in rectangular regions. It has an arbitrary order accuracy and low CPU requirements which makes it practical for treating large-scale problems.The method is based on a pseudospectral Fourier approximation and a polynomial subtraction technique. Fast convergence of the Fourier series is achieved by removing the discontinuities at the corner points using polynomial subtraction functions. These functions have the same discontinuities at the corner points as the sought solution. In addition to this, they satisfy the Laplace equation so that the subtraction procedure does not generate nonperiodic, nonhomogeneous terms.The solution of a boundary value problem is obtained in a series form in O(N log N ) floating point operations, where N 2 is the number of grid nodes. Evaluating the solution at all N 2 interior points requires O(N 2 log N ) operations.Key words. boundary value problem, Poisson equation, rectangular region, spectral method, corner discontinuities, polynomial subtraction AMS subject classifications. 35J05, 45L10, 65P05PII. S1064827595288589 1.Introduction. An important step in the development of fast numerical solvers for elliptic equations in complicated domains is an algorithm for the solution of boundary value problems for constant coefficient elliptic equations in rectangular regions. After solving a nonhomogeneous equation with some "convenient" boundary conditions, one must solve, in a correction step, the homogeneous problem with specified (Dirichlet or Neumann) boundary conditions. As a resulting boundary distribution may not satisfy the required compatibility conditions at the corner points, singularities may arise in the process of the solution, thus destroying the convergence rate even if the final solution is smooth.When the computational region is discretized on a N × N grid using some loworder (finite difference or finite element) scheme, the resulting system of linear algebraic equations is represented by a sparse matrix. Such a matrix can be inverted in O(N 2 ) (or O(N 2 log 2 N ) arithmetic operations using a "fast solver" [4]. However, since the method is of low order, the resolution N must be very large if a high accuracy is desired.Application of high-order (pseudo) spectral methods, based on global expansions into orthogonal polynomials, e.g., Chebyshev polynomials, to the solution of elliptic equations, results in a full matrix problem. The cost of inverting such a matrix using the best current algorithms is O(N 3 ) operations [2]. Besides, the accuracy decreases considerably as the dimension N 2 of the system grows due to accumulation of roundoff errors.When using the Chebyshev method, the fast computation of the expansion coefficients requires that the problem be discretized on a nonuniform grid. For timedependent problems, when elliptic equations arise due to a time-discretization pro-

A Comparison of Splittings and Integral Equation Solvers for a Nonseparable Elliptic Equation

“…Numerical examples, including comparison with previous results for simple setups as well as more challenging problems, are presented in Section 6. Our main conclusions are that -integral methods for (1,2) can be stably implemented on rectangular domains with high-order accuracy also for reformulations that involve the convolution with the gradient of the logarithmic kernel;…”

“…Constructing a fast direct solver for the entire problem (1,2), that is, a rapidly computable inverse or Greens' function, is difficult when σ is nonseparable. Fast schemes for this problem typically rely on a splitting of the operator into a supposedly dominant part which is used as a left-or right preconditioner and can be inverted with a fast solver, and a remaining part which can be viewed as a perturbation.…”

“…Possible splittings of (1), where the dominant operator appears on the left hand side, include ∆u = − ∇σ · ∇u σ + s σ ,…”

“…In the numerical examples we restrict our attention to the special case of Ω being a square. The paper is organized as follows: Section 2 gives a review of previous work on iterative schemes and fast solvers for (1,2). Interestingly, integral methods have been discussed to some extent but, seemingly, never been tested.…”

“…-different reformulations of (1,2) can lead to algorithms with very different convergence properties for rapidly varying σ. This seems not to have been noticed by earlier investigators.…”