## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# New Laplace and Helmholtz solvers

Edited by David L. Donoho, Stanford University, Stanford, CA, and approved April 18, 2019 (received for review March 11, 2019)

## Abstract

Numerical algorithms based on rational functions are introduced that solve the Laplace and Helmholtz equations on 2D domains with corners quickly and accurately, despite the corner singularities.

The Laplace and Helmholtz equations are the basic partial differential equations (PDEs) of potential theory and acoustics, respectively. Suppose a domain Ω bounded by a polygon P is given and (to begin with the Laplace case) we seek the unique function

The standard techniques for solving such a problem numerically are the finite element method (FEM) (1) and boundary integral equations (2). However, these methods face a challenge in calculating accurate solutions because of singularities at the corners (3). The mathematical basis of an algorithm for meeting this challenge, whose computational implications have not been noticed before, is a theorem in the field of approximation theory published by D. J. Newman in 1964 (4). Newman considered the approximation of

Our algorithm achieves root-exponential convergence for solving PDEs by approximating **1**] with root-exponential convergence. To find such functions computationally, one exploits the fact that since the poles are prescribed the problem is linear. Expansion coefficients are found by least-squares fitting in sample points on the boundary, a routine problem of linear algebra involving a matrix of dimensions about

Our method can be regarded as a variant of the method of fundamental solutions (MFS), in which solutions are also approximated via finite sums (6). The MFS differs from our approach in that each term is normally a monopole or point charge rather than a dipole (an exception in the context of Maxwell equations is ref. 7), it is not normally applied with exponential clustering (an exception is ref. 8), and it is not founded upon results of rational approximation theory (although rational functions are used in ref. 5).

Fig. 1 shows a computed solution in an L-shaped domain. For a problem that is generic in the sense of having singularities at the corners, the simplest choice of boundary data is

The numerical solution of PDEs has been at the heart of scientific computing since computers were invented, and the Laplace equation in a planar domain is as fundamental a problem in this area as any. There are two main classes of methods for solving such problems: FEM and boundary integral equations. To gather information on how our approach compares with these, in November 2018 we posed the problem of Fig. 1 as a challenge to the international numerical analysis community via the email list NA Digest (9). Specifically, we asked for a computation of

The other well-known methods are boundary integral equations, which are advantageous because of good conditioning and because the solution is represented in one dimension, along the boundary, rather than two dimensions in the domain (2). Here one begins by solving an integral equation to determine a density function

Whereas integral equations make use of a continuous distribution of dipoles **1**] regardless of z, requiring no special calculations when z is near a boundary arc or a corner. One may regard this as a zero-dimensional representation of the solution, whose complexity is determined only by the complexity of the function being represented, not by that of Ω or P.

Our method may be generalized to other elliptic PDEs, notably the Helmholtz equation **1**] to

There is a historical context that may shed light on the paucity of previous literature on solving PDEs via rational functions and their generalizations. The mathematicians of the 19th century concentrated on functions that were analytic or piecewise analytic, that is, smooth and representable by convergent Taylor series. Most physical applications are of this kind. In the 20th century, however, mathematicians turned to new challenges of less-smooth functions, developing advanced tools for analyzing fine distinctions of regularity (i.e., smoothness). Overwhelmingly, such tools became the standard framework for computational mathematics, too, and in particular FEM experts almost invariably derive and analyze their algorithms in the language of Sobolev spaces (1), in which precise distinctions are made, for example, between a function with one derivative of smoothness and a function with one and a half derivatives. This kind of analysis entails a bias toward low-accuracy methods tuned to problems with limited smoothness. The simplicity and speed of the methods proposed here are a reminder that there is also a place for numerical analysis based on less pessimistic smoothness assumptions.

In closing, let us summarize some pros and cons of the method. It is limited to problems of small or medium size in two dimensions (at least as currently developed) and to certain constant-coefficient PDEs with known special solutions. It faces challenges of ill-conditioning that are not yet understood and do not arise with boundary integral equations, and since the convergence is only root-exponential, it would not be competitive at very high accuracies (say,

## Acknowledgments

We thank Alex Barnett, Timo Betcke, Leslie Greengard, Vladimir Rokhlin, and Kirill Serkh for advice.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: trefethen{at}maths.ox.ac.uk.

Author contributions: A.G. and L.N.T. designed research; A.G. and L.N.T. performed research; and L.N.T. wrote the paper.

The authors declare no conflict of interest.

- Copyright © 2019 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution License 4.0 (CC BY).

## References

- ↵
- Schwab C

- ↵
- Serkh K,
- Rokhlin V

- ↵
- ↵
- Newman DJ

- ↵
- Hochman A,
- Leviatan Y,
- White JK

- ↵
- Barnett AH,
- Betcke T

- ↵
- Eremin Yu A,
- Sveshnikov AG

- ↵
- Liu Y

- ↵
- Trefethen LN

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Physical Sciences
- Applied Mathematics