An approach to building explicit time-marching stencil computation schemes for thetransient2Dacousticwaveequationwithoutusingﬁnite-diﬀerenceapproxima- tionsisproposedandimplemented. Itisbasedonusingtheintegralrepresentation formula (Poisson’s formula) that provides the exact solution of the initial-value problem for the transient 2D scalar wave equation at any time point through the initial conditions. For the purpose of constructing a two-step time-marching al- gorithm, a modiﬁed integral representation formula involving three time levels is also employed. It is shown that integrals in the two representation formulas are exactly calculated if the initial conditions and the sought solution at each time level as functions of spatial coordinates are approximated by stencil interpolation polynomialsintheneighborhoodofanypointina2DCartesiangrid. Asaresult, ifauniformtimegridischosen,theproposedtime-marchingalgorithmconsistsof twonumericalprocedures: 1)thesolutioncalculationattheﬁrsttime-stepthrough theinitialconditions;2)thesolutioncalculationatthesecondandnexttime-steps using a generated two-step numerical scheme. Three particular explicit stencil schemes (with ﬁve, nine and 13 space points) are built using the proposed approach. Their stability regions are presented. The obtainedstencilexpressionsarecomparedwiththecorrespondingﬁnite-diﬀerence schemes available in the literature. Their novelty features are discussed. Simula- tion results with new and conventional schemes are presented for two benchmark problemsthathaveexactsolutions. Itisdemonstratedthatusingthenewﬁrsttime- stepcalculationprocedureinsteadoftheconventionalonecanprovideasigniﬁcant improvement of accuracy even for later time steps. Keywords: Transient 2D acoustic wave equation, Initial-boundary value Email address: khutoryn@drexel.edu (Naum Khutoryansky) Preprint submitted to Elsevier May 1, 2019 problems, Poisson’s formula, Stencil computation, Explicit two-step numerical scheme, First time-step expression 1. Introduction Stencil computations are widely implemented in many numerical algorithms that involve structured grids, including ﬁnite-diﬀerence techniques in such di- verse areas as acoustics, ﬂuid dynamics, solid mechanics and electromagnetism. Classical methods to obtain ﬁnite-diﬀerence approximations of their governing partialdiﬀerentialequationsaremostlybasedonTaylorexpansions. However,for somestencilswithmultiplenodes,coeﬃcientsoftheseapproximationscannotbe uniquely determined without additional conditions (see, e.g., [8]) which creates competing ﬁnite-diﬀerence schemes for the same problem of interest and for the same set of stencil nodes. Inanexplicitﬁnite-diﬀerencescheme,asolutionvalueateachpointinatime- spacegridiscalculatedusingalinearcombinationofvaluesatitsspatialneighbors from previous time steps. Among the explicit schemes for the transient 2D wave equation the main attention in the literature has been given to two-step schemes (whichoperateoverthreetimelevels tk+1= (k+1)t,tk=ktandtk1= (k1)t where tis a ﬁxed time increment). They have been intensively studied and reviewed in many articles and books (see, e.g., [4, 8, 2, 3, 13, 10]). Duetotheabovementioneduncertaintiesthatbuildingﬁnite-diﬀerenceschemes forthetransient2Dwaveequationcanencounter,itisreasonabletoconsiderbuild- ingexplicitstencilcomputationschemesforthisequationbasedonprinciplesdif- ferentfromtheﬁnite-diﬀerenceapproximationofdiﬀerentialoperators. Thereare veryfewpapersinthisdirection. Anintegralevolutionformulathatcanbeuseful for building two-step schemes was derived in [1]. This formula was implemented in [11] for the 2D case using a least square piecewise polynomial approximation andquadraturealgorithmbutforthesakeofsimplicityignoringtheﬁrsttime-step issues ("startup issues"). The two-step schemes pose some challenges when imposing the initial con- ditions. To calculate the value of a sought solution uat the ﬁrst time level t1, one needs values of ufrom time levels t0andt1. The initial condition for u provides the required values at t0. The needed information at the earlier time t1 must be inferred. The initial condition for v=¶u/¶tcan be used to do this in diﬀerentways[13]. Theconventionalapproachdescribedinmanytextbooks(see, e.g., [10]) uses the central diﬀerence for approximating that initial condition. For 2 any point xinR2, the value of u(x,t1)is inferred as u(x,t1)2tv(x,0). Even thoughthisadd-hocapproachisattractiveduetoitssimplicity,itisworthwhileto consider the numerical schemes for both the ﬁrst time step and the other steps in the framework of a uniﬁed approach. This paper focuses on building explicit time-stepping stencil computation schemes for the transient 2D acoustic (scalar) wave equation without using ﬁnite- diﬀerence methods. The approach implemented here is based on using a famous integralrepresentationformula(Poisson’sformula)thatprovidestheexactsolution of the initial-value problem for the transient 2D scalar wave equation at any time point through the initial conditions. For the purpose of constructing a two-step time-marching numerical scheme, a modiﬁcation of Poisson’s formula involving three time levels is also employed. A general form of the proposed time-stepping algorithm is presented in section 3 where separate integral expressions are given for the ﬁrst time-step and for the further steps. A main part of the paper deals with deriving explicit stencil computation schemesfromtheseintegralexpressionsusingpolynomialinterpolationonstencils andexactintegrationwhichisdescribedinsections4and5. Particularnumerical schemes for 5-, 9- and 13-point stencils presented as linear combinations of node values are obtained in section 6 where exact coeﬃcients of the linear combina- tions are expressed through the Courant number. For each scheme, two separate expressions are derived: 1) a one-step expression for the ﬁrst time-step through the stencil interpolated initial conditions; 2) a two-step expression for the second and next time-steps. All the obtained expressions for the ﬁrst time-step have not been previously presented in the literature. A generated two-step expression for the 9-point stencil is also new. In contrast, the obtained two-step expressions for 5- and 13-point stencils coincide with the corresponding ﬁnite-diﬀerence sten- cil expressions. It is shown by simulation that the derived numerical schemes can signiﬁcantly improve accuracy of stencil computations in comparison with conventional approaches. 2. Representation formulas Consider a Cauchy problem for the transient 2D scalar wave equation ¶2u ¶tc2¶2u ¶x2 1+¶2u ¶x2 2 =0,u=u(x,t),x= (x1,x2)inR2,(1) ujt=0=u0(x),¶u ¶tjt=0=v0(x). (2) 3 Its solution is given by the representation formula which is often named Pois- son’s formula for the 2D wave equation [5]: u(x,t) =1 2pct2Z jyxj2<t2tu0(y) +tÑu0(y)(yx) +t2v0(y)p c2t2jyxj2dy1dy2(3) where Ñdenotes the gradient operator and y= (y1,y2)2R2is a variable of integration. Rewriting the previous integral on the unit disk, we get u(x,t) =1 2pZ jzj<1u0(x+ctz) +ctÑu0(x+ctz)zp 1jzj2dz1dz2 +t 2pZ jzj<1v0(x+ctz)p 1jzj2dz1dz2(4) where z= (z1,z2)2R2is a new variable of integration. Formulas(3)and(4)arealsovalidfornegative twhichcanbeprovedbasedon thetimereversalpropertyofthewaveequation. Bysubstituting tfortinformula (4) and changing ztozinside its integrals, we obtain the following analog of that formula for negative time: u(x,t) =1 2pZ jzj<1u0(x+ctz) +ctÑu0(x+ctz)zp 1jzj2dz1dz2 t 2pZ jzj<1v0(x+ctz)p 1jzj2dz1dz2(5) The only diﬀerence between the right-hand sides of (4) and (5) is the opposite signsofthesecondterm. So,onecaneliminatethistermbysumming(4)and(5). Shiftingtheinitialmomentintheresultingformulafrom t=0tot=t,weobtain the following expression involving three time points with a time increment t: u(x,t+t) +u(x,tt) =1 pZ jzj<1u(x+ctz,t) +ctÑu(x+ctz,t)zp 1jzj2dz1dz2 (6) where the right-hand side does not include the time derivative v. A similar representation formula involving three time levels has been derived (in a diﬀerent form) in [1]. Both formulas (4) and (6) will be used below to build a time-marching stencil computation algorithm for the wave equation. 4 3. The proposed time-stepping algorithm in the integral operator form Considerauniformtimegrid ft0=0,t1=t,:::,tk=kt,:::gwhere tisaﬁxed time-step. Denote by uk(x)the restriction of u(x,t)to a moment t=tk. Next, denoteby A(x,t)andB(x,t)thefollowingintegraloperatorsactingoncontinuous functions deﬁned in R2: A(x,t)f() =1 2pZ jzj<1f(x+ctz) +ctÑf(x+ctz)zp 1jzj2dz1dz2(7) B(x,t)f() =t 2pZ jzj<1f(x+ctz)p 1jzj2dz1dz2 (8) Thetime-steppingalgorithmproposedinthispaperconsistsoftwoprocedures based on the representation formulas (4) and (6): 1)The procedure for the ﬁrst time-step which, according to (4), calculates u1(x)as u1(x) =A(x,t)u0() +B(x,t)v0(),x2R2(9) 2)Theprocedureforthesecondandnexttime-steps which,using(6)for t=tk, calculates uk+1(x)as uk+1(x) =2A(x,t)uk()uk1(x),x2R2,k=1,2,... (10) Whileformula(10)involvesthreetimelevels(twotime-steps),formula(9)for the ﬁrst time-step incorporates only two time levels (one time-step) without any ﬁnite diﬀerence approximation of the time derivative. 4. Using polynomial interpolation on stencils Consideratwo-dimensionaluniformCartesiangrid f(i1h,i2h)gwhere i1andi2 areintegers,and histhegridspacinginbothdirections x1andx2. Supposethatthe evaluationpoint xinformulas(9)and(10)isagridpoint xi j= (ih,jh). Ourinten- tionistochooseastencilintheCartesiangridandreducetheintegralformulas(9) and(10)tolinearcombinationsofthestencilnodevaluesof uk(x)(k=0,1,2, :::) andv0(x). Such a reduction will be done by using polynomial interpolation. Assume that a particular stencil with mnodes is chosen for polynomial inter- polation. The corresponding index set f(q1,q2)gis denoted by Qm. The stencil indexcomponents q1andq2arenumberedrelativetothereferencingpointlocated 5 at the evaluation point xi j. So, polynomial interpolation in a neighborhood of the evaluation point xi jwill be carried out using interpolation points xi,j+xq1,q2=xi+q1,j+q2,(q1,q2)2Qm. (11) Following [12] we associate with Qma set of mdistinct bivariate monomials Mm=fxa(q1) 1xa(q2) 2,(q1,q2)2Qmg (12) where a(q) =( 2jqj1ifq<0 2qifq0.(13) The above function has a unique inverse function q(a) = (1)aa+1 2 (14) where []isthewholepartfunction. So,foreachindexvalue (q1,q2))2Qmthere exists a unique monomial from Mmand vice versa. Consider a polynomial space Pmspanned by Mm. We will use only those stencils for which the Lagrange interpolation problem is unisolvent in Pm(see [6]). Inthiscase,thereexistsaLagrangebasisfor Pmthatcanbebuiltasdescribed below. Suppose that there is an ordering imposed on the monomials in Mm fm1(x),:::,mm(x)g. (15) The corresponding stencil nodes are numbered accordingly using (12)-(14): fx(1),...,x(m)g. (16) Thus, one can compute the following matrix: D=h ms(x(r))i mm(17) If this matrix is non-singular, i.e., det(D)6=0, which means that the Lagrange basis exists, then the inverse matrix C=[csr]mm=D1(18) 6 can be calculated. Its components are instrumental in expressing the Lagrange basis functions through the chosen set of monomials: Ls(x) =m å r=1csrmr(x),s=1,..., m. (19) Once the Lagrange basis is obtained, we again need two indexes to denote the Lagrange basis functions in accordance with the two index notation for grid points. Thesequence(15)deﬁnesarelationship s=g(a1,a2)betweentheordinal number sand monomial exponents (a1,a2). Then, the resulting relationship betweenanindexpair (q1,q2)andthecorrespondingordinalnumber sisgivenby the following expression: s=g(q1,q2) =g(a(q1),a(q2)). (20) where a()is speciﬁed in (13). Therefore, by introducing a new (two index) notation for the Lagrange basis functions fq1,q2(x) =Lg(q1,q2)(x), (21) we get the interpolation formula for a continuous function f(x)in the form ˜f(x) =å (q1,q2)2Qmfi+q1,j+q2fq1q2(x) (22) Eventhoughdiﬀerentsetsofmonomialsandtheirsequencescanbeemployed forbuildingLagrangebases,wewilluseaparticularmethodofmonomialordering that is described below. Denote by Mthe set of all monomials xa1 1xa2 2where a1anda2are natu- ral numbers including 0. For each monomial xa1 1xa2 2, the corresponding ordinal number swill be assigned using the following function: s=g(a1,a2) =(a1+a2)(a1+a2+1) 2+( a1a2ifa2<a1, a2a1+1ifa2a1.(23) It is easy to prove that the function (23) provides a one-to-one correspondence between Mand the set of all positive natural numbers with the usual ordering. The order induced by (23) uses the total degree as the ﬁrst sorting parameter (similarly to the more common graded lexicographic order) while the diﬀerence 7 between individual degrees with the same total degree is used as the next sorting variable. The setMendowed with the order induced by (23) will be denoted by M. Theinitialsegmentof Mwith mmembersisdenotedbelowby M m. Itwillbe shown in section 6 that such ordered monomial sets play a useful role in building particular numerical schemes. Example1. ConsiderbuildingtheLagrangebasisforthespaceofcompletesecond degreepolynomials. Inthiscase,themonomialbasissequenceorderedaccording to (23) is as follows: M 6=f1,x1,x2,x1x2,x2 1,x2 2g (24) The corresponding interpolation stencil node sequence will be written according to (14) as f(0,0),(h,0),(0,h),(h,h),(h,0),(0,h)g The index set Q6for this stencil is shown in Figure 1. 21 0 1 221012 21 0 1 221012 Figure 1: Index set for the 6-point interpolation stencil 8 Calculating matrix Dwe get D=2 66666641 1 1 1 1 1 0h0h h 0 0 0hh0h 0 0 0 h20 0 0h20h2h20 0 0 h2h20h23 7777775 Thismatrixisnon-singularwith det(D) =4h8. Asaresult,wegettheLagrange basis as f0,0=L1=1+x1x2 h2x2 2 h2x2 1 h2,f1,0=L2=x1 2hx1x2 h2+x2 1 2h2, f0,1=L3=x2 2hx1x2 h2+x2 2 2h2,f1,1=L4=x1x2 h2, f1,0=L5=x1 2h+x2 1 2h2,f0,1=L6=x2 2h+x2 2 2h2.(25) OnecanseethattheLagrangebasis(25)includesmonomialsinthescaledvariables x1/handx2/hwith coeﬃcients independent of h. Itisworthtomentionherethatnode (h,h)willdisappearinthecorresponding numerical scheme (see section 6). 5. Calculating the integrals Letx=xi jin the integral operators (7) and (8). Without loss of generality, assume that the origin of the 2D coordinate system is located at xi jwhich can be achieved by a parallel translation of coordinates. Let f(y)appearing in (7) and (8) be a monomial in variables y1/handy2/h m(y) =y1 ha1y2 ha2,(a1,a2)2N2. (26) Then the integrals in (7) and (8) can be exactly calculated and expressed through the Courant number l=ct h(27) Indeed, using the table of integrals of Gradshteyn & Ryzhik [7], expressions (7) and (8) are reduced to the following exact values: A(xi j,t)m() =0,B(xi j,t)m() =0ifa1ora2are non-negative odd integers , (28) 9 A(xi j,t)m() =(a11)!!(a21)!! (a1+a21)!!la1+a2, B(xi j,t)m() =t a1+a2+1A(xi j,t)m() ifa1anda2are both non-negative even integers(29) where ()!!is the double factorial. It is assumed that (1)!!=1,0!!=1. The above formulas allow one to exactly calculate integrals (7) and (8) when f()is a polynomial from the Lagrange basis (see section 4). 6. Particular explicit two-step schemes Now we can start building some numerical schemes by transforming the pro- cedures (9) and (10) into algebraic expressions. All the functions u0(x),v0(x) anduk(x),k=1,2,...included in these procedures will be interpolated in a sten- cil’s center neighborhood using the same stencil nodes. The following standard notations for grid values of the solution and initial conditions will be used in the computation schemes: u0 i j=u(xi j,0),v0 i j=v(xi j,0), (30) uk i j=u(xi j,kt),k=1,2,:::. (31) 6.1. The 5-point stencil Consider the monomial basis (24) used in example 1 and the corresponding Lagrange basis (25). Using formulas from section 5 and notation (21) we get A(xi j,t)f00() =12l2,A(xi j,t)f1,1() =0, A(xi j,t)f1,0() =A(xi j,t)f0,1() =1 2l2, (32) B(xi j,t)f0,0() =t(12 3l2),B(xi j,t)f1,1() =0, B(xi j,t)f1,0() =B(xi j,t)f0,1() =t 6l2, (33) So, all coeﬃcients for node (h,h)disappear and the corresponding computa- tional scheme contains only 5 spatial points which is shown in Figure 2. 10 21 0 1 221012 21 0 1 221012 Figure 2: Index set for the 5-point numerical scheme 6.1.1. A new ﬁrst time-step expression As a result, the proposed 5-point numerical scheme is as follows: 1)for the ﬁrst time-step u1 i j=u0 i j+tv0 i j+l2 2 u0 i1,j+u0 i+1,j+u0 i,j1+u0 i,j+14u0 i j +tl2 6 v0 i1,j+v0 i+1,j+v0 i,j1+v0 i,j+14v0 i j ;(34) 2)for the second and next time-steps uk+1 i j=2uk i juk1 i j +l2 uk i1,j+uk i+1,j+uk i,j1+uk i,j+14uk i j ,k=1,2,...(35) Theconventionalapproachforthe5-pointstencil(describedinmanytextbooks) uses the same procedure (35) for the second and next time-steps. However, for the ﬁrst time-step, the conventional approach uses another formula (rather than (34)) based on the central diﬀerence for approximating the initial condition for v combined with (35) for k=0(see, e.g., [10]): u1 i j=u0 i j+tv0 i j+l2 2 u0 i1,j+u0 i+1,j+u0 i,j1+u0 i,j+14u0 i j .(36) 11 Comparing (34) and (36) one can see that the ﬁrst three terms of the right-hand partsoftheseequationscoincide,buttheforthtermpresentintheright-handside of (34) is absent in (36). Thus, the diﬀerence between these two ﬁrst time-step expressions depends on properties of v0(x). If, for example, v0(x)is a linear function of spatial coordinate, then there is no diﬀerence between (34) and (36). However,formoregeneralcasesencounteredinpractice,thediﬀerencemayexist and can be inﬂuential from the accuracy point of view which is shown in the two examples below (subsections 6.1.2 and 6.1.3). Thestabilityconditionforbothnumericalschemesutilizing(35)iswellknown from von Neumann stability analysis [13]: llmax=p 2 2. (37) 6.1.2. Standing wave simulation To make a numerical comparison of both approaches, we use the initial and boundary conditions corresponding to a standing wave exact solution: ue(x1,x2,t) =sin(2px1)sin(2px2)sin(2p 2pct). (38) This solution creates the following pair of initial conditions for the numerical simulation: u0(x1,x2) =0,v0(x1,x2) =2p 2pcsin(2px1)sin(2px2).(39) Weconsidertheunitsquare W= [0,1]2asaspaceregionforthenumericalsolution and apply boundary condition u=0on¶Wgenerated by (38). In addition, we assume that c=1. Let ntbe the number of time-steps and n2be the space discretization number related to hby equation nh=1. The above numerical schemes have been employed over a spatial index set f(i,j),0<i,j<ngusing grid boundary conditions uk 0,j=uk n,j=uk i,0=uk i,n=0, 0i,jn,k0. The accuracy of both numerical schemes have been estimated using the relative L2error deﬁned as E(n,nt) =0 B@ånt k=1ån i=0ån j=0 uk i jue(ih,jh,kt)2 ånt k=1ån i=0ån j=0(ue(ih,jh,kt))21 CA1/2 (40) 12 Therelative L2errorsfortheproposednew5-pointscheme(34) &(35)(basedon Poisson’s formula) and the conventional one (36) &(35) are denoted by EP5(n,nt) andEC5(n,nt), respectively. Calculated values of EP5(n,nt),EC5(n,nt)and their ratio for diﬀerent combinations of n,ntandl=0.707are presented in table 1 below. For every value of nused, three values of ntare considered: nt=1(to check errors after the ﬁrst step), nt=nandnt=2n(to demonstrate the error accumulation process). Table 1: Wave (38) simulations using the 5-point stencil. n n t l EP5 EC5 EC5/EP5 10 1 0.707 9.08431046.893810275.9 10 10 0.707 9.15401046.894510275.3 10 20 0.707 9.16041046.894510275.3 20 1 0.707 5.47671051.6636102304 20 20 0.707 5.68001051.6638102293 20 40 0.707 5.73721051.6638102290 40 1 0.707 3.39241064.12301031215 40 40 0.707 4.03311064.12341031022 40 80 0.707 4.49281064.1234103918 80 1 0.707 2.11581071.02851034861 80 80 0.707 4.38201071.02861032347 80 160 0.707 6.58241071.02861031563 Table1clearlyshowsthatthedemonstratedaccuracyofthenewschemeismuch higher than that of the conventional one with the error ratio EC5(n,n)/EP5(n,n) exceeding 103for more dense grids. So, using the new ﬁrst time-step calculation procedure instead of the conventional one is able to provide a signiﬁcant increase of accuracy even for later time-steps. Observed increases of the relative L2errors from nt=ntont=2nare negligible (less than 106). 13 6.1.3. Simulation of a trailing wave produced by a point impulse source Another benchmark example that we have simulated in the same space region W= [0,1]2is a trailing wave solution for t0corresponding to a point impulse source wave excited from the center of region W= [0,1]2at a negative moment t=d: ue(x1,x2,t) =1 2pcH(c(t+d)r)p c2(t+d)2r2,r=q (x10.5)2+ (x20.5)2,(41) where H()is the Heaviside function and d>0is large enough to allow the corresponding wavefront to move outside of the region Wbefore t=0. This conditionissatisﬁedifweassumethat cd>1/p 2. Forthenumericalsimulation, itisassumedthatthesolution(41)speciﬁesboundaryvaluesof u(x,t)on¶Wand the following pair of initial conditions (setting c=1for the simplicity sake): u(x1,x2,0) =1 2pp d2r2,v(x1,x2,0) =d 2p(d2r2)3/2inW.(42) The resulting relative L2errors (40) for d=0.75are presented in Table 2. They show that the accuracy of the new scheme is approximately 10 times higher than that of the standard one. Table 2: Wave (41) simulations using the 5-point stencil. n n t d l EP5 EC5 EC5/EP5 10 10 0.75 0.707 1.14781031.11241029.69 20 20 0.75 0.707 3.82211043.822810310.0 40 40 0.75 0.707 1.07991041.128110310.4 80 80 0,75 0.707 2.83791053.000010410.6 6.2. A 9-point square stencil Considertheinitialsegment M 11ofMwiththelastmemberequalto x2 1x2 2: f1,x1,x2,x1x2,x2 1,x2 2,x2 1x2,x1x2 2,x3 1,x3 2,x2 1x2 2g (43) Thissetincludesallthemonomialsoftotaldegree 3andonemonomialoftotal degree =4, and is the minimal initial segment of Mthat includes the bivariate 14 tensor-product of the monomial bases for the second degree polynomials in each spatial coordinate. The corresponding index set Q11is shown in Figure 3 where both solid dots and empty circles denote interpolation points. The corresponding Lagrangebasiscanbeeasilycalculatedusingmatrix(17)andpresentedsimilarto (25). However, it is preferable to avoid presenting long expressions that include 11 monomials and display only terms that will be used later. Denote by fe q1,q2a part of fq1,q2that includes all the monomial terms with even exponents in both coordinates. Then we get fe 0,0=1x2 1 h2x2 2 h2+x2 1x2 2 h4,fe 2,0=fe 0,2=0, fe 1,0=1 2x2 1 h2x2 1x2 2 h4 ,fe 0,1=1 2x2 2 h2x2 1x2 2 h4 ,fe 1,1=1 4x2 1x2 2 h4.(44) 6.2.1. A new 9-point time-stepping scheme Itfollowsfromtheresults(28)and(29)ofsection5thatfunctions fe q1,q2rather than the complete Lagrange basis will be used in building numerical schemes. Therefore, two points (-2,0) and (0.-2) (presented by empty circles in Figure 3) willdisappearinthecorrespondingnumericalscheme. Theremainingnodes(solid circlesinFigure3)createthe9-pointsquare-shapedcomputationalstencil. Using 321 0 1 2321012 321 0 1 2321012 Figure 3: Index set for the 11-point interpolation stencil 15 (28)-(29) we obtain: A(xi j,t)f00() =12l2+l4 3,A(xi j,t)f1,1() =l4 12, A(xi j,t)f1,0() =A(xi j,t)f0,1() =l2 2l4 6,(45) B(xi j,t)f00() =t 12l2 3+l4 15 ,B(xi j,t)f1,1() =tl4 60, B(xi j,t)f1,0() =B(xi j,t)f0,1() =tl2 6 1l2 5 .(46) Now we can use (9) and (10) to build a new 9-point numerical scheme similar to (34)-(35). However, to avoid writing long expressions, some additional notations will be needed: dk i j(q1,q2) =uk i+q1,j+q2+uk iq2,j+q1+uk i1q1,jq2+uk i+q2,jq14uk i,j, e0 i j(q1,q2) =v0 i+q1,j+q2+v0 iq2,j+q1+v0 iq1,jq2+v0 i+q2,jq14v0 i,j, k=0,1,2, :::.(47) Thus, the following time-stepping numerical scheme is derived: 1)for the ﬁrst time-step u1 i j=u0 i j+tv0 i j+l2 2" 1l2 3 d0 i j(1,0) +l2 6d0 i j(1,1)# +tl2 6" 1l2 5 e0 i j(1,0) +l2 10e0 i j(1,1)# ;(48) 2)for the second and next time-steps uk+1 i j=2uk i juk1 i j+l2" 1l2 3 dk i j(1,0) +l2 6dk i j(1,1)# ,k=1,2,:::. (49) AvonNeumannstabilityanalysis(see,e.g.,[13])providesthestabilitycondi- tion for the above numerical scheme as llmax=p 3p 3 20.796. (50) 16 Toourknowledgethisnumericalschemehasnotbeenpresentedpreviouslyin the literature. A conventional explicit nine-point square-shaped scheme (dubbed theisotropicscheme(see,e.g.,[14])hasadiﬀerentformobtainedusinga9-point ﬁnite diﬀerence approximation of the Laplace operator in the two-dimensional space [9]: uk+1 i j=2uk i juk1 i j+l2" 2 3dk i j(1,0) +1 6dk i j(1,1)# ,k=1,2,:::.(51) Comparingthecoeﬃcientsin(49)and(51)itiseasytoseethattheterm dk i j(1,1) islessinﬂuentialinthenew9-pointschemethanintheconventional9-pointscheme (with a coeﬃcient ratio equal to l2). The stability condition for the conventional scheme is less restrictive than for the new one: llmax=p 3/20.866. On the other hand, the new scheme has some accuracy advantages in comparison with the conventional 9-point approach for l0.796which is shown in the next subsection. 6.2.2. Numerical simulation Weconsidersimulationresultsfortheproposed9-pointscheme(48)-(49)using initial and boundary conditions generated by the exact solutions of the previous subsection. Let c=1. We assume that the unit square W= [0,1]2is used as a space grid region, the time grid interval is [0,l]with nt=n. A comparison is made with simulated results on the same grids for the conventional scheme (51). Even though no expression for the ﬁrst time-step corresponding to the conventional9-pointschemeispresentedin[14],theusualapproachbasedonthe centraldiﬀerenceforapproximatingtheinitialconditionfor vcombinedwith(51) provides the corresponding expression u1 i j=u0 i j+tv0 i j+l2 2" 2 3d0 i j(1,0) +1 6d0 i j(1,1)# (52) which is used for numerical simulation. In addition, simulation results are provided for the new scheme (49) but us- ing instead of (48) the conventional approach based on the central diﬀerence for approximating the initial condition for v. The corresponding ﬁrst time-step expression has the form u1 i j=u0 i j+tv0 i j+l2 2" 1l2 3 d0 i j(1,0) +l2 6d0 i j(1,1)# (53) 17 which coincides with (48) if the fourth term in its right-hand side is omitted. The simulation results are presented in the next two tables and cover both standing wave and trailing wave simulations. We have used two values for the simulation: l=0.707, as in subsection 6.1, and l=0.796, according to the stability condition (50). Relative L2errors for the schemes (48) &(49), (53) &(49) and (52) &(51) are denoted by EP9,ECP9andEC9, respectively. According to the simulated data, the new scheme (48) &(49) appears to be more accurate than the other two schemes for both l=0.707andl=0.796. A demonstrated higher accuracy of the scheme (48) &(49) in comparison with (53) &(49) emphasizes the importance of using the new ﬁrst time-step procedure instead of the conventional one. Table 3: Wave (38) simulations using 9-point stencils n n t l EP9 ECP9 EC9 10 10 0.707 3.70581028.82351021.1741101 10 10 0.796 2.95871029.60961021.1241101 20 20 0.707 8.93331032.13041022.8002102 20 20 0.796 8.06971032.34101022.7523102 40 40 0.707 2.37231035.26831036.8821103 40 40 0.796 2.57371035.82961036.8668103 80 80 0.707 7.55731041.31191031.7084103 80 80 0.796 1.02741031.45781031.7187103 ComparisonsoftheobtainedresultswithTables1and2showthataccuracyof theboth9-pointschemesinthesimulatedcasesfor l=0.707ismuchlowerthan that of the new 5-point scheme. 6.3. A 13-point stencil The next numerical stencil is based on complete bivariate interpolation poly- nomials of the fourth degree. Consider the initial segment M 15ofM: f1,x1,x2,x1x2,x2 1,x2 2,x2 1x2,x1x2 2,x3 1,x3 2,x2 1x2 2,x3 1x2,x1x3 2,x4 1,x4 2g(54) 18 Table 4: Wave (41) simulations using 9-point stencils n n t d l EP9 ECP9 EC9 10 10 0.75 0.707 2.60211031.31561021.5301102 10 10 0.75 0.796 2.89061031.71171021.8658102 20 20 0.75 0.707 8.98891044.35691034.9178103 20 20 0.75 0.796 9.48731045.53391035.9647103 40 40 0.75 0.707 2.92001041.25571031.3892103 40 40 0.75 0.796 3.48211041.56731031.6650103 80 80 0.75 0.707 1.04491043.30191043.6165104 80 80 0.75 0.796 1.49011044.07101044.2991104 Thissetincludesallthemonomialsoftotaldegree 4. Thecorrespondingindex setQ15is presented in Figure 4 where both thirteen solid dots and two empty circles denote interpolation points. After calculating the corresponding Lagrange basis (see section 4) and using (28)-(29),onecandeterminethatallthecoeﬃcientsinexpressions(9)-(10)related to the two empty circles will disappear. As a result, we get a numerical scheme that involves the 13-point stencil (the solid dots in Figure 4). The notations (47) will be used in presenting the scheme to make the expressions more compact: 1)for the ﬁrst time-step u1 i j=u0 i j+tv0 i j+l2 2" 42l2 3d0 i j(1,0) +l2 6d0 i j(1,1) +l21 12d0 i j(2,0)# +tl2 6"4 32l2 5 e0 i j(1,0) +l2 10e0 i j(1,1) +l2 201 12 e0 i j(2,0)# ;(55) 2)for the second and next time-steps uk+1 i j=2uk i juk1 i j+l2" 42l2 3d0 i j(1,0) +l2 6d0 i j(1,1) +l21 12d0 i j(2,0)# , k=1,2,:::. (56) 19 321 0 1 2 3202 321 0 1 2 33210123 Figure 4: Index set for the 15-point interpolation stencil Formula(56)forthesecondandnexttime-stepscompletelycoincideswiththat obtained previously [2, 3] using the ﬁnite-diﬀerence method. However, expression (55) for the ﬁrst time-step has not been presented in the literaturebefore. Anadvantageofusingthisexpressionratherthantheconventional one (based on the central diﬀerence for approximating the initial condition for v) is shown in Table 5 where simulation results for the standing wave (38) solution are presented. The relative L2errors for the conventional and new approaches are denotedby EC13andEP13,respectively. Forsimplicitysake,theperiodicboundary conditions have been incorporated in the simulation taking into account that the solution (38) is periodic in both spatial directions. Since the maximal Courant number needed for stability of this scheme is 1/p 2, a value of l=0.707was used. The large values of the error ratio EC13/EP13demonstrated in Table 5 are explained by a higher accuracy of the new ﬁrst time-step expression (55) in com- parison with the conventional one. On the other hand, it is worth to notice that the relative L2error EC13(corresponding to the conventional ﬁrst time-step ap- proach) in Table 5 has almost the same values as the error EC5in Table 1 despite usingahigherdegreeinterpolationstencilinTable5. Thatis,anerrorintroduced at the ﬁrst time-step probably suppresses advantages of using a higher degree interpolation at later time-steps. 20 Table 5: Wave (38) simulations using the 13-point stencil. n n t l EP13 EC13 EC13/EP13 10 10 0.707 4.21461056.89381021.6357103 20 20 0.707 6.60041071.66361022.5205104 40 40 0.707 1.14711084.12301034.0537105 80 80 0.707 2.888410101.02851033.5608106 7. Conclusion Anewmethodisproposedtobuildexplicittime-steppingstencilcomputation schemesforthetransient2Dacousticwaveequation. ThemainideaistousePois- son’sformulathatprovidestheexactsolutionoftheinitial-valueproblemthrough the initial conditions. Implementation of this formula together with polynomial stencilinterpolationoftheinitialconditions(andexactintegration)foraparticular 2D stencil creates a uniﬁed time-marching scheme that combines computation expressions for the ﬁrst time-step and for the next steps. Particularexplicitstencilcomputationschemes(with5,9and13spacepoints) are derived in the paper. Their stability conditions are determined. All of the obtained ﬁrst time-step computation expressions are diﬀerent from those used in conventional ﬁnite-diﬀerence methods. Simulation comparison results are pre- sented for two benchmark problems that have exact solutions. It is demonstrated by simulation that the proposed computation approach maintains an accuracy ad- vantageincomparisonwithconventionalﬁnitediﬀerenceschemeswhichismostly attributed to the new ﬁrst time-step computation expressions. Acknowledgments This research did not receive any speciﬁc grant from funding agencies in the public, commercial, or not-for-proﬁt sectors.