A High Order Numerical Method for Scattering from Locally Perturbed Periodic Surfaces
Ruming Zhang








A High Order Numerical Method for Scattering from Locally Perturbed Periodic Surfaces
    Ruming ZhangCenter for Industrial Mathematics, University of Bremen
; 
    
=====================================================================================




In this paper, we will introduce a high order numerical method to solve the scattering problems with non-periodic incident fields and (locally perturbed) periodic surfaces. For the problems we are considering,  the classical methods to treat quasi-periodic scattering problems no longer work, while a Bloch transform based numerical method was proposed in <cit.>. This numerical method, on one hand, is able to solve this kind of problems convergently; on the other hand, it takes up a lot of time and memory during the computation. The motivation of this paper is to improve this numerical method, from the regularity results of the Bloch transform of the total field, which have been studied in <cit.>. As the set of the singularities of the total field is discrete in , and finite in one periodic cell, we are able to improve the numerical method by designing a proper integration contour with special conditions at the singularities. With a good choice of the transformation, we can prove that the new numerical method could possess a super algebraic convergence rate. This new method improves the efficient significantly. At the end of this paper, several numerical results will be provided to show the fast convergence of the new method. The method also provides a possibility to solve more complicated problems efficiently, e.g., three dimensional problems, or electromagnetic scattering problems.




§ INTRODUCTION

  
In this paper, we will propose an efficient numerical method of solving scattering problems from locally perturbed periodic surfaces with non-periodic incident fields. For the special case that a quasi-periodic incident field scattered by a smooth enough periodic surface,  a classical way is to reduce the problem into one periodic cell, see <cit.>. Then the reduced problem could be solved by numerical methods, such as the integral equation method (see <cit.>), and the finite element method (see <cit.>). However, when the incident field or the surface is no longer (quasi-)periodic, the classical method fails,  new numerical methods are required for these more difficult problems. In <cit.>, a Floquet-Bloch based numerical method was applied to solve this kind of problems. In this paper, the method is improved by reducing the computational  complexity significantly.


The Floquet-Bloch based methods were used in <cit.> for the numerical solutions in (locally perturbed) periodic scattering problems with non-periodic incident fields.  For theoretical results, see <cit.> and <cit.>. There is also another way to solve the locally perturbed periodic problems,  by studying the operators defined on the transfparent edge of a periodic cell. In <cit.>, a method that approximates the Dirichlet-to-Neumann map by solving an operator equation is applied, while in 
<cit.>, the problems were solved by approximating the Scattered-to-Scattered map by a so-called doubling recursive procedure. 


 The method in this paper follows the results in <cit.>. To solve the scattering problems with a locally perturbed surface, firstly we have to transform the problem into one defined on a periodic domain. With the Bloch transform, the problem is then transformed into an equivalent  coupled system of quasi-periodic problems, which is a three dimensional problem defined in a bounded domain. The error comes from two parts in the algorithm – the finite element method and the numerical inverse Bloch transform. As the error brought by the finite element method could be estimated by classical methods, in this paper, we will improve the numerical scheme by the study of the numerical inverse Bloch transform. To study this, a detailed behaviour of the Bloch transformed solutions is required, which has been shown in <cit.>.
 
Denote by _Ω the Bloch transform defined on the periodic domain Ω, then the total field u is transformed into w(α,x)=(_Ω u)(α,x). From <cit.>, the singularities in α locate in a discrete set in , with a square root like singularity in the neighbourhoods.  As the inverse Bloch transform has the representation of the integral on  a bounded interval, in this paper, we will pick up a new curve and let the inverse Bloch transform be defined on the curve such that the newly defined integrand satisfies some smoothness and periodic conditions. This technique to treat integrations with singularities is quite classical, and it has been described in detail and  applied to the numerical solution of the scattering problems from bounded obstacles with corners in Chap 3.5, <cit.> (pp. 83-85). From the well known result of the convergent rule for periodic functions, see <cit.>, we can prove that with the new integral curve, the numerical method will converge much faster, depending on the choice of the new integration curve. 

The high order numerical method is available when the incident fields satisfy certain conditions. As far as we know, both the point sources and the Herglotz wave functions satisfy the conditions, thus the scattering problems with these incident fields  could be solved numerically by the high order method introduced in this paper. Numerical examples for these cases will be presented at the end of this paper. However, the plane waves do not satisfy the conditions, so the numerical method fails in this case. 



The rest of the paper is organized as follows. In Section 2, we will review some results for Bloch transform and the locally perturbed periodic problems. In Section 3, we recall the singularity results for the Bloch transform of the total field. In Section 4, we will introduce the inverse Bloch transform defined on a new contour. In Section 5, we will study the finite element method with the new inverse Bloch transform. In Section 6, some numerical examples will be shown. 








§ BLOCH TRANSFORM AND THE LOCALLY PERTURBED PERIODIC PROBLEMS






 §.§ Introduction to the Scattering problems

Suppose ζ is a Liptschiz periodic function with period Λ, and the periodic surface is defined by

    Γ:={(x_1,ζ(x_1)): x_1∈}.
 
Let ζ_p be a function such that supp(ζ_p-ζ) is compact in . For simplicity, assume that supp(ζ_p-ζ)⊂(-Λ/2,Λ/2], then define the locally perturbed surface by

    Γ_p:={(x_1,ζ_p(x_1)): x_1∈}.

Define Λ^*=2π/Λ, and

    =(-Λ/2,Λ/2],  =(-Λ^*/2,Λ^*/2]=(-π/Λ,π/Λ].

We can define the space above the surfaces Γ and Γ_p by

    Ω:={(x_1,x_2): x_1∈, x_2>ζ(x_1)}
          Ω_p:={(x_1,x_2): x_1∈, x_2>ζ_p(x_1)}.

For simplicity, we assume that Γ and Γ_p are both above the straight line ×{0}, and H be a positive number such that H>max{sup{ζ},sup{ζ_p}}. Define the line Γ_H=×{H} and the spaces with finite height

    Ω^H:={(x_1,x_2): x_1∈, ζ(x_1)<x_2<H}
          Ω_p^H:={(x_1,x_2): x_1∈, ζ_p(x_1)<x_2<H}.

We can also define the curves and  domains in one periodic cell:

    Γ^Λ=Γ∩,  Γ^Λ_H=Γ_H∩;
          Ω^Λ=Ω∩×,  Ω^Λ_H=Ω^H∩×.




Given an incident field u^i, the scattering problem is to find a total field u such that it satisfies the Helmholtz equation in the space Ω_p 

    Δ u+ k^2 u=0   in Ω_p

and the scattered field u^s:=u-u^i propagates upward, i.e.,

    u^s(x)=1/2π∫_ e^ı x_1 ξ+ı√(k^2-|ξ|^2)(x_2-H)u^s(ξ,H)ξ̣   x_2>H,

where u^s(ξ,H) is the Fourier transform of u^s(·,H). The upward propagating condition is equivalent to 

    ∂ u/∂ x_2=T^+ u|_Γ_H+(∂ u^i/∂ x_2-T^+[u^i|_Γ_H])   on Γ_H,

where T^+ is the Dirichlet-to-Neumann map (see <cit.>) defined by

    T^+ ϕ(x_1)=ı/2π∫_√(k^2-|ξ|^2) e^ı x_1 ξ+ı√(k^2-|ξ|^2)(x_2-H)ϕ(ξ)ξ̣

that is bounded from H^1/2_r(Γ_H) to H^-1/2_r(Γ_H) for all |r|<1. Define the total field u=u^i+u^s, then we can formulate the variational problem, i.e., to find u∈H^1_r(Ω^H_p) such that for any ϕ∈ H^1(Ω^H_p) with a compact support such that

    ∫_Ω^H_p[∇ u·∇ϕ-k^2uϕ]x̣-∫_Γ_HT^+(u|_Γ_H) ϕṣ=∫_Γ_H[∂ u^i/∂ x_2-T^+(u^i|_Γ_H)]ϕṣ,
 
where H_r^s(Ω^H_p) is a weighted space defined by

    H_r^s(Ω^H_p):={ϕ∈𝒟'(Ω^H_p):  (1+|x_1|^2)^r/2ϕ(x)∈ H^s(Ω^H_p)},

and H_r^s(Ω^H_p) is a subspace defined by

    H_r^s(Ω^H_p)={ϕ∈ H_r^s(Ω^H_p): ϕ|_Γ_p=0}.


In <cit.>, the unique solvability of the variational form was proved.

For any |r|<1, given any u^i∈ H^1_r(Ω^H_p), there is a unique solution u∈H^1_r(Ω_p^H) of the variational problem (<ref>).






 §.§ Bloch transform and well-posedness of the scattering problems


Firstly, we have to introduce the Bloch transform defined on the periodic space Ω^H. For any function ϕ∈ C_0^∞(Ω^H), define

    (_Ωϕ)(α,x)=∑_j∈ϕ(x_1+Λ j,x_2)e^-ıαΛ j,  α∈; x∈Ω^H.


In this paper, we denote the one dimensional Bloch transform defined on any periodic domain in ^2 by _Ω.

For any fixed α∈, the function (_Ωϕ) is α-quasi-periodic with period Λ, i.e.,

    (_Ωϕ)(α,(
    x_1+Λ j
    x_2
    ))=e^ıαΛ(_Ωϕ).

For the positive integer m, let the space H^m(;H^s(Ω^Λ_H)) be defined with the norm

    ϕ_H^m(;H^s(Ω^Λ_H))=[∑_ℓ=0^m∫_∂_α^ℓϕ(α,·)^2_H^s(Ω^Λ_H)]^1/2,

then extend the definition to any real number r by interpolation and duality arguments. Define the subspace H_0^r(;H^s_α(Ω^Λ_H)) with a further condition that ϕ∈ H^r(;H^s(Ω^Λ_H)) is Λ^*-periodic in α and ϕ(α,·)∈ H^s_α(Ω^Λ_H). The mapping property of the Bloch transform was proved in <cit.>.

The Bloch transform extends to an isomorphism between H^s_r(Ω) and H_0^r(;H^s_α(Ω^Λ_H)) for any s, r∈. When s=r=0, _Ω is an isometry with the inverse operator 

    (^-1_Ω^Hw)(
    x_1+Λ j
    x_2
    )=[Λ/2π]^1/2∫_ w(α,x)e^ıΛ jαα̣,   x∈Ω^Λ_H, j∈.

Moreover, the L^2-adjoint operator J^*_Ω equals to its inverse.


As the Bloch transform is defined on periodic domains, we have to transform the problem defined on a locally perturbed domain into one defined on a periodic domain. Following <cit.>, define a  diffeomorphism Φ_p from Ω^H into Ω_p^H by

    Φ_p(x_1,x_2)=(x_1,x_2+(x_2-H)^3/(ζ(x_1)-H)^3(ζ_p(x_1)-ζ(x_1)),

which has a support contained in Ω^Λ_H. Then u_T:=u∘Φ_p satisfies the following variational problem

    ∫_Ω_H[A_p∇ u_T·∇ϕ_T-k^2 c_p u_Tϕ_T]x̣   -∫_Γ_HT^+(u_T|_Γ_H)ϕ_Tṣ
       =∫_Γ_H[∂ u^i/∂ x_2-T^+(u^i|_Γ_H)]ϕ_Tṣ

for all ϕ_T∈ H^1(Ω^H) with compact support, where

    A_p(x)=|∇Φ_p(x)|[(∇Φ_p(x))^-1((∇Φ_p(x))^-1)^T]∈ L^∞(Ω^H,^2× 2),
          c_p(x)=|∇Φ_p(x)|∈ L^∞(Ω^H).


Let w=_Ω u_T then w∈ L^2(;H^1_α(Ω^Λ_H)), and satisfies the following variational problem for any ϕ∈ L^2(;H^1_α(Ω^Λ_H))

    ∫_ a_α(w(α,·),ϕ(α,·))α̣+b(^-1_Ω w,^-1_Ωϕ)=∫_∫_Γ^Λ_Hf(α,·)ϕ(α,·)ṣα̣,

where 

    a_α(w(α,·),ϕ(α,·))   =∫_Ω^Λ_H[∇ w(α,·)·∇ϕ(α,·)-k^2w(α,·)ϕ(α,·)]x̣
       -∫_Γ_HT^+_α(w(α,·)|_Γ_H)ϕ(α,·)
    
          b(ξ,ψ)=[Λ/2π]^1/2∫_Ω^Λ_H(A_p-I_2)∇ξ·∇ψx̣-k^2[Λ/2π]^1/2∫_Ω^Λ_H(c_p-1)ξψx̣,
          f(α,·)=∂_Ω u^i(α,·)/∂ x_2-T^+_α[(_Ω u^i)(α,·)|_Γ^Λ_H],

with T^+_α defined by

    T^+_α(ψ)=ı∑_j∈√(k^2-|Λ^* j-α|^2)ψ(j)e^ı(Λ^*j-α)x_1  ψ=∑_j∈ψ(j)e^ı(Λ^*j-α)x_1.




The system is decoupled when the surface is purely periodic. In this case, Φ_p is the identity operator, thus A_p=I_2 and c_p=1. From the definition of the term b, it disappears in this case, thus the system has the form of

    ∫_ a_α(w(α,·),ϕ(α,·))α̣=∫_∫_Γ^Λ_Hf(α,·)ϕ(α,·)ṣα̣.





At the end of this subsection, we will list some of  the useful results in <cit.>.


Assume u^i∈ H^1_r(Ω^H_p) for some r∈[0,1), then u_T∈H^1_r(Ω^H) satisfies (<ref>) if and only if w=_Ω u_T∈ H_0^r(;H^1_α(Ω^Λ_H)) satisfies (<ref>). 
If Γ_p is Lipschitz, then (<ref>) is uniquely solvable in H^r_0(;H^1_α(Ω^Λ_H)) for all u^i∈ H^1_r(Ω^H_p), r∈[0,1).



If u^i∈ H^2_r(Ω^H_p) for r∈[0,1) and ζ,ζ_p∈ C^2,1()， then w(α,·)∈H^2_α(Ω^Λ_H) and u_T=^-1_Ω w∈H^2(Ω^H).



If Γ_p is Lipschitz continuous and u^i∈ H^1_r(Ω^H_p) for r∈(1/2,1), then w∈ L^2(;H^1_α(Ω^Λ_H)) equivalently satisfies for all α∈ and ϕ_α∈H^1_α(Ω^Λ_H) that 

    a_α(w(α,·),ϕ_α)+b(^-1_Ω w,ϕ_α)=∫_Γ^Λ_Hf(α,·)ϕ_αṣ,






 §.§ Numerical method


Assume that _h is a family of regular and quasi-uniform meshes (see <cit.>) of Ω^Λ_H, such that the width of each mesh is not larger than h. Suppose {ϕ^(l)_M}_l=1^M is the piecewise linear nodal function that equals to one at the l-th nodal x^(ℓ)=(x_1^(ℓ),x_2^(ℓ))^⊤ while equals to zero at other nodal points. Then the space spanned by {ϕ^(l)_M}_l=1^M is a subset of H^1_0(Ω^Λ_H). 

For the domain , define the uniformly located grid points 

    α_N^(j)=-π/Λ+2π j/NΛ,   j=1,2,…,N,

and {ψ_N^(j)}_j=1^N be the basis of the functions that are piecewise constant on [α^(j)_N-π/(NΛ),α^(j)_N+π/(NΛ)]. The finite element space X_N,h is formulated by

    X_N,h={v_N,h(α,x)=∑_j=1^N∑_l=1^M v^(j,l)_N,hexp(-ıα_N^(j)x_1) ψ^(j)_N(α)ϕ^(l)_M(x): v^(j,l)_N,h∈}.

Then the numerical scheme is to seek for a finite element solution w_N,h∈ X_N,h such that 

    ∫_ a_α(w_N,h,ϕ_N,h)α̣+b(^-1_Ω w_N,h,^-1_Ωϕ_N,h)=∫_∫_Γ^Λ_Hf(α,·)ϕ_N,hṣα̣

for all ϕ_N,h∈ X_N,h. The unique solvability and  convergence of the finite element scheme of the variational form (<ref>).


Assume that u^i∈ H^2_r(Ω^H_p) for r≥ 1/2, ζ,ζ_p∈ C^2,1(), then is uniquely solvable in X_N,h. The error between the numerical solution w_N,h and exact one w is bounded be

    w_N,h-w_L^2(;H^l_α(Ω^Λ_H))≤ C h^1-l(N^-r+h)f_H^r_0(;H^1/2_α(Γ^Λ_H)).





For the special case that the surface is non-perturbed, the same error estimate of the numerical result is also obtained. Moreover, when u^i satisfies that _Ω u^i∈ W_p^1,p(;H_α^1(Ω^Λ_H)) for some p∈[1,2), then 

    w_N,h-w_L^2(;H_α^ℓ(Ω^Λ_H))≤ C(h^2-ℓ+N^-1)f_L^2(;H^1/2_α(Γ^Λ_H)).

The space W_p^1,p(;H_α^1(Ω^Λ_H)) is defined by

    w^p_W_p^1,p(;H_α^1(Ω^Λ_H))=∫_[w(α,·)^p_H^1_α(Ω^Λ_H)+∂/∂αw(α,·)^p_H^1_α(Ω^Λ_H)]α̣.

For details see <cit.>.



We will introduce the numerical method in <cit.> briefly in the rest of this section. Let w^s(α,x)=w(α,x)-(_Ω u^i)(α,x) and w_0(α,x)=e^ıα x_1w^s(α,x), then w_0∈ L^2(;H_0^1(Ω^Λ_H))).

Recall that for the choice of X_N,h, piecewise constant basis functions are chosen for α in , we can define the  discrete space as follows. Still denote by {ϕ_M^(ℓ)}_ℓ=1^M the piecewise linear nodal functions in Ω^Λ_H, and {ψ_N^(ℓ)}_ℓ=1^N the piecewise constant basic function in . Suppose there is a positive integer M'<M such that ϕ_ℓ|_Γ^Λ=0 for any ℓ=1,…,M' and ϕ_ℓ|_Γ^Λ≠0 for any ℓ=M'+1,…,M.. Thus  span{ϕ_M'^(ℓ)}_ℓ=1^M'⊂H_0^1(Ω^Λ_H). 

Define the discrete space by

    Y^(j)_N,h={v^(j)_N,h=∑_ℓ=1^Mv_N,h^(j,l)ψ_N^(j)(α)ϕ_M^(l)(x): v_N,h^(j,l)∈}⊂ L^2(;H_0^1(Ω^Λ_H))

for j=1,…,N and also the subspace of functions that vanish on Γ^Λ

    Y^(j)_N,h={v^(j)_N,h=∑_ℓ=1^M'v_N,h^(j,l)ψ_N^(j)(α)ϕ_M^(l)(x): v_N,h^(j,l)∈}⊂ L^2(;H_0^1(Ω^Λ_H)).

Then setting

    Y^0_N,h=Y^(1)_N,h⊕⋯⊕ Y^(N)_N,h,  Y^0_N,h=Y^(1)_N,h⊕⋯⊕Y^(N)_N,h,

then for any w_0∈ X_N,h, there is a unique vector (w_0^(j))_j=1^N∈ Y^0_N,h such that w_0=∑_j=1^N exp(-ıα_N^(j)x_1) w^(j)_0.  Then the inverse Bloch transform ^-1_Ω applied to the function w_0 equals to the numerical inverse Bloch transform ^-1_Ω,N applied to {e^-ıα_n^(j)x_1 w_0^(j)(α_N^(j),x)}_j=1^N

    _Ω^-1 w_0=_Ω,N^-1({e^-ıα_n^(j)x_1 w_0^(j)(α_N^(j),x)}_j=1^N)=Λ^*/N∑_j=1^N e^-ıα_N^(j)x_1 w_0^(j)(α_N^(j),x).


Suppose w_0=∑_j=1^N exp(-ıα_N^(j)x_1) w^(j)_0 where w_0^(j)=∑_m=1^M w_N,h^(j,m)ψ_N^(j)ϕ_M^(m), omitting the details in <cit.>, the final linear system is set up as follows.

    (
    A_1     0    ⋯    0     C_1
    
    0       A_2    ⋯    0     C_2
    ⋮   ⋮   ⋮   ⋮   ⋮
    
    0     0   ·    A_N     C_N
    
    B_1     B_2    ·    B_N     I_M
    )(
    W_1
     W_2
    ⋮
     W_N
     U
    )=(
    
    F_1
     F_2
    ⋮
     F_N
     0
    ),

where A_j, B_j, C_j, I_m are M× M matrices and F_j are M× 1 vertices, W_j=(w_N,h^(j,1),…,w_N,h^(j,m))^⊤ and U=(u_h^(1),…,u_h^(M))^⊤, where

    u_h^(ℓ)=Λ^*/N∑_j=1^N e^-ıα_N^(j)x_1^(ℓ)w_N,h^(j,ℓ).


The matrix linear system (<ref>) is a M(N+1)× M(N+1)-sparse matrix with 3N+1 non-zero M× M blocks. In <cit.>,  a GMRES iteration method with an incomplete LU-decomposition based pre-conditioner was applied to the solution of the linear system. The convergent rate with respect to N is N^-r where 1/2≤ r<1, which means that if the mesh size h is small enough, a small error requires a relatively large N. However, if N gets larger, there are several disadvantages to solve the linear system (<ref>).


  * Setting up the matrix takes up a lot of time and memory. 

  * The time of solving the system is long.

To solve the scattering problems more efficiently, we would like to look for a new method that converges faster in N.






§ REGULARITY RESULTS FOR THE SCATTERING PROBLEMS


A natural way to look for a new numerical method that converges faster as N tends to ∞ is to go deeper into the regularity of the field w(α,·) with respect to α∈. Theorem a in <cit.> gave a result for non-perturbed surfaces with incident plane waves. It could be a good inspiration for the problems in this paper.


Suppose the incident field u^i(α,x)=e^ıα x_1-ı√(k^2-α^2)x_2 for α∈(-k,k), and it is scattered by a smooth enough sound soft periodic surface.
 a) The total field, denoted by u(α,x), depends continuously on α∈(-k,k) and analytically on α in the set

    {α∈(-k,k):  (nΛ^*-α)^2≠ k^2 for every  n∈}.


 b) Let α_0∈(-k,k) such that (n_0Λ^*-α_0)^2=k_0^2 for some n_0∈. Then there is a neighbourhood U of α_0 and quasi-periodic functions v(α,·), w(α,·)∈H^1_α(Ω^Λ_H) that depend analytically on α∈ U and satisfies u(α,x)=v(α,x)+√(k^2-|n_0Λ^*-α|^2)  w(α,x).



 



From the proof of Theorem a in <cit.>, the singularity of the total fields with respect to incident plane waves comes from the DtN map. The singularity occurs at the points where at least one term of the DtN map vanishes (so called Wood anomalies, see <cit.>). However, for the scattering problems in this paper, the singularities may also come from the slow decay of the incident fields. In this section, we will only focus on the singularity comes from the DtN map. We will firstly consider the case that the incident field with compact support, and then extend the results to more generalized and commonly used incident fields. We will give a brief description of these singularity results, for details see <cit.>.




 §.§ Notations and  function spaces




 Define the set of Wood anomalies by

    :={α∈: ∃  n∈,  such that  |Λ^*n-α|=k},

which is a countable set in . Let {α_j}_j∈ be an ascending sequence of all the points in . To define the continuous or analytical "dependence" rigorously, we have to define the following function spaces first.


Let ⊂ be an interval (= included), W⊂^2 be a bounded domain, and define the space of functions defined in the  × W that depend analytically on the first variable by

    C^ω(;S(W)):={f∈'(× W): ∀ α_0∈, ∃ δ>0, s.t., ∀ α∈(α_0-δ,α_0+δ)∩,.
    .∃  C>0, f_n∈ S(W), s.t., f(α,x)=∑_n=0^∞ (α-α_0)^n f_n(x), f_n_S(W)≤ C^n},

where S(W) is a Sobolev space defined on the domain W. In this paper, W=Ω^Λ_H, the function space S(W) could be either H^n(Ω^Λ_H) or H^n(Ω^Λ_H) for n=0,1,2. Define the space of functions that depend C^n-continuously on the first variable

    C^n(;S(W)):={f∈'(× W): ∀α∈, j=0,…,n,∂^j f(α,·)/∂α^j∈ S(W),.
    .
    moreover, ∂^j f(α,·)/∂α^j_S(W) is uniformly bounded for α∈},

thus we can define the space C^∞(;S(W)) in the same way

    C^∞(;S(W)):={f∈'(× W): ∀α∈, j=0,…,∞,∂^j f(α,·)/∂α^j∈ S(W),.
    .
    moreover, ∂^j f(α,·)/∂α^j_S(W) is uniformly bounded for α∈},

We can also define the subspace of  C^n(;S(W)) for finite interval =[A_0,A_1] by

    C_p^n(;S(W)):={f∈ C^n(; S(W)): ∀ j=0,…,n,.
                           .
    lim_t→ A_0^+∂^j f/∂ t^j(t,·)_S(W)=lim_t→ A_1^-∂^j f/∂ t^j(t,·)_S(W)=0}

Similarly, we can also define C_p^∞(;S(W)). Define the subspace C^ω(;S_α(W)) by

    C^ω(;S_α(W)):={f(α,·)∈ C^ω(;S(W)); f
     is α-quasi-periodic for any  α∈},

and also define C^n(;S_α(W)) and C^∞(;S_α(W)) in the same way.





 §.§ Regularity properties for Bloch transformed solutions


Firstly, assume that the incident field u^i is compactly supported in Ω^H_p, thus u^i|_Γ_H∈ H^1_r(Γ_H)  for any  r∈, especially, for any r>1/2. Thus f(α,·) in (<ref>) belongs to the space H_0^r(;H^-1/2_α(Γ^Λ_H))∩ C^ω(;H^-1/2_α(Γ^Λ_H)). From  Theorem <ref>, the Bloch transformed scattered field w∈ H_0^r(;H^1_α(Ω^Λ_H)). Moreover,  for any fixed α∈ and ϕ_α∈H^1_α(Ω^Λ_H),

    a_α(w(α,·),ϕ_α)+b(u_T,ϕ_α)=∫_Γ^Λ_H f(α,·)ϕ_αṣ.


Set w_α(x):=exp(ıα x_1)w(α,x), then w_α is a periodic function for any fixed α∈. Define ϕ_α:=exp(-ıα x_1)ϕ(x) for any periodic function ϕ. Then replace w and ϕ_α by w_α and ϕ, for each w_α, it satisfies the following variational equation

    ∫_Ω^Λ_H[_α w_α·_αϕ-k^2w_αϕ]-∫_Γ^Λ_HT^+_α[w_α|_Γ^Λ_H]ϕṣ=∫_Γ^Λ_Hg(α,·)ϕṣ-b_α(u_T,ϕ),
 
where _α=-ıα(1,0)^⊤, T^+_α is defined by

    T^+_α(ψ)=ı∑_j∈√(k^2-|Λ^*j-α|^2)ψ(j)e^ıΛ^* x_1,  ψ=∑_j∈ψ(j)e^ıΛ^* j x_1,

the sysquilinear form b_α(·,·) is defined by

    b_α(ξ,ψ)=[Λ/2π]^1/2∫_Ω^Λ_He^ıα x_1[(A_p-I)ξ·_αψ-k^2(c_p-1)ξψ]x̣,

and 

    g(α,x)=e^ıα x_1f(α,x).


As f∈ C^ω(;H^-1/2_α(Γ^Λ_H)), the function g(α,·)∈ C^ω(;H^-1/2_0(Ω^Λ_H)). So the right hand side of equation (<ref>) depends analytically on α, for any g∈ C^ω(;H^-1/2_0(Ω^Λ_H)) and ϕ∈H^1_0(Ω^Λ_H). From Riesz representation theorem, there is an operator F_α that maps C^ω(;H^-1/2_0(Ω^Λ_H)) to the space H^1_0(Ω^Λ_H), such that

    ∫_Γ^Λ_Hg(α,·)ϕṣ-b_α(u_T,ϕ)=∫_Ω^Λ_HF_α(g)ϕx̣.
 
Moreover, the operator F_α depends analytically on α, see Chapter 7, Section 1 in <cit.>. Following the proof of Theorem a in <cit.>, we can obtain the following property of w_α from the singularity with respect to α in the DtN map T^+_α.


Suppose the incident field u^i∈ H^1(Ω_H) has a compact support. Then  w_α(x) that solves (<ref>) depends continuously in α∈ and analytically on ∖. If α_0∈, i.e., there is an n_0∈ such that |n_0Λ^*-α_0|=k. Then there is a neighbourhood U of α_0 and two functions v_1^α, v_2^α∈H^1_0(Ω^Λ_H) that depend analytically on α, such that w_α=v_1^α+√(k^2-|n_0Λ^*-α|^2) v_2^α.


The singularity result could have a simplified form.


If α_0∈, then there is a neighbourhood U of α_0 and two functions w_1, w_2∈ C^ω(U;H^1_α(Ω^Λ_H))  such that w(α,·)=w_1(α,·)+√(α-α_0) w_2(α,·).



For any α_0∈, there is either 1) an n_0∈ such that |Λ^*n-α_0|=k or 2) two n_1, n_2∈ such that |Λ^*n_1-α_0|=|Λ^*n_2-α_0|=k.

Suppose |Λ^* n_j-α_0|=k, then Λ^* n_j- k=α_0 or Λ^* n_j+ k=α_0. If Λ^* n_j-k=α_0, 

    √(k^2-|Λ^* n_j-α|^2)=√(α-α_0)√(2k+α_0-α),

where the second term is analytic in a small neighbourhood of α_0.  Thus the total field w_α could be written into the form w_α=w_1^α+√(α-α_0) w_2^α, where w_1^α=v_1^α and w_2^α=√(2k+α_0-α) v_2^α. Let w_1(α,·)=e^-ıα x_1w_1^α and w_2(α,·)=e^-ıα x_1w_2^α, then w(α,·)=w_1(α,·)+√(α-α_0) w_2(α,·). The case that Λ^*n_j+k=α_0 is similar.  The proof is finished.




We have just investigated the decomposition of the Bloch transformed field w(α,·), when  u^i is compactly supported. However, the result for this special case could be  extended to more generalized cases. Firstly, define the following two conditions for a function u∈ C^0(;S(W)).


  * For any subinterval _0⊂∖, u∈ C^ω(_0;S(W)).

  * For any α_j∈∩, there is a small enough δ>0 with two pairs v_1, w_1∈ C^ω((α_j-δ,α_j];S(W)) and v_2, w_2∈ C^ω([α_j,α_j+δ);S(W)) such that

    u=v_1+√(α-α_j) w_1 for α∈(α_j-δ,α_j],
          u=v_2+√(α-α_j) w_2 for α∈[α_j,α_j+δ).


Then the function spaces are defined as follows.

    ^ω_c(;S(W);):={ u∈ C^0(;S(W)): u satisfies conditions 1 and 2},

and denote its subspace of all α-quasi-periodic functions by ^ω_c(;S_α(W);).



The incident field u^i∈ H_r^2(Ω^H_p) for r>1/2 satisfies that its Bloch transform _Ω u^i∈_c^ω(;H^2_α(Ω^Λ_H);).





From <cit.>, the two types of incident fields, i.e., the half space Green's function and the Herglotz wave functions, satisfy Assumption <ref>. Set Λ=2π for simplicity,  then Λ^*=1.

The first step is to investigate the function e^ı√(k^2-|n-α|^2) x_2. From the Taylor's expansion,

    e^ı√(k^2-|n-α|^2) x_2=∑_ℓ∈(ı√(k^2-|n-α|^2)x_1)^ℓ/ℓ!

For any α∈∖, e^ı√(k^2-|n-α|^2) x_2 is an analytic function with respect to α. For any α_0∈, there is an integer n such that either k=n-α_0 or -k=n-α_0. Take  k=n-α_0 for example,

    e^ı√(k^2-|n-α|^2) x_2
    =∑_ℓ∈(ı√(k^2-|n-α|^2)x_2)^2ℓ/(2ℓ)!+∑_ℓ∈(ı√(k^2-|n-α|^2)x_2)^2ℓ+1/(2ℓ+1)!
    
    =   ∑_ℓ∈(|n-α|^2-k^2)^ℓx_2^2ℓ/(2ℓ)!+ı√(α-α_0)√(2k+α_0-α)∑_ℓ∈(|n-α|^2-k^2)^ℓx_2^2ℓ+1/(2ℓ+1)!.

As the two series above converges absolutely, they are both analytic functions with respect to α. In a small neighbourhood of α_0, √(2k+α_0-α) is also analytic, thus there are two functions ψ_1(α,x_1) and ψ_2(α,x_1) that depend analytically in a small neighbourhood of α_0 such that

    e^ı√(k^2-|n-α|^2) x_1=ψ_1(α,x_1)+√(α-α_0)ψ_2(α,x_1).


1) The half space Green's function, defined by

    G(x,y)=ı/4[H_0^(1)(k|x-y|)-H_0^(1)(k|x-y'|)]

with the point y=(y_1,y_2)^⊤ be above the (locally perturbed) periodic surface and y'=(y_1,-y_2)^⊤. It is easy to check that the Bloch transform has the following representation

    (_Ω G(·,y))(α,x)=1/2π∑_j∈e^ı(j-α)(x_1-y_1)+ı√(k^2-|j-α|^2)y_2 sinc(√(k^2-|j-α|^2)x_2)x_2,

where sinc is an analytic function defined by sin t/t for t≠ 0 with sinc(0)=1. As the Taylor's series sinc=∑_j=0^∞(-t^2)^n/(2n+1)!, the function sinc(√(k^2-|j-α|^2)x_2) depends analytically on α∈. Thus the only factor that may cause singularity is e^ı√(k^2-|j-α|^2) for finite number of j's. From (<ref>), the function _Ω G(·,y) also have the similar form of decomposition as (<ref>).

2) The upward/downward propagating Herglotz wave function defined by

    (H^∓ g)(x)=∫__±e^ı k x· dg(d)ṣ(d)=∫_-π/2^π/2e^ı k(cos t x_1±sin t x_1)ϕ(t)ṭ.

where _+ (_-) is the upper (lower) unit circle. (H^∓ g) satisfies Assumption (<ref>)  when ϕ(t)=h(cos(t))cos(t) where h is an analytic function defined in [0,1] with h(0)=0.

Take H^-g for example. From <cit.>, the Bloch transform has the following representation

    (_Ω H^-g)(α,x)   =∑_|j-α|<ke^ı(j-α)x_1-ı√(k^2-|j-α|^2)x_2ϕ(arcsin(j-α)/k)/√(k^2-|j-α|^2)
       =1/k∑_|j-α|<ke^ı(j-α)x_1-ı√(k^2-|j-α|^2)x_2h[√(1-|j-α|^2/k^2)]

This is a finite series, we only need to investigate it by the decompositions of the exponential terms and the property of the  analytic function h. The calculation process is omitted.




The plane waves do not satisfy Assumption <ref>, for the plane waves belong to H_r^2(Ω^H_p) for r<-1/2.


With the incident fields satisfying Assumption <ref>, the following results hold, see <cit.>.



Suppose u^i  satisfies Assumption  <ref>, then w∈ H_0^r(;H_α^2(Ω^Λ_H)) satisfies (<ref>) belongs to the space _c^ω(;H^2_α(Ω^Λ_H);).



When α_0∈∖, as _Ω u^i∈_c^ω(;H_α^2(Ω^Λ_H)), it depends analytically on α in a neighbourhood U⊂∖. From Corollary <ref>, w(α,·) also depends analytically on α∈ U.

When α_0∈∩. Consider the case that α∈(α_0-δ,α_0]. Suppose there are v_1(α,·), v_2(α,·)∈ C^ω((α_0-δ];H^2_α(Ω^Λ_H)) such that _Ω u^i=v_1(α,·)+√(α-α_0) v_2(α,·), then the Bloch transformed fields w_1(α,·), w_2(α,·) with respect to v_1(α,·), v_2(α,·) have the following decompositions

    w_1(α,·)=w_1^1(α,·)+√(α-α_0) w_1^2(α,·);   w_2^α=w_2^1(α,·)+√(α-α_0) w_2^2(α,·),

where w_j^ℓ(α,·)∈ C^ω((α_0-δ];H^2_α(Ω^Λ_H)) for j, ℓ=1,2. Then the Bloch transformed field w with respect to the field _Ω u^i could be written into

    w(α,·)   =w_1(α,·)+√(α-α_0)w_2(α,·)
       =[w_1^1(α,·)+(α-α_0)w_2^2(α,·)]+√(α-α_0) [w_1^2(α,·)+w_2^1(α,·)].

Similar result could also be obtained for α∈[α_0,α_0+δ). 

Thus w(α,·)∈_c^ω(;H^2_α(Ω^Λ_H);)), the proof is finished.


The decomposition of the Bloch transformed total field could also be written by smooth functions.


Suppose u^i satisfies Assumption <ref>, then the Bloch transform w(α,·)=(_Ω u)(α,·) depends continuously in α∈. For each j∈, for α∈[α_j,α_j+1], there are three  functions w_ℓ^j(α,·)∈ C^∞([α_j,α_j+1];H^2(Ω^Λ_H)) where ℓ=0,1,2 such that

    w(α,·)=w_0^j(α,·)+√(α-α_j) w_1^j(α,·)+√(α-α_j+1) w_2^j(α,·).






§ THE INVERSE BLOCH TRANSFORM


In <cit.>, the accuracy of the numerical inverse Bloch transform plays an important role in the convergence of the numerical schemes. In this section, we will discuss a new numerical method for the inverse Bloch transform ^-1_Ω w.  For references we refer to <cit.>, Chap 3.5.



 From the previous section, the Bloch transformed solution w(α,·) only has square-root like singularities at the points in , i.e., any α∈ such that |Λ^* n-α|=k for some n∈. Thus there is an n∈, such that α=Λ^*n-k or α=Λ^*n+k. Thus we can define the real number $̨ by 

    :̨=min{|Λ^*n-k|: n∈},

 then≤̨Λ^*/2andcould also be represented by$̨
 
    ={Λ^*n±:̨ n∈}.
 
 There are two different kinds of representations of points in  by $̨.

  *  Case 1, =̨mΛ^*/2 where m=0,1, the set of singularities ={+̨nΛ^*: n∈}.

  *  Case 2, ≠̨mΛ^*/2 for m=0,1, ={+̨mΛ^*, -+̨nΛ^*: m,n∈}.



For simplicity, we redefine one periodic cellby=(-,̨Λ^*-]. Define the singular points inas follows.

  *  For Case 1, let a_0=-$̨,a_1=Λ^*-$̨.

  *  For Case 2, let a_0=-$̨,a_1=$̨, a_2=Λ^*-$̨.

The  redefinition is proper as the fieldw(α,·)isΛ^*-periodic inαfor any fixedx,

    ∫_-Λ^*/2^Λ^*/2w(α,·)α̣=∫_-^Λ^*-w(α,·)α̣=^-1_Ω w.


 §.§ Change of variables



To design a new algorithm with a greater convergence rate, we intend to use a "better" contour to take place of the straight line in the integral of the inverse Bloch transform. 
The strategy adopted in this paper is illustrated and utilized in Chap 3.5, <cit.>. In this section, we will recall the method and apply it to the numerical approximation of the inverse Bloch transform. Define a  monotonic functiong∈ C^∞[a_j,a_j+1]such that 

    g(a_j)=a_j, g(a_j+1)=a_j+1;  g'>0 in (a_j,a_j+1).

 Moreover, assume that the following conditions are also satisfied.

Assume that there is a δ>0, such that

    g(t)-a_j=C_j h^2(t-a_j),   t∈[a_j,a_j+δ);
          a_j+1-g(t)=C_j h^2(a_j+1-t),   t∈(a_j+1-δ,a_j+1],
where the function h∈ C^∞[0,δ)  satisfies either of the following two conditions.
  *  There is a positive integer m∈ such that h(t)=O(t^m+1) as t→ 0^+.

  * h(t)=o(t^n) for any n∈ as t→ 0^+.

The contours  in one periodic cell(Λ=2π)are shown in Figure <ref>. For Case 1, we takek=1as an example, and the contour is plotted in (a); for Case 2,k=1.2is taken as an example, see (b).




























With the transform ofα=g(t), the inverse Bloch transform has the representation

    (^-1w)(x)=∫_w(α,·)α̣=∫_ w(g(t),·)g'(t)ṭ.
 
For Case 1,=(a_0,a_1],

    (^-1w)(x)=∫_a_0^a_1 w(g(t),·)g'(t)ṭ;
 
for Case 2,=(a_0,a_2],

    (^-1w)(x)=∫_a_0^a_1 w(g(t),·)g'(t)ṭ+∫_a_1^a_2 w(g(t),·)g'(t)ṭ.
 

Let the interval[A_0,A_1]be either[a_0,a_1](Case 1 & 2)  or[a_1,a_2](Case 2). Define 

    v(t,·):=w(g(t),·)g'(t),   t∈[A_0,A_1],

then the inverse Bloch transform

    (^-1_Ω w)(x)=∫_ v(t,x)ṭ,   x∈Ω^Λ_H.



We will study the Bloch transform whenwsatisfies certain conditions.
w∈ C^0([A_0,A_1];S_α(W)) has a form of 

    w=w_2+√(α-A_0)  w_0+√(α-A_1) w_1,

 where w_0, w_1, w_2∈ C^∞([A_0,A_1];S(W)).

If there is an m∈ such that h(t)=O(t^m+1) as t→ 0^+, for any w satisfies Assumption <ref>, v defined by (<ref>) belongs to the  function space C^2m_p([A_0,A_1];S_g(t)(W))∩ C^∞([A_0,A_1];S_g(t)(W)).
First, consider the case that S(W)=L^2(Ω^Λ_H).
Consider the function with α near the point A_0. From Assumption <ref>, for any 0<δ<A_1-A_0, in the interval   [A_0,A_0+δ], w(α,x) has the form of 

    w(α,x)=w_0(α,x)+√(α-A_0) w_0(α,x),

where w_0(α,x)=w_2(α,x)+√(α-A_1)  w_1(α,x) belongs to C^∞([A_0,A_0+δ];S(W)). 
Then for any t∈[A_0,A_0+δ],

    v(t,x)   =w_0(g(t),x)g'(t)+w_0(g(t),x)√(g(t)-A_0)  g'(t)
       =2 C_0 w_0(A_0+h^2(t-A_0),x)h(t-A_0)h'(t-A_0)
       +2 C_0 w_0(A_0+h^2(t-A_0),x)h^2(t-A_0)h'(t-A_0).
From the representation of w_0 and v, 

    v(t,·)_S(W)≤ C(g)[w_0_S(W)+w_1_S(W)+w_2_S(W)].

Take the first term for example. As h∈ C^∞[0,δ) and g∈ C^∞[A_0,A_1], the function v∈ C^∞([A_0,A_1];S_g(t)(W)).  As h(t)=O(t^m+1), h(t-A_0)h'(t-A_0)=O((t-A_0)^2m+1), then 

    lim_t→ A_0^+[h(t-A_0)h'(t-A_0)]^(ℓ)(t)=0,     for any ℓ=0,1,…,2m.
 
From the ℓ-th derivative w.r.t. t
    ∂^ℓ/∂ t^ℓ[w_0(A_0+h^2(t-A_0),x)h(t-A_0)h'(t-A_0)]
    
    =
       ∑_j=0^ℓ(ℓ
    j
    )∂^j/∂ t^j[w_0(A_0+h^2(t-A_0),x)]∂^ℓ-j/∂ t^ℓ-j[h(t-A_0)h'(t-A_0)],

in each term, the first element ∂^j/∂ t^j[w_0(A_0+h^2(t-A_0),x)] is a finite sum of the terms if form of ∂^ℓ/∂ t^ℓ[w_0(A_0+h^2(t-A_0),x)]Π_m=0^N [h^(m)(t-A_0)]^n(m), where N is some positive integer and n(m)≥ 0 in . Thus the norm ∂^j/∂ t^j[w_0(A_0+h^2(t-A_0),x)]_S(W) is uniformly bounded for t→ A_0^+. From lim_t→ A_0^+[h(t-A_0)h'(t-A_0)]^ℓ(t)=0, 

    ∂^ℓ v(t,x)/∂ t^ℓ_S(W)→ 0  as t→ A_0^+.

Similarly, we can also prove that ∂^ℓ v(t,x)/∂ t^ℓ_S(W)→ 0 as t→ A_1^-, ℓ=0,1,…,2m, thus v∈ C_p^2m([A_0,A_1];S_g(t)(W)). 

The case that S(W)=H^n(Ω^Λ_H) for any n∈ could be proved similarly, thus we omit the proof here.

A direct corollary shows a similar result for the case thath(t)=o(t^n)ast→ 0^+,∀ n.

If h(t)=o(t^n) as t→ 0^+ for any positive integer n, for any w satisfies Assumption <ref>, v∈ C_p^∞([A_0,A_1];S_g(t)(W)).
 

From the smoothness of the functionv(t,·), the bounds of the derivatives ofvdepend on the functionswandg, and the order of the derivative.



Suppose w satisfies Assumption <ref>, v is defined by (<ref>). If h(t)=O(t^m+1) for some m∈, then v belongs to the functions space C_p^2m([A_0,A_1];S_g(t)(W))∩ C^∞([A_0,A_1];S_g(t)(W)) with

    ∂^j/∂ t^jv(t,·)_S(W)≤ C(j,g)w_S(W),   j=0,1,…,2m;

if h(t)=o(t^n) for any n∈, then v∈ C_p^∞([A_0,A_1];S_g(t)(W)) with 

    ∂^j/∂ t^jv(t,·)_S(W)≤ C(j,g)w_S(W),   j=0,1,…,∞,

where C(j,g) is a constant that depends on j and g, where 

    w_S(W):=w_0_S(W)+w_1_S(W)+w_2_S(W).


 §.§ Numerical method for the inverse Bloch transform


In this subsection, we will analyse the new numerical method of the inverse Bloch transform. Consider the definite integral

    ∫_A_0^A_1 v(t,·)ṭ=∫_A_0^A_1w(α,·)α̣,
 
wherewsatisfies Assumption <ref> andvis defined by (<ref>). We will approximate the function v by trigonometrical interpolation (for references see <cit.>), and then study the convergence of the numerical integration based on the approximation. Let[A_0,A_1]be divided uniformly intoNsubintervals, whereNis assumed to be even in this paper. Let the uniformly located grid points

    t_j=A_0+A_1-A_0/Nj,    j=1,…,N,

and define the basic functions by

    ϕ_N^(m)(t)=1/N∑_l=-N/2+1^N/2exp(ı l (t-t_m)2π/A_1-A_0), j=1,2,…,N,

thenϕ_N^(m)(t_ℓ)=δ_ℓ,m, where

    δ_ℓ,m=
    1,  ℓ=m;
    
    0,  ℓ≠ m.

The basic functions {ϕ_N^(m)}_m=1^N  are orthogonal.

For any j,j'∈{1,…,N}, the inner product shows

    ∫_A_0^A_1ϕ^(j)_N(t)ϕ^(j')_N(t)ṭ
    
    =   1/N^2∫_A_0^A_1∑_ℓ=-N/2+1^N/2∑_ℓ'=-N/2+1^N/2exp(ıℓ (t-t_j)2π/A_1-A_0)exp(-ıℓ' (t-t_j')2π/A_1-A_0)ṭ
    
    =   A_1-A_0/N^2∑_ℓ=-N/2+1^N/2exp(ıℓ(t_j'-t_j)2π/A_1-A_0)=A_1-A_0/N^2∑_ℓ=-N/2+1^N/2exp(ıℓ (j'-j)2π/N)
    ,
 
then the integral equals to 0 if and only if j≠ j'. Thus

    ∫_A_0^A_1ϕ^(j)_N(t)ϕ^(j')_N(t)ṭ=A_1-A_0/Nδ_j,j'.

The basic functions {ϕ_N^(m)}_m=1^N  are orthogonal. The proof is finished.


Letv_Nbe the interpolation ofvin the space span{ϕ_N^(m)}_m=1^N, with the representation

    v_N(t,·)=∑_ℓ=1^N ϕ^(ℓ)_N(t)v(t_ℓ,·),

thenv_N(t_ℓ,·)=v(t_ℓ,·)∈ S(W)forℓ=1,2,…,N. Then the definite integration 

    ∫_A_0^A_1v_N(t,·)ṭ=∑_ℓ=1^N v(t_ℓ,·)∫_A_0^A_1ϕ^(ℓ)_N(t)ṭ=A_1-A_0/N∑_j=1^N v(t_j,·),

which coincides with the formula of the trapezoidal rule. 

We will study the error estimate of the trapezoidal rule of the integration ofv. Let's begin with the Fourier series of the functionv(t,·).  A classical Minkowski integral inequality is needed, see Theorem 202, <cit.>.


Suppose (S_1,μ_1) and (S_2,μ_2) are two measure spaces and F: S_1× S_2→ is measurable. Then the following inequality holds for any p≥ 1
    [∫_S_2|∫_S_1F(y,z)μ̣_1(y)|^pμ̣_2(z)]^1/p≤∫_S_1(∫_S_2|F(y,z)|^pμ̣_2(z))^1/pμ̣_1(y)
Suppose v belongs to ∈ C_p^2m([A_0,A_1];S_g(t)(W))∩ C^∞([A_0,A_1];S_g(t)(W)). If h=O(t^m+1) as t→ 0^+ for some positive integer m, then v has the Fourier series with respect to t has the form
 
    v(t,x)=∑_j∈v_j(x)e^ı j t2π/A_1-A_0,

where v_j∈ S(W) and its norm for j∈^+
    v_j_S(W)≤|A_1-A_0/2π j|^2m+1v_C^2m+1([A_0,A_1];S(W)).
First, consider the case that S(W)=L^2(Ω^Λ_H).
As v(t,·)∈ C_p^2m([A_0,A_1];S_g(t)(W))∩ C^∞([A_0,A_1];S_g(t)(W)),  for any j=0,…,2m, ∂^j/∂ t^jv(t,x)∈ L^2([A_0,A_1];S(W)). So the Fourier series of v(t,x) with respect to t exists.  From Theorem <ref>, the function v∈ C_p^2m([A_0,A_1];S(W)). By integration by part, for j≠ 0
    v_j(x)   =1/A_1-A_0∫_A_0^A_1 v(t,x) e^-ı j t2π/A_1-A_0ṭ
       =1/A_1-A_0(A_1-A_0/2ıπ j )^2m+1∫_A_0^A_1d^2m+1/d t^2m+1v(t,x)e^-ı j t2π/A_1-A_0ṭ,
 
then from Lemma <ref>,

    ∫_A_0^A_1d^2m+1/d t^2m+1v(t,x)e^-ı j t2π/A_1-A_0ṭ_S(W)
    
    =   (∫_Ω^Λ_H|∫_A_0^A_1d^2m+1/d t^2m+1v(t,x)e^-ı j t2π/A_1-A_0ṭ|^2x̣)^1/2
    ≤   ∫_A_0^A_1(∫_Ω^Λ_H|d^2m+1/d t^2m+1v(t,x)e^-ı j t2π/A_1-A_0|^2x̣)^1/2ṭ
    ≤   ∫_A_0^A_1d^2m+1/dt^2m+1v(t,·)_S(W)ṭ≤ (A_1-A_0)v_C^2m+1([A_0,A_1];S(W)),

then

    v_j(x)_S(W)≤|A_1-A_0/2π j|^2m+1v_C^2m+1([A_0,A_1];S(W)).
The case that S(W)=H^n(Ω^Λ_H) for any n∈ could be proved similarly, thus we omit the proof here. The proof is finished.


A direct corollary of the lemma shows the case thath(t)=o(t^n)for anyn∈. 

If h=o(t^n) as t→ 0^+ for any positive integer n, then v has the Fourier series with respect to t has the form
 
    v(t,x)=∑_j∈v_j(x)e^ı j t2π/A_1-A_0,

where v_j∈ S(W) and its norm

    v_j_S(W)=o(j^-2n).


Then we can get the error between the exact functionvand its Fourier interpolationv_Ndefined in (<ref>) forα∈, whenNis large enough.


If v∈∈ C_p^2n([A_0,A_1];S_g(t)(W))∩ C^∞([A_0,A_1];S_g(t)(W)) for some n∈, then the  difference between v(t,·) and its approximation v_N(t,·) is bounded by

    v-v_N_L^2([A_0,A_1];S(W))≤CN^-2n-1/2v_C_p^2m([A_0,A_1];S(W)).
First, let S(W)=L^2(Ω^Λ_H).
From the definition of v_N(t,·), i.e., the projection of v into the subspace span{c_n(x)exp[ı n2π/A_1-A_0], n=-N/2+1,⋯,N/2, c_n∈ S(W)},

    v- v_N^2_L^2([A_0,A_1];S(W))
    =   ∫_A_0^A_1∫_Ω^Λ_H|v-v_N|^2x̣ṭ
    
    =   ∫_A_0^A_1∫_Ω^Λ_H|∑_j∈∖[-N/2+1,N/2]v_l e^ı j t2π/A_1-A_0|^2x̣ṭ
    ≤    (A_1-A_0)∑_j∈∖[-N/2+1,N/2]v_j_S(W)^2.

With the result of Lemma <ref>,

    v-v_N^2_L^2([A_0,A_1],S(W))   ≤ (A_1-A_0)∑_j∈∖[-N/2+1,N/2]|A_1-A_0/2π j|^4n+2v^2_C^2n+1([A_0,A_1];S(W))
       ≤ CN^-4n-1v^2_C_p^2m([A_0,A_1];S(W)).
The case that S(W)=H^n(Ω^Λ_H) for any n∈ could be proved similarly, thus we omit the proof here.
 The proof is finished.





The error estimation of the trapezoidal rule will be shown in the following theorem.



 If v∈  C_p^2n([A_0,A_1];S_g(t)(W))∩ C^∞([A_0,A_1];S_g(t)(W)) for some n∈, then the error of the numerical integration of S is bounded by

    ∫_A_0^A_1(v-v_N)(t,·)ṭ_S(W)≤CN^-2n-1/2v_C^2n+1([A_0,A_1];S(W)).
First, let S(W)=L^2(Ω^Λ_H).
With the result of Theorem <ref> and Theorem <ref>,

    ∫_A_0^A_1(v- v_N)(t,·)ṭ^2_S(W)
    =   ∫_Ω^Λ_H|∫_A_0^A_1(v-v_N)(t,·)ṭ|^2x̣
    ≤   ∫_A_0^A_1v-v_N^2_L^2(Ω^Λ_H)ṭ
    ≤   v-v_N^2_L^2([A_0,A_1];S(W)).

Then 

    ∫_A_0^A_1(v-v_N)(t,·)ṭ_S(W)   ≤v-v_N_L^2([A_0,A_1],S(W))
       ≤CN^-2n-1/2v_C^2n+1([A_0,A_1];S(W)).
The case that S(W)=H^n(Ω^Λ_H) for any n∈ could be proved similarly, thus we omit the proof here.
then the proof is finished.



The result could be applied to the error estimation of  the numerical  inverse Bloch transform (<ref>).  

Given a positive integer N,   define the nodal points and inverse Bloch trasnform for the two conditions.

  *  If =[a_0,a_1], then let t_j=a_0+j(a_1-a_0)/N, then the inverse Bloch transform

    (^-1_Ω.Nw)(x)=a_1-a_0/N∑_j=1^N v(g(t_j),x)g'(t_j)

  *  if =[a_0,a_2], then let t_j=a_0+j(a_1-a_0)/N, t_j+N=a_1+j(a_2-a_1)/N, then the inverse Bloch transform

    (^-1_Ω.Nw)(x)=a_1-a_0/N∑_j=1^N v(g(t_j),x)g'(t_j)+a_2-a_1/N∑_j=N+1^2N v(g(t_j),x)g'(t_j).

If h(t)=O(h^m+1) for some m∈,  the error of the numerical inverse Bloch transform is bounded by

    ^-1_Ω,Nw-^-1_Ω w_S(W)≤CN^-2m-1/2v_C^2m+1([A_0,A_1];S(W)),  ℓ=1,2;

if h(t)=o(h^n+1) for any n∈,   the error of the numerical inverse Bloch transform is bounded by

    ^-1_Ω,Nw-^-1_Ω w_S(W)≤CN^-2n-1/2v_C^2n+1([A_0,A_1];S(W)),  ℓ=1,2,

for any n∈.


§ FINITE ELEMENT METHOD FOR THE SCATTERING PROBLEMS


In this section, we assume that the incident fieldu^isatisfies Assumption <ref>. With the transform thatα=g(t)in, the transformed inverse Bloch transform 

    (^-1_Ω v)(
    x_1+Λ j
    
    x_2
    )=[Λ/2π]^1/2∫_ v(t,x)e^-ı g(t)Λ jṭ   x∈Ω^Λ_H.
Replace α by g(t) in (<ref>), then the function v(t,x) satisfies the following equation for any ϕ∈ L^2(;H^1_g(t)(Ω^Λ_H))
    ∫_a_g(t)(v(t,·),ϕ(t,·))ṭ+ b(^-1_Ω v, ^-1_Ω[ϕ(t,·)g'(t)])=∫_∫_Γ^Λ_HF(t,·)ϕ(t,·)ṭ

whereF(t,·)=f(g(t),·)g'(t). 

Define the  space L^2(;H^s_g(t)(Ω^Λ_H);g) by 

    L^2(;H^s_g(t)(Ω^Λ_H);g):={ϕ(t,·)=ψ(g(t),·)g'(t); ψ(α,·)∈ L^2(;H^s_α(Ω^Λ_H))}.

Then define the norm in L^2(;H^s_g(t)(Ω^Λ_H);g) by

    ϕ^2_L^2(;H^s_g(t)(Ω^Λ_H);g)=ψ^2_L^2(;H^s_α(Ω^Λ_H))=∫_ψ(g(t),·)^2_H^s_g(t)(Ω^Λ_H)g'(t)ṭ.

We can also define the function space H_0^r(;H^s_g(t)(Ω^Λ_H);g) in the same way. 
By replacing ϕ(t,·)g'(t) with ψ(t,·)∈ L^2(;H^s_g(t)(Ω^Λ_H);g), the variational problem (<ref>) is equivalent to

    ∫_a_g(t)(v(t,·),ψ(t,·))[g'(t)]^-1ṭ+ b(^-1_Ω v, ^-1_Ωψ)=∫_∫_Γ^Λ_HF(t,·)ψ(t,·)[g'(t)]^-1ṭ.
 
The transformed inverse Bloch transform J^-1_Ω is also defined in H_0^r(;H^s_g(t)(Ω^Λ_H);g), we will show that the operator is bounded. Take v∈ L^2(;L^2(Ω^Λ_H);g),

    _Ω^-1v^2_L^2(Ω^H)=
    _Ω^-1 w^2_L^2(Ω^H)≤w^2_L^2(;L^2(Ω^Λ_H))=v^2_L^2(;L^2(Ω^Λ_H);g).

We can also prove that ^-1_Ω is a bounded operator from H_0^r(;H^s_g(t)(Ω^Λ_H);g) to H_r^s(Ω^H).


From the procedure in the arguments above, the following equivalence result is obtained by simply change of variables.

Suppose the incident field u^i∈ H_r^1(Ω^H_p) for some r∈[0,1) and the functions ζ, ζ_p are Lipschitz continuous functions. Then w=_Ω u_T∈ H_0^r(;H^1_α(Ω_H^Λ)) satisfies  variational problem (<ref>) if and only if v∈ H_0^r(;H^1_α(Ω_H^Λ);g) is the solution to (<ref>).



From the equivalence between (<ref>) and (<ref>), the following result is a simple corollary of Theorem <ref> and Theorem <ref>.

The variational problem (<ref>) is uniquely solvable in H_0^r(;H^1_α(Ω_H^Λ);g).  Further more, if u^i satisfies Assumption <ref>,  ζ, ζ_p∈ C^2,1(), when h(t)=O(t^m+1), the unique solution v to (<ref>) belongs to the space C_p^2n(;H^2_g(t)(Ω^Λ_H)). 

We only need to prove the unique solvability. Suppose v∈ H_0^r(;H^s_g(t)(Ω^Λ_H);g) is the solution to (<ref>) for F=0, from the definition of the space H_0^r(;H^s_g(t)(Ω^Λ_H);g), there is a w∈ H_0^r(;H^1_α(Ω^Λ_H)) such that v(t,·)=w(g(t),·)g'(t). From Lemma <ref>, w(α,·) is the unique solution to (<ref>) with f=0. From Theorem <ref>, the variational problem (<ref>) has at most one solution, so w=0, then v=0. Then the injectivity is proved. From the Bounded Inverse Theorem, the problem is uniquely solvable. 

When u^i satisfies Assumption <ref> and ζ, ζ_p∈ C^2,1(), the solution w to (<ref>) belongs to the space H_0^r(;H^2_α(Ω^Λ_H)) and belongs to 𝒜_c^ω(;H^2_α(Ω^Λ_H);), then from Corollary <ref>w satisfies Assumption <ref>. From Theorem <ref>, when h=O(t^m+1), v belongs to the space C_p^2n([a_0,a_1];H^2_g(t)(Ω^Λ_H)) for any a_0, a_1∈ such that (a_0,a_1)∩. From the definition of the space C_p^2n(;H^2_g(t)(Ω^Λ_H)) and the new definition of ,  v∈ C_p^2n(;H^2_g(t)(Ω^Λ_H))⊂ H_0^2n(;H^2_g(t)(Ω^Λ_H)). 


We will build up the finite element space, take=[a_0,a_1]for example. Let the uniformly located grid points be{t_j}_j=1^Nand Fourier basic functionsϕ^(j)_Nbe defined as in (<ref>)-(<ref>). For notations defined inΩ^Λ_H, recall the triangular meshℳ_H, the piecewise linear basic functions{ψ^(ℓ)_M}_ℓ=1^Mand the finite dimensional subspaceV_0= span{ϕ^(ℓ)_M}_ℓ=1^M⊂ H_0^1(Ω^Λ_H). Define the finite element spaceX_N,hby

    X_N,h={v_N,h(t,x)=∑_j=1^N∑_l=1^M v^(j,ℓ)_N,hexp(-ı g(t_j)x_1)ϕ^j_N(t)ψ^ℓ_M(x)}.

If v∈  C_p^2n([a_0,a_1];H^2_g(t)(Ω^Λ_H)), then the error between the approximation of v in X_N,h, i.e., v_N, and v, is bounded by

    min_v_N,h∈X_N,hv-v_N,h_L^2([A_0,A_1];H^1_g(t)(Ω^Λ_H))≤C(N^-2n-1/2+h)v_C^2n+1([A_0,A_1];H^2_g(t)(Ω^Λ_H)),

for large enough N and small enough h>0, where C depends on n and g.

From Theorem <ref>, v_N(t,x)=∑_j=1^Nϕ_N^(j)(t)v(t_j,x) is an approximation of v with the error 

    v-v_N_L^2(;H_g(t)^1(Ω^Λ_H))≤CN^-2m-1/2v_C^2m+1([A_0,A_1];H_g(t)^1(Ω^Λ_H)).
 

For each t_j, the function exp(ı g(t_j)x_1)v(t_j,x) is approximated by functions in the subspace V_0= span{ψ^(ℓ)_M, ℓ=1,…,M}, i.e., there is a series of coefficients v^(j,ℓ)_N,h∈ such that

    v(t_j,x)=∑_ℓ=1^M v^(j,ℓ)_N,hexp(-ı g(t_j)x_1)ψ^ℓ_M(x)

for any x in the grid points of ℳ_h, with the error estimate

    v(t_j,·)-∑_ℓ=1^M v^(j,ℓ)_N,hexp(-ı g(t_j)x_1)ψ^ℓ_M(·)_H^1_g(t)(Ω^Λ_H)   ≤ C hv(t_j,·)_H^2_g(t)(Ω^Λ_H)
       ≤ C h v_C_p^2n([a_0,a_1];H^2_g(t)(Ω^Λ_H)).

For v_N,h=∑_j=1^N∑_ℓ=1^M v^(j,ℓ)_N,hexp(-ı g(t_j)x_1)ϕ_N^(j)(t)ψ^ℓ_M(x), and let v_j^h(x):=v(t_j,x)-∑_ℓ=1^M v_N,h^(j,ℓ)exp(-ı g(t_j)x_1)ψ^ℓ_M(x), then v_j^h_H^1_g(t)(Ω^Λ_H)≤ C h v_H_0^2n([a_0,a_1];H^2_g(t)(Ω^Λ_H)). With the result of Lemma <ref>,

    v_N-v_N,h^2_L^2(;L^2(Ω^Λ_H))
    
    =    ∫_∫_Ω^Λ_H|∑_j=1^Nϕ_N^(j)(t)v(t_j,x)-∑_j=1^N∑_ℓ=1^M v^(j,ℓ)_N,hexp(-ı g(t_j)x_1)ϕ_N^(j)(t)ψ^ℓ_M(x)|^2x̣ṭ
    
    =   ∫_∫_Ω^Λ_H|∑_j=1^Nϕ_N^(j)(t)v^h_j(x)|^2x̣ṭ
    ≤   ∫_∫_Ω^Λ_H∑_j,j'=1^N ϕ_N^(j)(t)ϕ_N^(j')(t)v^h_j(x)v^h_j'(x)x̣ṭ
    
    =   a_1-a_0/N∑_j=1^Nv^h_j^2_L^2(Ω^Λ_H).

Similarly, we can have the H^1-estimate 

    v_N-v_N,h^2_L^2([a_0,a_1];H^1_g(t)(Ω^Λ_H))   ≤a_1-a_0/N∑_j=1^Nv^h_j^2_H^1_g(t)(Ω^Λ_H)
       ≤ C^2[a_0-a_1]h^2v^2_C_p^2n([a_0,a_1];H^2_g(t)(Ω^Λ_H))

With the inequalities (<ref>) and (<ref>), 

    v-v_N,h_L^2(;H^1_g(t)(Ω^Λ_H))   ≤v-v_N_L^2(;H^1_g(t)(Ω^Λ_H))+v_N-v_N,h_L^2(;H^1_g(t)(Ω^Λ_H))
       ≤CN^-2n-1/2v_C^2n+1([A_0,A_1];H^1(Ω^Λ_H))+Ch v_C_p^2n([A_0,A_1];H^2(Ω^Λ_H))
       ≤ C(N^-2n-1/2+h)v_C^2n+1([A_0,A_1];H^2_g(t)(Ω^Λ_H)).

where C depends on n and g.


With this results, we can estimate the error of the approximation of the Bloch transformed total fieldwby functions inX_N,h, when=[a_0,a_1].

Suppose u^i satisfies Assumption <ref>, g satisfies Assumption <ref>, ζ, ζ_p∈ C^2,1(). u is the total field with the incident field u^i and surface Γ_p, and  u_T:=u∘Φ_p, and define w is the Bloch transform of u_T, i.e., w:=_Ω u_T. Let v be defined by (<ref>), then v satisfies the variational equation (<ref>), for any ϕ∈ L^2(;H^1_α(Ω^Λ_H)). Moreover, if h=O(t^m+1) for some m∈ as t→0^+, the approximation of v in X_N,h is bounded by

    min_v_N,h∈X_N,hv-v_N,h_L^2(;H^1(Ω^Λ_H))≤ C(N^-2m-1/2+h)v_C^2m+1([A_0,A_1];H^2(Ω^Λ_H));

if h=o(t^n+1) for any n∈ as t→0^+,  the approximation of v in X_N,h is bounded by

    min_v_N,h∈X_N,hv-v_N,h_L^2(;H^1(Ω^Λ_H))≤ C(N^-2n-1/2+h)v_C^2n+1([A_0,A_1];H^2(Ω^Λ_H)),

for any n∈.

As u^i satisfies Assumption <ref> and ζ, ζ_p∈ C^2,1(), from Corollary <ref>, for the interval =[a_0,a_1], there are three functions w_0, w_1, w_2∈ C^∞(;H^2_α(Ω^Λ_H)) such that

    w=w_2+√(α-A_0)  w_0+√(α-A_1) w_1.

So w also satisfies the conditions in Assumption <ref>. From the definition

    v(t,·)=w(g(t),·)g'(t).

From Theorem <ref>, there are two different cases. 

1) If h(t)=O(t^m+1) as t→ 0^+ for some m∈, v∈ C_p^2m(;S(W))∩ C^2m+1(;S(W)), then from Theorem <ref>, the approximation of v in X_N,h is bounded by

    min_v_N,h∈X_N,hv-v_N,h_;H^1(Ω^Λ_H))≤ C(N^-2m-1/2+h)v_C^2m+1([A_0,A_1];H^2(Ω^Λ_H)).



2) If h(t)=o(t^n+1) as t→ 0^+ for any n∈, v∈ C_p^∞(;S(W)), then from Theorem <ref>, the approximation of v in X_N,h is bounded by

    min_v_N,h∈X_N,hv-v_N,h_L^2(;H^1(Ω^Λ_H))≤ C(N^-2n-1/2+h)v_C^2n+1([A_0,A_1];H^2(Ω^Λ_H))

for any n∈.


Similar settings and results could be obtained for=[a_0,a_2], with nodal pointst_jare chosen in the same way as Theorem <ref>. Following the proof of Theorem 9 in <cit.>, the error estimate of the finite element method is obtained.

Assume that u^i∈ H_r^2(Ω_H) for r≥ 1/2 satisfies Assumption <ref>, ζ, ζ_p∈ C^2,1(), then the linear system (<ref>) is uniquely solvable in X_N,h for large enough N and small enough h>0. 

1) If h(t)=O(t^m+1) for some positive integer m, then the solution w_N,h∈X_N,h satisfies the error estimate
 
    v-v_N,h_L^2(;H^ℓ(Ω^Λ_H))≤ C h^1-ℓ(N^-2m-1/2+h)v_C^2m+1([A_0,A_1];H^2(Ω^Λ_H)),  ℓ=0,1,

where C depends on m,u,g.

2)If h(t)=o(t^n+1) for any positive integer n, then the solution w_N,h∈X_N,h satisfies the error estimate
 
    v-v_N,h_L^2(;H^ℓ(Ω^Λ_H))≤ C h^1-ℓ(N^-2n-1/2+h)v_C^2n+1([A_0,A_1];H^2(Ω^Λ_H)),  ℓ=0,1,

for any n∈, where C depends on n,u,g.


§ NUMERICAL RESULTS


In this section, we will show some numerical results with the finite element method introduced in this paper. The periodic function is defined by

    ζ(t)=1.9+sin (t)/3-cos(2t)/4

with the local perturbation defined by

    ζ_p(t)-ζ(t)=exp(1/t^2-1)sin[π(t+1)],   where t∈[-1,1];
    
    0,   otherwise.

In all of the numerical examples, we fixH=4, and the periodΛ=2π, thusΛ^*=1.





We will show the results for two groups. For each group, we consider two different choices of the integration contour, i.e., the functiong. 

 Choice 1.g(t)is defined in[A_0,A_1]by  

    g_1(t)=A_0+A_1-A_0/∫_A_0^A_1 (s-A_0)^3(s-A_1)^3ṣ[∫_A_0^t (s-A_0)^3(s-A_1)^3ṣ],

thush(t)=O(t^2)ast→ 0^+. 

 Choice 2.g(t)is defined in[A_0,A_1]by 

    g_2(t)=A_0+A_1-A_0/∫_A_0^A_1exp(1/(s-A_0)(s-A_1))ṣ[∫_A_0^t exp(1/(s-A_0)(s-A_1))ṣ],

thush(t)=o(t^n)ast→ 0^+for anyn∈. 





Each group contains four different examples. In the four examples, two of them are with non-perturbed surfaceΓand another two are with locally perturbed surfaceΓ_p. The incident field is defined by a half-space Green's function, that is defined as

    u^i_G(x,y)=ı/4[H_0^(1)(k|x-y|)-H_0^(1)(k|x-y'|)]

wherey=(y_1,y_2)^⊤andy'=(y_1,-y_2)^⊤, or the upward propagating Herglotz wave function

    u^i_+(x)=∫_-π/2^π/2exp[ı k(sinθ x_1+cosθ x_2)]cos^2(θ)θ̣,

or the downward propagating Herglotz wave function

    u^i_-(x)=∫_-π/2^π/2exp[ı k(sinθ x_1-cosθ x_2)]cos^2(θ)θ̣.
 Group 1. In this group, we show four upward propagating incident fields scattered by both perturbed and non-perturbed surface. As the upward propagating incident fields satisfy the Helmholtz equation above the surface and also satisfy the radiation condition, from the boundary condition that u^s=-u^i and the unique solvability of the direct scattering problem, u^s=-u^i above the (perturbed) surface. The parameterhis fixed to be0.01. TheL^2-relative errors, defined byu^s_N,h-u^s_L^2(Γ^Λ_H)/u^s_L^2(Γ^Λ_H), forN=4,8,16,32,64in Table <ref> forg=g_1and in Table <ref> forg=g_2.

  *  The surface is perturbed. The incident field is the half-space Green's function u^i_G(·,y) with y=(0.5,0.4)^⊤,  k=1.

  *  The surface is not perturbed. The incident field is the half-space Green's function u^i_G(·,y) withy=(-1,0.4)^⊤, k=1.

  *  The surface is perturbed. The incident field is the  upward propagating Herglotz wave function u^i_+, k=√(2)
  *   The surface is not perturbed. The incident field is the upward propagating Herglotz wave function u^i_+, k=1.5.
 Group 2. In this group, we show four downward propagating incident fields scattered by both perturbed and non-perturbed surface. As the exact value of the scattered field could not be written out, we set the solution with N=256 be the "exact value", and show the relative errors forN=4,8,16,32,64in Table <ref> forg=g_1and in Table <ref> forg=g_2. Figure <ref> plots in logarithmic scale of the relativeL^2-errors for the examples, see (a) forg=g_1and (b) forg=g_2. The mesh sizehis fixed to be0.02.

  *  The surface is perturbed. The incident field is the half-space Green's function u^i_G(·,y) with y=(0.5,3)^⊤, k=1.

  *  The surface is non-perturbed. The incident field is the half-space Green's function u^i_G(·,y) with y=(-1,3)^⊤, k=1.5.

  *  The surface is perturbed. The incident field is the downward propagating Herglotz wave function u^i_-, k=√(2).

  *  The surface is non-perturbed.  The incident field is the downward propagating Herglotz wave function u^i_-, k=2.01.


First let's go to the examples in Choice 1. 
In Table <ref> and <ref>, the relative error decays fast at first, and then stays at the level of10^-5asNgets larger (N≥ 16). As we have shown, theL^2-relative error is bounded byC(N^-2n-1.2+h)h, whenhis small enough, the error is brought byNis comparatively larger, then whenNis relatively small, the error decays fast whenNincreases. However, whenhis not small engough, the dominant part of error is brought byh, then the error will stay at a certain level whenNincreases. Note that the first two examples Group 1 were shown in <cit.> and <cit.>, as we list in Table <ref>.


Compared to the first and second columns in Table <ref> and Table <ref>, even the result produced by our new method in this  paper atN=16is much better than that produced by the old method atN=320. This means during the computational procession, we can save a lot of  time and memory space in setting up the matrix and solving the linear system <ref>.


Now let's go to the examples in Choice 2. 
In Table <ref> and <ref>, the error decays asNgets larger, and the logarithmic scale of theL^2-errors have been plotted in Figure <ref>. 
Figure <ref> shows the convergences of the new method for bothg=g_1andg_2. In (a), the error decreases at the rate that is aboutO(N^-4). Form=1, the expected error isO(N^-2m-1/2=O(N^-2-1/2), thus the error decays faster than expected. In (b), the errors decrease much more faster, thus it shows the super algebraic convergence rate of the numerical method.  The results in the two tables are good illustrations for the result in Theorem <ref>. With these results shown above, we are confident to say that the new method is very efficient to solve the (locally-perturbed) periodic scattering problems numerically.




§ DEDICATION.
 This paper is dedicated to  Prof. Dr. Armin Lechleiter, who passed away in January 2018 at a very young age. He was a talented mathematician, a helpful colleague, and a nice friend. His mathematical theory and numerical analysis for the
scattering problems from locally perturbed periodic surfaces, which have been cited several times, are essential for the work in this paper. The author is very grateful for his scientific contribution to this area and his advise during their collaboration. He will always be missed.




§ ACKNOWLEDGEMENTS.
 The author was supported by the University of Bremen and the European Union FP7 COFUND under grant agreement n^∘ 600411.


alpha