{'fraction_non_alphanumeric': 0.03439373750832778, 'fraction_numerical': 0.010409726848767488, 'mean_word_length': 4.732219570405728, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Solving Laplace problems with the AAA algorithm

We present a novel application of the recently developed AAA algorithm to the solution of Laplace 2D problems; an application to conformal mapping is also shown as a particular case. These classes of problems have also been addressed by means of dedicated software implementations based on rational function approximation, with unprecedented speed and accuracy. Still we show how the AAA algorithm, conceived to be a flexible and domain independent method, is capable of addressing these situations within a unified framework.


Introduction
Recently, problems of the Laplace type and conformal mapping, which can be shown to reduce to a particular form of a Laplace problem [1], have been addressed by the strikingly efficient software implementation based on rational approximations called the "lightning Laplace solver" (LLS) in domain with corners [2], and the "Vandermonde with Arnoldi" algorithm for the stabilization of matrix computation associated with power sums in regions with smooth boundaries [3]. 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 root-exponential convergence [4,2]. 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 Antoulas-Anderson (AAA) algorithm [5], 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.

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 u(z) = ( z) 2 = x 2 on an L-shaped 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. 1e-13, half their number should achieve about 1e-6, 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 u(z) 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 [2]: given a real-valued function u(z) defined on a contour γ in the complex plane, and dividing it into an interior region Γ ⊕ and an exterior one Γ , find a function w(z) analytic in Γ ⊕ for which w(z) ≈ u(z) and w(z) ≈ v(z): often in applications these are referred to as the scalar potential and the stream function respectively, of the complex potential w(z). c is the centre of a power sum describing the "smooth" part of the problem, placed in about the middle of the geometry, whilst {p k } 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 {p k } ∈ Γ ⊕ and N with a negative value in formula (1), c resulting in an additional arbitrary pole placed in the interior region. The magnitude of N is discussed later. From now on, we pursue a different path from [2]: the idea is to build the two partsw(z) andŵ(z) in sequence. First, a point c is chosen which lies about in the middle of the region, and {a n = α n + iβ n } are determined by solving the least-squares (LS) problem where {z j } is the given set of points describing γ, and {u(z j )} the corresponding function values. The real singular partû(z) of the problem is pulled out by simply subtractinĝ Thenû(z) is approximated by the AAA algorithm, unwanted poles are dropped, and a second LS problem is solved in order to determine {b k = δ k + i k }: The reason for there being two constants in (1) is that b 0 adjusts the value of a 0 in the smooth part, in order to approximate {u(z j )} better with the second LS fitting. The parameter N 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 u(z) and its first LS fittingũ(z) along the boundary would be a function severely affected by the Gibbs phenomenon, not restricted to corner neighbourhoods. Figure 2 illustrates such situation, whereû(z) can not be approximated correctly by the AAA algorithm. Therefore, N 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.  (1), far from being beneficial, typically causesû(z) to behave as shown: widely oscillating and "noisy". In such situations the AAA algorithm is unable to return sensible approximations, even if a very large degree is allowed.
We are now ready to give an overview of the algorithm. We exploit the recently announced Arnoldi stabilization for Vandermonde matrices [3], which provides superior stability to the ill-conditioned computation of the smooth part in the first step: 1. Approximate the smooth partũ(z) by means of a power sum, analytic in the desired region; 2. Pull the singular partû(z) out of u(z), and approximate it by means of the AAA algorithm; 3. Use the poles returned by the AAA approximation in the other region to solve the second LS problem forû(z); 4. Build a function handle for w(z), and estimate the maximum error |u − w| ∞ on the boundary.
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 cut-and-paste 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 1e-13. 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 |u − w| ∞ 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 [6], carrying the most recent version of the AAA algorithm; The polyfitA.m and polyvalA.m functions described in [3], for stable computations of power sums; these must be found in the environment path; The SC-Toolbox package [7] for the polygon() and isinpoly() functions, to build a polygon and check whether points fall within it. A comparison throughout the domain with the values of u(z) computed by the LLS, called with a tolerance for maximal absolute error of 1e-13, 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 |u − w| ∞ is very low. As a testbed [8] consider the computation of u(0.99 + 0.99 i), the exact value being 1.0267919261073... We get w(0.99 + 0.99 i) = 1.025 (-0.16%). However, consider that, as reported in [9], experts using specialized FEM codes gave solutions to 2-4 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 w(0.99+0.99 i) = 1.0267918, and |u− w| ∞ = 4.82e-5 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 1e-6 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.

Discussion
The method described is conceived to be general-purpose, 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: Riemann-Hilbert problems is one, where sectionally analytic functions are to be determined from boundary conditions.