Bernstein-Polynomials-Based Highly Accurate Methods for One-Dimensional Interface Problems

Abstract

A new numerical method based on Bernstein polynomials expansion is proposed for solving one-dimensional elliptic interface problems. Both Galerkin formulation and collocation formulation are constructed to determine the expansion coefficients. In Galerkin formulation, the flux jump condition can be imposed by the weak formulation naturally. In collocation formulation, the results obtained by B-polynomials expansion are compared with that obtained by Lagrange basis expansion. Numerical experiments show that B-polynomials expansion is superior to Lagrange expansion in both condition number and accuracy. Both methods can yield high accuracy even with small value of N.

1. Introduction

In this paper, we consider the following two-point boundary value problem: 𝛽𝑢_𝑥_𝑥+𝑢=𝑓+𝑣𝛿(𝑥−𝛼)+12𝛽⁻+𝛽⁺𝑤𝛿′(𝑥−𝛼),𝑥∈Ω=(𝑎,𝑏),(1.1) with boundary conditions 𝑢(𝑎)=𝑢_𝑎,𝑢(𝑏)=𝑢_𝑏,(1.2) where 𝑎<𝛼<𝑏,𝛿(𝑥) is the Dirac delta function and 𝛿′(𝑥) is the dipole source term. The function 𝛽(𝑥) is allowed to be discontinuous at 𝑥=𝛼. For simplicity here we assume that 𝑓(𝑥) is smooth function. Due to the presence of singular source, the solution 𝑢(𝑥) possesses the following jump relations [1]: [𝑢]=𝑤,𝛽𝑢_𝑥=𝑣,(1.3) where the jump notation [⋅] is defined as [𝑢]=lim_𝑥→𝛼+𝑢(𝑥)−lim_{𝑥→𝛼−}𝑢(𝑥)=𝑢⁺−𝑢⁻.(1.4) This problem is referred to as the interface problem and is used in various applications of physics, engineering, and biological sciences, see [2–4] and the references therein.

For interface problems, since sharp interfaces or local jumps exist within the solution domain, any high-order method, such as the spectral method, suffers from the Gibbs phenomenon [5]. Here we evoke items we care most for solving interface problems. In 1993, the immersed interface method (IIM) was proposed for interface problems [1]. It is a two-order finite difference method based on Cartesian grids by incorporating the jump relations into difference schemes. The authors have constructed the IIM-based ADI finite difference scheme for 2D nonlinear convection diffusion interface problems [6]. In [7], a high-order method was developed for both discontinuous coefficients and singular source based on finite element method. To enhance the accuracy, the modified Hermite polynomials are used for the basic functions in each element. The matched interface and boundary method was proposed in [8] for elliptic interface problems with discontinuous coefficients and singular source. In [9], the coupling interface method was developed for elliptic interface problems. Recently, Shin and Jung [5] developed the spectral collocation method for one-dimensional interface problems (1.1) with 𝑤=0, in which Lagrange basis functions are chosen as the trial functions.

In this work, we also consider this problem and take Bernstein polynomials basis as the trial functions. Bernstein polynomials are useful polynomials in computer-aided geometric design because of their excellent properties [10]. Recently there are some works that used Bernestein polynomials as basis for numerically solving differential equations [11, 12], integral equations [13, 14], and so on, but none of them is about interface problems. Our method is different from Shin and Jung’s in three aspects. (i)The trial functions are chosen as Bernstein polynomials due to its nice properties. These polynomials defined on an interval form a complete basis over the interval. Each of these polynomials is positive and their sum is unity. (ii)The Galerkin formulation is constructed for this problem. Since the Bernstein polynomials are algebraic polynomials, the mass matrix and stiff matrices can be evaluated efficiently. (iii)The B-polynomial-based collocation formulation is given collocated with both equidistant points and spectral points, respectively. Unlike Lagrange basis functions, the B-polynomial differential matrix can be computed easily and the condition number of resulting linear system is smaller than that of Lagrange basis, which is shown in the numerical examples.

The rest of this paper is organized as follows. In Section 2, the Bernstein polynomials on interval [𝑎,𝑏] and its derivatives are given. The Galerkin formulation and B-polynomial-based collocation formulation are shown in Section 3. Two numerical examples are given and analyzed in Section 4 before the conclusions are made in Section 5.

2. Bernstein Polynomials Basis

The general form of Bernstein polynomials of 𝑛th degree on interval [𝑎,𝑏] is defined as [12] 𝐵_𝑖,𝑛(𝑥)=⎛⎜⎜⎝𝑛𝑖⎞⎟⎟⎠(𝑥−𝑎)^𝑖(𝑏−𝑥)^𝑛−𝑖(𝑏−𝑎)^𝑛,𝑖=0,1,…,𝑛,(2.1) where the binomial coefficients are given by (^𝑛𝑖)=𝑛!/𝑖!(𝑛−𝑖)!, with 𝑛!=1×2×⋯×𝑛 for 𝑛≥1 and 0!=1. These (𝑛+1) B-polynomials of degree 𝑛 form a complete basis over the interval [𝑎,𝑏]. It is easy to show that any given polynomial of degree 𝑛 can be expressed in terms of linear combination of the B-polynomials basis functions. The B-polynomials can be generated by a recursive definition: 𝐵_𝑖,𝑛(𝑥)=𝑏−𝑥𝑏−𝑎𝐵^{𝑖,𝑛−1}(𝑥)+𝑥𝑏−𝑎𝐵^{𝑖−1,𝑛−1}(𝑥).(2.2) The derivatives of the 𝑛th degree B-polynomials are combinations of B-polynomials of degree 𝑛−1, which can be formulated as 𝐵^𝑖,𝑛=𝑛𝑏−𝑎𝐵^{𝑖−1,𝑛−1}𝐵_{𝑖,𝑛}−𝐵_{𝑖,𝑛−1},=𝑛(𝑛−1)(𝑏−𝑎)²𝐵_{𝑖,𝑛}𝐵_{𝑖−2,𝑛−2}−2𝐵_{𝑖−1,𝑛−2}+𝐵_{𝑖,𝑛−2},=𝑛(𝑛−1)(𝑛−2)(𝑏−𝑎)³𝐵_{𝑖−3,𝑛−3}−3𝐵_{𝑖−2,𝑛−3}+3𝐵_{𝑖−1,𝑛−3}−𝐵_{𝑖,𝑛−3},(2.3) where we set 𝐵_𝑖,𝑛=0 if 𝑖<0 or 𝑖>𝑛. In favor of these recursive relations, the differentiation matrix of Bernstein basis can be evaluated conveniently, while the computation of differentiation matrix of high-order Lagrange basis function may suffer from certain difficulties [15].

3. Numerical Methods

In this section, we give the numerical method for solving interface problem (1.1)–(1.3). It includes Galerkin formulation and collocation formulation. Without loss of generality, we assume that only one interface 𝑥=𝛼 exists in Ω and 𝛽 is piecewise constant in each subdomain. The multiple interface can be handled analogously. Thus the entire domain is divided by 𝛼 into two parts Ω₁=(𝑎,𝛼) and Ω₂=(𝛼,𝑏).

Suppose that the solutions in Ω₁ and Ω₂ are 𝑢₁ and 𝑢₂, respectively, then the problem (1.1)–(1.3) can be regarded as the following two smooth problems: 𝛽𝛽₁₂𝑢^1^+𝑢₁=𝑓,𝑥∈Ω₁,with𝑢₁(𝑎)=𝑢_𝑎,𝑢^2^+𝑢₂=𝑓,𝑥∈Ω₂,with𝑢₂(𝑏)=𝑢_𝑏,(3.1) together with jump conditions (1.3), which make (3.1) closed.

In each subdomain Ω_𝑘,𝑘=1,2, the solution 𝑢_𝑘(𝑥) can be approximated by 𝑈_𝑘(𝑥)=^𝑁𝑘_𝑖=0𝐶_𝑘,𝑖𝐵_{𝑖,𝑁𝑘}(𝑥),𝑥∈Ω_𝑘,𝑘=1,2,(3.2) where 𝑁_𝑘 is the degree of B-polynomials on each subdomain. According to the interpolation property of B-polynomials at two endpoints, we can easily get 𝐶_1,0=𝑢_𝑎 and 𝐶_2,𝑁2=𝑢_𝑏.

3.1. Galerkin Formulations

The variational formulation of (3.1) reads −−^𝛼𝑎_𝑏𝛼𝛽₁𝑢^1𝑣^1d𝑥+^𝛼𝑎𝑢₁𝑣₁d𝑥=^𝛼𝑎𝑓𝑣₁d𝑥−𝛽₁𝑢^1𝑣₁||^𝛼𝑎,𝛽₂𝑢^2𝑣^2d𝑥+^𝑏𝛼𝑢₂𝑣₂d𝑥=^𝑏𝛼𝑓𝑣₂d𝑥−𝛽₂𝑢^2𝑣₁||^𝑏𝛼,(3.3) where 𝑣_𝑘∈𝐻¹(Ω_𝑘) is arbitrary and satisfies certain boundary conditions.

Plugging (3.2) into (3.3) and replacing 𝑣_𝑘 with B-polynomials basis produce −^𝑁1_𝑖=0^𝛼𝑎𝛽₁𝐵^𝑖,𝑁₁−𝐵^𝑗,𝑁_1𝑁2d𝑥+^𝛼𝑎𝐵_𝑖,𝑁1𝐵_𝑗,𝑁1d𝑥𝐶_1,𝑖=^𝛼𝑎𝑓(𝑥)𝐵_𝑗,𝑁1d𝑥−𝛽₁𝑢^1𝐵_𝑗,𝑁1||^𝛼𝑎,_𝑖=0^𝑏𝛼𝛽₂𝐵^𝑖,𝑁₂𝐵^𝑚,𝑁₂d𝑥+^𝑏𝛼𝐵_𝑖,𝑁2𝐵_𝑚,𝑁2d𝑥𝐶_2,𝑖=^𝑏𝛼𝑓(𝑥)𝐵_𝑚,𝑁2𝑑𝑥−𝛽₂𝑢^2𝐵_𝑚,𝑁2||^𝑏𝛼,(3.4) with 𝑗=0,1,…,𝑁₁,𝑚=0,1,…,𝑁₂.

Observing that 𝐶_1,0=𝑢_𝑎, 𝐶_2,𝑁2=𝑢_𝑏 and 𝐵_{𝑗,𝑁𝑘}(𝑥), 𝑗=1,2,…,𝑁_𝑘−1, have zeros at the endpoints, (3.4) can be reformulated as ^𝑁1_𝑖=1𝐶_1,𝑖^𝛼𝑎−𝛽¹=𝐵^𝑖,𝑁₁𝐵^𝑗,𝑁₁+𝐵_𝑖,𝑁1𝐵_𝑗,𝑁1d𝑥_𝛼𝑎𝑓(𝑥)𝐵^𝑗,𝑁₁+𝛽₁𝑢_𝑎𝐵^0,𝑁₁𝐵^𝑗,𝑁₁−𝑢_𝑎𝐵_0,𝑁1𝐵_𝑗,𝑁1d𝑥,𝑗=1,…,𝑁¹−1,(3.5)_𝑁2−1_𝑖=0𝐶_2,𝑖^𝑏𝛼−𝛽²=𝐵^𝑖,𝑁₂𝐵^𝑚,𝑁₂+𝐵_𝑖,𝑁2𝐵_𝑚,𝑁2d𝑥_𝑏𝛼𝑓(𝑥)𝐵^𝑚,𝑁₂+𝛽₂𝑢_𝑏𝐵^𝑁_2,𝑁2𝐵^𝑚,𝑁₂−𝑢_𝑏𝐵_𝑁2,𝑁2𝐵_𝑚,𝑁2d𝑥,𝑚=1,…,𝑁²−1.(3.6) Let 𝑗=𝑁₁ in (3.5) and 𝑚=0 in (3.6), we get ^𝑁1_𝑖=1𝐶_1,𝑖^𝛼𝑎−𝛽¹=𝐵^𝑖,𝑁₁𝐵^𝑁_1,𝑁1+𝐵_𝑖,𝑁1𝐵_𝑁1,𝑁1d𝑥_𝛼𝑎𝑓(𝑥)𝐵^𝑁₁^,𝑁₁+𝛽₁𝑢_𝑎𝐵^0,𝑁₁𝐵^𝑁_1,𝑁1−𝑢_𝑎𝐵_0,𝑁1𝐵_𝑁1,𝑁1d𝑥−𝛽¹𝑢^−𝑥,(3.7)_𝑁2−1_𝑖=0𝐶_2,𝑖^𝑏𝛼−𝛽²=𝐵^𝑖,𝑁₂𝐵^0,𝑁₂+𝐵_𝑖,𝑁2𝐵_0,𝑁2d𝑥_𝑏𝛼𝑓(𝑥)𝐵^0,𝑁₂+𝛽₂𝑢_𝑏𝐵^𝑁_2,𝑁2𝐵^0,𝑁₂−𝑢_𝑏𝐵_𝑁2,𝑁2𝐵_0,𝑁2d𝑥+𝛽²𝑢^+𝑥.(3.8) Adding (3.7) and (3.8) together, combined with jump condition [𝛽𝑢_𝑥]=𝛽₂𝑢^+𝑥−𝛽₁𝑢^−𝑥=𝑣, yields ^𝑁1_𝑖=1𝐶_1,𝑖^𝛼𝑎−𝛽¹+𝐵^𝑖,𝑁₁𝐵^𝑁_1,𝑁1+𝐵_𝑖,𝑁1𝐵_𝑁1,𝑁1d𝑥_𝑁2−1_𝑖=0𝐶_2,𝑖^𝑏𝛼−𝛽²=𝐵^𝑖,𝑁₂𝐵^0,𝑁₂+𝐵_𝑖,𝑁2𝐵_0,𝑁2d𝑥_𝛼𝑎+𝑓(𝑥)𝐵^𝑁₁^,𝑁₁+𝛽₁𝑢_𝑎𝐵^0,𝑁₁𝐵^𝑁_1,𝑁1−𝑢_𝑎𝐵_0,𝑁1𝐵_𝑁1,𝑁1d𝑥_𝑏𝛼𝑓(𝑥)𝐵^0,𝑁₂+𝛽₂𝑢_𝑏𝐵^𝑁_2,𝑁2𝐵^0,𝑁₂−𝑢_𝑏𝐵_𝑁2,𝑁2𝐵_0,𝑁2d𝑥+𝑣.(3.9) Another jump condition [𝑢]=𝑤 implies −𝐶_1,𝑁1+𝐶_2,0=𝑤.(3.10) Equations (3.5), (3.6), (3.9), and (3.10) form a linear system with 𝑁₁+𝑁₂ equations and 𝑁₁+𝑁₂ unknowns: 𝐶_1,1,𝐶_1,2,…,𝐶_1,𝑁1,𝐶_2,0,𝐶_2,1,…,𝐶_2,𝑁2−1.(3.11)

3.2. Bernstein Collocation Methods

Substitution of (3.2) into (3.1) produces residuals 𝑅_𝑘=^𝑁𝑘_𝑖=0𝐶_𝑘,𝑖𝛽_𝑘𝐵^{𝑖,𝑁}_𝑘+^𝑁𝑘_𝑖=0𝐶_𝑘,𝑖𝐵_{𝑖,𝑁𝑘}−𝑓,𝑘=1,2.(3.12) In each interval [𝑎,𝛼] and [𝛼,𝑏], define the collocation points: 𝑋𝑋₁=𝑥_1,𝑖∣𝑎=𝑥_1,0<𝑥_1,1<⋯<𝑥_1,𝑁1−1<𝑥_1,𝑁1=𝛼₂,=𝑥_2,𝑖∣𝛼=𝑥_2,0<𝑥_2,1<⋯<𝑥_2,𝑁2−1<𝑥_2,𝑁2=𝑏.(3.13) Here the points in 𝑋_𝑘 can be equidistant points, Legendre-Gauss-Lobatto (L-G-L) points, or Chebyshev-Gauss-Lobatto (C-G-L) points.

Note that 𝐶_1,0 and 𝐶_2,𝑁2 in (3.2) are known. Collocation of (3.12) at points 𝑋_𝑘 yields ^𝑁1_𝑖=1𝐶=_1,𝑖𝑓𝛽₁𝑥𝐵^{𝑖,𝑁}_11,𝑗𝑥_1,𝑗+𝐵_𝑖,𝑁1𝑥_1,𝑗−𝑢_𝑎𝛽¹𝐵^0,𝑁₁𝑥_1,𝑗+𝐵_0,𝑁1𝑥_1,𝑗,𝑗=1,…,𝑁₁−1,_𝑁2−1_𝑖=0𝐶=𝑓_2,𝑖𝛽𝑥₂𝐵^{𝑖,𝑁}_22,𝑚𝑥_2,𝑚+𝐵_𝑖,𝑁2𝑥_2,𝑚−𝑢_𝑏𝛽²𝐵^𝑁_2,𝑁2𝑥_2,𝑚+𝐵_𝑁2,𝑁2𝑥_2,𝑚,𝑚=1,…,𝑁₂−1.(3.14) From (3.2) we can easily get 𝑈^𝑘(𝑥)=^𝑁𝑘_𝑖=0𝐶_𝑘,𝑖𝐵^𝑖,𝑁_𝑘(𝑥),𝑥∈Ω_𝑘,𝑘=1,2.(3.15) Jump condition [𝑢]=𝑤 and [𝛽𝑢_𝑥]=𝑣 imply −−𝐶_1,𝑁1+𝐶_2,0=𝑤,_𝑁1_𝑖=1𝐶_1,𝑖𝛽₁𝐵^𝑖,𝑁₁(𝛼)+^𝑁2−1_𝑖=0𝐶_2,𝑖𝛽₂𝐵^𝑖,𝑁₂(𝛼)=𝑣+𝛽₁𝑢_𝑎𝐵^0,𝑁₁(𝛼)−𝛽₂𝑢_𝑏𝐵^𝑁_2,𝑁2(𝛼).(3.16) Equations (3.14) and (3.16) produce the linear systems of 𝑁₁+𝑁₂ equations with 𝑁₁+𝑁₂ unknowns 𝐶_1,1,𝐶_1,2,…,𝐶_1,𝑁1,𝐶_2,0,𝐶_2,1,…,𝐶_2,𝑁2−1.(3.17)

4. Numerical Experiments

In this section, we give two examples to verify the accuracy of proposed numerical method. The first example is the one in which 𝑣≠0 and 𝑤=0, while in the second example both 𝑤 and 𝑣 are not zeros. We compare our results (B-polynomials-based) with that of Shin and Jung’s (Lagrange-polynomials-based) [5]. In all cases, we take 𝑁₁=𝑁₂=𝑁 and the resulting linear systems are almost block diagonal and solved by BiCGStab algorithm.

Example 4.1. Consider the interface problem given in [5] 𝛽𝑢_𝑥_𝑥+𝑢=𝑢1(𝑎+)=𝑣𝛿𝑢(𝑥(𝑏−)=𝛼)0,,𝑥∈(𝑎,𝑏),(4.1) with the following exact solution: 𝑢(𝑥)=⎧⎪⎪⎨⎪⎪⎩𝐶¹cos𝑥√𝛽₁+𝐶²sin𝑥𝐶√𝛽₁+1,𝑥∈(𝑎,𝛼),₃cos𝑥√𝛽₂+𝐶⁴sin𝑥√𝛽₂+1,𝑥∈(𝛼,𝑏),(4.2) where 𝑎=0, 𝑏=5, 𝛼=5/3, and 𝑣=10. The jump conditions read [𝑢]=0, [𝛽𝑢_𝑥]=𝑣. 𝐶_𝑖s can be determined by the boundary and jump conditions.
Both Galerkin formulation and collocation method are used to solve this problem numerically. The condition number (Cond) of the resulting linear system is given in each computation. Table 1 gives the convergence analysis of Galerkin formulation with different coefficients 𝛽₁ and 𝛽₂, in which the 𝐿₂ and 𝐻₁ norms are shown. It can be seen that the error reduces rapidly as the order of the B-polynomials increases.
The convergence analysis of collocation formulation is shown in Tables 2–5, in which the B-polynomials basis and Lagrange basis are compared. In each table, we use two types of collocation points to compare the results. It shows that the results of spectral collocation points (L-G-L or C-G-L points) are more accurate than the equidistant collocation points. Comparing Table 2 with Table 3, we can conclude that the condition number of linear systems derived from B-polynomials is much smaller than that derived from Lagrange polynomials. And the error from the former is smaller than the latter. Similar analysis can be got by comparing Table 4 with Table 5.

Table 1 

Convergence analysis of Example 4.1 by Galerkin formulation.

Table 2 

Convergence analysis of Example 4.1 by Bernstein collocation method (𝛽₁=100, 𝛽₂=10).

Table 3 

Convergence analysis of Example 4.1 by Lagrange collocation method (𝛽₁=100, 𝛽₂=10).

Table 4 

Convergence analysis of Example 4.1 by Bernstein collocation method (𝛽₁=10, 𝛽₂=100).

Table 5 

Convergence analysis of Example 4.1 by Lagrange collocation method (𝛽₁=10, 𝛽₂=100).

Example 4.2. In this example, the solution 𝑢 has nonzero jump across 𝛼. The following interface problem is considered 𝛽𝑢_𝑥_𝑥+𝑢=1+𝑣𝛿(𝑥−𝛼𝑢)(𝑎+)=21𝑢𝛽(₁𝑏+)=𝛽₂0.𝑤𝛿^(𝑥−𝛼),𝑥∈(𝑎,𝑏),(4.3) The jump conditions are [𝑢]=𝑤, [𝛽𝑢_𝑥]=𝑣. The exact solution is 𝑢(𝑥)=⎧⎪⎪⎨⎪⎪⎩𝐶¹cos𝑥√𝛽₁+𝐶²sin𝑥𝐶√𝛽₁+1,𝑥∈(𝑎,𝛼),₃cos𝑥√𝛽₂+𝐶⁴sin𝑥√𝛽₂+1,𝑥∈(𝛼,𝑏),(4.4) where 𝑎=0, 𝑏=5, 𝛼=5/3, 𝑤=10, 𝑣=10, and 𝐶_𝑖s can be determined by boundary and jump conditions.
The similar convergence analysis results can be obtained compared with Example 4.1. Since the nonzero jump 𝑤 just affects the right-hand side of the resulting linear systems, the condition numbers in Tables 6, 7, 8, 9, and 10 are unchanged compared with that in Tables 1–5, while the error in Tables 6–10 is much larger than that in Tables 1–5. The regularity of the solution can affect the accuracy of the numerical algorithm enormously.

Table 6 

Convergence analysis of Example 4.2 by Galerkin formulation.

Table 7 

Convergence analysis of Example 4.2 by Bernstein collocation method (𝛽₁=100, 𝛽₂=10).

Table 8 

Convergence analysis of Example 4.2 by Lagrange collocation method (𝛽₁=100, 𝛽₂=10).

Table 9 

Convergence analysis of Example 4.2 by Bernstein collocation method (𝛽₁=10, 𝛽₂=100).

Table 10 

Convergence analysis of Example 4.2 by Lagrange collocation method (𝛽₁=10, 𝛽₂=100).

5. Conclusions

In this paper, a new numerical method based on B-polynomials expansion is proposed for solving one-dimensional interface problems. We give two methods to evaluate the expansion coefficients, the Galerkin formulation, and the collocation formulation. Both methods can yield highly accurate results with small number of B-polynomials. In collocation method, the Lagrange polynomials are used to compare with B-polynomials. It is shown by numerical examples that B-polynomials are superior to Lagrange polynomials in both condition number and accuracy, especially when collocated with equidistant points. In theoretical aspect, since the B-polynomials basis is equivalent to power basis or Lagrange basis under certain invertible transformations, theoretical analysis of the proposed method may be done similarly, which is a part of our future research plan. The method can be extended to problems with multiple interfaces easily.

Acknowledgments

The authors thank the supports from the Institute of Applied Software of Central South University. The authors also thank the anonymous reviewers and the editor for their valuable comments. This work was supported partially by National Natural Science Foundation of China (no. 51174236), National Basic Research Program of China (no. 2011CB606306), and Hunan Provincial Innovation Foundation For Postgraduate (no. CX2011B080).