eu SOLVING THE DIRICHLET PROBLEM WITH DISCONTINUOUS COEFFICIENTS IN BOUNDED MULTIPLY CONNECTED REGIONS USING A BOUNDARY INTEGRAL EQUATION WITH THE GENERALIZED NEUMANN KERNEL

This paper presents a numerical method for solving the Dirichlet problem with discontinuous coefficients in bounded multiply connected regions. The method is based on reducing the problem of solving the Dirichlet problem with discontinuous coefficients to a problem of solving Dirichlet problem with continuous coefficients. The Dirichlet problem with continuous coefficients is then solved using a combination of a uniquely solvable boundary integral equation with generalized Neumann kernel and the Fast Multipole Method. Numerical results and comparison are given to illustrate the efficiency of the suggested method. Received: March 5, 2014 c © 2014 Academic Publications, Ltd. url: www.acadpubl.eu Correspondence author 448 M. Aghaeiboorkheili, A.H.M. Murid AMS Subject Classification: 30E25, 31B10

AMS Subject Classification: 30E25, 31B10 Key Words: Dirichlet problem with discontinuous coefficients, boundary integral equations, generalized Neumann kernel, simply connected region, multiply connected region

Introduction
The boundary integral equation method is a classical method for solving the Dirichlet and Neumann boundary value problems.The classical boundary integral method for the Dirichlet problem and the Neumann problem with continuous coefficients are in the form of second kind Fredholm integral equations with the Neumann kernel.These integral equations are derived by representing the solutions as the potential of double layer for the Dirichlet problem and single layer for the Neumann problem [1].
Recently, the interplay of the Riemann-Hilbert problem and integral equation with the generalized Neumann kernel has been investigated in [18,19].It has been shown that the problem of conformal mapping, Dirichlet problem, Neumann problem, and mixed Dirichlet-Neumann problem can all be treated as Riemann-Hilbert problems [12,15,16,19,20].Hence they can be solved efficiently using integral equations with the generalized Neumann kernel.
Dirichlet problem with discontinuous coefficients in rectangular regions have been solved by separation of variable method, Fourier series and Green function [6,7].Also, Dirichlet problem with discontinuous coefficients in smooth boundary regions have been solved by conformal mapping and Poisson's integral formula [2,4,8,9].
In this paper, we solve Laplace's equation on bounded simply connected region with Dirichlet condition with discontinuous coefficients via an integral equation with the generalized Neumann kernel.The results of this paper extend the results presented in [15] where the Dirichlet problem involved continuous coefficients only.

Dirichlet Problem with Discontinuous Coefficients in Bounded
Simply Connected Region

Preliminaries
Let G be a bounded simply connected region in the extended complex plane C with boundary Γ in counterclockwise direction(see Figure 1).Suppose that α is a fixed point in G.The curve Γ is assumed to has a 2π-periodic twice continuously differentiable complex function η(t) with non-vanishing first derivative, i.e., Assume that H be the space of all real Hölder continuous functions on the boundary Γ.From now on, for complex-valued or real-valued function ψ ∈ H defined on the boundary Γ and for t ∈ J, we will not distinguish between ψ(η(t)) and ψ(t).
We define the function A by with a fixed point α ∈ G.The generalized Neumann kernel formed with A is defined by [18] The kernel N is continuous with Define also a real kernel M by [18] The kernel M can be represented by with a continuous kernel M 1 which takes on the diagonal the values Hence, the operator [18] is a Fredholm integral operator and the operator is a singular integral operator.
Theorem 1. [10,18] If F is a solution of the Dirichlet boundary value problem Re [F (η(t))] = γ(t) with boundary values and Im F (0) = 0, then the imaginary part of F satisfies the integral equation

Dirichlet Problem with Discontinuous Coefficients
Suppose that γ is a given function with finite number of jump discontinuities at the points t k ∈ [0, 2π) for k = 1, 2, . . ., l. Denote the limiting values of γ(t) at the points of discontinuities by Find a function u harmonic in G + such that its boundary values on smooth boundary Γ satisfy [4,8,9] u(η(t)) = γ(t).
The unique solution u of the Dirichlet problem (13) can be regarded as a real part of an analytic function in G, F (z).We allow the boundary Γ to have corner points at η(t k ) for k = 1, 2, . . ., l.Let where each value of the argument arg t k is chosen such that the branch cut from the point η(t k ) to the point at ∞ is outside of the domain G.
represents the interior angle between the left and right tangents to the boundary Γ at η(t k ).We assume that 0 The representation (36) is different from the one used in [15].The branch cut of log t k is chosen to be outside of the domain G from the point η(t k ) to the point at where and It is clear that the function φ(t) is continuous on [0, 2π] − {t 1 , . . ., t l }.The function φ(t) is also continuous at the each point t j for j = 1, 2, . . ., m.To show that, we should show that Since γ k (t) is continuous at t j for each j = k.Hence, we have Thus we have succeeded in reducing the problem of solving the Dirichlet problem with discontinuous coefficients (13) to a problem of solving Dirichlet problem with continuous coefficients (37).Without lose the generality, we assume that Im then by Theorem (1) we have By solving (41), we obtain ψ, then we obtain f (z) from ( 40), so we can obtain F (z) from (36) and our final answer u(z) will be the real part of the F (z).

Dirichlet Problem with Discontinuous Coefficients in Bounded
Multiply Connected Region

Notations and Auxiliary Material
Let G be a Bounded multiply connected region, of connectivity m + 1 ≥ 1, in the extended complex plane C whit boundary Γ = ∪ m j=0 Γ j consisting of m + 1 smooth closed Jordan curves Γ j , j = 0, 1, 2, . . ., m.The curves Γ 1 , . . ., Γ m are inside the curveΓ 0 .The complement G − is defined as G − := C \ G (see Figure 2).Suppose that α is a fixed point in G.The curve Γ 0 has counterclockwise orientations and the curves Γ 1 , . . ., Γ m always have clockwise orientations(in fact the orientation of boundary Γ is such that G be on the left of Γ ).A 2π-periodic twice continuously differentiable complex function η j (t) with non-vanishing first derivative ηj (t) = dη j (t)/dt = 0, t ∈ J j := [0, 2π], (21) j = 0, 1, 2, . . ., m is a parameterization for the curve Γ j .we will represent the disjoint union of the intervals J j by J as a total parameter domain .Automatically, the whole boundary Γ will be parameterized by complex function η defined on J by Assume that H be the space of all real Hölder continuous functions on the boundary Γ.Since η is smooth, a function φ ∈ H can be explained by φ(t) := φ(η(t)), t ∈ J, as a real Hölder continuous 2π-periodic functions φ(t) of the parameter t ∈ J, i.e., φ(t) := with real Hölder continuous 2π-periodic functions φj defined on J j ; and vice versa.
From now on, for complex-valued or real-valued functions ψ defined on the boundary Γ and for t ∈ J, we will not differentiating between ψ(η(t)) and ψ(t).For t ∈ J k , the values ψ(t) will be indicated by ψ k (t).
Suppose that γ ∈ H is a given function , the Dirichlet problem can be described as follows (see e.g.[1, p.307], [8, p. 145] and [9, p. 164]): Suppose that χ [j] be the piecewise constant function defined on J by for j = 0, 1, . . ., m.Then, we define the space S by S = span{χ [0] , χ [1] , . . ., χ [m] }. (25) We define also the space S by It follows from the definition of the space S that a function h ∈ S if and only if h can be written as with real constants h 0 , h 1 , . . ., h m .We define an operator R : S → S by Theorem 2 ([15]).Let γ be a given function.Then, there exists a unique function h ∈ S and a unique function µ such that are boundary values of a single-valued analytic function f in G with Im f (α) = 0.The function µ is the unique solution of the integral equation

Dirichlet Problem with Discontinuous Coefficients
Suppose that γ(t) is a given function in the following form where for at least one j = 0, 1, . . ., m the function γ j (t) has finite number of first kind discontinuities at t j k ∈ J i , k = 1, . . ., l j .Find a function u continuous on the closure G(except probably at t j k ), harmonic in G, such that its boundary values satisfy on Γ u(η(t)) = γ(t).
In the other hand, the discontinuity points of the functions γ j (t) on Γ are in the following form: it means we have l 0 discontinuity points on Γ 0 , l 1 discontinuity points on Γ 1 , . . ., l m discontinuity points on Γ m .
The unique solution u of the Dirichlet problem (33) can be regarded as a real part of an analytic function F (z) in G. Suppose that the discontinuity points of the function γ j (t) are at the points t j k ∈ [0, 2π) for k = 1, 2, . . ., l j .Denote the limiting values of γ j (t) at the points of discontinuities by where each value of the argument arg is chosen such that the branch cut form the point η j (t j k ) is outside of the domain G. Then b + j,k − b − j,k represents the interior angle between the left and right tangents to the boundary Γ at η j (t j k ).We assume that 0 where and It is clear that the function φ(t) is continuous on [0, 2π]−{t j k : j = 0 . . ., m & k = 1 . . ., l j }.The function φ(t) is also continuous at the each point t j k .To prove that, we should show that Since γ j k (t) is continuous at t s i for each i = k and s = j.Hence, we have Thus we have succeeded in reducing the problem of solving the Dirichlet problem with discontinuous coefficients (33) to a problem of solving Dirichlet problem with continuous coefficients (37).Without lose the generality, we assume that Im[f ( 0 then by Theorem (2) we have By solving (41), we obtain ψ, then we obtain f (z) from ( 40), so we can obtain F (z) from (36) and our final answer u(z) will be the real part of the F (z).

Numerical Examples
Since the functions A j and η j are 2π-periodic, a reliable procedure for solving the integral equation (41) numerically is by using the Nyström method with the trapezoidal rule [1].Using the Nyström method with trapezoidal rule reduces solving the integral equation ( 41) to solving a linear system [1] where B is the n × n discretizing matrix of the operator N, C is the n × n discretizing matrix of the operator M, I is the identity matrix, and x is an n × 1 vector of the approximate values of the unknown function ψ at the nodes and y is an n × 1 vector of the values of the function φ at the nodes.Since the integral equation ( 41) is uniquely solvable, then for sufficiently large n, the linear system ( 42) is uniquely solvable [1].See [11,12,13,14,15,17] for more details.For some ideas on how to solve numerically boundary integral equations with the generalized Neumann kernel on regions with corners, see [17].
In this paper, the linear system (42) will be solved by the MATLAB function gmres.The function gmres can be used with a matrix-vector product function, i.e., it is not necessary to form the matrices B and C explicitly.The matrix-vector product functions for the matrices B and C are defined using the function zfmm2dpart in the MATLAB toolbox FMMLIB2D developed by Greengard and Gimbutas [5].Using the Fast Multipole Method (FMM), the product of the matrix B by a vector can be computed in O(n) operations.Computing the matrix-vector product Cy in the right-hand side of the linear system (42) requires using both the FMM and the FFT which requires O(n ln n) operations.Thus, the linear systems (42) will be solved in O(n ln n) operations.See [14,13] for more details.For comparisons, we use MATLAB toolbox pdetool to obtain the solutions of examples by finite element method.In the following numerical examples, u(z) is the exact solution of the Dirichlet boundary value problem with discontinuous coefficients, u n (z) is the approximate solution obtained with n node points on the boundary obtained by the presented method, and u f em (z) is the approximate solution obtained by the finite element method.
Example 1. Find a function u harmonic in the unit disk having the boundary value 1 in the right half plane and 0 in the left half plane.This example has been considered in [3, Example 4, p. 248] using the conformal mapping and Poisson's integral formula method.The exact solution is given by This function u has two points of discontinuities at t 1 = π 2 (z 1 = A) and t 2 = 3π 2 (z 2 = B).By using the explanations in Section 2.2 we have: So we have, = u(t) − 1 π arg t 1 (η(t) − i) + 1 π arg t 2 (η(t) + i).The surface plot of the solution for Example 1 obtained by our method is shown in Figure 4 and the absolute error |u n (z) − u(z)| for the entire region obtained with n = 1024 by our method is shown in Figure 5.In Figure 6 we can see the absolute error |u f em (z) − u(z)| for the entire region obtained by finite element method with 130561 points.As it is clear in the Figure 5 and Figure 6, our method is more accurate than finite element method, especially at the points near the discontinuity points.

Figure 1 .
Figure 1.A bounded simply connected region in the extended complex plane

Figure 2 .
Figure 2. A bounded multiply connected region in the extended complex plane

Figure 3 .
Figure 3.The region for Example 1

Figure 4 .
Figure 4.The surface plot of the solution for Example 1 that is obtained by our method