HUNGARIAN JOURNAL OF INDUSTRIAL CHEMISTRY VESZPRÉM Vol. 33(1-2). pp. 119-124. (2005) NEUMANN BOUNDARY VALUE PROBLEMS WITH BEM AND COLLOCATION N. HERRMANN1 AND V. KARNJANATAWEE2 1Institut für Angewandte Mathematik, Universität Hannover, Hannover, Germany 2Dept. of Mathematics, King Mongkut University of Technology Thonburi, Bangkok, Thailand This paper is an introduction into the Boundary Element Method (BEM) with collocation to find the numerical solution of two different types of Neumann problems. At first we start with the Laplace equation and continue later with the heat equation on a bounded convex domain with smooth boundary in two dimensions. We will show how to transform the governing problem into a boundary integral equation which can be solved by dividing the boundary into a finite number of segments and applying the collocation method. We finish presenting an example of the heat equation. Keywords: Boundary element method, Laplace equation, heat equation, Neumann boundary value problem, Collocation Introduction Let Ω be a convex open domain in 2R with smooth boundary . Consider the Laplace equation: Ω∂ in Ω, ( ) 0=∆ y,xu ( ) ( y,xg n y,xu = ∂ ∂ ) on Ω∂ , where 2 2 2 2 yx ∂ ∂ + ∂ ∂ =∆ and the heat equation ( ) ( ) ( ) , y t,y,xu x t,y,xu t t,y,xu 2 2 2 2 ∂ ∂ + ∂ ∂ = ∂ ∂ where ( ) ( ]T,t,y,x 0∈Ω∈ with initial and boundary conditions, ( ) ( ,y,xf,y,xu =0 ) ( ) Ω∈y,x ( ) ( ) ( ) ( ]T,t,y,x,t,y,xg n t,y,xu 0∈Ω∂∈= ∂ ∂ By using the fundamental solution we will develop the Boundary Integral Equation (BIE) from the governing boundary value problem and will create a Fredholm integral equation of the second kind. According to the fact that the solution of a Neumann value problem is not unique the same is true for the approximate solution which we get by using BEM and collocation. Piecewise constant and piecewise linear functions will be used to form the ansatz functions which lead to a simple linear system of equations. Neumann Boundary Value Problem for the Laplace Equation Let Ω be a convex open domain in 2R with smooth boundary Ω∂ . Consider the Neumann boundary value problem (NBVP): ( ) 0=∆ y,xu in Ω, (1) ( ) ( y,xg n y,xu = ) ∂ ∂ on Ω∂ . (2) It is well known that for two dimensions ( ) π ξ ξ 2 − −= xlog ,xE (3) is the Green's function or the fundamental solution for the operator ∆, that means it is the solution of the problem ( ) ( ) 2R∈∀−−=∆ ξξδξξ ,xx,xE . (4) It can be seen that E is not defined at ξ=x , where E is singular. 120 Boundary Integral Equation For we follow the definition of the distribution Ω∈x δ, and with (1) and (4) ( ) ( )( )uxxu ξδ −= ( ) ( ) ( )( ) ( )(∫Ω ∆−∆= ξξξξξ ξ du,xE,xEu ) The second Green-Gauß formula gives ( ) ( ) ( ) ( ) ( ) Ω∈∀⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ − ∂ ∂ = ∫ Ω∂ xdun ,xE ,xE n u xu ξσξ ξ ξ ξ .(5) There is no singularity for . Therefore for Ω∈x Ω∈x we can compute u(x) by knowing the boundary data u(x) and n u ∂ ∂ on . We will refer to this formula as the representation of the solution Ω∂ u(x) for . Ω∈x Our Neumann problem from above gives n u ∂ ∂ on . We have to find Ω∂ u(x) = g(x) on . Ω∂ Both integrals in (5) ( ) ( ) ( ) ( )∫∫ Ω∂Ω∂ ∂ ∂ ∂ ∂ ξ ξ ξ σξ ξ σξ ξ du n ,xE d,xE n u and (6) are well defined for . But for we have a jump of magnitude Ω∈x Ω∂∈x ( ) 2 xu for the second integral : ( ) ( ) ( ) ( ) ( )∫∫ Ω∂Ω∂Ω∂∈→ ∂ ∂ += ∂ ∂ ξ ξ ξ ξ σξ ξ σξ ξ du n ,xExu du n ,xE lim xx 00 20 and so we obtain the boundary integral equation ( ) ( ) ( ) ( ) ( ) ( )∫ Ω∂ ⎟⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∂ ∂ − ∂ ∂ += ξ ξ σξ ξξ ξ du n ,xE n u ,xE xu xu 00 0 0 2 (7) Ω∂∈∀ 0x Then we get from (7) a Fredholm integral equation of the second kind ( ) ( ) ( ) ( ) ( ) ξξ ξ σ ξ ξσξ ξ d n u ,xEdu n ,xExu ∫∫ Ω∂Ω∂ ∂ ∂ = ∂ ∂ + 2 (8) Ω∂∈∀x where u on the left hand side is unknown and ( ) ( )xg n xu = ∂ ∂ on is given. Ω∂ Collocation Method The integral equation (9) does generally not admit a solution in closed form. We will show how to solve such a problem numerically using the collocation method. Given ( ) ( ) ,x,xg n xu Ω∂∈= ∂ ∂ solve for the unknown Dirichlet data u(x), Ω∂∈x from ( ) ( ) ( ) ( ) ( ) ,dg,xEdu n ,xExu ξξ ξ σξξσξ ξ ∫∫ Ω∂Ω∂ =∂ ∂ + 2 (9) Ω∂∈∀x The solution of (9) is not unique since if we replace u(x) by ( ) ( ) cxuxu += where c is constant then u is a solution, too. We need a compatibility condition so that the Neumann boundary value problem is well posed: g(x) should satisfy ( ) ( ) ( )∫∫∫ Ω∂Ω∂Ω =∂ ∂ =∆= σσ dxgd n xu dxxu0 (10) We use now the collocation method to discretize the problem. We subdivide the boundary into n arcs and call the midpoints of the arcs with . We take the ansatz: Ω∂ n,,, ΓΓΓ K21 nx,,x,x K21 1xxn = (11) ( ) ( )∑ = = n i ii xaxu ~ 1 χ where ( )xiχ is the characteristic function on and are unknown parameters. That means iΓ ia ( )xu~ is a piecewise constant approximation to u(x) on . As collocation points we use the midpoints. This leads to Ω∂ ( ) ( ) ( ) ,dg,xEd n ,xE a a j n i j i j i ∫∑ ∫ Ω∂ = Γ = ∂ ∂ + ξξ ξ σξξσ ξ 12 (12) nj ≤≤1 These are n linear equations with the unknowns , so that we have na,,a K1 n × n linear system of equations Neumann Boundary Value Problem for the Heat Equation Let Ω be a bounded convex domain with smooth boundary Γ=Ω∂ in 2R . Consider the heat equation ( ) ( ) ( ) , y t,y,xu x t,y,xu t t,y,xu 2 2 2 2 ∂ ∂ + ∂ ∂ = ∂ ∂ (13) ( ) ( ]T,t,y,x 0∈Ω∈ with initial and boundary conditions, ( ) ( ),y,xf,y,xu =0 ( ) Ω∈y,x (14) ( ) ( ) ( ) ( ]T,t,y,x,t,y,xg n t,y,xu 0∈Γ∈= ∂ ∂ . (15) The problem will be transformed into an integral equation by using the fundamental solution and will be solved by applying the collocation method in the same manner as for the Laplace equation (1) and (2). Boundary integral equation The well-known Green’s function or fundamental solution for the heat equation 121 ( ) ( ) ( ) ⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ ≤ > −= − −− τ τ τπτξ τ ξ t te tt,x;,E t x if0 if 4 1 4 2 , (16) is used as weight function to generate the integral equation ( ) ( ) ( )∫ ∫ Ω =⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ −∆ T ddt,;,E t,u ,u 0ττ τ τ ξxξ ξ ξξ where ( t,u x ) is the unknown solution of the heat equation, ( ),, 21 ξξ=ξ ( ),y,x=x 2 2 2 2 1 2 ξξ ∂ ∂ + ∂ ∂ =∆ ξ . Let Ω be a bounded convex domain in 2R with smooth boundary Ω∂=Γ . By using the Gauss-Green formula the integral becomes: ( ) ( ) ( ) ( ) ( ) ( ) ( )∫∫ ∫ ∫ ∫ ΩΓ Γ + ∂ ∂ − ⋅ ∂ ∂ = ξξξτσ τξ τξ τστξ τξ ξ ξ ξ ξ dt,x;,Efdd n t,x;,E ,u ddt,x;,E n ,u t,xu t t 0 0 0 for ( ) ( ]T,t,x 0×Ω∈ The left hand side of the formula has to be replaced by the famous jump relation if x is on the boundary. So we are led to a boundary integral equation for Γ∈x and . Tt ≤<0 ( ) ( ) ( ) ( ) ( ) ( ) ( )∫ ∫ ∫ ∫ ∫ Ω Γ Γ + ∂ ∂ − ⋅ ∂ ∂ = ξξξ τσ τξ τξ τστξ τξ ξ ξ ξ ξ dt,x;,Ef dd n t,x;,E ,u ddt,x;,E n ,u t,xu t t 0 2 1 0 0 Substitute the given boundary condition into ( ) ξ τξ n ,u ∂ ∂ the boundary integral equation becomes ( ) ( ) ( ) ( ) ( ) ( ) ( )∫∫ ∫ ∫ ∫ ΩΓ Γ +⋅= ∂ ∂ + ξξξτστξτξ τσ τξ τξ ξ ξ ξ dt,x;,Efddt,x;,E,g dd n t,x;,E ,ut,xu t t 0 2 1 0 0 Collocation Method We divide the boundary into n segments as shown. Ω∂=Γ Let . Γ∈=+ 1121 xx,x,,x,x nnK Define ( ){ }101 12 ≤≤−+=∈=Γ ′ + λλλ ,xxx:x iii R where n,,i K1= . Then is a polygon in Ω since Ω is convex. nΓ′Γ′ UKU1 Let jn′ v be the outer unit normal vector on iΓ′ , i n,,i maxh Γ′= = K0 be the step width in space and for N∈m m T k = be the step width in time and we get the time steps ( ,m). jkt j = K,,j 10= Using collocation method the piecewise linear spline functions ( ) n,,i, otherwise x, xx xx x, xx xx x i ii i i ii i i K1 0 1 1 1 1 1 = ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ Γ′∈′ − ′− Γ′∈′ − −′ =′ + + − − − ϕ are applied to space and the piecewise constant characteristic functions ( ) ( ] ( ] m,,j, t,tt t,tt t jj jj j K1 1 0 1 1 = ⎪⎩ ⎪ ⎨ ⎧ ∈ ∉ = − − χ are applied to time. With both we form the ansatz function ( ) ( ) (∑ ∑ = = ′=′ m j n i ji j i txut,xu ~ 1 1 χϕ ) where u are the unknown coefficients which we have to find. j i At the time step p of segment i the boundary integral equation becomes xn+1=x1 x2 xixi+1 xn Γ Γ' ni 122 ( ) ( ) ( ) ( ) ( ) ( )∫ ∫ ∫ ∫ ∫ Ω′ Γ′ ′ Γ′ ′ ′ ′′′+ ′′⋅′= ′ ∂ ′∂ ′+ ξξξ τστξτξ τσ τξ τξ ξ ξ ξ dt,x;,Ef ddt,x;,E,g dd n t,x;,E ,u~u pi t t pi t pip i p p 0 2 1 0 0 This is a system of linear equations. The coefficients can be written in short form as ( ) ( ) ∫ ∫ − Γ′ ′ ′ ′ ∂ ′∂ ′= q q t t qi j pq ij ddn t,x;,E b 1 τσ τξ ξϕ ξ ξ ( ) ( ) ( ) ( )∫ ∫ ∫ Ω′ Γ′ ′ ′′′+ ′′⋅′= ξξξ τστξτξ ξ dt,x;,Ef ddt,x;,E,gc pi t t pi p i p 0 0 The boundary integral equation at time step p of segment i is rewritten as: [ ] [ ] ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ + ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ −= ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ + ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∑ − = p n p q n q p q pq ij p n p pp ij c c u u b u u b MMMO 111 1 1 2 1 2 1 which should be solved to get the solution. Example As example we consider the heat equation ( ) ( ) ( ) ( ) ( ]10 2 2 2 2 ,t,y,x, y t,y,xu x t,y,xu t t,y,xu ∈Ω∈ ∂ ∂ + ∂ ∂ = ∂ ∂ with the following initial and Neumann boundary conditions, ( ) ,,y,xu 10 = ( ) Ω∈y,x (17) ( ) ( ) ( ]100 ,t,y,x, n t,y,xu ∈Γ∈= ∂ ∂ (18) where Ω is the unit circle. Forming the boundary integral equation and substituting the data given in (17) and (18) we get ( ) ( ) ( ) ( )∫∫ ∫ ΩΓ = ∂ ∂ + ξτσ τ τ ξ ξ dt,;,Edd n t,;,E ,ut,u t xξ xξ ξx 0 2 1 0 , ( ) ( ]10,t,y,x ∈Γ∈=x To apply the collocation method we divide the perimeter of the unit circle into 6 segments and discretize the time into 10 time steps. With piecewise linear functions used in space and piecewise constant functions ( 61 ,,i,i K=ϕ ) ( )101 ,,j,j K=χ used in time, the ansatz function is ( ) ( ) (∑ ∑ = = ′=′ 10 1 6 1j i ji j i txut,xu ~ χϕ ) where u are the unknown constant collocation points. ji Then the boundary integral equation at time step p becomes ( ) ( ) ( ) 610 2 1 0 ,,i,dt,;,E dd n t,;,E ,u~u pi t pip i K=′′= ∂ ′∂ ′+ ∫ ∫ ∫ Ω Γ′ ′ ′ ξ τσ τ τ ξ ξ xξ xξ ξ So at the first time step we have the equation ( ) ( )∫ ∫ ∫ Ω Γ′ ′ ′ ′′= ∂ ∂ +++ ξξ τσϕϕ ξ ξ dt,x,,E dd n E uuu i t i i 1 0 1 6 1 61 1 1 1 0 2 1 K , where i = 1, …, 6 and ( ) ξξ τξ ′′ ∂ ′∂ = ∂ ∂ n t,x;,E n E ii 1 1 . This is a system of linear equations ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ∫ ∫∫ ∫ ∫ ∫∫ ∫ ∫ ∫∫ ∫ Γ ′∪Γ′ ′Γ′∪Γ′ ′ Γ ′∪Γ′ ′Γ′∪Γ′ ′ Γ ′∪Γ′ ′Γ′∪Γ′ ′ 1 6 1 2 1 1 0 2 1 6 0 1 1 6 0 2 1 2 0 1 1 2 0 6 1 1 0 1 1 1 1 21 1 16 1 21 1 16 1 65 1 16 2 1 00 0 2 1 0 00 2 1 u u u dd n E dd n E dd n E dd n E dd n E dd n E tt tt tt M M L MM L L OMM L L τσϕτσϕ τσϕτσϕ τσϕτσϕ ξξ ξξ ξξ ( ) ( ) ( )⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ′ ′ ′ = ∫ ∫ ∫ Ω′ Ω′ Ω′ 16 12 11 ;0 ;0 ;0 t,,E t,,E t,,E x x x ξ ξ ξ M . The quantities ξ ′∂ ∂ n Ei 1 with ( ) ( ) ( ) ⎪ ⎪ ⎩ ⎪⎪ ⎨ ⎧ ≤ > −=′ − −′− τ τ τπτξ τ ξ 1 1 4 11 if0 if 4 1 1 2 t te tt,x;, t x i i E are ( ) ( ) ( )τξ ξ ξ τπξ − −′ − ′ ′ ′ − ⋅−′ −=⋅ ′∂ ∂ = ∂ ∂ 1 2 4 2 1 11 8 tiii i e t E n E xξ nxξ n . The term ( ) ξξ ′⋅−′ nx i can be shown in the picture below ix y jx 1+jxξ ′ ξ ′n 123 It is clear that ( ) ( ) ( )( ) ξξ ξξ ′′ ⋅−′+−=⋅−′ nn yxyx ii ji dxy =−= . jd are constant. Then the elements of the coefficient matrix of above system are ( ) ( ) ( ) ( ) ∫ ∫ ∫ ∫∫ ∫ Γ′ ′ − −′ − Γ′ ′ − −′ − − Γ′∪Γ′ ′ − − − −= ∂ ∂ −− j i j i jj t j t x j t j t x j t j i dde t d dde t d dd n E 1 1 2 1 1 1 2 1 1 0 4 1 0 4 1 1 0 1 8 8 ξ τ ξ ξ τ ξ ξ στϕ τπ στϕ τπ τσϕ . We integrate analytically with respect to the time variable and get ∫ ∫∫ ∫ Γ′ ′ −′ − Γ′ ′ −′ − − Γ′∪Γ′ ′ −′ − −′ −= ∂ ∂ −− j i j i jj de x d de x d dd n E j t x i j j t x i j t j i ξ ξ ξ ξ ξ σϕ ξπ σϕ ξπ τσϕ 1 2 1 1 2 1 1 4 2 4 2 1 0 1 2 2 . We calculate each term numerically and obtain the system of linear equations ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ 3282980 3282980 3282980 3282980 3282980 3282980 1 6 1 5 1 4 1 3 1 2 1 1 . . . . . . u u u u u u A ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = 0.50020200.000490000200.00049002020 0020200.50020200.000490000200.00049 0.000490020200.50020200.00049000020 0000200.000490020200.50020200.00049 0.000490000200.0004900202050002020 0020200.000490000200.000490020200.5 ... ... ... ... .... ... A The solution of the first time step is ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ = ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ 650050 650050 650050 650050 650050 650050 1 6 1 5 1 4 1 3 1 2 1 1 . . . . . . u u u u u u Conclusions As shown in this paper the Boundary Element Method (BEM) can be applied to elliptic and parabolic differential equations. The benefit of the method compared to conventional methods such as finite different and finite element method is the reduction of dimension of the problem by one since the BEM deals only with the corresponding boundary integral equation. Moreover for applying the collocation method the unknown trial or ansatz function can be formed by simple functions i.e. piecewise constant and piecewise linear functions. However the disadvantage of the method is the knowledge of the fundamental solution. So its application is limited to linear differential equations SYMBOLS pq ij pq ij pq ij c,b,a parameters h space step k time step i, j index parameters m number of time intervals n number elements jn v unit normal vector p, q time parameters t time variable ( )t,y,xu solution n iu collocation parameters x, y space variable x space vector Greek letters α a constant χ, ϕ trial function Γ boundary of the domain Γ' polygonal boundary ξ, τ integral variable Ω domain Ω' polygonal domain REFERENCES 1. Brebria C.A. and Walker S.: Boundary Element Techniques in Engineering, Newnes - Butterworths, London, 1980 2. Costabel M.: Boundary Integral Operators for the heat equation, Integral Equations OperatorTheory, 1990, vol.13, 498 - 552 3. Costabel M., Stephan E.P. :., On the Convergence of Collocation Methods for Boundary Integral Equations on Polygons, Mathematics of Computation, 1987, vol.49, 461 -478 4. Costabel M., Onishi K., Wendland W. L.: Boundary Element Collocation Method for the Neumann Problem of the Heat Equation, Academic Press Inc., 1987 5. Dyck U.: Randelement-Lösungen für die Wärmeleitungsgleichung, Master thesis in Mathematics, Hannover University, Germany, 1992 6. Friedman A.: Partial Differential Equations of Parabolic Type, Prentice - Hall. Inc., 1964 7. Herrmann N.: BEM with Collocation for the Heat Equation with Neumann and mixed boundary 124 values, AMS Contemporary Mathematics, 2002, vol.295, 265 - 277 8. Herrmann N.: Improved Method for solving the Heat Equation with BEM and Collocation, AMS Contemporary Mathematics, 2003, vol.329, 165 - 174 9. Herrmann N.: Time discretization of linear parabolic problems, Hungarian Journal of Industrial Chemistry, 1991, vol.19, 275 - 281 10. Herrmann N.: Numerical Problems in Determining Pore-Size Distribution in Porous Material, Workshop at the University of Budapest, Invited lecture, 1995 11. Herrmann N., Siefer J., Stephan E.P., and Wagner R., Mathematik und Umwelt, Edition Univ. Hannover, Theodor Oppermann Verlag, Hannover, 1994 12. Herrmann N. and Stephan E.P., FEM und BEM Einführung, Eigendruck Inst. f. Angew. Math., Univ. Hannover, 1991 13. Iso Y.: Convergence of Boundary Element Solutions for the Heat Equation, Journal of Computational and Applied Mathematics, 1991, vol.38, 201 – 209