On the numerical solution of minimal energy problems

We study the constructive and numerical solution of minimizing the energy for the Gauss variational problem involving the Newtonian potential. As a special case, we also treat the corresponding condenser problem. These problems are considered for two two-dimensional compact, disjoint Lipschitz manifolds Γ j ⊂ ℝ3, j = 1, 2, charged with measures of opposite sign. Since this minimizing problem over an affine cone of Borel measures with finite Newtonian energy can also be formulated as the minimum problem over an affine cone of surface distributions belonging to the Sobolev–Slobodetski space , which allows the application of simple layer boundary integral operators on Γ, a penalty approximation for the Gauss variational problem can be used. The numerical approximation is based on a Galerkin–Bubnov discretization with piecewise constant boundary elements. To the discretized problem, the projection-iteration is applied where the matrix times vector operations are executed with the fast multipole method. For the condenser problem, we solve the dual problem which reduces in our case to solving two linear boundary integral equations. Here the fast multipole method provides an efficient solution algorithm. We finally present some convergence studies and error estimates.


Introduction
Carl Friedrich Gauss investigated in [1] the variational problem of minimizing the energy of the Newtonian potential over nonnegative charges on the boundary surface of a given domain. The sign condition was given up in connection with boundary integral equation methods where distributional boundary charges have been introduced for solving boundary value problems (see e.g. Costabel's article [2] and also [3]).
A different generalization of the original Gauss variational problem maintaining the sign restriction but employing measures as charges has independently grown into an eminent branch of modern potential theory (see e.g. [4] and the extensive works [5][6][7][8][9][10][11] and for two dimensions [12]). In this article, we consider this so-called generalized Gauss variational problem (from now on just Gauss problem) on two two-dimensional nonintersecting compact Lipschitz surfaces À 1 , À 2 & R 3 which are loaded by measures ¼ 1 À 2 , j 2 M þ (À j ), j ¼ 1, 2, and try to determine the corresponding equilibrium state also numerically by appropriate boundary element approximation. To this end, we first show that every Borel measure with finite Newtonian energy defines uniquely a distribution on À ¼ À 1 [ À 2 belonging to the Sobolev-Sobodetski space H À 1 2 ðÀÞ such that the linear functional defined by the measure and the one defined by the L 2 -duality with coincide on C 1 (À). Moreover, the minimizing Gauss problem over the measures is equivalent to minimizing the energy over a corresponding affine cone in H À 1 2 ðÀÞ which can be expressed in terms of the simple layer boundary integral operator on À. This allows one to approximate the Gauss problem by employing a penalty formulation. The latter corresponds to a nonlinear variational problem on the whole space H À 1 2 ðÀÞ whose solutions provide error estimates depending on the penalizing parameter.
For the penalized variational problem, the use of piecewise constant boundary elements on triangulations on À and corresponding Galerkin-Bubnov discretization results in a convex, finite-dimensional minimization problem with a symmetric, positive definite system matrix. The latter can be solved by the gradient-projection method applied to the discrete system [13,Chapter I]. For the corresponding matrix times vector multiplications, the fast multipole method [14] is used. The convergence of the gradient-projection method depends on the penalty parameter as well as on the boundary element's mesh width and becomes very slow for large penalty and small mesh width. In particular, the charges converge faster for smaller penalty whereas the total charge values are better approximated for larger penalty. Therefore we are using a cascading approach.
In [10], the dual problem to the Gauss problem is defined and it is shown that its solution coincides with the one to the original Gauss problem if the latter exists. If the exterior source of energy is zero then the Gauss problem becomes a condenser problem for which the dual formulation is easier to treat [9,11]. If the condenser consists of just two plates without irregular points as considered here, and for the special weight function g ¼ 1, the dual problem significantly simplifies since then the potential values of the minimizing solution are constant on À j , j ¼ 1, 2 and sum up to 1. In this case, the dual problem can be reduced to the solution of two linear boundary integral equations whose solution provides the desired equilibrium state. The latter can be approximated also by piecewise constant charges, and the fast multipole boundary element method provides here an efficient solution algorithm.
In applications, the numerical solution of the Gauss as well as of the condenser problem is of great interest if for practical reasons in electrical engineering on À 1 only nonnegative and on À 2 only nonpositive charges are allowed (see 'capacitors' in [15]). The numerical solution of the condenser problem also has applications in approximation theory and the development of efficient numerical integration [16].
We present numerical experiments for the Gauss as well as for the condenser problem. In particular, we show that the cascading approach for the penalty method can be used for the Gauss variational problem. For the condenser problem, two types of examples, (a) and (b), are presented where for À 2 surfaces are used which approximate an infinitely long surface. The corresponding numerical results are in agreement with theoretical results in [5,8], where under the geometry of example (b) it was shown that for the infinitely long surface an equilibrium minimizing measure does not exist, in contrast to example (a).

The generalized Gauss variational problem
We are going to consider the Gauss problem in R 3 , i.e. the problem of minimizing the Newtonian energy of signed Borel measures in the presence of an external field f. The corresponding class of admissible measures (or charges) is supposed to be associated with a condenser, which is treated here as an ordered pair A ¼ (À 1 , À 2 ) of two-dimensional compact, Lipschitz manifolds À 1 , À 2 & R 3 such that To formulate the problem, we need the following notation.
Let M ¼ MðR 3 Þ be the -algebra of Borel measures on R 3 , equipped with the vague topology, i.e. the topology of pointwise convergence on the class C 0 (R 3 ) of all real-valued continuous functions on R 3 with compact support (see e.g. [17,18]). The mutual Newtonian energy of 1 , 2 2 M is given by the formula (certainly, if the integral on the right is well defined -as a finite number or AE1). For 1 ¼ 2 we get the energy I( 2 ) :¼ I( 2 , 2 ) of 2 , and for 1 ¼ x -the value of the potential of 2 at x 2 R 3 , denoted by U 2 (x); namely, dð yÞ jx À yj , where 2 M and x 2 R 3 : consist of all 2 MðR 3 Þ with finite energy. Since the Newtonian kernel jx À yj À1 is strictly positive definite (see e.g. [18] and the references given there), the bilinear form I( 1 , 2 ) defines on E a scalar product. Hence, E can be treated as a pre-Hilbert space with the norm Given a closed set F & R 3 , let MðF Þ consist of all 2 MðR 3 Þ supported by F, and let M þ ðF Þ be the convex cone of all nonnegative 2 MðF Þ. We also write E þ ðF Þ :¼ M þ ðF Þ \ E: The condenser A ¼ (À 1 , À 2 ) is supposed to be loaded by charges The collection of all those will be denoted by E(A); it is a convex cone in E.
Further, let g be a given continuous, positive function on À :¼ À 1 [ À 2 and a ¼ (a 1 , a 2 ) 2 R 2 a given vector with a j 40, j ¼ 1, 2. Then the set of admissible charges for the Gauss problem is given by EðA, a, gÞ :¼ Observe that E(A, a, g) is an affine, convex cone in the pre-Hilbert space E.
In addition, let f denote a given continuous function on À, characterizing an exterior source of energy. Then The Gauss problem now reads: A minimizing measure 0 is unique (if exists). This follows from the strict positive definiteness of the Newtonian kernel and the convexity of the class of admissible measures; see [7,Lemma 6]. But what about the existence of 0 ?
Assume for a moment that at least one of À 1 and À 2 is noncompact. Then it is not clear at all that the equilibrium state in the Gauss variational problem can be attained. Moreover, it has recently been shown by the third author that, in this case, a minimizing measure 0 in general does not exist; necessary and sufficient conditions for 0 to exist were given in [5,6,8].
However, in the case under consideration, where both À 1 and À 2 are assumed to be compact, the Gauss variational problem has a (unique) solution in the cone E(A, a, g). Indeed, this follows from the vague compactness of E(A, a, g) when combined with the fact that the Gauss functional G f is vaguely lower semicontinuous on E(A) (cf [4]).
Remark If f ¼ 0, g ¼ 1 and À is only one connected manifold, then the potential U 0 (x) is constant for x 2 À, and d 0 is given by the so-called natural layer or Robin density [3].

The variational formulation in a trace space
In this section, we show that for the compact plates À 1 , À 2 given as two closed Lipschitz surfaces in R 3 , the Gauss problem (2.1) is equivalent to a variational problem in a particular Sobolev-Slobodetski space on À ¼ À 1 [ À 2 . Let & R 3 be the domain (bounded or unbounded) with the boundary @ R 3 ¼ À, and let H 1 2 ðÀÞ be the space of traces of the Sobolev space H 1 () onto À [19].
Note that À is Lipschitz. Let C 1 (À) be the trace space on À of C 1 (R 3 ), and define for ' 2 C 1 (À) Then C 1 (À) is dense in the trace space H 1 2 ðÀÞ, its closure with respect to the norm is given by (3.1) [19].
Moreover, the surface measure ds on À is well defined and generates on C 1 (À) the L 2 scalar product, We denote that extension by the same symbol (Á, Á). Since À is Lipschitz, the function space C 1 (À) is dense in each of the spaces H 1 2 ðÀÞn L 2 (À), and H À 1 2 ðÀÞ. For the Gauss variational problem (2.1) we shall find an equivalent formulation which is based on special distributions concentrated on À with densities ' 2 H À 1 2 ðÀÞ. These define bounded linear functionals on H In H À 1 2 ðÀÞ, let us define the affine cone KðA, a, gÞ :¼ where a 1 , a 2 are the two given positive constants. From now on, for the given positive function g we require g 2 CðÀÞ \ H 1 2 ðÀÞ. As we shall see, the solution of the Gauss problem (2.1) can be obtained with the help of the simple layer potential operator which for 2 H À 1 2 ðÀÞ is well defined almost everywhere on the Lipschitz surface À [3,20,21].
Let us collect some of the properties of V, for Lipschitz surfaces À proved by Costabel [20] and Verchota [21], in the following theorem.
where ' denotes equivalence. We shall call the distribution in H À 1 2 ðÀÞ associated with AE.
Proof Let 2 C 1 (À), then : ð yÞdsð yÞ dAEðxÞ and with Fubini's theorem, where É, and the energy scalar product in (3.11) is thus well defined. Therefore, for all 2 C 1 ðÀÞ, where c AE and c 0 are not depending on . Since C 1 (À) is dense in H 1 2 ðÀÞ, AE defines a bounded linear functional on H 1 2 ðÀÞ, and there exists a uniquely defined element 2 H À 1 2 ðÀÞ such that (3.9) holds. In order to show (3.10) we consider (y) :¼ jx À yj À1 for any fixed x 2 R 3 n À as a function of y. Then 2 C 1 (À) and ð yÞdsð yÞ ¼ VðxÞ because of (3.9), and V(x) here denotes the simple layer potential for x = 2 À. This relation holds for every x = 2 À, therefore the potentials U AE and V coincide everywhere in as well as in c ¼ R 3 n . As is well known (see e.g. [18, Theorem 1.20]), Hence, U AE j 2 H 1 () and U AE j c 2 H 1 ( c ). Moreover, U AE solves the following transmission problem (see [3, Section 5.6.3] and [22]): In (3.13), almost everywhere means with respect to the surface measure on À and the limits are supposed to be nontangential, while in (3.14) n @ and n @ c denote the interior unit normal vectors with respect to and c , respectively, and the equality holds in the weak sense. Then it follows [20] that which is equivalent to kk 2 H À 1 2 and, furthermore, equals kAEk 2 E due to (3.12). This completes the proof of Theorem 3.2. g As a consequence of Theorem 3.2, we conclude for the solution of the Gauss problem the following.
Proof Since A is compact, there exists a unique solution 0 of the Gauss problem in E(A, a, g) (Section 1). Hence, I( 0 )51. Because of Theorem 3.2, there exists In order to find a solution of the Gauss problem (2.1) with a suitable algorithm we replace the affine cone K(A, a, g) by a cone K(À) with vertex at 0 by employing Lagrange multipliers for the side conditions R À j g' j ds ¼ a j , j ¼ 1, 2. Then we arrive at the following problem: which is the minimizer of Namely,

gÞ:
Here %40 is a penalty parameter to be chosen later on appropriately.
For ' 0 we have the same estimate due to 2ðÀÞ : g THEOREM 3.5 Let ' 0 be the minimizer of (2.1) and ' % be the minimizer of (3.17).
Then for we find the following estimates: , and a 0 :¼ minfa 1 , a 2 g: Here, kVk denotes the norm of the operator V : H À 1 2 ðÀÞ ! H 1 2 ðÀÞ. Moreover, The proof of Theorem 3.5 is given in the appendix. Theorem 3.5 shows that the solution of the Gauss problem (2.1) can be approximated by the solutions of the minimizing problem (3.17) by choosing % in (3.17) large enough. The advantage of (3.17) lies in the fact that for (3.17) the cone K(À) is a cone with vertex 0 and not an affine cone as in (2.1). For solving (3.17) we shall use the gradient projection method to be described in Section 5.

The condenser problem
This section deals with the case where f ¼ 0, g ¼ 1, and a 1 ¼ a 2 ¼ 1. Observe that then G f (A, a, g), being equal to the minimum of the energy over the class of all 2 E(A) with (À j ) ¼ 1 for j ¼ 1, 2, is actually equal to [cap A] À1 , where cap A denotes the Newtonian capacity of the condenser A. Hence, is actually the Newtonian equilibrium measure on A.
On the other hand, it follows from results of the third author [5,9,11] that Ã A , the Newtonian equilibrium measure on A, can also be obtained as the (unique) solution to the (dual) minimal energy problem over the class of all Borel measures of finite energy supported by À such that U ðxÞ ¼ CðÞ for all x 2 À 1 , where C() is a constant depending on . What is important is that the measures , admissible in the dual variational problem, are no longer required to belong to E(A); that is, the measures j À j , j ¼ 1, 2, are allowed to be signed.
In view of Theorem 3.2, this leads us to the following formulation in H À 1 2 ðÀÞ. THEOREM 4.1 Let ¼ 1 À 2 2 H À 1 2 ðÀÞ with supp j & À j , j ¼ 1, 2, be the distribution associated with the Newtonian equilibrium measure Ã ¼ Ã A . Then is also the minimizer of the energy over the set H À 1 2 ðÀÞ consisting of all ¼ 1 À 2 2 H À 1 2 ðÀÞ with supp j & À j satisfying the system of boundary integral equations where c 2 R is a parameter.

THEOREM 4.2
The distribution ¼ 1 À 2 associated with the Newtonian equilibrium measure Ã ¼ Ã A can be obtained with the following procedure: Solve the system of boundary integral equations solve the following system of boundary integral equations for :
Proof Because of (4.3), the energy in (4.2) can also be written as where 1 , 2 and the energy kk 2 V depend on c, the constant value of the potential V on À 1 . The energy kk 2 V becomes minimal if @ @c where _ j :¼ @ j =@c. Since the solution of the minimum problem in Theorem 4.1 is the distribution associated with the Newtonian equilibrium measure Ã A , satisfies the conditions ð j , 1Þ ¼ cap A for j ¼ 1, 2, the last two terms in (4.7) cancel each other and the equation (4.7) now reads Àcd þ ð _ 2 , 1Þ ¼ 0, which gives the required value of C.
Taking at C the derivatives @/@c of the solution of the minimum problem, we find i.e. system (4.4). Since the simple layer potential operator V on À 1 Â À 2 is strongly elliptic [3,20], the system (4.4) has a unique solution _ 2 H À 1 2 ðÀÞ. Then C is uniquely determined provided d 6 ¼ 0, and with c :¼ C we obtain as the (unique) solution of the strongly elliptic boundary integral equation system (4.6). g Proof Without loss of generality, we can assume À 2 to contain boundary limit points of the unbounded component 1 of R 3 n À, for if not, we replace the roles of À 1 and À 2 . Let us consider the simple layer potential U _ in R 3 nÀ, which has boundary values almost everywhere on À. Then we conclude from (4.4) that U _ ¼ 1 a.e. on À: Since U _ 2 H À 1 2 ðÀÞ, the normal derivative @U _ @n exists almost everywhere on À and belongs to H À 1 2 ðÀÞ [21]. Here, n denotes the interior normal vector with respect to 1 . With the jump relation of the normal derivative of the simple layer potential U _ across À [20] and since U _ ðxÞ ¼ 1 for all x in the bounded components of R 3 nÀ, we find Since U _ ðxÞ ¼ Oðjxj À1 Þ for jxj ! 1, Green's theorem can be applied to U _ in 1 , which yields

The approximation of the Gauss problem by the use of piecewise constant charges
Since the minimizing problem (3.17) is quadratic, one could compute the minimizer by using the gradient-projection method: where P K is the H À1/2 (À)-orthogonal projection onto K(À) [13, I(3.13)], and V 0 f,% the Frechet derivative of the functional V f,% in (3.17). But a numerical realization of P K is as difficult as solving the original problem (3.17) itself. Therefore we shall apply a gradient-projection iteration to the discretized formulations instead.
For solving (3.17) numerically, let us use a quasiregular family of triangulations T h of the surface À ¼ À 1 [ À 2 [24,Chapter 10], where h denotes the maximum mesh width of the elements in T h . On the triangulation, we introduce piecewise constant functions S 0 h ðÀÞ & L 2 ðÀÞ & H À1=2 ðÀÞ as the trial as well as the test space.
Correspondingly, we define The finite-dimensional approximation to (3.4) then reads: Find the minimizer ' %h 2 K h (À) of the quadratic functional If M denotes the number of elements in T h , which here is also the dimension of S 0 h , then minimizing (5.3) on K h (À) defines a quadratic programming problem with M linear constraints given by (5.2). Since the bilinear form As M, in general, is rather large, standard methods for finding the discrete minimizer of (5.3) requires long computing times (see [16], where just few hundred degrees of freedom were used). On the other hand, (5.3) is of the type (3.8) in [13, Chapter I]. We therefore write the Frechet derivative of (5.3) operating on ' %h 2 S 0 h in variational form: This Frechet derivative can be realized by using a multipole approximation [14] and is equivalent to a linear system of M equations on M variables collected in the vector x 2 R M , the coefficients of ' 1 0h À ' 2 0h with respect to the corresponding basis defined by the characteristic functions of ' Inserting the representations of ' 1 h À ' 2 h with the coefficients x into (5.3), the discretized minimizing problem now reads: Minimize ð5:6Þ But this problem is a problem in R M with a positive definite symmetric matrix A and, hence, of the type (3.8) in [13]. Therefore, we now may use the gradient projection method [13,Chapter I (3.12)] in the Euclidean space R M :

ð5:7Þ
where is appropriately chosen, 0 The projection P R M þ is given by Obviously, kP R M þ k L 2 ðR M Þ,L 2 ðR M Þ ¼ 1 and, therefore, (5.7) converges in R M as the geometric series P 1 k¼0 ð1 À Þ k . Although all matrix times vector multiplications are executed with the fast multipole approximation [14], because of (5.9) and as will be seen from the computations, the convergence turns out to be rather slow.

Some computations with the projection iteration
As computational examples, we consider a sphere for À 1 and for À 2 a rotational symmetric surface with radius R(x 1 ) decaying to zero for x 1 ! 1. The surfaces are given by RðXÞ for x 1 ¼ Xg:

ð5:11Þ
In particular, we choose The infinite rotational body is cut off at X ¼ 4.
For first tests of the gradient-projection method, we consider the Gauss problem (2.1) with g ¼ 0.1, f (x) ¼ 1000jx À x 0 j À1 where x 0 ¼ (À4.0, 0, 0) > , and use an approximation of À 1 and À 2 (X ¼ 4) with 2048 and 2088 triangles, respectively, see Figure 1. We computed the charges a 1 ¼ 848.8032 and a 2 ¼ 564.4113 for % ¼ 0 and then tried to recover the charges for several values of % 2 {1, 10, 100, 1000} by employing the gradient-projection iteration (5.7). We observed that decreases rapidly for large % and hence, a very slow convergence. In addition, the residual is still large after millions of iteration steps for large %, but at the same time we observed a bad approximation of the charge distribution.
Next, we checked the convergence of the approximate total chargesã 1 andã 2 for % 2 {0, 10, 100, 1000} in the case of active constraints for chosen values of a 1 ¼ 400 and a 2 ¼ 300. In Table 1, we observe convergence of the computed approximations to a 1 and a 2 for increasing %, which was proposed in (3.24). A relative accuracy of 10 À4 was used as stopping criterion in the gradient-projection method. For % ¼ 1000, however, the iteration was stopped after 3 Â 10 6 steps.
Since the charges converge very slowly for large %, we now use a cascading approach. We start with small values of % to compute a good approximation of the charge distribution with a few iteration steps. Then this approximation serves as an initial guess for a new run of the gradient-projection method with a larger % to improve the approximation of the total charges. In each run, we limit the maximal number of iteration steps. Repeating this approach several times, we get a good approximation with a significantly smaller total number of iteration steps.
For our example, we choose % 2 {10, 50, 100, 300, 500, 700, 800, 900, 950, 1000}, limit the number of iteration steps in each run to 30,000 and compare the results with the results of the gradient-projection method for % ¼ 1000 and 3 Â 10 6 iterations. In Table 2, we see a good agreement of the valuesã 1 andã 2 of the gradient-projection method and the cascading approach with a total of 3 Â 10 5 iterations. In addition, the functional of the approximate solution of the cascading approach is smaller. This is reflected even better in Figure 1. Altogether we get a much better approximation of the charge distribution by the cascading approach than with the pure gradient-projection method.  6. The approximation of the condenser problem For the numerical solution of the two systems (4.4) and (4.6) with (4.5) we use the same boundary triangulations T h of À ¼ À 1 [ À 2 as in Section 5. But now as test and trial functions in the Galerkin approximation we employ the full space S 0 h and set with h' the characteristic function to the triangle ' 2 T h \ À j where I 1 ¼ {1, . . . , M 1 } for the triangulation of À 1 and I 2 ¼ {M 1 þ 1, . . . , M} for the triangulation of À 2 . With the Galerkin matrix A 0 of V on À 1 Â À 2 , the Galerkin equations to (4.4) now are equivalent to the linear system for the coefficients _ x of Then the approximations of d in (4.5) and C are obtained by solving the Galerkin equations of (4.6), i.e. the additional linear system which is the desired approximation. The solutions of the linear systems (6.2) and (6.3) are obtained with a preconditioned conjugate gradient method as in [14] where the preconditioner to (6.2) as well as to (6.3) is based on an artificial multilevel approach due to Steinbach in [25].
The corresponding matrix times vector multiplications are executed by using the fast multipole method [14] for the simple layer potential operator.
Since for this condenser problem the two linear systems (4.4) and (4.6) in H À 1 2 ðÀÞ are approximately solved by Galerkin's method, standard a priori error estimates for h ! 0 are available [26]. LEMMA 6.1 Let _ , 2 H s ðÀÞ with À 1 2 s 1, and À1 À s t s 1, and let the family of triangulations be quasiregular. Then there hold the asymptotic a priori If the family of triangulations is not quasiregular then (6.4) holds only for À1 À s t À 1 2 s 1.
Since the constant right-hand sides in (4.4) and (4.6) are in C 1 (À), the regularity, i.e. s depends only on the regularity of À.
In Table 3, we compare the solutions of the condenser problem for several levels of discretization of the two bodies, which were considered in Section 5.1 and are given by (5.11) and (5.12) for X ¼ 4. The meshes are created by uniform refinement of the coarsest mesh with 308 triangles, where the new vertices are projected onto À 1 and À 2 . We consider the approximate solution 7 on the finest refinement level with about five million triangles as a reference solution and compute an indicated error as error :¼ k L À 7 k L 2 ðÀÞ =k 7 k L 2 ðÀÞ : We observe a reduction of the error but a low order of convergence eoc as predicted in (6.4). But the convergence is a lot better for the approximation of the charges at the tip of À 2 , which is of further interest in this example.
The charges at the tip of the rotational body À 2 are defined by the integral Z for increasing values of the length X. For these computations, we use highly adaptive meshes with 50,240 up to 44,6372 triangles for the approximation of À 1 and À 2 . In Table 4, the charges at the tips of À 2 are divided by the areas at the tips for several X. In addition, the computed densities, which are given by the quotients of the charges and the areas of the tips, indicate the existence of a limit, which was proven to exist in [5]. With these computations, however, we reach the limits of the currently available software based on the fast and data-sparse realization of the single layer potential by the fast multipole method. For larger X, additional techniques will be necessary such as, e.g. appropriate domain decomposition methods.
As a second example, we consider again two bodies given by (5.11), where now RðrÞ ¼ e Àr for r ! 0: ð6:6Þ The solution for a discretization with 82,364 boundary elements for X ¼ 4 is given in Figure 2. The solution is positive on À 1 and negative on À 2 as required. Next, we investigate the asymptotic behaviour of the charges at the tips for X 2 {4, 6, 8}. About 10 5 boundary elements were used for the discretization of À 2 for each X. In Table 5, the extremal values of these computations are given. The fast decrease in the minimum observed at the tip, indicates already the divergence of the charges at the tips for X ! 1, as predicted in [5].  Table 4. Charges (6.5) and their densities at the tips for a ¼ X À 10 À8 . The plots of the solution on the sphere À 1 for X 2 {4, 6, 8} show no visible difference. We only can observe a slight increase in the maximal values, see Table 5. On the disc at x 1 ¼ 0, which is opposite to the sphere, the solutions do not change for X 2 {4, 6, 8}. We plotted the solution for X ¼ 6 in Figure 3 to indicate that the charge distributions have the correct signs.
Since the radius R in (6.6) decreases exponentially, we had to adapt the mesh size of the boundary elements also exponentially along the body. Otherwise, it was not possible to run the simulations for X ¼ 6 and X ¼ 8. For larger X, the number of shape regular elements increases so rapidly that the computations were not possible anymore.
As in our first example, we are interested in the charges (6.5) at the tips for several cutoffs X. These values are divided by the areas of the tips, see Table 6, for X 2 {4, 6, 8}. The computed densities at the tips indicate their divergence for X ! 1 as predicted in [5].

Convergence for h ! 0
For the convergence behaviour for h ! 0 and X fixed in Lemma 6.1, the regularity of À and, hence, the value if s for , _ 2 H s ðÀÞ is needed. In both examples À 1 2 C 1 whereas À 2 has an edge À ¼ {x ¼ (0, cos ', sin ')}, ' 2 [0, 2), where the tangent vectors on À 2 perpendicular to À perform a jump of angle 7 8 (with respect to 1 ). Moreover, our example is rotationally symmetric. Therefore and _ , the normal derivatives of U and U _ , respectively, have the form r À4/7 (r) where r is the distance from to y ¼ (x 1 , R(x 1 ) cos ', R(x 1 )sin ') 2 À 2 for x 1 40, and to x ¼ (0, (1 À r)cos ', Figure 3. Charge density at the ball and the rotational body for X ¼ 6. Table 6. Charges at the tips and densities for X 2 {4, 6, 8} and a ¼ X À 10 À8 . (1 À r)sin ') 2 À 2 for x 1 ¼ 0, and where is a smooth function of r. At the second edge of À 2 at x 1 ¼ X, the corresponding tangent vector's jump is smaller and the solutions are less singular. Using Fourier transform in R 2 , it then follows that , _ 2 H s ðÀÞ for s 5 1 14 . Consequently, we find from Lemma 6.1 that k _ h À _ k H t ðÀÞ þ k h À k H t ðÀÞ ch 1=14Àt for À 15 14 t s 5 1 14 and jC h À Cj þ jðcapAÞ h À capAj ch s 0 with s 0 5 8 7 :