Solving Laplace problems with the AAA algorithm

01/26/2020 ∙ by Stefano Costa, et al. ∙ IEEE 0

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.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 5

page 6

page 7

page 8

page 9

page 11

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 root-exponential 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 Antoulas-Anderson (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 L-shaped boundary discretized by evenly sampled points.

Figure 1: Poles returned by the AAA approximation of on the L-shaped boundary (left) and closeup of the reentrant corner (right); note the pole clustering notwithstanding the even boundary sampling.

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 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 real-valued 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 least-squares (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.

Figure 2: Choosing a large value for in (1), far from being beneficial, typically causes 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 Tr19b, which provides superior stability to the ill-conditioned computation of the smooth part in the first step:

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

  2. Pull the singular part out of , 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 ;

  4. Build a function handle for

    , and estimate the maximum error

    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 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 SC-Toolbox 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 L-shaped domain interior

1stp = 0.01;                                             % boundary step discretization
2c = 0.5+0.5*i;                                          % power sum centre
3
4z = [0:stp:2]; z = [z 2+i*[0:stp:1]];                   % L-shaped boundary definition
5z = [z [2:-stp:1]+i]; z = [z 1+i*[1:stp:2]];
6z = [z [1:-stp:0]+2*i]; z = [z i*[2:-stp:0]].’;
7p = polygon([0 2 2+i 1+i 1+2*i 2*i]);
8
9u = real(z).^2;                                         % Dirichlet boundary condition
10
11N = 10+ceil(log(length(z)));                            % smooth part, nr. of terms
12[a,H] = polyfitA(z-c,u,N);                              % Vandermonde with Arnoldi
13wz = @(x) polyvalA(a,H,x(:)-c);                         % function handle, smooth part
14kz = imag(wz(c));
15
16up = u-real(wz(z));                                     % singular part of b.c.
17[r,pol] = aaa(up,z,’mmax’,1000,’lawson’,0);             % AAA approximation
18m = find(isinpoly(pol,p,1E-16)==0);                     % find poles outside
19polm = pol(m);
20B = 1./(z-polm.’);                                      % coefficient matrix
21B(:,end+1) = 1;
22B = [real(B) -imag(B)];
23b = reshape(B\up,[],2);                                 % LS approximation
24b = b(:,1)+i*b(:,2);
25wp = @(x) [1./(x(:)-polm.’) ones(size(x(:)))]*b;        % function handle, singular part
26kp = imag(wp(c));
27
28w = @(x) wp(x)+wz(x(:))-i*(kz+kp);                      % function handle, solution
29maxerr = norm(u-real(w(z)),’inf’)                       % maximum error at the boundary

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.47e-7 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 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 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 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 = 1.0267918, and

= 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.

3.2 Analytic solution in the blade-shaped domain interior,

1h = linspace(0,2*pi,500)’;                              % discretization, angle theta
2c = 0;                                                  % power sum centre
3z = 2*cos(h)+i*(sin(h)+2*cos(h).^3);                    % blade-shaped boundary definition
4p = polygon(z);
5
6u = real(z).^2;                                         % Dirichlet boundary condition
7
8N = 10+ceil(log(length(z)));                            % smooth part, nr. of terms
9[a,H] = polyfitA(z-c,u,N);                              % Vandermonde with Arnoldi
10wz = @(x) polyvalA(a,H,x(:)-c);                         % function handle, smooth part
11kz = imag(wz(c));
12
13up = u-real(wz(z));                                     % singular part of b.c.
14[r,pol] = aaa(up,z,’mmax’,1000,’lawson’,0);             % AAA approximation
15m = find(isinpoly(pol,p,1E-16)==0);                     % find poles outside
16polm = pol(m);
17B = 1./(z-polm.’);                                      % coefficient matrix
18B(:,end+1) = 1;
19B = [real(B) -imag(B)];
20b = reshape(B\up,[],2);                                 % LS approximation
21b = b(:,1)+i*b(:,2);
22wp = @(x) [1./(x(:)-polm.’) ones(size(x(:)))]*b;        % function handle, singular part
23kp = imag(wp(c));
24
25w = @(x) wp(x)+wz(x(:))-i*(kz+kp);                      % function handle, solution
26maxerr = norm(u-real(w(z)),’inf’)                       % maximum error at the boundary

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.10e-7. The listing is the same as in the previous case except in the boundary definition.

3.3 Analytic solution in the L-shaped domain exterior

1stp = 0.01;                                             % boundary step discretization
2c = 0.5+0.5*i;                                          % power sum centre
3
4z = [0:stp:2]; z = [z 2+i*[0:stp:1]];                   % L-shaped boundary definition
5z = [z [2:-stp:1]+i]; z = [z 1+i*[1:stp:2]];
6z = [z [1:-stp:0]+2*i]; z = [z i*[2:-stp:0]].’;
7p = polygon([0 2 2+i 1+i 1+2*i 2*i]);
8
9u = real(z).^2;                                         % Dirichlet boundary condition
10
11N = 10+ceil(log(length(z)));                            % smooth part, nr. of terms
12[a,H]=polyfitA(1./(z-c),u,N);                           % Vandermonde with Arnoldi
13wz = @(x) polyvalA(a,H,1./(x(:)-c));                    % function handle, smooth part
14kz = imag(wz(inf));
15
16up = u-real(wz(z));                                     % singular part of b.c.
17[r,pol] = aaa(up,z,’mmax’,1000,’lawson’,0);             % AAA approximation
18m = find(isinpoly(pol,p,1E-16)==1);                     % find poles inside
19polm = pol(m);
20B = 1./(z-polm.’);                                      % coefficient matrix
21B(:,end+1) = 1;
22B = [real(B) -imag(B)];
23b = reshape(B\up,[],2);                                 % LS approximation
24b = b(:,1)+i*b(:,2);
25wp = @(x) [1./(x(:)-polm.’) ones(size(x(:)))]*b;        % function handle, singular part
26kp = imag(wp(inf));
27
28w = @(x) wp(x)+wz(x(:))-i*(kz+kp);                      % function handle, solution
29maxerr = norm(u-real(w(z)),’inf’)                       % maximum error at the boundary

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.07e-5. 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 L-shaped domain interior onto/from disk interior

1stp = 0.01;                                             % boundary step discretization
2c = 0.5+0.5*i;                                          % power sum centre
3
4z = [0:stp:2]; z = [z 2+i*[0:stp:1]];                   % L-shaped boundary definition
5z = [z [2:-stp:1]+i]; z = [z 1+i*[1:stp:2]];
6z = [z [1:-stp:0]+2*i]; z = [z i*[2:-stp:0]].’;
7p = polygon([0 2 2+i 1+i 1+2*i 2*i]);
8
9u = -log(abs(z-c));                                     % Dirichlet boundary condition
10
11N = 10+ceil(log(length(z)));                            % smooth part, nr. of terms
12[a,H] = polyfitA(z-c,u,N);                              % Vandermonde with Arnoldi
13wz = polyvalA(a,H,z-c);                                 % smooth part
14
15up = u-real(wz);                                        % singular part of b.c.
16[r,pol] = aaa(up,z,’mmax’,1000,’lawson’,0);             % AAA approximation
17m = find(isinpoly(pol,p,1E-16)==0);                     % find poles outside
18B = 1./(z-pol(m).’);                                    % coefficient matrix
19B(:,end+1) = 1;
20V = [real(B) -imag(B)];
21b = reshape(V\up,[],2);                                 % LS approximation
22b = b(:,1)+i*b(:,2);
23wp = B*b;                                               % singular part
24
25w = wp+wz;                                              % solution
26maxerr = norm(u-real(w),’inf’)                          % maximum error at the boundary
27
28g = (z-c).*exp(w);                                      % conformal mapping function
29f = aaa(g,z);                                           % polygon onto disk, AAA approx.
30finv = aaa(z,g);                                        % disk onto polygon, AAA approx.

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.62e-8. 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 L-shaped boundary match those computed with the SC-Toolbox within some 1e-3.

4.2 L-shaped domain exterior onto/from disk interior

1stp = 0.01;                                             % boundary step discretization
2c = 0.5+0.5*i;                                          % power sum centre
3
4z = [0:stp:2]; z = [z 2+i*[0:stp:1]];                   % L-shaped boundary definition
5z = [z [2:-stp:1]+i]; z = [z 1+i*[1:stp:2]];
6z = [z [1:-stp:0]+2*i]; z = [z i*[2:-stp:0]].’;
7p = polygon([0 2 2+i 1+i 1+2*i 2*i]);
8
9u = log(abs(z-c));                                      % Dirichlet boundary condition
10
11N = 10+ceil(log(length(z)));                            % smooth part, nr. of terms
12[a,H] = polyfitA(1./(z-c),u,N);                         % Vandermonde with Arnoldi
13wz = polyvalA(a,H,1./(z-c));                            % smooth part
14
15up = u-real(wz);                                        % singular part of b.c.
16[r,pol] = aaa(up,z,’mmax’,1000,’lawson’,0);             % AAA approximation
17m = find(isinpoly(pol,p,1E-16)==1);                     % find poles inside
18B = 1./(z-pol(m).’);                                    % coefficient matrix
19B(:,end+1) = 1;
20V = [real(B) -imag(B)];
21b = reshape(V\up,[],2);                                 % LS approximation
22b = b(:,1)+i*b(:,2);
23wp = B*b;                                               % singular part
24
25w = wp+wz;                                              % solution
26maxerr = norm(u-real(w),’inf’)                          % maximum error at the boundary
27
28g = (1./(z-c)).*exp(w);                                 % conformal mapping function
29f = aaa(g,z);                                           % polygon onto disk, AAA approx.
30finv = aaa(z,g);                                        % disk onto polygon, AAA approx.

Conformal mapping with = 17 in the smooth part and 204 poles, of which 124 fall in the interior region. = 7.61e-7.

4.3 Doubly connected domain onto/from annulus

1stp = 0.01;                                             % boundary step discretization
2c = -0.25-0.25*i;                                       % power sum centre;
3
4s1 = [0:stp:2];                                         % outer Jordan curve
5z1 = [i*sqrt(2)+s1*exp(i*5/4*pi)]; z1 = [z1 -sqrt(2)+s1*exp(i*7/4*pi)];
6z1 = [z1 -i*sqrt(2)+s1*exp(i*1/4*pi)]; z1 = [z1 sqrt(2)+s1*exp(i*3/4*pi)].’;
7p1 = polygon([i*sqrt(2) -sqrt(2) -i*sqrt(2) sqrt(2)]);
8
9s2 = [0:stp:0.5];                                       % inner Jordan curve
10z2 = s2-0.5-0.5*i; z2 = [z2 s2*exp(i*pi/2)+0.0-0.5*i];
11z2 = [z2 s2*exp(i*pi)+0.0+0.0*i]; z2 = [z2 s2*exp(i*3/2*pi)-0.5+0.0*i].’;
12p2 = polygon([-0.5-0.5*i 0.0-0.5*i 0.0+0.0*i -0.5+0.0*i]);
13
14z = [z1; z2];
15u = -log(abs(z-c));                                     % Dirichlet boundary condition
16
17N = 10+ceil(log(length(z)));                            % smooth part, nr. of terms
18[a,H] = polyfitA(z-c,u,N);                              % Vandermonde with Arnoldi
19wz = polyvalA(a,H,z-c);                                 % smooth part
20
21up = u-real(wz);                                        % singular part of b.c.
22[r,pol] = aaa(up,z,’mmax’,1000,’lawson’,0);             % AAA approximation
23m = find((isinpoly(pol,p1,1E-16)==0) | ...              % find poles outside
24(isinpoly(pol,p2,1E-16)==1));
25B = 1./(z-pol(m).’);                                    % coefficient matrix
26B(:,end+1) = 1;
27B(:,end+1) = 0; B(length(z1)+1:end,end) = 1;            % add column for modulus
28V = [real(B) -imag(B)];
29b = reshape(V\up,[],2);                                 % LS approximation
30b = b(:,1)+i*b(:,2);
31modulus = exp(-b(end))
32maxerr = norm(u-real(B*b+wz),’inf’)                     % maximum error at the boundary
33b(end) = 0;
34wp = B*b;                                               % singular part
35
36w = wp+wz;                                              % solution
37
38g = (z-c).*exp(w);                                      % conformal mapping function
39f = aaa(g,z);                                           % polygon onto disk, AAA approx.
40finv = aaa(z,g);                                        % disk onto polygon, AAA approx.

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.78e-5.

5 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.

Acknowledgements

I am grateful to Lloyd N. Trefethen for his valuable suggestions that helped to improve this paper.

References