A SECOND-ORDER ABSORBING BOUNDARY CONDITION FOR TWO-DIMENSIONAL PERIDYNAMICS

. The aim of this paper is to develop numerical analysis for the two-dimensional peridynamics which depicts nonlocal phenomena with artificial boundary conditions (ABCs). To this end, the artificial boundary conditions for the fully discretized peridynamics are proposed. Then, the numerical analysis of the fully discretized scheme is developed such that the ABCs solve the corner reflection problem with second-order accuracy. Finally numerical examples are given to verify theoretical results.

The fractional Laplacian (− 2  ) /2 is a special case of the integro-differential operator ℒ  .We must point that the displacement  for the real peridynamics is a vector.But for convenience we deal with scalar peridynamics here.The study for the scalar cases can be extended to the vector cases just as lattice dynamics.Therefore, we will generally refer to scalar peridynamics that we treat here as peridynamics from now on.
The artificial boundary conditions (ABCs) are always applied to avoid difficulty caused by unboundedness of the spatial domain while solving problems defined on an unbounded domain.To this end, one needs to reformulate the unbounded problem into a problem defined on a bounded computational domain by applying appropriate ABCs to the domain's boundary.Inappropriate boundary conditions will lead to numerical reflection that disrupt the numerical solution in the computational domain while the waves touch the artificial boundaries.Weckner's work has indicated the difficulties in designing boundary condition [39] while studying peridynamics in an unbounded medium.An ideal boundary condition can efficiently avoid numerical reflection by annihilating the waves that touch the artificial boundaries.In this paper the ABCs are applied to the boundary of a rectangular computational domain to reformulate the unbounded problem (1) into an initial-boundary-value problem.
These years, many contributions of ABCs (see e.g.[2,3,6,14,15,[17][18][19]35]) have been devoted to the local PDEs, such as wave equation and Schrödinger equation.Comparing to local PDEs, much less works on ABCs of nonlocal PDEs are available in the literature.For one-dimensional cases, we refer to [42,43] for the boundary treatments of the nonlocal heat equations and to [30,41] for the boundary treatments of nonlocal Schrödinger equations.As for one-dimensional nonlocal wave equations, we refer to [12].In [5,21] the PML approach is applied to study the one-dimensional heat and Schrödinger equations involving fractional operators.Furthermore, a general PML approach is introduced to study fractional PDEs in [1,4].Some local in-time ABCs are also introduced to treat the one-dimensional peridynamics in [34,37].However, for two and three dimensional cases, these algorithms have difficulties in dealing with the corner problem.In [11], ABCs for two-dimensional nonlocal PDEs are proposed, but without numerical analysis.In general, there are much less works devoted to ABCs for two-dimensional nonlocal models.A crucial difficulty for designing ABCs of two-dimensional nonlocal models is to compute the Green's functions accurately [20,29,31,36].In the present work, a numerical analysis is developed for the ABC proposed in [32].It can be proved that the numerical ABC is of second order by computing the Green's functions accurately [32].Thus, the ABC solves the problem of corner reflection efficiently (inappropriate boundary condition can be considered physically as an artificial obstacle setting on the boundary, which will make the waves that should have gone out of the computational domain to be reflected back at the boundary.This kind of reflection is more obvious at the corner).
The aim of the present paper is (i) to develop a numerical analysis for the numerical scheme of (1) with ABCs derived in [32] and (ii) to verify the numerical accuracy of the approach.In Section 2, a full discretization of ( 1) is introduced by a Crank-Nicolson time discretization and a spatial asymptotically compatible scheme.Then, the exact ABC proposed in [32] for the fully discretized version is stated.In addition, some properties of the ABCs are stated for numerical analysis.In Section 3, the numerical analysis of the scheme with ABCs is developed.In Section 4, we consider the wave equation, namely, ⎧ ⎪ ⎨ ⎪ ⎩   (, , ) = △(, , ) +  (, , ), (, ) ∈ R 2 ,  ∈ (0,  ], (, , 0) =  0 (, ),   (, , 0) =  1 (, ), (, ) ∈ R 2 , lim  2 + 2 →+∞ (, , ) = 0,  ∈ (0,  ]. ( In this section, we also develop the numerical analysis for the semi-discretized version of (3) with the exact ABCs by the same manner.In Section 5 we report some numerical examples.Finally, we draw a conclusion in Section 6.

Absorbing boundary conditions for the fully discretized peridynamics
In this section we give the result of the accurate absorbing boundary conditions and the full scheme for the two-dimensional peridynamics (1).
In the same manner, we use the symbol () to denote the location of the -th exterior layer node (black node in Fig. 1).We also define: In the same manner, the exterior layer vector is defined by ( We have left the expressions of () and () in the Appendix A.
The -transform of a complex-valued sequence (  ) ∈N is defined by In the same manner, the -transforms for a sequence of  ×  matrices {A  } ∈N and for a sequence of vectors {v  } ∈N ∈ C  , are defined by: respectively.In addition, the convolution ⋆ of two sequences {  }  and {  }  is defined by: while the convolution of two matrices sequences {A  }  and {B  }  is defined by: Finally, we introduce the function such that By the compact support assumption of initial data  0 and  1 , the computational domain ℛ  (0, 0) can be chosen such that the initial data is zero outside the selected computational domain, including its interior layer (i.e.compactly supported in the domain consisting of the white points in Fig. 1).In this setting, the unknowns in the exterior layer (black nodes in Fig. 1) are required to complete the discretized system (7), the complete relation is discretized boundary condition.To this end, the discretized ABC (9) in Theorem 1 is proposed [32].The discretized ABC (9) depicts the relation between the unknowns   () in the exterior layer (black nodes in Fig. 1) and the values   () (1 ≤  ≤ |ℬ Int |) in the interior layer (yellow nodes in Fig. 1).
We remarked the main computing cost comes from (9).For total computational step  , the computing cost for every convolution in ( 9) is ( 2 ).Thus, for fixed , the computing cost of ( 9) is Then, we know that the total computing cost of the boundary condition is

Errors analysis for the two dimensional peridynamics
According to the assumption in Introduction, we recall that   is homogeneous in the exterior domain.It may be a priori inhomogeneous in the computational domain.However, from now on we assume that   is also homogeneous inside the computational domain to simplify the presentation of the proofs.The   has a singularity behaviour   () ∼ || −2−2 with  < 1 around || = 0.In addition, if   is non homogeneous in the domain of computation, the proofs can be adapted but at the price of more complexity with the similar singularity types.

Estimates for the solution 𝑢 of system (1)
We first propose some useful estimates for the solution  to the continuous problem (1).Firstly we introduce Let us assume that ,  and ℓ are all nonnegative integers.The initial data  0 and  1 belong to  ++ℓ+2 (R 2 ) and the source term  belongs to   ([0,  ],  ++4 (R 2 )).Then for the solution  to the continuous problem (1), the following bounds hold, where  > 0 is a constant depends on  , , , ℓ,   and initial data  0 ,  1 .
Proof.It is clearly that The term in the last line of ( 20) This term tends to 0 as  → ∞.Thus, by taking  → ∞ in (20) one has and By the above two equalities one has Thus, taking  by   in (21) we have

Now let us consider the governing equation
By Summing up (22) multiplied by   , and integrating the equality on R 2 , by (24) and dominated convergence theorem we obtain By integral with respect to variable , the above inequality will give us Then, by Gronwell's inequality one has This gives By the same manner, taking  ℓ−1        on the both sides of ( 22) one has Multiplying ( 26) by  ℓ         and integrating it on R 2 , one has Further more, taking  ℓ−2        on the both sides of ( 22) for ℓ ≥ 2, we derive By induction one has Taking the above estimate into (27), one derives which gives (18).
In addition, we have the following equality, where ̂︀ ( 1 ,  2 , ) is the Fourier transform of (, , ) with respect to the variables  and  such that This gives us )︁ , which proves (19).

Truncated error estimates
Let us introduce   ,ℓ such that for (, ℓ) ∈ Z 2 with  ,ℓ () = (ℎ, ℓℎ, ) where (, , ) designates the solution to (1).The operator  ,ℎ is introduced in (7 where  depends only on  ,   and initial data  0 ,  +  4  6 Applying (33) with  =   =  , we have with In the same way by using (34) one has with Recalling [13] |ℒ  (, , ) we have For any function (, , ) smooth enough, there exists (x, ỹ, t) On the right hand side of (38) Thus, summing up the square of the first two terms in the right hand side of ( 38) with (, ℓ) ∈ Z 2 we obtain In the same manner, the square summation of the third term in the right hand side of ( 38) with (, ℓ) ∈ Z 2 satisfies the following estimate, and the square summation of the fourth term in the right hand side of ( 38) with (, ℓ) ∈ Z 2 satisfies the following estimate, Combining the estimates ( 39)-( 41) one derives By ( 37) and ( 42), finally we get Due to Lemma 3, above estimate implies (29).By the same procedure, one has Therefore, we obtain which leads to ( 30)-( 32) by Lemma 3.
The above equality is equivalent to As a consequence, we have while the term (54) can be written as By the above equality, equations ( 53) and (54), we deduce It is clear that (49) can be written as . Thus, we can rewrite (48) as Multiplying (55) by ℎ 2 and summing up it from  = 1 to  − 1, one has The first term for R.H.S of (56) can be evaluated as Firstly we estimate the terms on the right hand side of (57) with respect to the numerical solution   ,ℓ , namely, Consequently, by ( 16) of Theorem 1, as  → ∞, above estimate satisfies On the other hand, as  → ∞, the terms on the right hand side of (57) with respect to the exact solution  ,ℓ () satisfy Recalling (17) with  0 ,ℓ =  0 (ℎ, ℓℎ) and  0 ,ℓ =  1 (ℎ, ℓℎ) we derive that (ℛ  (ℎ,ℓℎ)×[0, ]) , Therefore we have By (62), above equality gives By Lemma 4, equation ( 61) and above estimate one gets which is equivalent to by which a Gronwall inequality [8] yields: Multiplying the above equality by ℎ 2 ( +  +  −  ) +    ,ℓ , and summing up all the terms for (, ℓ) ∈ ℛ  (0, 0), we get )︀ /(2 ).
After the same manipulation as previous, as  → ∞ one has )︀ , which leads to Thus, we derive that Finally, for any (, ℓ) ∈ Z 2 we have which gives (65).
We remark that the main difficulty in numerical analysis is the complicated estimate on the boundary, but the exact discretized boundary conditions help us to avoid the difficulty.
We also point that frog-leap scheme will be more suitable for wave equation due to the good CFL condition.It will be interesting to analysis the stability and error by using frog-leap scheme.

Accurate boundary conditions for wave equation on rectangular domain
The initial data  0 (, ),  1 (, ) and the source  (, , ) of peridynamics (3) are assumed to be compactly supported in [−  ,   ] × [−  ,   ].Under the settings of Section 2.2, based on CN time scheme at time   =  we obtain the fully discretized version for (3), namely, with which is exactly a special case of semi-discretized peridynamics (7).Thus, the boundary conditions can be obtained by Theorem 1 for the special  ,ℓ proposed in (64).
In the same manner, by refining  and fixed ℎ = 2 −7 , for  = 5 (with   =  ) the numerical errors are ploted in the Figure 4 with  = 0.125 and  = 0.25.We observe that the scheme is around two-order in time.
In addition, we show the behaviour of the numerical solution with initial data  0 and  1 = 0 at different times for  = 0.25 in Figure 5.We observe that the wave packet goes through the corners without reflection.

Numerical example of wave equation
The wave equation ( 3) is considered with the initial data  0 and  1 as the followings,   We employ the 4-th spatial order scheme in Remark 1 and the corresponding ABCs to implement the computation.The computational domain is set as [−2.5, 2.5] 2 and the total computational time is set as  = 4.We still compute the reference solutions  ref on a larger spatial domain [−16, 16] 2 with ℎ = 2 −8 and  = 2 −16 .The  ∞ -norm error  ℎ, ∞, is defined by the same way as it was in Section 5.1 by refining ℎ, fixing  = 2 −16 and  = 4 (with   =  ).The numerical errors are ploted in Figure 6.We observe that the scheme is around four-order in space.
By refining  and fixing ℎ = 2 −8 , the numerical errors are ploted in Figure 7 with  = 4 (  =  ).We observe that the scheme is around two-order in time.
In the Figure 8 we shows the behaviour of the numerical solution with initial data  0 (and  1 = 0) at different times.We observe that the wave packet goes through the corners without reflection.