1 Introduction
Recently, problems of the Laplace type and conformal mapping, which can be shown to reduce to a particular form of a Laplace problem Tr19, have been addressed by the strikingly efficient software implementation based on rational approximations called the ”lightning Laplace solver” (LLS) in domain with corners GT19b, and the ”Vandermonde with Arnoldi” algorithm for the stabilization of matrix computation associated with power sums in regions with smooth boundaries Tr19b. The relevant literature also shows how it is crucial that when singularities like corners are present, sample points be exponentially clustered in their neighbourhood to achieve rootexponential convergence GT19a; GT19b. Our method starts from the aim of giving a solution also for those cases, frequently found in physics and engineering, where only a set of values associated to a generic ”reasonably discretized” contour is given, and results with an accuracy of a fraction of a percent are sufficient. Being based on the flexible Adaptive AntoulasAnderson (AAA) algorithm NOT19, it can help whenever parametric boundaries or explicit function definitions are not available.
This paper is focused on coding. After a brief presentation of the method, illustrating the underlying idea and its mathematical definition, we give a collection of selected problems along with the listings for their solution.
2 Method and Algorithm
The AAA algorithm computes rational barycentric approximations of scalar functions defined on subsets of the complex plane, also returning their poles, residues and zeros: this allows representations in terms of partial fractions. Figure 1 shows the poles returned after the approximation of on an Lshaped boundary discretized by evenly sampled points.
We first notice that poles are arranged in a somewhat exponentially clustered fashion on both sides of the boundary to represent corners, irrespective of sampling frequency; the trouble is that they are spread everywhere with no chance to exclude them from a particular region. Roughly speaking, if the whole set provides a tolerance of e.g. 1e13, half their number should achieve about 1e6, which is still decent indeed, so the idea is as simple as this: if we want to solve a Laplace problem with Dirichlet boundary condition in the inside region, we must drop the inner poles; if we want to solve outside, then we must drop the outer ones: in other words, we drop the subset impairing analyticity. Clearly this plain concept alone is not sufficient, as trivial examples can show, therefore a replacement for the excluded subset need be found.
Our method, from a formal viewpoint, has much in common with the LLS fully described in GT19b: given a realvalued function defined on a contour in the complex plane, and dividing it into an interior region and an exterior one , find a function analytic in
(1) 
for which and : often in applications these are referred to as the scalar potential and the stream function respectively, of the complex potential . is the centre of a power sum describing the ”smooth” part of the problem, placed in about the middle of the geometry, whilst is the subset of poles in shaping the ”singular” part, namely contour corners and irregular function behaviours. The interior region is defined with respect to : we let be the region located on the left of the boundary when traversing it counterclockwise. Conversely, if we want a solution in we take and with a negative value in formula (1), resulting in an additional arbitrary pole placed in the interior region. The magnitude of is discussed later.
From now on, we pursue a different path from GT19b: the idea is to build the two parts and in sequence. First, a point is chosen which lies about in the middle of the region, and are determined by solving the leastsquares (LS) problem
(2)  
(3)  
(4) 
where is the given set of points describing , and the corresponding function values. The real singular part of the problem is pulled out by simply subtracting
(5) 
Then is approximated by the AAA algorithm, unwanted poles are dropped, and a second LS problem is solved in order to determine :
(6)  
(7)  
(8) 
The reason for there being two constants in (1) is that adjusts the value of in the smooth part, in order to approximate better with the second LS fitting. The parameter in (1) deserves particular attention. One usually works with sample points (almost) regularly spaced in applications, but picking a value with the same order of magnitude as sample number would imply an overwhelming amount of terms with unacceptable computing time. Even worse, too large a number would compromise the subsequent approximation of the singular part: the difference between and its first LS fitting along the boundary would be a function severely affected by the Gibbs phenomenon, not restricted to corner neighbourhoods. Figure 2 illustrates such situation, where can not be approximated correctly by the AAA algorithm. Therefore, is determined by taking the logarithm of sample number: it is an engineering choice, typically giving rise to some 10 terms in the smooth part, and we don’t claim it to be optimal though in a number of cases it has appeared to work well. In cases of concave shapes close to a circle (a square, a regular polygon, an oval, or the like) an even smaller value need be used.
We are now ready to give an overview of the algorithm. We exploit the recently announced Arnoldi stabilization for Vandermonde matrices Tr19b, which provides superior stability to the illconditioned computation of the smooth part in the first step:

Approximate the smooth part by means of a power sum, analytic in the desired region;

Pull the singular part out of , and approximate it by means of the AAA algorithm;

Use the poles returned by the AAA approximation in the other region to solve the second LS problem for ;
In the next sections, detailed solutions to some problems of the Laplace type and conformal mapping are given. Some general remarks:

Large portions of the code are repeated in different listings: we always give the full solution ready for cutandpaste into the computing environment. This widespread uniformity and standardization, we believe, will make the reader appreciate the generality of the method.

Boundary is always represented starting from a uniform discretization step, to give evidence of the AAA algorithm’s capability of clustering poles as necessary, thus partially compensating any lack of exponential sampling close to critical points.

The AAA algorithm is called with a high maximum degree of 1000 as the default value of 100 is seldom sufficient to reach the default tolerance of 1e13. In all the proposed examples an approximation with about 200 poles is obtained;

Lawson iterations for minimax approximation are very likely to cause poles to reposition badly, so they must be turned off. If doubts arise about the possibility of getting a number of zeros greater than the number of poles, one iteration only should be used.

The maximum error on the boundary is computed for available samples, and accuracy can not be confirmed by checking on a twice or four times finer mesh as the LLS does, therefore in general it can not be considered an accurate upper bound. Experiments show that a much bigger error is to be expected throughout the domain if no particular care is taken in the neighbourhood of critical points.
The following are needed to run the examples:

A MATLAB 6 or Octave 4.2 environment, or newer;

The Chebfun package DHT14, carrying the most recent version of the AAA algorithm;

The polyfitA.m and polyvalA.m functions described in Tr19b, for stable computations of power sums; these must be found in the environment path;

The SCToolbox package SCTB for the polygon() and isinpoly() functions, to build a polygon and check whether points fall within it.
3 Laplace problems with Dirichlet boundary condition
3.1 Analytic solution in the Lshaped domain interior
Field solution with uniform discretization step of 0.01 at the boundary, = 17 in the smooth part and 186 poles, of which 66 fall in the exterior region. = 6.47e7 close to that returned by the LLS when called with default parameters, even though this can not be considered a true upper bound for the error.
A comparison throughout the domain with the values of computed by the LLS, called with a tolerance for maximal absolute error of 1e13, and therefore considered ”exact”, shows clearly on the left a possible drawback of picking evenly sampled points along the boundary. In this case an excellent pointwise boundary approximation does not override the high local error at the reentrant corner, even though in this very point is very low. As a testbed Tr18 consider the computation of , the exact value being 1.0267919261073… We get = 1.025 (0.16%). However, consider that, as reported in GT19c
, experts using specialized FEM codes gave solutions to 24 digits of accuracy, and it is still much better than values obtained from many commercial and open source FEM packages, usually affected by errors between 1% and 10% if no particular attention is paid to mesh generation near critical points.
On the right, the solution after adding 50 clustered points per side for each corner (i.e. 100 points/corner) to the evenly spaced set. This time we get = 1.0267918, and
= 4.82e5 represents a reliable upper bound. Additional points are obtained from interpolation of the given samples very easily, by means of calling
logspace(6,0,50) to determine their displacement from corners along the boundary. In other words they start from a distance of 1e6 from corners and proceed logarithmically spaced. The degree of the AAA approximation is limited to 200 to gain in speed with no drawbacks, as in this case clustered points make the difference. Whenever precise interpolations can be computed from given sets of values, this simple technique is recommended to improve final results greatly.3.2 Analytic solution in the bladeshaped domain interior,
Field solution with = 17 in the smooth part and 167 poles, of which 34 fall in the exterior region. The smooth contour is represented by a polygonal union of tiny segments, getting smaller at the blade ends. = 4.10e7. The listing is the same as in the previous case except in the boundary definition.
3.3 Analytic solution in the Lshaped domain exterior
Field solution with = 17 in the smooth part, meaning = 17 in the first power sum of (1), and 205 poles, of which 124 fall in the interior region. = 3.07e5. The listing differs from the first case (solution in the interior region) in the definition of the smooth part, and in the isinpoly() function asking for interior poles.
4 Conformal mapping
4.1 Lshaped domain interior onto/from disk interior
Conformal mapping (see Tr19; Tr19b) with uniform discretization step of 0.01 at the boundary: from disk onto polygon (mapping of concentric circles, left) and from polygon onto disk (mapping of square grids of different size, centre and right); this = 17 in the smooth part and 192 poles, of which 71 fall in the exterior region. = 5.62e8. Even though sample points are not exponentially clustered near corners, an acceptable accuracy is maintained: images of grid lines at a distance of 0.01 from the Lshaped boundary match those computed with the SCToolbox within some 1e3.
4.2 Lshaped domain exterior onto/from disk interior
Conformal mapping with = 17 in the smooth part and 204 poles, of which 124 fall in the interior region. = 7.61e7.
4.3 Doubly connected domain onto/from annulus
Conformal mapping (see Tr19) with uniform discretization step of 0.01 at both Jordan curves, from annulus of modulus 0.314 onto polygon with hole (left), and from polygon onto annulus (the image lines of a square grid, right). = 17 in the smooth part and 246 poles, of which 108 fall outside the polygon, i.e. outside the outer square or inside the inner one. = 4.78e5.
5 Discussion
The method described is conceived to be generalpurpose, and clearly can not compete in speed nor accuracy with the specialized lightning Laplace solver, even so all the problems addressed, here and in other unreported tests, are solved in a few seconds on a reasonably modern laptop PC with good accuracy, nearly all of the time being consumed by the AAA approximation process. We emphasize again the fact that none of the solutions presented makes use of clustered samples in the neighbourhood of corners: while this clearly remains a possibility for the user, the AAA algorithm has always proven to be so flexible as to return decent approximations for any reasonably discretized contour. The same happens in cases of functions, not presented in this paper, that have singularities on smooth boundaries: the AAA algorithm clusters poles close to critical points following a similar pattern.
Besides the Laplace solver presented in this paper, and its possible applications to engineering computations, the idea of approximating analytic functions on adjacent regions by means of set of poles exterior to each other could suggest novel approaches in other topics: RiemannHilbert problems is one, where sectionally analytic functions are to be determined from boundary conditions.
Acknowledgements
I am grateful to Lloyd N. Trefethen for his valuable suggestions that helped to improve this paper.
Comments
There are no comments yet.