A Least Squares Approach for Determining the Coefficients of the CVBEM Approximation Function

The Complex Variable Boundary Element Method (CVBEM) is a technique for modeling boundary value problems (BVPs) of the Laplace type [1-3]. As part of the CVBEM modeling process, it is necessary to make several implementation decisions, including the following:

1. A scheme for determining the initial locations of candidate nodes and candidate collocation points (or boundary data, in general).

2. Determination of the basis functions to comprise the approximation function.

3. Selection of an algorithm for choosing which subset of the candidate nodes and boundary data should be used in the formulation of the CVBEM model.

4. Selection of a metric for measuring approximation error.

5. Selection of a technique for determining the coefficients of the CVBEM approximation function (in this case, collocation or least squares).

Several of these topics have recently been reviewed. For example, a scheme for determining the initial locations of candidate nodes was proposed in [4]. Discussion regarding the selection of basis functions for use with the CVBEM can be found in [5,6]. Algorithms for selecting the computational nodes and collocation points to be used in the CVBEM modeling process have been developed in $[7,8]$. The fourth implementation decision refers to the selection of a norm for measuring the error of the CVBEM approximation function. Commonly-used norms for this purpose include the $\ell_1, \ell_2$, and $\ell_{\infty}$; however, in principle, any norm of the form $\ell_p, p \in \mathbb{N}$, can be used.

The subject of this paper pertains to the last of the implementation decisions listed above namely, the selection of a technique for determining the coefficients of the CVBEM approximation function. In particular, this work investigates a collocation approach and, separately, a least squares approach, for determining the coefficients of the CVBEM approximation function.

In most of the recent works related to the CVBEM, collocation has been used to determine the coefficients of the approximation function [4-11]. Since the collocation approach is well-documented in the aforementioned references, the focus of this paper is on developing the least squares approach. Figure 1 illustrates the logic of the CVBEM algorithm in the context of the least squares approach.

The difference between the collocation and least squares approaches has to do with how the boundary data are used when computing the coefficients of the CVBEM approximation function. In particular, in the collocation approach, two collocation points are selected for each term in the CVBEM approximation function. Hence, for $n$ terms in the CVBEM approximation function, $2 n$ collocation points are selected at which to apply the boundary conditions. Thus, in the collocation approach, only $2 n$ of the available boundary data points are used when computing the coefficients of the CVBEM approximation function. This is in contrast to the least squares approach in which all of the available boundary data are utilized in each iteration of the algorithm when computing the coefficients of the CVBEM approximation function.

In the collocation approach, the CVBEM approximation function is guaranteed to satisfy the boundary conditions at each of the collocation points. However, in the least squares approach, there is no guarantee that the CVBEM approximation function satisfies the boundary conditions at any location on the boundary. Rather, the least squares approach minimizes the $\ell_2$ norm of the difference between the values of the CVBEM approximation function and the target boundary conditions along the problem boundary.

The computational outcomes for both the collocation and least squares approaches demonstrate the measured error decreasing in the $\ell_{\infty}$ norm as the number of linearly independent terms used in the CVBEM approximation function increases. These terms each correspond to a source point (i.e. computational node) located in the exterior of the problem domain, and the placement of these source points is determined according to the node-positioning algorithms (NPAs) discussed in [7, 8]. Furthermore, it is noted that in this implementation of the CVBEM, the candidate nodes are located strictly in the exterior of the problem domain. This placement scheme is contrary to the usual boundary element numerical models in which the model nodes are placed on the problem boundary and are also used to describe the problem boundary geometry. In the present CVBEM implementation, the disconnection between the locations of computational nodes and the problem boundary geometry provides a significant increase in degrees of freedom that can be used to further improve the modeling success with regard to satisfying the problem boundary conditions.

The CVBEM approximation function is as follows:

where $\left(c_j=a_j+i b_j\right) \in \mathbb{C}$, and $g_j(\mathrm{z})=\lambda_j(x, y)+i \mu_j(x, y)$ is a complex variable basis function that is analytic in the problem domain. The coefficients of eqn (1) are complex numbers of the form $c_j=a_j+i b_j$, where each $a_j$ and $b_j$ is an unknown value that needs to be determined. These values are determined using either collocation or least squares so as to satisfy the given boundary conditions as closely as possible in a given norm.

In the collocation approach, the $a_j$ and $b_j$ values are determined by solving a square matrix equation of the form $A x=b$, where $A \in \mathbb{R}^{2 n \times 2 n}, x \in \mathbb{R}^{2 n \times 1}$, and $b \in \mathbb{R}^{2 n s 1}$. More specifically,

where the $i^{\text {th }}$ rows of $\Lambda$ and $M$, denoted $\Lambda_{i,:}$ and $M_{i, .,}$ respectively, are as follows:

In eqns (2) and (3), ( $x_i, y_i$ ) denotes the location of the $i^{\text {th }}$ collocation point on the problem boundary, and $\phi\left(x_i, y_i\right)$ denotes the value of the Dirichlet boundary condition at that collocation point. The matrix equation for handling mixed boundary conditions is similar and described in [6].

The advantage of the collocation approach is that the resulting CVBEM approximation function will satisfy the given boundary conditions at each of the collocation points. That is, theoretically, the error of the CVBEM approximation function is 0 at each of the collocation points. However, in practice, the error can technically be nonzero due to numerical considerations such as the effects of truncation and finite precision arithmetic, particularly if the matrix, $A$, in eqn (2) is ill-conditioned.

Figure 1 depicts a typical situation when evaluating the error along the problem boundary of a CVBEM model developed using the collocation approach. Importantly, at the collocation points, depicted as black dots, the absolute value of the error function is numerically 0 .

Using the least squares approach, the CVBEM approximation function has the same form as in eqn (1). In particular, there are still $2 n$ coefficients of the CVBEM approximation function to determine. The primary difference between the collocation approach and the least squares approach has to do with how these coefficients are calculated. In the least squares approach, the $2 n$ coefficient values are determined by solving an over-determined matrix equation of the form $A x=b$, where $A \in \mathbb{R}^{N_{\mathrm{B}} \times 2 \mathrm{n}}, x \in \mathbb{R}^{2 n \times 1}$, and $b \in \mathbb{R}^{N_{\mathrm{B}} \times 1}$, where $N_{\mathrm{B}}>2 n$ denotes the number of known boundary condition data points. Using the least squares approach, the objective is to obtain the coefficient vector, $x^*$, that minimizes the following:

In particular, the matrix equation is as follows:

where the $i^{\text {th }}$ rows of $\Lambda$ and $M$, denoted $\Lambda_{i,:}$ and $M_{i,:}$, respectively, are as follows:

In eqns (5) and (6), ( $x_i, y_i$ ) denotes the location of the $i^{\text {th }}$ boundary condition data point and $\phi\left(x_i, y_i\right)$ denotes the value of the Dirichlet boundary condition at that location. The matrix equation for handling mixed boundary conditions is similar and may be adapted from the matrix equation described in [6].

The advantage of the least squares approach is that it incorporates all of the boundary data when determining the coefficients of the CVBEM approximation function. The least squares approach would be beneficial in a situation in which using only $2 n$ collocation points would be insufficient with regard to describing possibly extreme variations in the target potential function on the problem boundary.

Another aspect of the least squares approach is regularization, such as Tikhonov regularization [12, 13]. Regularization is used to address potential computational issues such as overfitting to the given boundary conditions. When Tikhonov regularization is applied, the objective is to find the coefficient vector, $x_{\text {REG}}^{\star}$, minimizing the following:

In eqn (7), $a \in \mathbb{R}$ is known as the regularization parameter. The difference between eqn (4) and eqn (7) is that eqn (7) employs a penalty on the norm of the coefficient vector, $x$, such that solution vectors with large norms are not preferred. Additionally, Tikhonov regularization can improve the condition of the matrix $A$, which may lead to more numerically stable solutions for the coefficients. The improved conditioning is a consequence of pre-pending a scalar multiple of the identity matrix on top of the standard least squares matrix, as indicated in eqn (8). That is, the lower $N_{\mathrm{B}}$ rows of the matrix equation in eqn (8) are exactly the same as in eqn (5), but the first $2 n$ rows are a scalar multiple of the identity matrix. Since Tikhonov regularization can be used to improve the condition number of a matrix, it may be an appropriate tool to remedy ill-conditioned CVBEM models of BVPs. When Tikhonov regularization is employed, the matrix equation is as follows:

where the $i^{\text {th }}$ rows of $\Lambda$ and $M$ are as described in eqn (6).

Let $\Omega \subset \mathbb{C}$ denote a simply-connected problem domain. The basis functions used to develop the CVBEM models in this work are of the form $\left(z-z_j\right) \ln {\alphaj}\left(z-z_j\right)$, and the resulting CVBEM approximation function is as follows:

where $z_j \in \mathbb{C} \backslash \Omega$, and $a_j \in \mathbb{R}$ denotes the angle of the rotated branch cut of the complex logarithm function. Typically, the branch cuts are rotated radially away from a point in the interior of the problem domain, such as the centroid of the problem domain.

The basis functions in eqn (9) are derived from the Cauchy integral equation, as described in [ $2,3,14$ ]. The points $z_j$ in eqn (9) are singularities of the basis functions and are the locations of the computational nodes that are evaluated by the NPA. The points $z_j$ are referred to as computational nodes because they do not have a physical meaning in the context of the CVBEM solution within the problem domain. Instead, they are consequences of the particular choice of basis functions, which have singularities. However, the existence of these singularities is computationally important because determining suitable locations for them introduces new degrees of freedom that can be optimized using a NPA, such as the ones described in $[7,8]$.

The NPAs used in this work are referred to as NPA1 (developed in [7]) and NPA2 (developed in [8]). These algorithms are heuristics for determining subsets of the candidate nodes to use in the CVBEM model. In both NPAs, the nodes are selected one-at-a-time. Specifically, the CVBEM model begins with 0 nodes. To select the first node for incorporation in the CVBEM model, each candidate node is assessed as a 1 -node CVBEM model. Then, the node that can be used to satisfy the given boundary conditions most accurately is added to the CVBEM model. To determine the second node, all of the candidate nodes are assessed as 2 -node models when added to the already-selected first node. The new node that results in a 2 -node model that can be used to satisfy the given boundary conditions most accurately is added to the CVBEM model. This process repeats until $n$ nodes have been selected. In general, when selecting the $k^{\text {th }}$ node, the previous ( $k-1$ )-selected nodes are considered fixed and all of the remaining candidate nodes are assessed with regard to their performance in a $k$-node CVBEM model when added to the already-selected (k-1) nodes.

The difference between NPA1 and NPA2 is that in NPA1, once a node is selected, it is permanently added to the CVBEM model. In NPA2, a nodal position refinement procedure is implemented, which allows for previously-selected nodes to be re-evaluated based on the subsequent selection of nodes to determine if a different node could be used in the model to more accurately satisfy the given boundary conditions. More details pertaining to the nodal position refinement procedure may be found in [8].

In the following demonstration problems, note the difference between the candidate nodes and boundary data, as depicted in Figs. 2 and 9 for problems 1 and 2, respectively. In particular, the black dots are the candidate nodes, whose locations correspond to the singularities of the basis functions defined in eqn (9), which are optimized using the NPA. The blue dots are the boundary data at which the boundary conditions are known. In the least squares approach, all of the boundary data are used when computing the coefficients of the CVBEM approximation function. In the collocation approach, only a subset of $2 n$ boundary data points, whose locations are also determined by the NPA, are used when computing the coefficients of the CVBEM approximation function.

A convenience of modeling with the CVBEM is that complex variable theory provides a simple technique for determining the maximum error of the CVBEM approximation function in the $\ell_{\infty}$ norm. Let $\phi \in \mathbb{R}^2$ denote a harmonic function that is the target solution of the BVP of interest, and let $\hat{\phi}=\Re[\hat{\omega}(z)] \in \mathbb{R}^2$ denote the CVBEM approximation of $\phi$. Error estimation is based on the following two key observations:

1. The target solution, $\phi$, satisfies Laplace's equation and is, therefore, harmonic in $\Omega$.

2. The CVBEM approximation function, denoted $\hat{\omega}$, is an analytic complex variable function. Consequently, $\hat{\phi}=\Re[\hat{\omega}]$, which is the CVBEM approximation of the target function $\phi$, is also harmonic in $\Omega$.

The error function used in this work is $|\varepsilon(x, y)|$, where $\varepsilon(x, y)$ is defined as follows:

Since $\varepsilon(x, y)$ is the difference between two harmonic functions, $\varepsilon(x, y)$ is itself harmonic in $\Omega$. Consequently, by the maximum modulus principle for harmonic functions, $|\varepsilon(x, y)|$ attains its maximum on $\partial \Omega$. Hence, the maximum error of the CVBEM approximation function occurs on $\partial \Omega$:

A reasonable approximation of $\max {(x, y) \in \partial \Omega}|\varepsilon(x, y)|$ can be obtained by computing the value of $|\varepsilon(x, y)|$ at many locations along the problem boundary. Generally, as $|\varepsilon(x, y)|$ is computed at more locations along the boundary, the estimation of $\max {\partial \Omega}|\varepsilon(x, y)|$ will improve provided that the error evaluation points are reasonably spaced.

The following demonstration problems contain stagnation points as features of their solutions. The stagnation points are among the most difficult areas of the target flow situations to model because of the relatively extreme curvature of the solutions in those areas. Given the computational difficulty of modeling these areas, it is of interest to compare the computational outcomes obtained using the collocation and least squares approaches specifically in these areas. Consequently, particular focus is given to the flownets developed by the CVBEM models at the stagnation points in Fig. 3c for problem 1 and Figs. 11 and 12 for problem 2.

The two methodologies are applied to modeling potential flow over a half-cylindrical obstacle. This demonstration problem leverages the difficult flow characteristics observed at the upstream and downstream stagnation points of the obstacle where, in the limiting situation, streamlines adjacent to the solid surface form a right angle. Additionally, the 'north pole' of the cylinder presents a difficult flow situation to model computationally because of the relatively extreme curvature of the solution in that area.

The analytic representation of the velocity potential for this problem is given in [15] as

For this demonstration, it is possible to take advantage of the symmetry of the target flow situation. In particular, the target potential flow is symmetric about the real axis. Therefore, the presented results are of the CVBEM approximation in the upper half-plane because the solution in the lower half-plane can be obtained by reflection. Of note, the target potential flow is also symmetric about the imaginary axis. However, this symmetry is not taken advantage of in this work due to the desire to model both the upstream and downstream stagnation points because of the computational difficulty of obtaining an accurate approximation of the potential flow in their vicinity.

A formal description of the test problem follows in Table 1. Since the exact solution is analytic everywhere except at $z=0$, the real and imaginary parts of $\omega$ are harmonic functions in $\mathbb{C} \backslash\{0\}$ and thus harmonic throughout $\Omega \subset \mathbb{C} \backslash\{0\}$.

Figure 3 depicts the locations of the candidate nodes and 5\% of the boundary data used in the formulation of the CVBEM models in this work. Recall, the least squares approach uses all of the available boundary data during the computation of the coefficients for the CVBEM approximation function, while in the collocation approach, only a subset of $2 n$ points, as determined by the NPA, are used to compute the coefficients.

Figure 4 depicts various flownet outcomes from the CVBEM model. In these figures, particular emphasis is given to the north pole of the half-cylinder, as well as to the stagnation points because these are the regions of the flow situation in which the curvature of the target solution is most extreme. Figure 3b depicts the left side of the obstacle, and Fig. 3c depicts a magnified version of the stagnation point at the left edge of the obstacle, $(-1,0)$. The stagnation points at $(-1,0)$ and $(-1,0)$ are difficult to model due to the extreme curvature of the target potential function in these areas. For this reason, it is of interest to ensure the CVBEM models yield high-fidelity approximations of these areas.

Figure 5 depicts the locations of the NPA2-selected nodes for the least-squares-based CVBEM model of the potential flow problem described in Table 1. Figure 5 depicts the locations of the NPA2-selected collocation points and NPA2-selected nodes for the collocation-based CVBEM model of the potential flow problem described in Table 1.

Using NPA2, 428 candidate nodes, 2,000 boundary data points, and $n=40$ terms in the CVBEM approximation function, the collocation and least squares approaches yield CVBEM models with maximum errors on the order of $10^{-7}$ and $10^{-8}$, respectively, as indicated in Table 2 . Since these models are highly accurate, superimposing their respective flownet outcomes would not help to compare the outcomes because they would be visually indistinguishable. Consequently, this work examines the accuracy of the CVBEM models with regard to approximating five specific potential lines of interest located near the left stagnation point, ( $-1,0$ ). The five potential lines of interest are depicted in Fig. 6. Each of the five potential lines of interest are defined by the value of their level curve. For this problem, potential lines can be defined for values $c \in[-10 / 3,10 / 3]$ as the set of points $(x, y)$ satisfying:

The five potential lines that are considered in this demonstration correspond specifically to the values $c \in\{-2.000000,-1.999999,-1.999996,-1.999991,-1.999984\}$.

Figure 8 shows the computational error for the collocation and least squares outcomes when approximating these five potential lines of interest. The collocation outcomes for each of the potential line approximations demonstrate the absolute error decreasing below $10^{-10.5}$. This occurs because of the presence of collocation points near the stagnation point, which locally shrink the computational error to 0 . Meanwhile, the least squares outcomes show that the error for each of the potential line approximations is consistently on the order of $10^{-10.5}$.

Figure 8 highlights one of the primary differences between the two approaches. Namely, in the collocation approach, the approximation error is numerically 0 at each of the collocation points (a subset of only $2 n$ points from the set of available boundary data points). On the other hand, the least squares approach seeks to minimize the approximation error at all of the boundary data points, which results in the approximation error usually never being exactly 0 anywhere, but tending to be small everywhere.

Although the computational results, given in Table 2, demonstrate consistent error reduction with increasing computational effort and model complexity (i.e. the number of linearly-independent terms in the CVBEM approximation function), the problem geometry itself is anticipated to have a direct influence on the level of accuracy that is achievable for the examined CVBEM models. For example, the prior paper [5] showed that simpler geometries are associated with further error reduction of the CVBEM model with comparable computational effort. Furthermore, as demonstrated in the second example problem of this paper, which has a more complicated geometry and more complicated boundary conditions, it is necessary to use considerably more terms in the CVBEM approximation function to achieve error levels comparable to the ones achieved while modeling this problem.

Figure 9 illustrates how the maximum error of each of the CVBEM models examined tended to decrease as the number of terms used in the CVBEM approximation function increased. This figure depicts the maximum error results as each new node is added to the CVBEM model using the collocation approach as well as the least squares approach. Results using NPAs 1 and 2 are shown for both of the approaches.

In this problem, two methodologies are applied to modeling potential flow in a corner and over two successive half-cylindrical obstacles. This demonstration problem has a more complicated geometry and more complicated boundary conditions than the previous example problem.

As discussed in the previous problem, the greatest computational difficulty when modeling these example problems has to do with modeling the curvature of the flow situation in the stagnation points. The previous example illustrated the CVBEM solution to a potential flow problem with two stagnation points. This example will illustrate the CVBEM solution to a potential flow problem with five stagnation points. Given the increased difficulty of this example problem, as indicated by the increased number of stagnation points to model, obtaining highly accurate CVBEM models will require using CVBEM models with more linearly independent terms in the approximation function. A formal description of the test problem follows in Table 3.

Figure 10 depicts the locations of the candidate nodes and 2.5% of the boundary data used in the formulation of the CVBEM models in this work. Recall, the least squares approach uses all of the available boundary data during the computation of the coefficients for the CVBEM approximation function, while in the collocation approach, only a subset of $2 n$ points, as determined by the NPA, are used to compute the coefficients. Figures 11-13 depict the flownet outcomes from modeling this problem with the CVBEM using the least squares approach. Special emphasis is given to the flownets developed in the five stagnation points of this flow situation.

Figure 14 depicts the locations of the NPA2-selected nodes for the least-squares-based CVBEM model of the potential flow problem described in Table 3. Figure 15 depicts the locations of the NPA2-selected collocation points and NPA2-selected nodes for the colloca-tion-based CVBEM model of the potential flow problem described in Table 3.

Figure 12 shows the flownets obtained from the CVBEM models in the vicinity of each of the four stagnation points associated with the two half-cylindrical obstacles. Specifically, there are upstream and downstream stagnation points at the edge of each obstacle where the obstacle intersects the horizontal axis. Special attention is paid to these points because of the relatively extreme curvature of the target potential function in these areas. This curvature makes the target potential function difficult to model with high accuracy.

Figure 16 illustrates how the maximum error of each of the CVBEM models examined tended to decrease as the number of terms used in the CVBEM approximation function increased. This figure depicts the maximum error results as each new node is added to the CVBEM model using the collocation approach as well as the least squares approach. Results using NPAs 1 and 2 are shown for both of the approaches. Interestingly, for this example problem, the use of NPA2 did not seem to result in much of an improvement over NPA1 for either the collocation or least squares approaches. This may be due to the general difficulty of modeling this problem given that the target flow situation incorporates five stagnation points.

Historically, the collocation approach has been the predominant technique for determining the coefficients of a CVBEM approximation function. This paper treats the technique for determining the CVBEM coefficients as an opportunity for customization within the CVBEM methodology and demonstrates the use of an alternative method - namely, the least squares approach. These two approaches are demonstrated and compared when applied to two benchmark problems in ideal fluid flow.

A primary advantage of the collocation approach is that it guarantees the resulting CVBEM approximation function will satisfy the given boundary conditions at no less than $2 n$ locations on the problem boundary, where $n$ is the number of linearly independent terms used in the CVBEM approximation function. On the other hand, a benefit of the least squares approach is that it incorporates all of the available boundary data when determining the coefficients of the CVBEM approximation function. This feature of the least squares approach should be considered for applications of the CVBEM in higher spatial dimensions in which there may be a lot of boundary data available, especially if the surface area of the problem domain is large. The collocation approach could suffer in this setting if the target potential function cannot be adequately described by the values of the boundary conditions at just $2 n$ locations.

The collocation approach has been used as the standard approach in CVBEM implementations since its initial description in [2]. However, this work demonstrates that the least squares approach can also be used to obtain highly-accurate CVBEM models of the benchmark problems of potential flow over a half-cylindrical obstacle (a Dirichlet BVP) and potential flow in a corner followed by two successive half-cylindrical obstacles (a mixed BVP). The success of the least squares approach in modeling these problems suggests the approach may be viable in other computational settings and should be considered during the formulation of CVBEM models.