Boundary integral equations on special surfaces and their applications

Cost-effective iterative processes are constructed for the solution of the boundary integral equations in the potential theory, defined on twoor three-dimensional surfaces of special type. Algorithms are discussed which are obtained by combining the method of boundary integral equations with grid methods. Examples are given of numerical solution of three-dimensional problems of magnetostatics for unbounded domains. The method of boundary integral equations, or, in a more general formulation, the boundary element method, is widely used in solving theoretical and applicationsoriented problems [2,8,9,16,50]. From the standpoint of numerical realization, it offers an alternative to the finite difference and finite element methods based on a differential formulation of the original problem. In general, nonsymmetric dense matrices are obtained when the boundary integral equations are discretized. Only some information on the spectrum localization is available for such matrices. In practical problems, the dimension of the system of linear algebraic equations may reach 10-10. The interest to iterative methods of solving such systems is therefore obvious. The rate of convergence when solving discrete counterparts of the integral equations of the second kind is pratically independent of the discretization step. Correspondingly, the solution of a system of order N to the accuracy ε = N~" (a > 0) requires O(N\nN) arithmetic operations even for methods like the stationary Richardson method. This review surveys efficient iterative methods for solving the integral equations of the potential theory on surfaces of special type. It is shown that these methods yield cost-effective combined algorithms which employ equations in the differential form and boundary integral equations. The analysis is carried out for twoand threedimensional problems of magnetostatics which arise when stationary Maxwell's equation is solved. In fact, a considerable part of the results formulated in the paper is valid for a wider range of problems in mathematical physics. Boundary integral relations for a number of differential operators are given, for example, in [2]. The current status of numerical algorithms for mathematical problems of magnetostatics is reflected in the monograph [16] and the references cited therein, and in the review papers [13,39,42,49]. A particular progress was made in the methods of solving grid equations [7,20,21,24,27,28,33,35-37] in bounded domains. Fast algorithms for equations with constant coefficients in rectangular, parallelepiped, cylindrical, etc. domains occupy a special place. This paper analyses efficient t Originally published in Russian as Preprint No. 137 of the Department of Numerical Mathematics of the USSR Academy of Sciences, Moscow, 1986. 464 E. P. Zhidkov and B. N. Khoromsky algorithms for solving integral equations at the boundaries of domains of these types. They can be used to develop cost-effective methods of solving boundary value problems in unbounded domains. Thus the integral equations of magnetostatics written for a domain with constant magnetic permeability are used to treat exactly the boundary conditions at infinity. A comparative analysis of the differential and integral formulations [8,13,16,39,42, 49,55,58] leads to the following conclusion. Combined methods, bringing together the advantages of each of these approaches [4,12,14,15,18,58], can be exploited to optimize the numerical algorithms for the boundary value problems of magnetostatics in unbounded twoor three-dimensional domains. In fact, one realizes then a version of the domain decomposition method (without overlapping) in which one of the subdomains is unbounded. The aspect of justifying such algorithms for elliptic problems is discussed in [1,4,11,14,18,26,30,31,52]. A generalization to the case of quasilinear elliptic equations was given in [25]. Methods of introducing the boundary condition at infinity in the differential formulation are discussed in [10,13,15,32]. Artificial boundary conditions [38] can be employed in two-dimensional problems in order to account for the asymptotic behavior of a function which is harmonic in the neighborhood of the infinitely removed point [38]. The prospects for the application of three-dimensional integral equations are outlined in [13,42,49]. Special algorithms are needed to carry out an efficient numerical realization of combined methods for problems in an unbounded domain, because integral equations must be solved a large number of times within some iterative process. In fact, each single solution of an integral equation in the general case is computationally more costly than the solution of a linear differential equation in the inner domain. For example, it will suffice to consider the surface of a cube in the case of the space R or the boundary of a square in the case of R. If the grid with n unknown variables is uniform in all directions or the approximation of integrals is piecewise-constant, the array necessary for the discretization of the boundary equation contains 72n and 32n elements, respectively. Direct methods of solution would require 72n and (64/3)w operations, respectively (we omit the cost of computing the matrix elements). These asymptotic estimates of the required computational cost remain valid for surfaces of more general type as well. The combined algorithms can be optimized with respect to computational cost if the known properties of integral operators are taken into account in choosing auxiliary surfaces. The main property is that the kernels of the integral operators for the potentials of the single and double layers are functions of only the difference of scalar arguments for the boundaries obtained either by translation or by rotation of a part of the boundary. If grid is uniform, the matrices of discretized integral equations for such surfaces contain blocks which are either Toeplitz or circulant matrices [46,47,51]. As a result, it becomes possible to reduce the memory required to store the matrices by an order of magnitude in n and to reduce the computational cost of solving the corresponding systems of equations [46,51,55,56]. Thus a numerical algorithm was constructed in [55] for an integral equation on the surface of a parallelepiped; it requires (asymptotically) even a smaller number of operations than the solution of a grid problem within the domain bounded by this surface. The aspects of computations with Toeplitz-type matrices are discussed, for example, in [46]. For a special surface one can choose in R surfaces of a ball, a cylinder, or a parallelepiped with a square base, and a circle or a rectangle in R. Cost-effective Boundary integral equations 465 algorithms for solving integral equations on these surfaces are given in Sections 2 and 3. In order to use this formulation at the inner points of the domain bounded by an auxiliary surface, one can approximate the original differential equation by the finite differences or finite elements. The sought function and its derivative along the normal to the boundary at the boundary points satisfy the discrete counterpart of the integral equation, with a dense coefficient matrix. The arising combined system of equations is solved by an alternating iterative procedure. A grid boundary value problem with the Dirichlet (Neumann) condition at the auxiliary boundary and an integral equation corresponding to the exterior Neumann (Dirichlet) problem are solved at each iteration, with subsequent relaxation. The use of a sequence of progressively finer grids is an efficient method of improving the convergence rate of the iterative process and of monitoring the accuracy [3,16,17,29,43,54]. Section 1 describes the properties of the boundaty integral operators and treats certain iterative methods of solving the corresponding integral equations in the general case. Sections 2 and 3 analyse the structure of matrices which approximate the integral equations on special surfaces, and also the algebraic properties of these matrices. Cost-effective iterative methods of solving the obtained systems of equations are then discussed, and the required numbers of memory locations and arithmetic operations for the realization of the algorithms suggested are estimated. Section 4 contains estimates of the rate of convergence of the alternating iterations without preliminary discretization of equations. The problem of choosing an optimal relaxation parameter is discussed. Section 5 gives the results of numerical experiments with a program which realizes a cost-effective method of solving the boundary integral equation on the surface of a parallelepiped. This section deals with the solution of three-dimensional boundary value problem in an unbounded domain, using an alternating iterative process on a sequence of grids. The above-mentioned program for solving the integral equation [56] was employed. The suggested method can be applied, for example, for the solution of exterior boundary value problems involving the Laplace operator, in twoand threedimensional domains of complex shape. 1. FORMULATION OF BOUNDARY INTEGRAL EQUATIONS. PROPERTIES OF THESE EQUATIONS Let us consider boundary integral equations obtained by a straight-forward application of Green's formulae for harmonic functions in the space R (n = 2,3). Let a closed Lyapunov surface [44] Γ divide the space R into two simply connected domains, ΩΙ and Ωβ, one interior and the other exterior, so that Ω^ΓυΩ«, = R, where Ωβ is an unbounded domain. Let us consider a problem (Problem A) of determining a function u(M) (MeR), which is harmonic in one of the domains Ω, or Ωβ, from the values it assumes at the boundary Γ (Dirichlet problem) or from the values assumed by its normal derivative on Γ (Neumann problem). In the case of the exterior problem (domain Ωβ) we assume M(OO) = 0. We denote by d/dn the derivative along the inward (with respect to Ω;) normal to the boundary Γ and set