Mechanical Quadrature Method and Splitting Extrapolation for Solving Dirichlet Boundary Integral Equation of Helmholtz Equation on Polygons

Abstract

We study the numerical solution of Helmholtz equation with Dirichlet boundary condition. Based on the potential theory, the problem can be converted into a boundary integral equation. We propose the mechanical quadrature method (MQM) using specific quadrature rule to deal with weakly singular integrals. Denote by ℎ_𝑚 the mesh width of a curved edge Γ_𝑚(𝑚=1,…,𝑑) of polygons. Then, the multivariate asymptotic error expansion of MQM accompanied with 𝑂(ℎ³_𝑚) for all mesh widths ℎ_𝑚 is obtained. Hence, once discrete equations with coarse meshes are solved in parallel, the higher accuracy order of numerical approximations can be at least 𝑂(ℎ⁵_max) by splitting extrapolation algorithm (SEA). A numerical example is provided to support our theoretical analysis.

1. Introduction

We consider Helmholtz equation with Dirichlet boundary condition: Δ𝑢(𝑥)+𝛼²𝑢(𝑥)𝑢=0,𝑥∈Ω,_𝑚(𝑥)=𝑔_𝑚(𝑥),𝑥∈Γ_𝑚Γ=^𝑑_𝑚=1, (1) where Ω⊂ℜ² is a polygonal domain with the curved boundary Γ=∪^𝑑_𝑚=1Γ_𝑚, 𝑑>1, and the angels between Γ_𝑚 and Γ_𝑚+1 are in (0,2𝜋], and the function 𝑔_𝑚(𝑥)=𝑔(𝑥)|_Γ𝑚 is known on Γ_𝑚.

By the potential theory, the solutions of (1) can be represented as a single-layer potential 𝑢𝑦=_ΓΦ𝑦,𝑥𝑣(𝑥)𝑑𝑠_𝑥,𝑦=𝑦₁,𝑦₂∈Ω, (2) where 𝑥=(𝑥₁,𝑥₂) and Φ(𝑦,𝑥) is the foundation solution of Helmholtz equation, which is given by Φ=𝑖𝑦,𝑥4𝐻₀⁽¹⁾𝛼𝑥−𝑦. (3) And 𝐻₀⁽¹⁾=𝐽₀+𝑖𝑁₀ is the Hankel functions of order zero and of the first kind, where 𝐽₀(𝑧)=^∞_𝑛=0(−1)^𝑛(𝑛!)²𝑧2^2𝑛 (4) for the Bessel function of order zero and 𝑁₀(𝑧)=2𝜋𝑧ln2+𝛾𝐽₀(𝑧)+2𝜋^∞_𝑛=1^𝑛_𝑚=11𝑚(−1)^𝑛+1(𝑛!)²𝑧2^2𝑛 (5) for the Neumann function of order zero. And 𝛾=0.57721… is Euler constant.

In what follows, in order to analyze properties of the kernel, we decompose the kernel Φ𝑥,𝑦=𝑘₁𝑥,𝑦+𝑘₂𝑥,𝑦, (6) where 𝑘₁(𝑥,𝑦)=−1/2𝜋ln|𝑥−𝑦| is logarithmic singular function and 𝑘₂(𝑥,𝑦)=𝑖/4−1/2𝜋(ln(𝛼/2)+𝛾)+𝑂(|𝑥−𝑦|ln|𝑥−𝑦|) is a smooth function. 𝑣(𝑥) is the solution of the following equation: 𝑔𝑦=_ΓΦ𝑦,𝑥𝑣(𝑥)𝑑𝑠_𝑥,𝑦=𝑦₁,𝑦₂∈Γ. (7) Equation (7) is weakly singular BIE system of the first kind, whose solution exists and is unique as 𝐶_𝑇≠1 [1], where 𝐶_𝑇 is the logarithmic capacity (i.e., the transfinite diameter). Once 𝑣(𝑥) is solved from (7), the function 𝑢(𝑦)(𝑦∈Ω) can be calculated by (2).

Galerkin and collocation methods [2, 3] are used to solve (7). However, the discrete matrix is full and each element has to calculate the weakly singular integral for collocation methods or the double weakly singular integral for Galerkin methods, which implies that the work calculating discrete matrix is so large as greatly to exceed to solve the discrete equations. When the numerical methods are applied, the accuracy of numerical solutions is lower at singular points [4] and the corresponding numerical results become unreliable, because the condition numbers are very large.

In the paper, MQM is proposed to calculate weakly singular integrals by Sidi quadrature rules [5], which makes the calculation of the discrete matrix becomes very simple and straightforward without any singular integrals. The convergence theory of approximations is given by estimating eigenvalues of the discrete matrix and using Anselone’s collectively compact convergent theory [6], which shows that the method retains the optimal convergence order 𝑂(ℎ³_max) and possesses the optimal condition number 𝑂(ℎ⁻¹_min). Since MQM possesses the multivariate asymptotic expansion of errors, we can construct SEA to obtain the convergence order 𝑂(ℎ⁵_max). Once discrete equations on some coarse meshes are solved in parallel, the accuracy of numerical solutions can be greatly improved by SEA.

This paper is organized as follows. Section 2 includes the singularity analysis of the integral kernels and the solution. In Section 3, the MQM is described. In Section 4, we can obtain multiparameter asymptotic expansion of errors and SEA is described. In Section 5, a numerical example is provided to verify the theoretical results.

2. The Integral Kernels and the Solution of Singularity Analysis

Define boundary integral operators on Γ_𝑚: 𝐾_𝑞𝑚𝑣_𝑚𝑦1=−2𝜋_Γ𝑚Φ𝑥,𝑦𝑣_𝑚(𝑥)𝑑𝑠_𝑥,𝑦∈Γ_𝑞,𝑚,𝑞=1,…,𝑑. (8) Then (7) can be converted into a matrix operator equation 𝐾𝑉=𝐺, (9) where 𝐾=[𝐾_𝑞𝑚]^𝑑_𝑞,𝑚=1, 𝑉=(𝑣₁(𝑥),…,𝑣_𝑑(𝑥))^𝑇, and 𝐺=(𝑔₁(𝑦),…,𝑔_𝑑(𝑦))^𝑇.

Assume that Γ_𝑚 can be described by the parameter mapping: 𝑥_𝑚(𝑠)=𝑥_1𝑚(𝑠),𝑥_2𝑚(𝑠)∶[0,1]⟶Γ_𝑚,0≤𝑠≤1, (10) with |𝑥^′_𝑚(𝑠)|=[(𝑥^′_1𝑚(𝑠))²+(𝑥^′_2𝑚(𝑠))²]^1/2>0, 𝑚=1,…,𝑑. Using Sidi periodic transformation [7] 𝜓_𝑝(𝜏)=𝜗_𝑝(𝜏)𝜗_𝑝(1)∶[0,1]⟶[0,1],𝑝∈𝑁, (11) where 𝜗_𝑝(𝜏)=∫^𝜏₀(sin(𝜋𝜌))^𝑝𝑑𝜌, then operator (8) is also converted into integral operator on [0,1]: 𝐵_𝑞0𝑣_𝑞(𝑡)1=−2𝜋¹₀ln2𝑒^−1/2sin(𝜋(𝑡−𝜏))𝑣_𝑞(𝜏)𝑑𝜏,𝐵_0𝑞𝑣_𝑞(𝑡)1=−2𝜋¹₀ln𝑥_𝑞(𝑡)−𝑥_𝑞(𝜏)2𝑒^−1/2(𝜋(−sin𝑡−𝜏))𝜋𝑖2𝛼+ln2+𝛾𝑣_𝑞(𝜏)𝑑𝜏,𝐾_𝑞𝑚𝑣_𝑚(𝑡)=¹₀Φ𝑥_𝑞(𝑡),𝑥_𝑚(𝜏)𝑣_𝑚(𝜏)𝑑𝜏,𝑞≠𝑚, (12) where 𝑥_𝑞(𝑡)=𝑥_𝑞(𝜓_𝑝(𝑡)) and 𝑣_𝑚(𝜏)=𝑣_𝑚(𝜓_𝑝(𝜏))|𝑥^′_𝑚(𝜓_𝑝(𝜏))|𝜓^′_𝑝. Then (9) can be rewritten as 𝐵₀+𝐵₁+𝐾𝑉=𝐺, (13) where 𝐵₀=diag(𝐵₁₀,…,𝐵_𝑑0), 𝐵₁=diag(𝐵₀₁,…,𝐵_0𝑑), 𝐾=[𝐾_𝑞𝑚]^𝑑_𝑞,𝑚=1, 𝑉=(𝑣₁(𝜏),…,𝑣_𝑑(𝜏))^𝑇, 𝐺=(𝐺₁(𝑡),…,𝐺_𝑑(𝑡))^𝑇, and 𝐺_𝑚(𝑡)=𝑔_𝑚(𝜓_𝑝(𝑡)).

Let 𝑏_𝑞0(𝑡,𝜏), 𝑏_0𝑞(𝑡,𝜏), and 𝑘_𝑞𝑚(𝑡,𝜏) be the kernels of the integral operators 𝐵_𝑞0, 𝐵_0𝑞, and 𝐾_𝑞𝑚, respectively. Then, the following results hold.(1)𝑏_𝑞0(𝑡,𝜏) is a logarithmic singular function on 𝑡∈[0,1] and 𝜏∈[0,1].(2)𝑏_0𝑞(𝑡,𝜏) is a continuous function on 𝑡∈[0,1] and 𝜏∈[0,1].(3)For Γ_𝑞∩Γ_𝑚=∅(|𝑞−𝑚|≠1or𝑑−1), 𝑘_𝑞𝑚(𝑡,𝜏) is a continuous function on 𝑡∈[0,1] and 𝜏∈[0,1].(4)The solution 𝑣_𝑞(𝑡) is a smooth function under (11).

Lemma 1. If Γ_𝑞∩Γ_𝑚≠∅(|𝑞−𝑚|=1 or 𝑑−1) and 𝑝≥2, then 𝑘_𝑞𝑚(𝑡,𝜏)(=sin^𝑝(𝜋𝑡)𝑘_𝑞𝑚(𝑡,𝜏)) and (𝜕𝑛/𝜕𝑡^𝑛)𝑘(𝑡,𝜏) (𝑛=1,2) are continuous on [0,1]².

Proof. Without loss of generality, we assume that the origin (0,0)=Γ_𝑚∩Γ_𝑞 of coordinates is a vertex with an interior angle 𝜃_𝑚+1. By (6), 𝑘_𝑞𝑚 can be expressed by 𝑘_𝑞𝑚()1𝑡,𝜏=−𝑎2𝜋ln²₀(𝑡)+𝑎₁(𝜏)²−2𝑎₀(𝑡)𝑎₁(𝜏)cos𝜃_𝑚+1+𝑘₂𝑥_𝑞,𝑥_𝑚, (14) where 𝑎₀(𝑡)=|𝑥_𝑞(𝜓_𝑝(𝑡))| and 𝑎₁(𝜏)=|𝑥_𝑚(𝜓_𝑝(𝜏))|. We assume that 𝑎₀(0)=𝑎₁(0)=0; namely, 𝑘_𝑞𝑚 has the logarithmic singularity at points (0,0) and is continuously differentiable in [0,1]²⧵(0,0). Consider 𝑘_𝑞𝑚()1𝑡,𝜏=−𝑎2𝜋ln²₀(𝑡)+𝑎₁(𝜏)²−12𝜋×ln1−2𝑎₀(𝑡)𝑎₁(𝜏)cos𝜃_𝑚+1𝑎²₀(𝑡)+𝑎₁(𝜏)²+𝑘₂𝑥_𝑞,𝑥_𝑚. (15) We have 2𝑎₀(𝑡)𝑎₁(𝜏)cos𝜃_𝑚+1𝑎²₀(𝑡)+𝑎₁(𝜏)²≤cos𝜃_𝑚+1<1. (16) So the singularity of 𝑘_𝑞𝑚(𝑡,𝜏) occurs in (1/2𝜋)ln(𝑎²₀(𝑡)+𝑎₁(𝜏)²). In order to remove the singularity, we define 𝜍(𝑡,𝜏)=sin^𝑝(𝜋𝜏)𝜓_𝑝(1)ln𝑎²₀(𝑡)+𝑎₁(𝜏)². (17) Since 𝜓^′_𝑝 have 𝑝 order zero, we have the error estimate |𝜍(𝑡,𝜏)|=𝑜𝜀^𝑝|ln𝜀|⟶0,as𝜀⟶0, (18) which shows that 𝜍(𝑡,𝜏) is bounded and continuous on (𝐶²[0,1])² when 𝜀/2≤𝑡, 𝜏≤𝜀 for any 𝜀>0.
Similarly, we can prove (𝜕/𝜕𝑡)𝜍(𝑡,𝜏) and (𝜕²/𝜕𝑡²)𝜍(𝑡,𝜏) are continuous on (𝐶[0,1])², which can be obtained by 𝜕𝜍(≤𝜕𝑡𝑡,𝜏)sin^𝑝(𝜋𝜏)𝜓_𝑝(1)2𝑎₀(𝑡)𝑥^′_𝑚+1𝜓_𝑝(𝑡)𝜓^′_𝑝(𝑡)𝑎²₀(𝑡)+𝑎²₁(𝜏)=𝑂𝜀^𝑝𝜀,∀𝑡,𝜏∈2,𝜀, (19) and (𝜕²/𝜕𝑡²)𝜍(𝑡,𝜏)=𝑂(𝜀^𝑝−1). As above, (𝜕^𝑛/𝜕𝑡^𝑛)𝜍(𝑡,𝜏)(𝑛=0,1,…,𝑝) are continuous on (𝐶[0,1])². We can prove that (𝜕^𝑛/𝜕𝑡^𝑛)𝑘_𝑞𝑚(𝑡,𝜏)(𝑛=0,1,…,𝑝) are also continuous on (𝐶[0,1])². The proof of Lemma 1 is completed.

3. Mechanical Quadrature Method

Let ℎ_𝑚=1/𝑛_𝑚(𝑛_𝑚∈𝑁,𝑚=1,…,𝑑) be mesh widths and let 𝑡_𝑚𝑗=(𝑗−1/2)ℎ_𝑚(𝑗=1,…,𝑛_𝑚) be nodes. Since an integral operator 𝑊 is continuous, by the trapezoidal or the midpoint rule [8], we can construct Nyström approximate operator 𝑊^ℎ for the integral operator 𝑊, defined by 𝑊^ℎ𝑣_𝑚(𝑡)=ℎ_𝑚^𝑛𝑚_𝑗=1𝑤𝑡,𝑡_𝑚𝑗𝑣_𝑚𝑡_𝑚𝑗, (20) and the error 𝑊^ℎ𝑣_𝑚(𝑡)−𝑊𝑣_𝑚(𝑡)=𝑂ℎ_𝑚^2𝑙,𝑙∈𝑁. (21) Since 𝐵_𝑞0 has the singularities on [0,1], by Sidi quadrature formula [5], we get the following approximations 𝐵^ℎ_𝑞0 of integral operator 𝐵_𝑞0: 𝐵^ℎ_𝑞0𝑣_𝑚(𝑡)1=−ℎ2𝜋^𝑚𝑛𝑚_{𝑗=1,𝑡≠𝑡𝑚𝑗}ln2𝑒^−1/2sin𝜋𝑡−𝑡_𝑚𝑗×𝑣_𝑚𝑡_𝑚𝑗−ℎ_𝑚2𝜋ln2𝜋𝑒^−1/2ℎ_𝑚(2𝜋)𝑣_𝑚(𝑡), (22) and the error 𝐵^ℎ_𝑞0𝑣_𝑚(𝑡)−𝐵_𝑞0𝑣_𝑚(𝑡)2=−𝜋^2𝑙−1_𝜇=1𝜁^′−2𝜇2𝜇!𝑣_𝑚(𝑡)^2𝜇ℎ_𝑚^2𝜇+1+𝑂ℎ_𝑚^2𝑙, (23) where 𝜁^′(𝑡) is the derivative of the Riemann zeta function.

Then, we can get the approximate equation of (13): 𝐵^ℎ₀+𝐵^ℎ₁+𝐾^ℎ𝑉^ℎ=𝐺^ℎ, (24) where 𝑉^ℎ=(𝑉^ℎ₁,…,𝑉^ℎ_𝑑), 𝑉^ℎ_𝑚=(𝑣_𝑚(𝑡₁),…,𝑣_𝑚(𝑡_𝑛𝑚)), 𝐵^ℎ₀=diag(𝐵^ℎ₁₀,…,𝐵^ℎ_𝑑0), 𝐵^ℎ_𝑞0=[𝑏_𝑞0(𝑡_𝑗,𝜏_𝑖)]^𝑛𝑞_𝑗,𝑖=1, 𝐵^ℎ₁=diag(𝐵^ℎ₀₁,…,𝐵^ℎ_0𝑑), 𝐵_0𝑞=[𝑏_0𝑞(𝑡_𝑗,𝜏_𝑖)]^𝑛𝑞_𝑗,𝑖, 𝐾^ℎ=[𝐾^ℎ_𝑞𝑚]^𝑑_𝑞,𝑚=1, 𝐾^ℎ_𝑞𝑚=[𝑘_𝑞𝑚(𝑡_𝑗,𝜏_𝑖)]^{𝑛𝑞,𝑛𝑚}_𝑗,𝑖=1, 𝐺^ℎ=(𝐺^ℎ₁,…,𝐺^ℎ_𝑑)^𝑇, and 𝐺^ℎ_𝑞=(𝑔_𝑞(𝑡₁),…,𝑔_𝑞(𝑡_𝑛𝑞)). Obviously, (24) is a system of linear equations with 𝑛(=∑^𝑑_𝑚=1𝑛_𝑚) unknowns. Once 𝑉^ℎ is solved by (24), the solution 𝑢(𝑦)(𝑦∈Ω) can be computed by 𝑢^ℎ𝑦=^𝑑_𝑚=1ℎ_𝑚^𝑛𝑚_𝑗=1Φ𝑦,𝑥𝑡_𝑚𝑗𝑣_𝑚𝑡_𝑚𝑗. (25)𝐵^ℎ_𝑞0 is symmetric circular matrix and has the form of 𝐵^ℎ_𝑞0−ℎ=circulate^𝑚2𝜋ln2𝑒^−1/2ℎ_𝑚𝜋−ℎ,…,_𝑚2𝜋ln2𝑒^−1/2sin𝜋𝑛_𝑚−1ℎ_𝑚. (26)

Lemma 2 (see [9]). (1) There exists a positive 𝑐₁>0 so that the eigenvalues 𝜆_𝛽(𝛽=1,…,𝑛_𝑚) of 𝐵^ℎ_𝑞0 satisfy 𝑐₁>𝜆_𝛽>1/(2𝜋𝑛_𝑚). (2) The condition number of 𝐵^ℎ_𝑞0 is 𝑂(𝑛_𝑚). (3) 𝐵_𝑞0 is invertible, and (𝐵^ℎ_𝑞0)⁻¹ is uniformly bounded with the spectral norm ‖(𝐵^ℎ_𝑞0)⁻¹‖=𝑂(𝑛_𝑚).

Based on Lemma 2, we immediately get the following corollary.

Corollary 3. (1) 𝐵^ℎ₀ is invertible, and (𝐵^ℎ₀)⁻¹ is uniformly bounded with the spectral norm ‖(𝐵^ℎ₀)⁻¹‖=𝑂(𝑛₀). (2) The condition number of 𝐵^ℎ₀ is 𝑂(𝑛₀), where 𝑛₀=max^𝑑_𝑚=1𝑛_𝑚.

From Corollary 3, we know that (24) is equivalent to 𝐸^ℎ+𝐵^ℎ₀⁻¹𝐵^ℎ₁+𝐾^ℎ𝑉^ℎ=𝐵^ℎ₀⁻¹𝐺^ℎ, (27) where 𝐸^ℎ denotes the unit matrix.

Now we give the following definitions to discuss the approximate convergence in (27).

Define a subspace 𝐶₀[0,1]={𝑣(𝑡)∈𝐶[0,1]∶𝑣(𝑡)/sin²(𝜋𝑡)∈𝐶[0,1]} of the space 𝐶[0,1] with a norm ‖𝑣‖^∗=max_0≤𝑡≤1|𝑣(𝑡)/sin²(𝜋𝑡)|. Let 𝑆^ℎ_𝑚⊂𝐶₀[0,1] be a piecewise linear function subspace with base points {𝑡_𝑖}^𝑛𝑚−1_𝑖=0, and let 𝑒_𝑗(𝑡) 𝑗=0,…,𝑛_𝑚−1 be base function satisfying 𝑒_𝑗(𝑡_𝑖)=𝛿_𝑗𝑖. Define a prolongation operator 𝑃^ℎ𝑚∶ℜ^ℎ𝑚→𝑆^ℎ𝑚 satisfying 𝑃^ℎ𝑚𝜔=^𝑛𝑚−1_𝑗=0𝜔_𝑚𝑒_𝑗(𝑡),∀𝜔=𝜔_𝑚1,…,𝜔_𝑚𝑛𝑚∈ℜ^𝑛𝑚. (28) Define a restricted operator 𝑅^ℎ𝑚∶𝐶₀[0,1]→ℜ^𝑛𝑚 satisfying 𝑅^ℎ𝑚𝑣_𝑚=𝑣_𝑚𝑡_𝑚1,…,𝑣_𝑚𝑡_𝑚𝑛𝑚∈ℜ^𝑛𝑚,∀𝑣_𝑚∈𝐶₀[0,1]. (29)

Replacing (𝐵^ℎ₀)⁻¹(𝐵^ℎ₁+𝐾^ℎ)=𝑃^ℎ(𝐵^ℎ₀)⁻¹𝑅^ℎ(𝐵^ℎ₁+𝐾^ℎ), 𝑉^ℎ=𝑃^ℎ𝑉^ℎ, and 𝐺=𝑃^ℎ𝐺^ℎ, we construct an operator equation 𝐸^ℎ+𝐵^ℎ₀⁻¹𝐵^ℎ₁+𝐾^ℎ𝑉^ℎ=𝐵^ℎ₀⁻¹𝐺^ℎ, (30) where 𝑃^ℎ=diag(𝑃₁,…,𝑃_𝑑) and 𝑅^ℎ=diag(𝑅₁,…,𝑅_𝑑). Obviously, if 𝑉^ℎ is the solution of (30), then 𝑅^ℎ𝑉^ℎ must be the solution of (27); conversely, if 𝑉^ℎ is the solution of (27), then 𝑃^ℎ𝑉^ℎ must be the solution of (30). In order to prove the convergence of MQM, we give the following lemma.

Lemma 4 (see [9]). The operator sequence {𝑃_𝑞(𝐵^ℎ_𝑞0)⁻¹𝑅_𝑞𝐵_𝑞0∶𝐶²[0,1]→𝐶[0,1]} is uniformly bounded and convergent to the embedding operator 𝐼.

Corollary 5. Let the Nyström approximation 𝑊^ℎ be defined by (20). (1) For Γ_𝑞∩Γ_𝑚=∅, one has 𝑃_𝑚𝐵_𝑚0⁻¹𝑅_𝑚𝑊−−→^ℎ𝑐.𝑐𝐵_𝑚0⁻¹𝑊,𝑖𝑛𝐶[0,1]⟶𝐶[0,1]. (31) (2) For Γ_𝑞∩Γ_𝑚≠∅, one has 𝑃_𝑚𝐵_𝑚0⁻¹𝑅_𝑚𝐾−−→^ℎ_𝑞𝑚^𝑐.𝑐𝐵_𝑚0⁻¹𝐾_𝑞𝑚,𝑖𝑛𝐶[0,1]⟶𝐶[0,1], (32) where 𝐾^ℎ_𝑞𝑚 is the Nyström approximation of integral operator 𝐾_𝑞𝑚 with the kernel 𝑘_𝑞𝑚(𝑡,𝜏)=sin^𝑝(𝜋𝑡)𝑘_𝑞𝑚(𝑡,𝜏) and −−→^𝑐.𝑐 denotes the collectively compact convergence.

Theorem 6. Let Γ=∪^𝑑_𝑚=1Γ_𝑚 with 𝐶_𝑇≠1 and Γ_𝑚(𝑚=1,…,𝑑) be smooth curve. The operator sequence {𝑃^ℎ(𝐵^ℎ₀)⁻¹𝑅^ℎ(𝐵^ℎ₁+𝐾^ℎ)} is collectively convergent to (𝐵₀)⁻¹(𝐵₁+𝐾) in 𝑉=(𝐶₀[0,1])^𝑑: 𝑃^ℎ𝐵^ℎ₀⁻¹𝑅^ℎ𝐵^ℎ₁+𝐾−−→^ℎ𝑐.𝑐𝐵₀⁻¹𝐵₁+𝐾. (33)

Proof. Let Θ={𝜐∶‖𝜐‖≤1,𝜐∈𝑉} be a unit ball. 𝐻={ℎ¹,ℎ²,…} is the grid step sequence, where ℎ^𝑛={ℎ₁^(𝑛),ℎ₂^(𝑛),…,ℎ_𝑑^(𝑛)} with 𝑛→∞, max_{1≤𝑚≤𝑑}ℎ_𝑚^(𝑛)→0. We arbitrarily take a sequence {𝑍_ℎ}_ℎ∈𝐻 in space Θ, where 𝑍_ℎ=(𝑍_ℎ1,𝑍_ℎ2,…,𝑍_ℎ𝑑)^𝑇 with ‖𝑍_ℎ𝑚/sin²(𝜋𝑡)‖_∞≤1, 𝑚=1,…,𝑑. We firstly consider the first component of 𝑃^ℎ(𝐵^ℎ₀)⁻¹𝑅^ℎ(𝐵^ℎ₁+𝐾^ℎ)𝑍_ℎ: ^𝑑_𝑚=1𝑃₁𝐵^ℎ₁₀⁻¹𝑅₁𝐵^ℎ₀₁+𝐾^ℎ_1𝑚𝑍_ℎ𝑚. (34) For Γ_𝑚∩Γ₁=∅, by Corollary 5, 𝑃₁(𝐵^ℎ₁₀)⁻¹𝑅₁(𝐵^ℎ₀₁+𝐾^ℎ_1𝑚)−−→^𝑐.𝑐𝐵⁻¹₁₀(𝐵₀₁+𝐾_1𝑚) in 𝐶[0,1], and there exists a convergent subsequence in {𝑃₁(𝐵^ℎ₁₀)⁻¹𝑅₁(𝐵^ℎ₀₁+𝐾^ℎ_1𝑚)𝑍_ℎ𝑚}. For Γ_𝑚∩Γ₁≠∅, we have 𝑃₁𝐵^ℎ₁₀⁻¹𝑅₁𝐾^ℎ_1𝑚𝑍_ℎ𝑚0,0=𝑃₁𝐵^ℎ₁₀⁻¹𝑅₁𝐾^ℎ_1𝑚𝑍_ℎ𝑚sin²(𝜋𝑡)_0,0≤𝑃₁𝐵^ℎ₁₀⁻¹𝑅₁𝐵_100,2𝐵⁻¹₁₀𝐾^ℎ_1𝑚2,0𝑍_ℎ𝑚sin²(𝜋𝑡)_∞. (35) By Lemma 4, there exists a convergent subsequence in {𝑃₁(𝐵^ℎ₁₀)⁻¹𝑅₁(𝐵^ℎ₀₁+𝐾^ℎ_1𝑚)𝑍_ℎ𝑚}. Based on the above two cases it is proved that there exists an infinite subsequence 𝐻⁽¹⁾⊂𝐻 such that the first component converges. Similarly, it can be concluded that there exists an infinite subsequence 𝐻^(𝑑)⊂𝐻^(𝑑−1)⊂⋯⊂𝐻⁽¹⁾⊂𝐻 such that 𝑃^ℎ(𝐵^ℎ₀)⁻¹𝑅^ℎ(𝐵^ℎ₁+𝐾^ℎ)𝑍_ℎ converges. Therefore, 𝑃^ℎ(𝐵^ℎ₀)⁻¹𝑅^ℎ(𝐵^ℎ₁+𝐾^ℎ) is collectively compact convergent sequence, and 𝑃^ℎ(𝐵^ℎ₀)⁻¹𝑅^ℎ(𝐵^ℎ₁+𝐾^ℎ)−→^𝑃(𝐵₀)⁻¹(𝐵₁+𝐾), where −→^𝑃 shows the point convergence. We complete the proof.

For the stability of MQM, we have the following corollary.

Corollary 7. Let Γ=∪^𝑑_𝑚=1Γ_𝑚 with 𝐶_Γ≠1, let Γ_𝑚 (𝑚=1,…,𝑑) be smooth curve, let 𝐵^ℎ₀, 𝐵^ℎ₁, and 𝐾^ℎ be the discrete matrices defined by (20) and (22), respectively, and let 𝜆_𝑖(𝑖=1,…,𝑛) be the eigenvalues of discrete matrix 𝑀^ℎ=𝐵^ℎ₀+𝐵^ℎ₁+𝐾^ℎ. Then, there exists the bound of condition number Cond𝑀^ℎ=max_{1≤𝑖≤𝑛}𝜆_𝑖𝑀^ℎmax_{1≤𝑖≤𝑛}𝜆_𝑖𝑀^ℎ=𝑂ℎ⁻¹_min, (36) where ℎ_min=min^𝑑_𝑚=1ℎ_𝑚, ℎ_𝑚=1/𝑛_𝑚, is the mesh step size of a curved edge Γ_𝑚.

4. Multiparameter Asymptotic Expansion of Errors and SEA

In this section, we derive the multivariate asymptotic expansion of solution errors and describe SEA. We first provide the main result.

Theorem 8. Let Γ=∪^𝑑_𝑚=1Γ_𝑚 with 𝐶_Γ≠1. There exists a vector function Φ=(𝜙₁,…,𝜙_𝑑)^𝑇 independent of ℎ=(ℎ₁,…,ℎ_𝑑) such that the following multiparameter asymptotic expansion holds at nodes: 𝑉^ℎ−𝑉=diagℎ³₁,…,ℎ^𝑑_𝑑Φ+𝑜ℎ³₀,ℎ₀=max_{1≤𝑚≤𝑑}ℎ_𝑚, (37) where 𝑉^ℎ∈𝑆^𝑑, 𝑉∈𝐶³[0,1].

Proof. By (20) and (22), there exists the asymptotic expansion 𝐺^ℎ−𝐺=diagℎ³₁,…,ℎ^𝑑_𝑑𝑃^ℎ𝑅^ℎ𝜔+𝑜ℎ³₀, (38) where ℎ=max^𝑑_𝑚=1ℎ_𝑚 and 𝜔=(𝜔₁,…,𝜔)^𝑇 with 𝜔_𝑚=−𝜁^′(−2)(𝐺_𝑚(𝑡))^′′/𝜋.
Using (13), (21), (23), and (24), we obtain 𝐵^ℎ₀+𝐵^ℎ₁+𝐾^ℎ𝑅^ℎ𝑉^ℎ−𝑉=𝑉^ℎ−𝑃^ℎ𝐵^ℎ₀+𝐵^ℎ₁+𝐾^ℎ𝑅^ℎ𝑉=𝑉^ℎ−𝐵₀+𝐵₁+𝐾𝑅^ℎ𝑉−diagℎ³₁,…,ℎ^𝑑_𝑑𝑃^ℎ𝑅^ℎ𝛾+𝑜ℎ³₀=𝐺^ℎ−𝐺+diagℎ³₁,…,ℎ^𝑑_𝑑𝑃^ℎ𝑅^ℎ𝛾+𝑜ℎ³₀=diagℎ³₁,…,ℎ^𝑑_𝑑𝑃^ℎ𝑅^ℎ𝜂+𝑜ℎ³₀, (39) where 𝛾=(𝛾₁,…,𝛾_𝑑)^𝑇 with 𝛾_𝑚=−𝜁^′(−2)(𝑥^′(𝑡)𝜓^′(𝑡)𝑣(𝑡))^′′/𝜋 and 𝜂=(𝜂₁,…,𝜂_𝑑)^𝑇 with 𝜂_𝑚=𝜔_𝑚+𝛾_𝑚. From Corollary 3, we can obtain 𝐸+𝑃^ℎ𝐵^ℎ₀⁻¹𝑅^ℎ𝐵^ℎ₁+𝐾^ℎ𝑉^ℎ−𝑉=diagℎ³₁,…,ℎ^𝑑_𝑑𝐵^ℎ₀⁻¹𝑃^ℎ𝑅^ℎ𝜂+𝑜ℎ³₀. (40) Constructing the auxiliary equation 𝐸+𝐵₀⁻¹𝐵₁+𝐾Φ=𝐵₀⁻¹𝜂 (41) and its approximate equation 𝐸+𝑃^ℎ𝐵^ℎ₀⁻¹𝑅^ℎ𝐵^ℎ₁+𝐾^ℎΦ^ℎ=𝐵^ℎ₀⁻¹𝑃^ℎ𝑅^ℎ𝜂 (42) and substituting (42) into (40), we obtain 𝐸+𝑃^ℎ𝐵^ℎ₀⁻¹𝑅^ℎ𝐵^ℎ₁+𝐾^ℎ×𝑉^ℎ−𝑉−diagℎ³₁,…,ℎ^𝑑_𝑑Φ^ℎ=𝑜ℎ³₀. (43) Since (𝐸+𝑃^ℎ(𝐵^ℎ₀)⁻¹𝑅^ℎ(𝐵^ℎ₁+𝐾^ℎ)) is uniformly bounded from Theorem 6, we get 𝑉^ℎ−𝑉−diagℎ³₁,…,ℎ^𝑑_𝑑Φ^ℎ=𝑜ℎ³₀. (44) Replacing Φ^ℎ with Φ, we can complete the proof.

The multiparameter asymptotic expansion (37) means that SEA can be applied to solve (7); that is, higher order accuracy 𝑜(ℎ³₀) at coarse grid points can be obtained by solving some discrete equations in parallel. The process of SEA is as follows [10].

Step 1. Take ℎ⁽⁰⁾=(ℎ₁⁽⁰⁾,…,ℎ_𝑑⁽⁰⁾) and ℎ^(𝑚)=(ℎ₁⁽⁰⁾,…,ℎ_𝑚⁽⁰⁾/2,…,ℎ_𝑑⁽⁰⁾), and solve (24) under mesh parameters ℎ^(𝑚) in parallel to get the numerical solutions 𝑉^ℎ(0)(𝑡_𝑚𝑗) and 𝑉^ℎ(𝑚)(𝑡_𝑚𝑗), 𝑚=1,…,𝑑, 𝑗=1,…,𝑛_𝑚.

Step 2. Compute 𝑢^ℎ(0)(𝑦) and 𝑢^ℎ(𝑚)(𝑦)(𝑦∈Ω), by (25), 𝑉^ℎ(0)(𝑡_𝑚𝑗), and 𝑉^ℎ(𝑚)(𝑡_𝑚𝑗).

Step 3. Compute the extrapolation on the coarse grids as follows: 𝑢^∗𝑦=87^𝑑_𝑚=1𝑢^ℎ(𝑚)𝑦−7𝑑−8𝑢^ℎ(0)𝑦. (45)

Step 4. Compute a posteriori error estimate on the coarse grids as follows: 𝑢𝑦−1𝑑^𝑑_𝑚=1𝑢^ℎ(𝑚)𝑦≤𝑢𝑦−^𝑑_𝑚=1𝑢^ℎ(𝑚)𝑦−7𝑑−8𝑢^ℎ(0)𝑦+871𝑑−1𝑑^𝑑_𝑚=1𝑢^ℎ(𝑚)𝑦−𝑢^ℎ(0)𝑦≤871𝑑−1𝑑^𝑑_𝑚=1𝑢^ℎ(𝑚)𝑦−𝑢^ℎ(0)𝑦+𝑜ℎ³₀. (46) In the actual calculation process, a posteriori error estimate is immediately used to verify the calculation accuracy.

5. Numerical Example

In this section, we carry out a numerical example for the Helmholtz equation by MQM and SEA, in order to verify the error and stability analysis in the previous sections. Let Re(𝑒_𝑝) be the real part of 𝑒_𝑃=|𝑢^ℎ(𝑃)−𝑢(𝑃)| and let Im(𝑒_𝑃) be the imaginary part of 𝑒_𝑃=|𝑢^ℎ(𝑃)−𝑢(𝑃)|, where 𝑃 denotes a point. Post-error and SEA-error denote the a posteriori error and the error after SEA once, respectively.

Example 1. Consider Helmholtz equation with 𝛼=√2 on a plate domain Ω. We describe the boundary Γ=∪⁴_𝑚Γ_𝑚 with Γ₁={(𝑥₁,𝑥₂)∣0≤𝑥₁≤1,𝑥₂=0}, Γ₂={(𝑥₁,𝑥₂)∣0≤𝑥₂≤1,𝑥₁=1}, Γ₃={(𝑥₁,𝑥₂)∣0≤𝑥₁≤1,𝑥₂=1}, and Γ₄={(𝑥₁,𝑥₂)∣0≤𝑥₂≤1,𝑥₁=0}. Dirichlet boundary conditions corresponding to the analytical solution 𝑒^{𝑖×(𝑥1+𝑥2)} were applied to the boundary. We compute the numerical solution 𝑢^ℎ with 𝑃=𝐴(0.5,0.5) and 𝑃=𝐵(0.6,0.6) using 𝜑₃(𝑡). The numerical results are listed in Tables 1 and 2. From Table 1, we can know that the convergence rates of 𝑢^ℎ are 𝑂(ℎ³_max) for MQM and that the convergence rates of 𝑢^ℎ are at least 𝑂(ℎ⁵_max) for SEA. From Table 2, we can see Cond|₍₂^𝑘+1_,2^𝑘+1_,2^𝑘+1_,2^𝑘+1₎/Cond|₍₂^𝑘_,2^𝑘_,2^𝑘_,2^𝑘₎≈2 (𝑘=3,4,5) to indicate Corollary 7. It verifies the stability of convergent theory for MQM. Those results coincide with the theoretical analysis made.

Table 1 

The numerical results.

Table 2 

The condition number.

6. Concluding Remarks

To close this paper, let us make a few concluding remarks.(1)Evaluation on entries of discrete matrices is very simple and straightforward without any singular integrals by MQM.(2)The numerical experiments show that MQM retains the optimal convergence order 𝑂(ℎ³_max) and possesses the optimal condition number 𝑂(ℎ⁻¹_min) which shows MQM owns the excellent stability. The approximate solutions accuracy order is at least 𝑂(ℎ⁵_max) after splitting extrapolation once, which is a great improvement in accuracy.

Conflict of Interests

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