Numerical method for scattering problems in periodic waveguides

In this paper, we propose a new numerical method for scattering problems in periodic waveguide, based on the newly established contour integral representation of solutions in a previous paper by the author (see [Zhadf]). For this kind of problems, solutions are obtained via the Limiting Absorption Principle and we all them LAP solutions. Based on the Floquet-Bloch transform and analytic Fredholm theory, an LAP solution could be written explicitly as an integral of quasi-periodic solutions on a contour, which depends on the periodic structure. Compared to previous numerical methods for this kind of problems, we do not need to adopt the LAP for approximation, thus a standard error estimation is easily carried out. Based on this method, we also develop a numerical solver for the halfguide problems. This method is also based on the result from [Zhadf], which shows that any LAP solution of a halfguide problem could be extended into a fullguide problem with a non-vanishing source term(not uniquely!). So first we approximate the source term from the boundary data by a regularization method, and then the LAP solution could be obtained from the corresbonding fullguide problem. At the end of this paper, we also show some numerical results to present the efficiency of our numerical methods.


Introduction
Numerical simulation for scattering problems in periodic waveguides is an interesting topic, due to its wide applications in optics, nanotechnology, etc. Due to the existence of eigenvalues, the unique solvability is not guaranteed thus the Limiting Absorption Principle is adopted. That is, the solution of the original problem is assumed to be the limit of a family of solutions with absorbing material, where the absorption parameter tends to 0. It is proved that the LAP holds for planar waveguides filled up with periodic material, and we refer to [Hoa11,FJ15,KL18a,KL18b] for different proofs.
In recent years, some numerical methods have been introduced to solve this kind of problems and most of the methods are based on the LAP. As the problem with positive absorption parameter is uniquely solvable and the unique solution decays exponentially at the infinity, the solution is also easily computed numerically. We refer to [EHZ09] based on a so-called recursive doubling process to approximate Robinto-Robin maps on left-and right-boundary of a periodic cell. Thus it is natural to approximate the original solution by the one when the positive absorption parameter is small enough. Based on this idea, the solution is approximated by the extrapolation with repsect to small absorption parameters in [ESZ09], and the method is also extended to scattering problems in locally perturbed periodic layers in [SZ09]. On the other hand, the LAP could also be applied to the choice of "proper" propagation modes. With the knowledge of the asymstotic behaviour of propagating modes when the absorption paramter tends to 0, it is possible to pick up certain modes that are propagating either to the left or to the right. With these propagating modes, the Dirichlet-to-Neumann maps on the boundaries of periodic cells could be approximated and the solution for the whole problem is then obtained. For details we refer to [JLF06,FJ09]. For other numerical methods we also refer to [LHL07,ALMBLA13,DS18].
In a recent paper by the author, the LAP solution is written in a simplified form, by the contour integral of quasi-periodic solutions with quasi-periodicities on a closed contour in a neighbourhood of the unit circle, which is decided only by the periodic structure. However, to guarantee that the simplified representation holds, it is required that the wave number should be any real number except for a countable subset. The simplified representation also provides an efficient method to solve the problems numerically. As the solutions for quasi-periodic problems are classic, the only difficulty is to find a proper contour for certain strucuture. For an observation of the dispersion map, the unit circle may or may not contain eigenvalues corresbonding to propagating modes. The eigenvalues are exactly poles of quasi-periodic solutions with respect to quasi-periodicities. When there is no eigenvalue on the unit circle, then the quasi-periodic problems are uniquely solvable and the curve is simply chosen as the unit circle. When there are eigenvalues on the unit circle, we may modify the it such that there is no eigenvalue on the new curve and the LAP holds for the integral representation on this curve. The curve could always be chosen as a piecewise analytic one as there are only finite number of eigenvalues on the unit circle. Finally, as the solution depends piecewise analytically on the curve, a high order numerical method is designed based on the contour integral.
The numerical method is then extended to halfguide problems. As any LAP solution of a halfguide problem is (not uniquely) extended to a solution of a fullguide problem with a non-vanishing source term, the numerical method contains two steps. The first step is to find out the source term that the fullguide problems approximates the halfguide problem, for given boundary data. We compute the boundary values for basis functions in the space that the source term lies in, and then find out the corresbonding coefficients by the least square method. As the choice of source term is not unique and the operator is compact, a Tikhonov regularization technique is adopted. With this source term, we can go to the second step, i.e., to solve the fullguide problem with the method developed in the same paper and approximate the solution of the halfguide problem.
The rest of this paper is organised as follows. We recall the integral representation for LAP solutions in the second section, and develop the high order numerical method on the curve in the third section. Then a finite element method for the fullguide problems is studied in Section 4. Then we extend the method to halfguid problem in Section 5. In Section 6, we present several numerical results of our proposed methods.

Solutions via limiting absorption principle
Let the waveguide Ω := R × [0, 1] be filled up with periodic material, and the upper and lower boundaries of Ω are denoted by Σ + = R × {1} and Σ − = R × {0}. We consider the solution of the following equation The source term f ∈ L 2 (Ω) is compactly supported and the function q is a real-valued periodic refractive index satisfies q(x 1 + n, x 2 ) = q(x 1 , x 2 ), x = (x 1 , x 2 ) ∈ Ω.
Then Ω = ∪ j∈Z Ω j . Define Γ j := {−1/2 + j} × Σ for any j ∈ Z, thus the left and right boundaries of the cell Ω j are Γ j and Γ j+1 . Let ∂ Ω j := ∂ Ω ∩ Ω j . For simplicity, we assume that Due to existence of eigenvalues, the uniqueness of the solution of (1) is not guaranteed, thus we adopt the Limiting Absorption Principle (LAP). That is, if k 2 is replaced by k 2 + iε for some ε > 0 in (1), and consider the solution u ε that satisfies From Lax-Milgram theorem, there is a unique solution u ε ∈ H 1 (Ω). As is proved in [Hoa11,FJ15,KL18b,Zhadf], the limit of u ε exists in H 1 loc (Ω) as ε → 0 + , and it is called the LAP solution through this paper.

Floquet-Bloch transform
A powerful tool to solve problems defined in periodic structures is the Floquet-Bloch transform. This Fourier-type transform is defined for the function ϕ ∈ C ∞ 0 (Ω) by the infinite series: The series is actually finite as ϕ has a compact support. However, the transform is extended to more generalized function spaces. First, we define some notations and spaces that will be useful in the rest of this paper. For any complex number z ∈ C \ {0}, define the z-quasi-periodic space H s z (Ω 0 ) (s ∈ R) by the subspace of H s (Ω) such that all the functions in this space can be extended into an H s loc (Ω) function with the help of the following extension operator: Especially, a function ϕ that belongs to H 1 z (Ω 0 ) should satisfy the following conditions on the boundary: To introduce the properties of the Floquet-Bloch transform more conveniently, we give the following definitions. A function ψ decays exponentially with the rate γ means that there is two constants γ > 0 and C > 0 such that |ψ(x 1 , x 2 )| ≤ C exp(−γ|x 1 |), ∀x 2 ∈ [0, 1].
Let S 1 = {z ∈ C : |z| = 1} be the unit circle, we can also define the Hilbert space L 2 (S 1 ; H s z (Ω 0 )) by is an open domain. If for any fixed z ∈ U , ϕ(z, ·) ∈ H s (Ω 0 ), then it depends analytically on z in U if for any z 0 ∈ U , there is a series {ϕ j : With these definitions, we can introduce the following property of the Floquet-Bloch transform, for details we refer to [LZ17b,Kuc93].
Theorem 1. The operator F has the following properties when z lies on (in the neighbourhood of ) S 1 : • For any s ∈ R, F is an isomorphism between H s (Ω) and L 2 (S 1 ; H s z (Ω 0 )); • (F ϕ)(z, ·) depends analytically on z in a neighbourhood of S 1 if and only if ϕ decays exponentially at the infinity.
Moreover, the inverse operator F −1 has the following representation: When k 2 is not an eigenvalue of the operator −q −1 ∆ defined in {ϕ ∈ L 2 (Ω) : ∆ϕ ∈ L 2 (Ω)}, then (1) is uniquely solvable in H 1 (Ω) and the solution decays exponentially when |x| → ∞. Especially, when k 2 is replaced by k 2 + iε for any k > 0 and ε > 0, the unique solution u ε ∈ H 1 (Ω) of (3) decays exponentially. For either of the two cases, we apply the Floquet-Bloch transform and let Then both w(z, ·) and w ε (z, ·) belongs to H 1 z (Ω 0 ) for any fixed z in its region of convergence, and depends analytically on z in this domain. Moreover, the functions satisfy the follow equations: Note that as f has a compact support in Ω 0 , F f = f for any z ∈ C \ {0}.
To consider the relationship between u and w, we have to introduce an important result from the Floquet theory.
Theorem 2. k 2 is an eigenvalue of −q −1 ∆ in {ϕ ∈ L 2 (Ω) : ∆ϕ ∈ L 2 (Ω)}, if and only if k 2 is an This result implies that when k 2 (or k 2 + iε) is not an eigenvalue, the problem (4) (or (5)) is uniquely solvable. Thus the problem is well-posed, and the solution u (or u ε ) could be obtained from the inverse Floquet-Bloch transform: The unit cirlce S 1 could also be replaced by any closed path encircling the origin that lies in the region of convergence of the Floquet-Bloch transform.

Quasi-periodic problems and Bloch wave solutions
From our previous paper, we have formulated the the LAP solution via the curve integral. First, when k 2 is not an eigenvalue, the representation (6) is already good enough thus we only need to consider the case when k 2 is an eigenvalue. In this case, from Theorem 2, there is at least one z ∈ S 1 such that (5) is not uniquely solvable. To study this problem, we have to introduce the following definition.
Definition 3. The complex nubmer z is a Floquet multiplier if (4) has a nontrivial solution with f = 0 in H 1 z (Ω 0 ), and the corresbonding solution is called a Bloch wave solution. Let the set of all Floquet multipliers be denoted by F, then it has the following properties.
• F is a discrete subset of C \ {0}; • z ∈ F if and only if z −1 ∈ F; • F ∩ S 1 is a finite set; • The only accumulation points in F are 0 and ∞.
Let v z (x) := z −x1 w(z, x) for any fixed z ∈ C \ {0}, then v z ∈ H 1 per (Ω 0 ) and satisfies the following equations Then it satisfies the following variational formulation For simplicity, we denote the left hand side by the sesquilinear form a z (v z , ϕ) defined in the space Let K z ∈ L H 1 per (Ω 0 ) be the bounded linear operator such that Thus it is compact and depends analytically on z ∈ C × . From the equivalence between the solution of (7) and (4), the problem (7) is uniquely solvable if and only if z / ∈ F. Then from analytic Fredholm theory (see [RS80]), the operator I − K z is invertible in C × except for the discrete set F. Moreover, (I − K z ) −1 is analytic in C × \ F and meromorphic in C × . With this result, we conclude the regularity property for w(z, x) with respect to z in the following theorem.
Proof. From the property of (I − K z ) −1 , v z is also analytic in C × \ F and meromorphic in C × . The only problem for w(z, x) is when z ∈ (−∞, 0). For any z 0 ∈ (−∞, 0) \ F, as w(z, x) depends analytically on z 0 in the neighbourhood of z 0 by choosing another proper branch. From the definition of F, as w(z, x) is a single-valued function for z in the small neighbourhood, it is also analytic. Thus it is analytic in C \ (F ∪ {0}). The proof is finished.
uniformly with respect to t. Thus the result also holds for w(a(t), ·).
Proof. As I − K z depends analytically on t and its inverse is bounded uniformly, v z H 2 (Ω0) ≤ C f L 2 (Ω0) for some C > 0 independent of t. Then we consider the derivative As a(t) = 0, the right hand side is bounded by For higher order derivatives, we can also have the same estimation. The proof is finished.
As F is discrete and symmetric, it is written as the union of the following sets: Moreover, z ∈ RS if and only if z ∈ LS. From Floquet theory, when the problem (1) is not uniquely solvable in H 1 (Ω), the set U S is not empty; while when the problem is uniquely solvable, U S is empty. Thus this set plays an important role in the solutions of (1). Actually, the set could also be decomposed due to the different property of corresbonding Bloch wave solutions.
For any z ∈ U S, it could be written as exp(iα) for a unique α ∈ (−π, π]. For any fixed α, there is a series of eigenvalues of −q −1 ∆ in the space H 1 exp(iα) (Ω 0 ) that tends to infinity. From [Kuc93], we can rearrange the eigenvalues in a proper way, such that they lie on a series of analytic functions µ n (α). Note that these analytic functions are not constants (see [SW02,SS02]). The graph of (α, µ n (α)) is called the dispersion diagram. We give two different examples of dispersion diagrams: 1. Example 1. q = 1 is a constant function in Ω, and its dispersion diagram is shown in Figure 7 (left).
The relationship between the quasi-periodic parameter α and the eigenvalue could be obtained analytically: 2. Example 2. q = 9 in a disk with center (0, 0.5) and radius 0.3 and q = 1 outside the disk. As the functions µ n could not be obtained analytically, the eigenvalues are computed by the finite element method based on regular triangle meshes in the periodic cell Ω 0 with the largest mesh size 0.005. We conclude the properties of the points in U S from the dispersion diagram.
Assumption 6. The set SU S is empty and the intersection RU S ∩ LU S is also empty.
From [Zhadf], there are only countable number of wavenumbers k > 0 such that neither SU S nor RU S ∩ LU S is empty. For simplicity, let

Integral formulation for LAP solutions
Let's return to the integral representation of u ε in (6) for any k > 0. To study the LAP solution, we have to consider the behaviour of poles of w ε (·, x) when ε → 0 + . For any z ± j ∈ F, let Z ± j (ε) be the corresbonding poles with k 2 + iε. First, we recal some results from [Zhadf].
Thus as for w ε (z, ·) that solves (5), it has poles at Z ± j (ε). Recall that when U S = ∅, the integrand has poles as ε → 0 + . To take the limit without having problems with poles, we choose a new curve C 0 (see Figure 3) such that it is the boundary of the following domain : . The parameter δ > 0 is chosen such that the following conditions are satisfied: • The intersection of any two balls B(z ± j , δ) and B(z ± ℓ , δ) when j = ℓ is empty.
. From the definition of the new curve, we can replace the unit circle by C 0 in (6), i.e., We take the limit ε → 0, as C 0 is away from any poles Z ± j (ε), it converges to w(z, x) in H 1 z (Ω 0 ) uniformly for any fixed z on C 0 , thus From the definition of C 0 and Corollary 5, the integrand w(z, x)z n−1 depends piecewisely analytically on z. Moreover, in each analytic curve segment, the derivatives of w(z, x) with respect to z is also uniformly bounded.
The well-posedness of LAP solutions is concluded in the following theorem. Theorem 8. Given any f ∈ L 2 (Ω 0 ) with compact support, there is a unique LAP solution in H 1 loc (Ω), and it could be represented by (8). There is a constant C = C(n) > 0 such that Moreover, from the regularity of the elliptic equation, u ∈ H 2 loc (Ω).

Approximation of curve integral 3.1 Parameterization of the integral curve
In this section, we consider the numerical approximation of the curve integral (6) or (8). w(z, ·) is supposed to be known for it is not difficult to be computed. We can always choose different n's such that the values are taken in different Ω n , but in this section n is fixed as 0 as an example. From the definition of C 0 , it is piecewise analytic thus w(z, x) is also piecewise analytic in z, and the corners of C 0 is also the singularities of w(z, x). Suppose the number of elements in RU S is P , from the symmetry of RU S and LU S, the number is the same for LU S. Thus U S = RU S ∪ LU S has 2P different elements when Assumption 6 holds. Especially, as 1 (or −1) in RU S if and only if 1 (or −1) in LU S, 1 (or −1) in RU S ∩ LU S. From Assumption 6, as RU S ∩ LU S = ∅, neither 1 nor −1 lies in U S. Then we have to figure out the exact parameter representation for the curve C 0 .
when j = 1, 2, . . . , 2P , where α − 2P +1 := α − 1 + 2π. For any j ∈ {1, 2, . . . , 2P }, C 0 ∩ ∂B exp(iα j ), δ j lies in the exterior or in the interior of S 1 , depending on either RU S or LU S that exp(iα j ) lies in. Let this curve segment be represented by exp(iα j ) + δ j exp(iθ). From direct calculation, there are two angles θ ± j such that exp(iα By choosing proper branches of the arccos function, θ ± j satisfy Thus the curve segement C 0 ∩ ∂B(exp(iα j ), δ j ) is parameterized as Thus the whole curve C 0 is now written as and is parameterized piecewisely. Thus the representation of u in (8) is written as where for any j = 1, 2, . . . , 2P , From the paramterization, the integrals could be written into the following forms: Both of the integrands are functions that depends analytically on t, and for any fixed t, it belongs to the space H 2 (Ω 0 ). Thus we only need to consider the numerical integration where a < b and g depends on x ∈ (a, b) and g(t, ·) ∈ S(W ), S(W ) is a Sobolev space defined on the open domain W .

then the integral is explicitly written as
Then we recall the result in [Zha18], and obtain the error estimate of the integral via (16).
Theorem 9 (Theorem 28, [Zha18]). Let s ∈ C N0−1 per ([−π, π]; S(W )) for some positive integer N 0 , then the error of the numerical integration is bounded by We apply the quadrature method (16) to the evaluation of I j and J j . For I j , let q j be the smooth function from [α + j , α − j+1 ] to [−π, π]. Then the approximation of I j , denoted by I N j , has the representation of We can estimate J j in the similar way. Let h j+2P be the smooth function from [θ − j , θ + j ] to [−π, π], then .
Thus u(x) is approximated by With the result of the theorem above, we can conclude that For any fixed z ∈ C 0 , the function w(z, x) could be approximated by w h (z, x). Then the LAP solution could be approximated by Then we conculde the error between u N,h and u in the domain Ω 0 .
Theorem 11. Let u N,h be the numerical solution comes from the finite element method and the integral approximation (22). Then the error is bounded by where C is a constant depnds on P and N 0 , but does not depend on N and h.
Proof. From the representations of u N,h and u, and also the results from (19) and The estimation of the H 1 −norm is also obtained: The proof is finished.

Half-guide problems
The method could also be extended to the numerical simulation of half-guide scattering problems. Let the half waveguide be defined by Ω + := ∪ ∞ n=1 Ω n (see Figure 4). Then we are looking for the LAP solution u + ∈ H 1 loc (Ω + ) such that it satisfies ∆u + + k 2 qu + = 0 in Ω + ; u + = ϕ on Γ 1 .
Let u + (ε) be the one associated with k 2 + iε for any ε > 0, then from the Lax-Milgram theorem again, there is a unique solution u + (ε) ∈ H 1 (Ω + ) that decays exponentially at the infinity. The solution u + which we are looking for, is the limit of u + (ε) as ε → 0. In the following, we also let Ω − := ∪ 0 n=−∞ Ω n for simplicity. Then Ω = Ω − ∪ Ω + . From the analysis in [Zhadf], the LAP solution is unique and there is always an f ∈ L 2 (Ω 0 ) such that Note that although the LAP solution u + is unique, the choice of f is multiple. This representation also implies that the LAP solution in Ω + could be extended to the LAP solution to a full guide problem, but the extension is not unique.
To solve this problem with the method based on (25), we have two steps. The first step is to use the property u + = ϕ on Γ 1 to determine a source term f . Then the next step is to use the source term to compute u + , with the same method as the full-guide problem. As the second step is exactly the same as the full-guide problem, we only discuss the first step in this section.
Let A be the operator defined by Thus A is a compact, surjective but non-injective operator with a dense range. We are looking for an f such that Af = ϕ on Γ 1 .
As A is compact and non-injective, we apply the Tikhonov regularization method to solve this equation. Given a small enough α > 0, we are looking for f α such that Suppose {µ n , u n , v n } is the singular system of the operator A, then µ 1 ≥ µ 2 ≥ · · · ≥ µ n ≥ · · · > 0. From the definition, Au n = µ n v n , A * v n = µ n u n .
As A is non-injective, N (A) = ∅. As A is a surjection, {v n } n∈N is a basis in H 1/2 (Γ 1 ). For any ϕ ∈ H 1/2 (Γ 1 ), there is a series {c n : n ∈ N} ⊂ C such that From direct calculation, Thus we estimate the error between Af α and ϕ.
We choose a family of basis functions for the space H 1 (Ω 0 ) and denote it by {ϕ ℓ } ∞ ℓ=1 . Then we fix a large enough integer M 0 > 0 and consider the approximation of f by f M0 Construct the matrix Φ := (Aϕ 1 , Aϕ 2 , . . . , Aϕ M0 ) , then we look for the matrix c(α) = (c 1 (α), c 2 (α), . . . , c M0 (α)) ⊤ such that We have approximated the function f by f M0 α with coefficients c(α), and then solve the equation (4) with the source term f M0 α . Finally, the LAP solution u is obtained by the integral representation (25). From the convergence of the numerical method to approximate f and solve (25), the full numerical scheme converges as M 0 → ∞, α → 0, h → 0 and N → ∞.

Numerical results
In this section, we present some numerical examples to show the efficiency of our numerical methods for both full-and half-guide problems. We take the example when the functions q and f are chosen as follows: The point a 0 = (0, 0.5), and ζ(t) is a C 8 -continuous function defined as First, we draw the dispersion diagram and the corresbonding z-space, see Figure 5. From the dispersion diagram, we can find out two stop bands, i.e., (2.956 ± 0.01, 7.574 ± 0.01) and (13.41 ± 0.01, 15.49 ± 0.01). Thus when k 2 = 5, it lies in the first stop band, thus there is no eigenvalue lies on S 1 in the z-space. In this case, u could be represented by (6) with the integral on the curve S 1 . However, when k 2 = 17, there are two points α 1 = 0.9576(±0.01) and α 2 = −0.9576(±0.01) on the dispersion curve with µ n (α) = 17, i.e., P = 1. Moreover, when α = α 1 (blue square on the left), the derivative of µ n is positive, thus exp(iα 1 ) ∈ RU S, i.e., it is related to the Bloch wave solution that propagates to the right. Similarly, when α = α 2 (green circle on the left), exp(iα 2 ) ∈ LU S, i.e., it is related to the Bloch wave solution that propagates to the left. Thus we can design the integral curve in (8) as the red curve on the right. Moreover, we also checked the condition number of the matrix obtained from the finite element method of (7), to find out a rough location of poles of the function w(z, ·) (see Figure 6).

Full-guide problems
In the first part of this section, we show some numerical examples for scattering problems in full-waveguide.
The source term f is defined as follows: With these settings, we can then calculate the value of u with different parameters. We fix N 0 = 6. For the finite element method, we choose h = 0.0025, 0.005, 0.01, 0.02 and N = 8, 16, 32, 64, 128, 256 for k 2 = 5, N = 4, 8, 16, 32, 64, 128 for k 2 = 17. We also obtain "exact solutions" for the corresbonding problem from the method introduced in [ESZ09, EHZ09] with h = 0.005 and the Lagrangian element is adopted. First we show the relative errors with different h and N for both cases, defined by where u N,h is the numerical solution with parameter N and H, and u exa is the "exact solution". The results for k 2 = 5 are shown in Figure 1 and for k 2 = 17 are shown in Figure 2. We can see that the relative erros decays as N gets larger and h gets smaller. Note that when N is large enough (e.g., N ≥ 32), the relative errors do not decay when N gets larger, this implies that the error brought by N is relatively smaller, compared with the error from h.  As the convergence rate with respect to h is classic, we are especially interested in the convergence with respect to N . We take h = 0.01 for both cases, and compute the relative error between u N,h and u 256,h for k 2 = 5, u 128,h for k 2 = 17. From the result in (17), the error is expected to decay at the rate of O(N −5.5 ). From the two pictures, the convergence is as fast as, or even faster than expected.

Half-guide problems
In the second part of this section, we also show some numerical examples for half guide problems. The boundary data is given by ϕ(x 1 , x 2 ) = sin x 2 + x 2 2 2 + exp(ix 1 ) cos(2x 2 ).
For all the examples, we approximate the source term by where α is the regularization paramter. As the results depend greatly on α, we show different results with corresbonding to differnt paramters, i.e., α = 10 −2 and 10 −5 . The "exact solutions" also comes from the method used for the full-guide problems, and the relative error err N,h is also defined in the same way. We show the results for k 2 = 5 with parameters N = 8, 16, 32 and h = 0.02, 0.01, 0.005 with α = 10 −2 in Table  3 and with α = 10 −5 in Table 4. For k 2 = 17, we show results for N = 4, 8, 16 and h = 0.02, 0.01, 0.005 with α = 10 −2 in Table 5, and with α = 10 −5 in Table 6. For all these cases, we can see that the error decays when N gets larger (especially when h is small) and h gets smaller (especially when N is large). However, the decay slows down significantly when the parameter N becomes sufficiently large or h becomes sufficiently small. The reason comes from the cut-off approximation of the series of f and the regularization process. We also notice that the relative errors corresbonds to α = 10 −2 is larger than that to α = 10 −5 , which is also as expected.