Numerical Solution for Elliptic Interface Problems Using Spectral Element Collocation Method

Abstract

The aim of this paper is to solve an elliptic interface problem with a discontinuous coefficient and a singular source term by the spectral collocation method. First, we develop an algorithm for the elliptic interface problem defined in a rectangular domain with a line interface. By using the Gordon-Hall transformation, we generalize it to a domain with a curve boundary and a curve interface. The spectral element collocation method is then employed to complex geometries; that is, we decompose the domain into some nonoverlaping subdomains and the spectral collocation solution is sought in each subdomain. We give some numerical experiments to show efficiency of our algorithm and its spectral convergence.

1. Introduction

Second order elliptic interface problems have many applications in engineering and science, in particular in fluid dynamics and material science [1–5]. Multiphase flow, crystal growth, modeling of the Stefan problem of solidification process, composite materials, cell and bubble deformation, and fiber suspension are few examples of the interface problems to be mentioned. Owing to discontinuity of coefficients of differential equations and singular source term, the solution and its derivative of interface problem are nonsmooth or even discontinuous across the interface. Hence usual numerical methods for partial differential equations do not work efficiently for these problems. Such kind of problems has attracted attention of scientists from both theoretical and computational approach over the past decades. Several high-order methods have been developed to deal with interface problem in terms of discontinuous coefficients and/or singular source terms. In [6] LeVeque and Li have investigated second order accurate finite difference method for one-dimensional interface problem. For one-dimensional problem, the standard finite element method can be used and if the interfaces lie on the grid point, second order accurate solution can be gained [7]. However by modifying basis function, Li obtained a second order accurate finite element method in the maximum norm in [7]. The finite difference and finite element method for one-dimensional problem are not much different. In [8], Shin and Jung used spectral collocation method and radial basis function method for one-dimensional interface problems and obtained very accurate solutions using the similar ideas given by Loubenets et al. in [9]. For two-dimensional elliptic interface problem, in [6], the immersed interface method (IIM) is developed for an elliptic interface problem with discontinuous coefficients and singular source term. By adding an extra node to the standard finite difference 5-point stencil, Huang and Li derived a second order accurate method in [10]; but the coefficient matrix is nonsymmetric and the method is stable for problems with piecewise constant coefficient. This problem is solved by a new version of IIM in [11]. In [12], Angelova and Vulkov presented second and fourth order finite difference scheme for elliptic problem with intersecting interface which has discontinuous variable coefficient and singular source term. However they assumed that interfaces are straight lines. In [13], Li et al. investigated finite element method based on Cartesian triangulation for interface problem involving discontinuous coefficients which is second order accurate in the maximum norm. Hanke and Loubenets proposed an immersed finite element method for elliptic interface problem with discontinuous coefficients and singular source term and present its second order convergence in 𝐿² norm in [14]. Kwak et al. have introduced an immersed finite element method based on the “broken” 𝑃₁-nonconforming piecewise linear polynomials on interface triangular element having edge average as degrees of freedom in [15]. These linear polynomials are broken to match the homogeneous jump condition along the interface which is allowed to cut through the element. They have proved optimal order of convergence in the 𝐻¹- and 𝐿²-norms. Latige et al. presented second order accurate finite volume method for interface problem with discontinuous coefficients in [16]. This paper is a generalization of one-dimensional work given in [8] to two-dimensional elliptic interface problems including curve boundaries.

In this paper, we first derive weak formulation of second order elliptic interface problem along with jump conditions. Then we give an algorithm to solve numerically these problems in the simple domain, namely, rectangle with straight line interface. We generalize the algorithm for domain with a curved boundary and interface by using a transformation so-called Gordon and Hall map. For a complicated domain, the spectral element collocation method is adopted by decomposing the domain into some nonoverlapping subdomains and then the spectral collocation method is used in each subdomain. Finally some numerical experiments are given to show efficiency and accuracy of our method.

2. Pseudo-Spectral Collocation Method

In this section we present an algorithm for solving an elliptic interface problem by the spectral collocation method. First we start with an elliptic interface problem with a discontinuous constant coefficient and a singular source term defined on a simple domain with line interface and then we extend it to an elliptic interface problem with any arbitrary domain and interface. We review some simple facts including function spaces, interpolation operators, and spectral matrices for the use of the spectral collocation method from now on. The standard notations and definitions are used for the weighted Sobolev spaces 𝐻^𝑠_𝜔(Ω) equipped with weighted inner product (⋅,⋅)_𝑠,𝜔 and corresponding weighted norms ‖⋅‖_𝑠,𝜔, 𝑠≥0, where 𝜔(𝑥,𝑦)=𝜔(𝑥)𝜔(𝑦) is the Legendre weight function when 𝜔(𝑡)=1 and Chebyshev weight function when 𝜔(𝑡)=1/√1−𝑡². Let 𝒫_𝑁 be the space of all polynomials of degree less than or equal to 𝑁 and let {𝜉_𝑖}^𝑁_𝑖=0 be the Legendre Gauss Lobatto (LGL) or Chebyshev Gauss Lobatto (CGL) points on [−1,1] such that −1=∶𝜉₀<𝜉₁<⋯<𝜉_𝑁−1<𝜉_𝑁∶=1. For Legendre case, {𝜉_𝑖}^𝑁_𝑖=0 is the set of zeros of (1−𝑡²)𝐿^′_𝑁(𝑡), where 𝐿_𝑁 is the 𝑁th Legendre polynomial and the corresponding quadrature weights {𝑤_𝑖}^𝑁_𝑖=0 are given by 𝑤₀=𝑤_𝑁=2𝑁(𝑁+1),𝑤_𝑗=2𝑁()1𝑁+1𝐿_𝑁𝜉_𝑗²,1≤𝑗≤𝑁−1. (1) For Chebyshev case, {𝜉_𝑖}^𝑁_𝑖=0 is the set of zeros of (1−𝑡²)𝑇^′_𝑁(𝑡), where 𝑇_𝑁 is the 𝑁th Chebyshev polynomial and the corresponding quadrature weights {𝑤_𝑖}^𝑁_𝑖=0 are given by 𝑤₀=𝑤_𝑁=𝜋2𝑁,𝑤_𝑗=𝜋𝑁,1≤𝑗≤𝑁−1. (2) For any continuous function 𝑢 on [−1,1], let 𝐼_𝑁𝑢 denote its Lagrange interpolation at collocation points {𝜉_𝑖}^𝑁_𝑖=0; that is, 𝐼_𝑁𝑢𝜉_𝑖=𝑢𝜉_𝑖,𝑖=0,1…,𝑁. (3) Let {𝜓_𝑗}^𝑁_𝑗=0⊂𝒫_𝑁 be the Lagrange basis functions of degree 𝑁 such that 𝜓_𝑗𝜉_𝑘=𝛿_𝑗𝑘∀𝑗,𝑘=0,1,…,𝑁. (4) Then 𝐼_𝑁𝑢(𝑥)=^𝑁_𝑗=0𝑢𝜉_𝑗𝜓_𝑗(𝑥). (5) The pseudo-spectral derivative 𝜕_𝑁𝑢 of a continuous function 𝑢 is defined to be the exact derivative of the interpolant of 𝑢; that is, 𝜕_𝑁𝑢𝜉_𝑖=^𝑁_𝑗=0𝜓^′_𝑗𝜉_𝑖𝑢𝜉_𝑗. (6) Then the pseudo-spectral derivative matrix 𝐷_𝑁 is 𝐷_𝑁𝑖,𝑗∶=𝜓^′_𝑗𝜉_𝑖. (7) Let 𝑈 be the vector valued function containing nodal values of 𝑢 at 𝜉_𝑗; that is, 𝑈=(𝑢(𝜉₀),…,𝑢(𝜉_𝑁))^𝑇; then the derivative vector valued function of 𝑢^′ is 𝐷_𝑁𝑈=(𝑢^′(𝜉₀),…,𝑢^′(𝜉_𝑁))^𝑇. If the interval [−1,1] is replaced by [𝑎,𝑏], we can use the following linear transformation: 𝑡=𝑏−𝑎2(𝑥+1)+𝑎∶[−1,1]⟶[𝑎,𝑏] (8) to find Gauss-points {𝜉̂_𝑗}^𝑁_𝑗=0 and the quadrature weights {𝑤_𝑗}^𝑁_𝑗=1𝜉̂_𝑗=𝑏−𝑎2𝜉_𝑗+1+𝑎,𝑤_𝑗=𝑏−𝑎2𝑤_𝑗. (9) Introducing 𝜓_𝑗(𝑡)∶=𝜓_𝑗(𝑥)=𝜓_𝑗2(),𝑏−𝑎𝑡−𝑎−1 (10) we easily obtain the following spectral matrix 𝐷_𝑁𝐷_𝑁=2𝐷𝑏−𝑎^𝑁by𝐷_𝑁𝑖,𝑗∶=𝜓^′_𝑗𝜉̂_𝑖=2𝜓𝑏−𝑎^′𝑗𝜉_𝑖. (11) The two-dimensional LGL and CGL nodes {𝐱_𝑖𝑗} and weights {𝑤_𝑖𝑗} are defined as 𝐱_𝑖𝑗=𝜉_𝑖,𝜉_𝑗,𝑤_𝑖𝑗=𝑤_𝑖𝑤_𝑗,𝑖,𝑗=0,1,…,𝑁. (12) Let 𝒬_𝑁 be the space of polynomials of degree less than or equal to 𝑁 with respect to each variable 𝑥 and 𝑦. The basis functions are also defined 𝜓_𝑖𝑗𝑥,𝑦=𝜓_𝑖(𝑥)𝜓_𝑗𝑦,𝑖,𝑗=0,1,…,𝑁. (13) We reorder the LGL and CGL points from bottom to top and then from left to right such that 𝐱_𝑗 for 𝑗=1,…,(𝑁+1)². Then pseudo-spectral derivative matrix on 2-dimensional space is defined via the Kronecker tensor product; that is, 𝑆_𝑥=𝐷_𝑁⊗𝐼_𝑁,𝑆_𝑦=𝐼_𝑁⊗𝐷_𝑁,𝑆_𝑥𝑥=𝐷²_𝑁⊗𝐼_𝑁,𝑆_𝑦𝑦=𝐼_𝑁⊗𝐷²_𝑁, (14) where 𝐼_𝑁 is identity matrix of the same order as 𝐷_𝑁.

Let Ω=(𝑎,𝑏)×(𝑐,𝑑) be a quadrilateral domain and let Γ={𝛼}×(𝑐,𝑑) be an interface separating the domain Ω into two subdomains Ω⁺ and Ω⁻, such that Ω=Ω⁺∪Ω⁻∪Γ. Here, Γ is referred to as an interface. The boundary of Ω is denoted by 𝜕Ω and 𝜕Ω⁺=Ω⁺∩𝜕Ω, 𝜕Ω⁻=Ω⁻∩𝜕Ω. Consider the following elliptic interface problem: −∇⋅𝛽∇𝑢+𝑘𝑢=𝑓+𝑔𝛿_ΓinΩ,𝑢=0on𝜕Ω, (15) where 𝑔∈𝐶(Γ),𝑓∈𝐿²(Ω) and 𝛽,𝑘∈𝐿^∞(Ω), with 𝛽=𝛽𝑥,𝑦⁺𝑥,𝑦𝑥,𝑦∈Ω⁺𝛽⁻𝑥,𝑦𝑥,𝑦∈Ω⁻,𝑘=𝜅𝑥,𝑦⁺𝑥,𝑦𝑥,𝑦∈Ω⁺𝜅⁻𝑥,𝑦𝑥,𝑦∈Ω⁻, (16) and 𝛿_Γ is a 2-dimensional delta function with the support along the interface Γ.

In order to solve this problem via the interface problem, we need to first derive the jump conditions. It is obvious from the given equation that there exists one jump discontinuity by the Dirac 𝛿-function. Using the standard finite element argument [6, 9], we have the following classical formulation: −∇⋅𝛽∇𝑢^±+𝑘𝑢^±=𝑓^±inΩ^±,𝑢^±=0on𝜕Ω^±, (17) along with the following interface jump conditions: [𝑢]|_Γ=0,[𝛽∇𝑢⋅𝐧]_Γ=𝑔, (18) where 𝐧=(𝑛₁,𝑛₂) denotes the outward normal vector on 𝜕Ω⁺. The interface jump is defined as follows: [𝑣]|_Γ=𝑣⁺−𝑣⁻, (19) where 𝑣⁺ and 𝑣⁻ are the traces of 𝑣|_Ω⁺ and 𝑣|_Ω⁻, respectively, on Γ. The existence and uniqueness of solution of (17) is studied in [17].

Here we consider the pseudo-spectral method for the problem (17)-(18) with 𝛽 and 𝜅 as piecewise continuous constants. Suppose that 𝑓 is a piecewise continuous function and 𝑔 is continuous on Γ. Suppose that 𝑢_𝑁∈𝒬_𝑁 is the pseudo-spectral approximation solution to problem (17)-(18). Then the approximation solution of the interface problem (17) can be expressed by 𝑢^±𝑥,𝑦=^𝑁𝑁_𝑖=0𝑗=0𝑢^±_𝑖𝑗𝜓^±_𝑖𝑗𝑥,𝑦. (20) Although it is possible to use different polynomial order approximation on each subdomain, we use the same polynomial order approximation, for the sake of simplicity. Hence we have the following equations: −𝛽^±Δ𝑢^±𝜉_𝑖,𝜉_𝑗+𝜅^±𝑢^±𝜉_𝑖,𝜉_𝑗=𝑓^±𝜉_𝑖,𝜉_𝑗∀𝜉_𝑖,𝜉_𝑗∈Ω^±, (21) which can be written in matrix-vector form as 𝐴₁𝑈=−𝛽⁺𝑆⁺_𝑥𝑥+𝑆⁺_𝑦𝑦+𝜅⁺𝐼𝟎𝟎−𝛽⁻𝑆⁻_𝑥𝑥+𝑆⁻_𝑦𝑦+𝜅⁻𝐼×𝑈⁺𝑈⁻=𝐹⁺𝐹⁻=𝐹₁, (22) where 𝑈^±=(𝑢^±(𝐱_𝑖)) and 𝐹^±=(𝑓^±(𝐱_𝑖)). Here 𝑆⁺_𝑥𝑥=2𝛼−𝑎²𝐷²_𝑁⊗𝐼_𝑁,𝑆⁺_𝑦𝑦=2𝑑−𝑐²𝐼_𝑁⊗𝐷²_𝑁,𝑆⁻_𝑥𝑥=2𝑏−𝛼²𝐷²_𝑁⊗𝐼_𝑁,𝑆⁻_𝑦𝑦=2𝑑−𝑐²𝐼_𝑁⊗𝐷²_𝑁. (23) Let 𝑈^±=[𝑈^±_bd𝑈^±_in𝑈^±_if]^𝑡, where 𝑈^±_if, 𝑈^±_bd, and 𝑈^±_in denote by the values of 𝑢^± at nodal points on the interface, boundaries, and interior of domain Ω^±, respectively. By jump conditions [𝑢]|_Γ=0 and [𝛽∇𝑢⋅𝐧]|_Γ=𝑔, we have 𝑢⁺−𝑢⁻=0 and 𝛽⁺∇𝑢⁺⋅𝐧−𝛽⁻∇𝑢⁻⋅𝐧=𝑔, respectively, which can be written in matrix-vector form as 𝑈⁺_if−𝑈⁺_if𝛽=0,₊𝑆⁺_𝑥⋅𝑛₁+𝛽⁺𝑆⁺_𝑦⋅𝑛₂𝑈⁺_if−𝛽⁻𝑆⁻_𝑥⋅𝑛₁+𝛽⁻𝑆⁻_𝑦⋅𝑛₂𝑈⁻_if=𝐺, (24) where 𝐺=(𝑔(𝐱_𝑖)) for 𝐱_𝑖 on Γ. The boundary conditions can be imposed as 𝑈^±_bd=0. (25) Now, the boundary and jump conditions can be represented as 𝐴₂𝑈=𝐼0000000000000𝐼00−𝐼000𝐼0000000000𝛽⁺𝑆⁺_𝑥𝑛₁+𝛽⁺𝑆⁺_𝑦𝑛₂00−𝛽⁻𝑆⁻_𝑥𝑛₁−𝛽⁻𝑆⁻_𝑦𝑛₂×𝑈⁺_bd𝑈⁺_in𝑈⁺_if𝑈⁻_bd𝑈⁻_in𝑈⁻_if=00000𝐺=𝐹₂. (26) We also note that 𝑆⁺_𝑥=2𝐷𝛼−𝑎^𝑁⊗𝐼_𝑁,𝑆⁺_𝑦=2𝐼𝑑−𝑐^𝑁⊗𝐷_𝑁,𝑆⁻_𝑥=2𝐷𝑏−𝛼^𝑁⊗𝐼_𝑁,𝑆⁻_𝑦=2𝐼𝑑−𝑐^𝑁⊗𝐷_𝑁. (27) Combining two systems, we have the following linear system: 𝐴𝑈=𝐹, (28) where 𝐴=𝐴₁+𝐴₂ and 𝐹=𝐹₁+𝐹₂. The resulting algebraic system can be efficiently solved by direct and iterative methods.

We now generalize the above algorithm for the case of discontinuous varying functions 𝛽 and 𝑘. We still assume that the domain Ω is rectangle and the interface is straight line as before. In this case, the first equation in (17) can be expanded as follows: −𝛽^±Δ𝑢^±+∇𝛽^±⋅∇𝑢^±+𝑘^±𝑢^±=𝑓inΩ^±. (29) Hence the pseudo-spectral method for (29) reads: −𝛽^±𝜉_𝑖,𝜉_𝑗Δ𝑢^±𝜉_𝑖,𝜉_𝑗+∇𝛽^±𝜉_𝑖,𝜉_𝑗⋅∇𝑢^±𝜉_𝑖,𝜉_𝑗+𝑘^±𝜉_𝑖,𝜉_𝑗𝑢^±𝜉_𝑖,𝜉_𝑗=𝑓𝜉_𝑖,𝜉_𝑗∀𝜉_𝑖,𝜉_𝑗∈Ω^±, (30) which can be written in matrix-vector form as 𝐴₁𝑈=𝐹₁ (31) with 𝐴₁=−𝜷⁺𝑆⁺_𝑥𝑥+𝑆⁺_𝑦𝑦+𝜷⁺_𝑥𝑆⁺_𝑥+𝜷⁺_𝑦𝑆⁺_𝑦+𝐤⁺𝐼𝟎𝟎−𝜷⁻𝑆⁻_𝑥𝑥+𝑆⁻_𝑦𝑦+𝜷⁻_𝑥𝑆⁻_𝑥+𝜷⁻_𝑦𝑆⁻_𝑦+𝐤⁻𝐼, (32)𝑈=[𝑈⁺𝑈⁻]^𝑡 and 𝐹₁=[𝐼_𝑁𝑓⁺𝐼_𝑁𝑓⁻]^𝑡, where 𝜷^±,𝜷^±_𝑥,𝜷^±_𝑦 and 𝐤^± are diagonal matrices of functions 𝛽,𝛽_𝑥,𝛽_𝑦, and 𝑘 evaluated at nodal points (𝜉_𝑖,𝜉_𝑗)∈Ω^±, respectively. In this case, the boundary and interface conditions remain unchanged.

The algorithm proposed before is to solve interface elliptic equations defined in a quadrilateral domain with straight line interface. One may generalize the algorithm to solve interface problems defined in a domain with curved boundaries using Gordon and Hall transformation [18, 19]. Here we briefly introduce the Gordon and Hall transformation on a simply connected domain Ω for reader’s convenience (see [18, 19] for more details.).

Let 𝐅 be a vector-valued function of two independent variables 𝑥 and 𝑦 over a domain 𝐒=[0,ℎ]×[0,ℎ] in the 𝑥𝑦-plane whose range is Ω in ℝ². We assume that 𝐅 is a continuous one-to-one transformation which maps 𝐒 onto the domain Ω such that 𝐅∶𝜕𝐒→𝜕Ω. Then we would like to construct a one-to-one function 𝐓∶𝐒→Ω which matches 𝐅 on the boundaries of 𝐒, so-called the boundary interpolant of 𝐹: 𝐓𝑦𝑦𝑦𝑦𝐓()()()()0,=𝐅0,,𝐓ℎ,=𝐅ℎ,,0≤𝑦≤ℎ,𝑥,0=𝐅𝑥,0,𝐓𝑥,ℎ=𝐅𝑥,ℎ,0≤𝑥≤ℎ. (33) Following the ideas given in [18, 19], one may choose the following simple transfinite bilinear Lagrange interpolant of 𝐅: 𝐓𝑦=𝑥𝑦𝑦𝑦𝑥𝑥,𝑥,𝑥,∶=1−ℎ𝐅𝑦+𝑥0,ℎ𝐅𝑦+𝑦ℎ,1−ℎ𝐅()+𝑦𝑥,0ℎ𝐅()−𝑥𝑥,ℎ1−ℎ𝑦1−ℎ𝐅()−𝑥0,01−ℎ𝑦ℎ𝐅()−𝑦0,ℎ1−ℎ𝑥ℎ𝐅()−𝑦ℎ,0ℎ𝑥ℎ𝐅(ℎ,ℎ). (34) Then we have the transformed function 𝑢(𝑥,𝑦)∶=𝑢(𝐓(𝑥,𝑦)) defined on the rectangular domain 𝐒 and the elliptic equation defined on the curved domain Ω is also transformed to an elliptic equation defined on the rectangular domain 𝐒 so that we can apply our proposed method to the transformed equation. The transformed equation can be represented as −𝑎^±Δ𝑢^±+𝐛^±⋅∇𝑢^±+𝑐^±𝑢^±=𝑓^±in𝐒^±, (35) whereas the coefficients 𝑎^±, 𝐛^±=(𝑏^±₁,𝑏^±₂), and 𝑐^± are functions in 𝑥 and 𝑦. Applying the spectral collocation method to the above equations we have the following matrix form: 𝐴₁=−𝐚⁺𝑆⁺_𝑥𝑥+𝑆⁺_𝑦𝑦+𝐛⁺₁𝑆_𝑥+𝐛⁺₂𝑆_𝑦+𝐜̂⁺𝐼𝟎𝟎−𝐚⁻𝑆⁻_𝑥𝑥+𝑆⁻_𝑦𝑦+𝐛⁻₁𝑆_𝑥+𝐛⁻₂𝑆_𝑦+𝐜̂⁻𝐼, (36) where 𝐚^±, 𝐛^±₁, 𝐛^±₂, and 𝐜̂^± are diagonal matrices of corresponding function values evaluated at nodal points 𝐱_𝑖∈𝐒^±, respectively. The jump and boundary conditions are similarly given as the previous 𝐴₂. For complex geometries we employ pseudo-spectral element method. First, we decompose the domain Ω into some nonoverlapping subdomains and then apply the spectral collocation method for each subdomain.

Finally, we give an example of the Gordon-Hall transformation for the domain given in Figure 2. The boundary interpolant is given by 𝐅()=𝑥−𝑥,0𝑥4()=𝑥,𝐅𝑥,ℎ𝑥4𝐅𝑦=23+2,0≤𝑥≤2,ℎ,21𝑦−2,𝐅𝑦=𝜋0,cos𝑦−12𝜋1+sin𝑦−12,0≤𝑦≤2, (37) and the explicit form of the transformation in (34) reduces to 𝐓𝑦=𝑥𝑥,1−2𝜋cos𝑦−12+𝑥𝑥1−2𝜋sin𝑦−123+1−43𝑥+4𝑦𝑥. (38)

3. Numerical Results

In this section, we first give an example defined on rectangle domain with straight line interface as in Figure 1. Afterward examples are defined on more complicated domain in which we have to use spectral element collocation method to solve them. Denote by 𝑢_𝑁 the discrete solution of interface problem and by 𝑒=𝑢−𝑢_𝑁, its errors. We present their 𝐿²_𝑤(Ω) and 𝐻¹_𝑤(Ω) discrete norm which is defined, respectively, as follows: ‖𝑒‖²_0,𝑤,𝑁=^𝑁_𝑖,𝑗=0𝑤_𝑖𝑗𝑒²𝐱_𝑖𝑗,‖𝑒‖²_1,𝑤,𝑁=‖∇𝑒‖²_0,𝑤,𝑁+‖𝑒‖²_0,𝑤,𝑁. (39)

Figure 1 

Interface domain.

Figure 2 

A curved domain.

Example 1 (straight line interface). Consider 𝛽𝑢_𝑥𝑥+𝑢_𝑦𝑦+𝑢=𝑓+𝑔𝛿_𝛼,inΩ=(0,𝐿)×(0,1),𝑢=0on𝜕Ω, (40) with the following exact solution: 𝑢⋅𝐶𝑥,𝑦=sin𝜋𝑦₁cos𝑥𝛾₋+𝐶₂sin𝑥𝛾₋+1,inΩ₋=()×(),𝐶0,𝛼0,1₃cos𝑥𝛾₊+𝐶₄sin𝑥𝛾₊+1,inΩ₊=(𝛼,𝐿)×(0,1), (41) where 𝛾₊=1√10,𝛾₋=1,√10𝛼=10𝜋6√,𝐿=10𝜋2,𝛽⁻=100,𝛽⁺=10, (42) and 𝐶_𝑖’s can be determined by the boundary and jump conditions. We have ,⋅𝑔=𝜈sin𝜋𝑦𝑓=sin𝜋𝑦1−𝜋²𝐶₁cos𝑥𝛾₋+𝐶₂sin𝑥𝛾₋+1,inΩ₋,1−𝜋²𝐶₃cos𝑥𝛾₊+𝐶₄sin𝑥𝛾₊+1,inΩ₊. (43)

Tables 1 and 2 show the spectral convergence in 𝐿²_𝜔(Ω)- and 𝐻¹_𝜔(Ω)-norm errors of Legendre and Chebyshev approximation for 𝜈=0 and 𝜈=10, respectively.

Table 1 

Error discretization for Example 1 with 𝜈=0.

Table 2 

Error discretization for Example 1 with 𝜈=10.

As we pointed out in Section 2, it is possible to use different polynomial order approximation on each subdomain. To do this for Example 1, we divide (0,𝛼) into (𝑁₁+1) intervals, (𝛼,𝐿) into (𝑁₂+1) intervals, and (0,1) into (𝑀+1) intervals corresponding to LGL or CGL points. For sake of simplicity we take 𝑁₁=𝑁₂=𝑁,𝑀=𝑁+2. The numerical result is given in Table 3 when 𝜈=0. We have similar results for 𝜈≠0. When the interface point 𝛼 is close to the endpoint, different polynomial order can be used to get accurate results.

Table 3 

Error discretization for Example 1 with 𝜈=0 and different polynomial order approximation.

In the following we present some numerical experiments in which the domain is Ω=(−11)² and interface Γ is the circle 𝑥²+𝑦²=1/4. To use spectral element collocation method, we decompose the domain Ω into twelve nonoverlapping subdomains as in Figure 3. It should be noted that there is no jump in solution 𝑢 and its normal derivative for the interface which is different from Γ; that is, [𝑢]|_Γ′=0,𝛽∇𝑢⋅𝐧_Γ′=0, (44) where Γ^′ is dotted line in Figure 3.

Figure 3 

Domain decomposition for the interface examples.

Example 2 (circular interface 1 [6]). Consider the elliptic equation of the form 𝑢_𝑥𝑥+𝑢_𝑦𝑦=2𝛿_Γ. (45) In this example, we have only singular source term. The exact solution is given by 𝑢=1𝑥,𝑦1,𝑟≤2()11+log2𝑟,𝑟>2, (46) where 𝑟=𝑥²+𝑦².

The interface conditions are given by [𝑢]|_Γ=0,[∇𝑢⋅𝐧]|_Γ=2. (47) Table 4 also shows the spectral convergence in 𝐿²_𝜔(Ω)- and 𝐻¹_𝜔(Ω)-norm errors of Legendre and Chebyshev approximation. The exact solution, approximated one, and its error for the Legendre case of this example are plotted in Figure 4.

Table 4 

Error discretization for Example 2.

(a)

(b)

(c)

Figure 4 

Exact solution, approximation solution, and its error for 𝑁=18.

Example 3 (circular interface 2 [6]). Consider the elliptic equation 𝛽𝑢_𝑥𝑥+𝛽𝑢_𝑦𝑦=𝑓+𝐶𝛿_Γ, (48) where 𝑓=8𝑥²+𝑦²𝛽=𝑥+4,𝑥,𝑦²+𝑦²1+1,𝑟≤21𝑏,𝑟>2. (49) The exact solution is given by 𝑟𝑢=²,1𝑟≤211−−18𝑏𝑏+𝑟⁴/2+𝑟²2+𝐶log(2𝑟)𝑏1𝑟>2. (50)

In this example, we have varying discontinuous coefficient and singular source term as well. Tables 4 and 5 show the numerical results for both Legendre and Chebyshev approximation for different values of 𝑏 and 𝐶. We observe that the errors in 𝐿²_𝜔(Ω)- and 𝐻¹_𝜔(Ω)-norm are decay exponentially regardless of discontinuous coefficients and singular source term. The exact solution, approximation, and its error for the Legendre case are plotted in Figures 5 and 6, for different values of 𝑏 and 𝐶.

Table 5 

Error discretization for 𝐶=0.1 and 𝑏=−3 for Example 3.

(a)

(b)

(c)

Figure 5 

Exact solution, approximation solution, and its error for 𝐶=0.1 and 𝑏=−3.

(a)

(b)

(c)

Figure 6 

Exact solution, approximation solution, and its error for 𝐶=0.1 and 𝑏=10.

Example 4 (circular interface 3 with jump condition [6]). Consider Laplace equation 𝑢_𝑥𝑥+𝑢_𝑦𝑦=0,inΩ,𝑢=0on𝜕Ω. (51) The exact solution is given 𝑢=𝑒𝑥,𝑦^𝑥𝑦1cos,𝑟≤2,10.𝑟>2. (52)

In this example we impose a jump in solution 𝑢 itself and also a jump in the normal derivative of 𝑢 as external constraints. Table 6 also shows the spectral convergence for the Legendre and Chebyshev pseudo-spectral element method although it has jump in solution and its normal derivative Table 7. The exact solution, approximated one, and its error for the Legendre case of this example are plotted in Figure 7.

Table 6 

Error discretization for 𝐶=0.1 and 𝑏=10 for Example 3.

Table 7 

Error discretization for 𝐶=0.1 and 𝑏=10 for Example 4.

(a)

(b)

(c)

Figure 7 

Exact solution, approximation solution, and its error for 𝑁=18.

4. Concluding Remarks

In this paper, we proposed the pseudo-spectral collocation method for second order elliptic interface problems with discontinuous coefficients and singular source term. First we derived interface conditions in variable and its normal derivative and two equations defined on each subdomain decomposing the whole domain into two subdomains based on the interface. Then we obtain very simple algorithm applying the spectral collocation method to each equation. It is shown that the proposed method can easily applied to interface problem defined on curved domain by using Gordon and Hall transformation. This means that the method can be easily extended to spectral element collocation method to solve interface problems defined on a complicated domain. The numerical experiments also showed that the method has the spectral convergence with high accuracy. Furthermore, the method can be adopted to solve Stokes or Navier-Stokes equations with interface which can be studied in a coming paper.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The second author was supported by Kyungpook National University fund 2013 and the third author was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (NRF2013R1A2009687).