A Schwarz-based domain decomposition method for the dispersion equation

We propose a Schwarz-based domain decomposition method for solving a dispersion equation consisting on the linearized KdV equation without the ad-vective term, using simple interface operators based on the exact transparent boundary conditions for this equation. An optimization process is performed for obtaining the approximation that provides the method with the fastest convergence to the solution of the monodomain problem


Introduction
The Korteweg -de Vries (KdV) equation, derived by [11] in 1895, models the propagation of waves with small amplitude and large wavelength, taking in account nonlinear and dispersive effects.In terms of dimensionless but unscaled variables, it can be written as [2] u t + u x + uu x + u xxx = 0 As done in [14] (and in [3] as a special case of their work), we will focus in this paper on the linearized KdV equation without the advective term : to which we will refer as dispersion equation.
The work developed here is inspired from [14] and [3].Nevertheless, our objectives are different from theirs.In this paper we propose an additive Schwarz method (ASM) for solving the dispersion equation (1) in a bounded domain, i.e., we decompose the computational domain in subdomains and solve the timedependent problem in each one of them.Our work focuses on the formulation of appropriate and optimized conditions on the interface between the subdomains, in order to minimize the error due to the domain decomposition method (DDM) and to accelerate the convergence of the method.
The interface boundary conditions (IBCs) proposed here are based on the exact transparent boundary conditions (TBCs) for the equation (1), derived by [14] and [3].The TBCs make the approximate solution in the computational domain coincide with the solution of the whole domain, but its exact computation is not doable in general [1].[14] and [3] proposed numerical approximations for these conditions, seeking to reduce the error created by the introduction of artificial boundaries.
In the work presented here, we do not propose approximate TBCs for reducing the error related to the finitude of the computational domain.In fact, we intend to reduce the error created by the decomposition of the domain and the introduction of an artificial interface boundary condition, in the context of a DDM.In other words, we study the effectiveness of the boundary conditions as IBCs, not as TBCs.As a consequence, our work shall not use the same reference solution as the one used by [14] and [3]: for validating their approaches, they compare their approximate solution with the exact solution in the whole domain.On the other hand, our reference solution is the approximate solution computed on the computational monodomain.Moreover, in order to isolate the error due to the DDM from that originated by time discretization, we study our method locally in time, i.e., along one time step.This paper is organized as follows : In Section 2, we describe the DDM used here and we recall the exact TBCs derived by [14] for equation (1).Then, we propose approximations for them, leading to very simple mixed-type conditions (avoiding, for example, integrations in time) to be used as IBCs in the DDM.Small modifications are proposed for these IBCs such that the solution of the DDM problem converges exactly to the reference solution (the solution of the monodomain problem).In Section 3, we perform a large set of numerical tests in order to optimize the IBCs, in the sense that we search the coefficients that provide the fastest convergence for the DDM iterative process.

Resolution of the dispersion equation using a domain decomposition method
We propose in this section a DDM for solving the problem We firstly present a brief review of the DDM considered here, the parallel or additive Schwarz method (ASM) and then we propose IBCs for applying it to the problem solved in this paper.

The Schwarz Method
Domain decomposition methods allow to decompose a domain Ω in multiple subdomains Ω i (that can possibly overlap) and solve the problem in each one of them.Therefore, one must find functions that satisfy the PDE in each subdomain and that match on the interfaces.
The first DDM developed was the alternating or multiplicative Schwarz method, in which the IBCs for computing the solution in a subdomain are function of the most updated solution in the neighbor subdomains.We consider here a modified version of this algorithm, introduced more recently by [12] and known as parallel or additive Schwarz method.In this algorithm, the IBCs for computing the solution u k i , in the subdomain Ω i and iteration k, are always constructed using the solution u k−1 j , j = i, of the previous iteration in the neighbor subdomains.
This modification originates an inherently parallel algorithm, which one naturally implements with parallel computing.The advantages obtained with the parallelism become more evident when the number of subdomains increases [12].
In the ASM, the boundary condition for the problem in Ω i , in each interface between the subdomains Ω i and Ω j , can be written as where (3), B i denotes the operator of the IBC.This operator allows the construction of more general Schwarz methods: in the original one, the IBCs are Dirichlet conditions (i.e., B i (u) = u ) [10,13].Without loss of generality, in the following we consider a domain Ω decomposed in two non-overlapping subdomains, Ω 1 and Ω 2 , with Γ = Ω 1 Ω 2 .
When implementing a Schwarz method, one must define appropriate operators B i such that: • There is a unique solution u i in each subdomain Ω i ; • The solution u i in each subdomain Ω i converges to u| Ωi , i.e., the solution u, restricted to Ω i , of the problem in the monodomain Ω; Moreover, one wants the method to show a fast convergence.
In fact, accordingly to [10], the optimal additive Schwarz method for solving the problem where A is a partial differential operator, is the one which uses as IBCs the exact TBCs for the problem, which are given by where ∂n i is the outward normal to Ω i on Γ , and the D2N (Dirichlet to Neumann) operator is defined by with α defined on Γ. v is solution of the following problem, solved in the complementary set of Ω i , denoted by The ASM using such exact TBCs is optimal in the sense that it converges in two iterations, and no other ASM can converge faster [10].Nevertheless, these TBCs, in general, are not simple to compute both analytically and numerically.More specifically, they are nonlocal in time, so they must be approximated for an efficient numerical implementation [1].These facts motivate us to look for simpler operators to use as IBCs in the ASM.We propose them based on the exact TBCs operators for the equation (1), as derived by [3].

Interface boundary condition operators based on the exact TBCs for the dispersion equation
In [3], TBCs are derived for the one-dimensional continuous linearized KdV equation (or Airy equation): where U 1 ∈ R, U 2 ∈ R + * and h is a source term, assumed to be compactly supported in a finite computational domain [a, b], a < b.
For the homogeneous initial boundary value problem the TBCs are given by [3, equations (2.17) -(2.18)] where L −1 denotes the inverse Laplace transform, * the convolution operator, s ∈ C, Re(s) > 0, is the Laplace frequency and λ 1 is, among the three roots of the cubic characteristic equation obtained when solving (4) in the Laplace space and in the complementary set of [a, b], the only one with negative real part.
In this paper, we focus on the special case U 1 = 0, U 2 = 1, which results on the dispersion equation (1).In this case, accordingly to [14], the only root with negative real part is The computation of the TBCs ( 5) is not simple due to the inverse Laplace transform, which makes these conditions nonlocal in time.Therefore, we propose approximations of the root (6) that avoid integrations in time, making the operators considerably simpler.
Obviously, we do not expect these operators to be as accurate as the approximate TBCs proposed by [3] (who derives TBCs for the discrete linearized KdV equation).Nevertheless, the objectives of our work and the work of [3] are very different: while they seek to minimize the error of the computed solution (compared to the analytical one) due to the boundary conditions, we want here to apply our operators as IBCs in a DDM.Therefore, our objective lays on the convergence of the DDM to the solution of the same problem in the monodomain, independently of the errors on the external boundaries.
We use the constant polynomial P 0 (s) = c for approximating λ 2 /s (which can be seen as a (0,0) order Padé approximation).Moreover, as a consequence of (6), we can approximate the other operands of the inverse Laplace transforms in (5) only in function of c : Replacing ( 7) in ( 5), using some well-know properties of the Laplace Transform (linearity and convolution) , we get the approximate transparent boundary conditions We notice that the approximation (8) has the same form as the exact TBCs for the equation (1) presented in [14] and [3], being the constant c an approximation for fractional integral operators.
We also remark that ( 8) are mixed-type boundary conditions (up to the second derivative of the solution), which, in the next section, we apply as IBCs in a DDM and we seek to optimize in order to accelerate the convergence of this method.The idea of using optimized boundary conditions in DDMs was already explored in [9] and [4], in the context of the Schrödinger equation.
Considering a discrete domain with mesh size ∆x and points x 0 , ..., x N and using some finite difference approximations, the operators (8) are discretized as

ASM with the proposed IBCs
With the operators Θ c i defined, we now apply them in a DDM, with two nonoverlapping subdomains Ω 1 = [a, 0] and Considering that we want to analyze and minimize the error due to the application of a DDM (isolating it from the error acumulated along the time steps, due to the temporal discretization), the reference solution u ref in our study is the solution of the monodomain problem (2) solved along onte time step (equation 10).Therefore, we implement a DDM to an evolution problem discretized in time (thus consisting in an ODE in space), an idea already explored by [5].
The external BCs Υ i , i = 1, 2, 3 (i.e., defined on ∂Ω i \Γ) are independent of the interface BCs.Here, we consider This choice was based on the simple form and implementation of these boundary conditions.Nevertheless, it does not have much importance in the study done here, as we want to analyze exclusively the behavior of the DDM.The only restriction for an appropriate study is that the external BCs for computing u ref must be the same Υ i , i = 1, 2, 3, used for each subdomain in the DDM, as we did in ( 11)-( 12) and (10).
The resolution of the problem (10) by the Additive Schwarz method and using the IBCs ( 9) is written as A simple analysis (for example in the Laplace domain) shows that the monodomain and DDM problems (10) and ( 11)-( 12) have an unique solution.
We also remark that our proposed DDM can be used for solving the problem (2), i.e., in a time window containing multiple time steps, by solving (11)- (12) in each time step, with the converged solution of the previous time step as initial data.
Remarks on the notation.In the following study of our proposed DDM, where we perform a spatial discretization, we introduce an extra subindex, so the solution is denoted as u k i,j , where i indicates the subdomain Ω i (or, in the case of the reference solution, i = ref , and in the convergence of the method, i = * ), j indicates the spatial discrete position and k indicates the iteration.

Spatial discretization of the problem
Concerning the spatial discretization, the monodomain Ω is divided in 2N +1 homogeneously distributed points, numbered from 0 to 2N .In all the analytical description, we consider that the two subdomains Ω 1 and Ω 2 have the same number of points, respectively x 0 , ..., x N and x N , ..., x 2N .The interface point x N is common to the two domains, having different computed solutions u k 1,N and u k 2,N in each one of them.Evidently, we expect, at the convergence of the method, that u ∞ As done in the initial numerical tests in the section 2.2, an implicit Finite Difference scheme is used here.For the interior points of each one of the domains, we consider a second order discretization for the third spatial derivative in equation (1): which is valid for j = 2, ..., N − 2 in the case i = 1; for j = N + 2, ..., 2N − 2 in the case i = 2; and for j = 2, ..., 2N − 2 in the case i = ref .
In the above expression, α i,j is a given data (for example, the exact or the converged solution in the previous time step).
For the points near the boundaries, we use second order uncentered discretizations or the appropriate boundary condition.Considering that one boundary condition is written for the left boundary and two for the right one, we have to impose an uncentered discretization only for the second leftmost point of the domain.For example, for the point x 1 : and similarly to the other points near the boundaries.
In the resolution of the problem in Ω 1 , two IBCs are imposed (corresponding to Θ 2 and Θ 3 ) to the discrete equations for the points x N −1 and x N .On the other hand, in the resolution of the problem in Ω 2 , only one interface boundary condition is used (corresponding to Θ 1 ), being imposed to the point x N .
Remark : modification of the reference solution.Even if the DDM with the proposed IBCs is compatible with the monodomain problem (which we will see that is not the case), the solution of the DDM does not converge exactly to u ref , for a reason that does not depend on the expression of the IBCs, but on the fact that for each domain we write two boundary conditions in the right boundary and only one on the left boundary.We are using a second order centered discretization for the third spatial derivative (which uses a stencil of two points in each side of the central point), implying that we must write an uncentered discretization for the point x N +1 when solving the problem in Ω 2 .Therefore, this point does not satisfy the same discrete equation as in the reference problem.In order to avoid this incompatibility and allow us to study the behavior of the DDM, we will modify the discretization for the point u N +1 in the monodomain problem, using the same second-order uncentered expression : Figure 1 resumes the discretizations imposed to each point in the monodomain and the DDM problems, as described above:

Corrections for the approximate IBCs
When using approximate IBCs in a Schwarz method, one should guarantee that the converged solutions u * satisfy the same equation as the solution u ref of the monodomain problem.Nevertheless, one can easily see that, in the convergence, the solution u * does not satisfy the discrete equation ( 13) on the points where the IBCs are imposed (the poins As pointed out by [8], a finite difference discretization of the IBCs requires a special treatment to be consistent with the monodomain discretization.Therefore, we will formulate modified IBCs in order to avoid this problem: with θ i , θ i given by It is straightforward to verify that the DDM problem with these modifications in the IBCs insure that the converged solution u * satisfies, in every point, the same discrete equations as the solution u ref of the monodomain problem (10).
In addition, we notice that all the modification terms θ i , θ i , i = 1, 2, 3, are of order O(∆x) (they are composed of discrete versions of time derivatives and second spatial derivatives multiplied by ∆x).It is essential to insure that these terms are small, for the consistency with the approximate IBCs Θ i to be fulfilled.

Numerical tests for optimizing the IBCs (speed of convergence)
Our objective now is to optimize the IBCs in the sense of minimizing the number of iterations of our method until the convergence.We perform a very large set of tests in order to find the coefficient c that provide the fastest convergence.To start with, we make this study with fixed time step and space step, in order to analyze exclusively the influence of the coefficient.
As we are interested in the speed with which the solution of the DDM method converges to the reference solution, the criteria of convergence used is The range of tested coefficients is [−10.0,20.0] (chosen after initial tests to identify a proper interval), with a step equal to 0.1 between them (or even smaller, up to 0.005, in the regions near the optimal coefficients), and the maximal number of iterations is set to 100.

Test varying the initial time step and the interface position
As said above, in the first set of tests we consider a fixed time step ∆t = 20/2560 = 0.0078125 and a fixed mesh size ∆x = 12/500 = 0.024.Moreover, we consider two subsets of tests, in order to study the speed of convergence with different initial conditions and different sizes of the subdomains: In all the cases, the reference solution u ref is the solution of the monodomain problem (10).
The results are summarized in Figure 2, with the number of iterations plotted as function of the coefficient c (for the positive coefficients).We can see a very similar behavior of all the curves, with two minima whose position do not depend on t 0 and α (approximately, c = 0.20 and c = 4.5).For c < 0, the curves are very similar, with two minima located at c = −0.10 and c = −1.35,approximately.Moreover, the minima closest to zero (c = −0.10 and c = 0.20) are both associated with very discontinuous peaks, while the other two minima are associated with smoother curves.A detail of the curves around each positive minima are shown in Figures 2c -2d and 2e -2f.Finally, we remark that, for some curves, the minimal number of iterations is associated with the coefficients closest to zero, and, for other ones, to the other minimum, but the minimal number of iterations are very similar (between 5 and 7).

Tests varying ∆t and ∆x
After verifying that the method behaves similarly for several initial conditions (i.e., for several values of t 0 ) and various positions of the interface, now we keep these parameters fixed (t 0 = 0 and α = 0.5) and make new tests with different values of ∆t (with fixed ∆x = 12/250) and different values of ∆x (with fixed ∆t = 0.02).
The number of iterations as functions of the coefficient, for some of the tests, are shown in Figure 4, in the case of positive coefficients.The results for negative coefficients are similar.
Figure 5 presents the optimal positive coefficient for each ∆t or ∆x (for one fixed value for the other coefficient).Considering the observation we did before about the similar results (i.e. the number of iterations until the convergence) for the four optimal coefficients, we only took into account, for the construction of this curve, the positive minimum farther from zero: it was done because, as shown in Figure 4, these minima have a strong dependency on ∆t or ∆x, and we will seek to study this relation.5 suggests a dependence of the optimal coefficient on (∆t) ν and (∆x) η , with 0 ≤ ν ≤ 1 and η < 0. In fact, performing some regressions with ∆t or ∆x fixed, we conclude that ν = 2 3 and η = −1 provide really well-fitted regression curves (with the coefficients of determination R 2 bigger than 0.99), both for the negative and the positive coefficients (although each one of these cases correspond to different curves).Therefore, we seek to model a function In order to validate the expressions (15) and ( 16), we used them to compute the optimal coefficients for several points (∆t, ∆x), with ∆t ∈ [0.0005, 0.3] and ∆x ∈ [12/5000, 12/50].For almost all the points in the considered domain, the computed optimal coefficient provides a fast convergence to the monodomain solution, with less than 20 iterations, what is also observed in the case of the negative coefficients.The numbers of iterations observed are not always the smallest ones that we could find (cf.Figures 2 to 4), because the expressions (15) and ( 16) are regressions constructed from optimal coefficients obtained among a discrete set of possible values.Nevertheless, they give a very good approximation for the optimal c for each (∆t, ∆x), and one could search around a small region around the computed c opt to obtain an even faster convergence.
The results presented in this section show that the DDM proposed here is able to provide a fast convergence toward the solution of the monodomain problem.Furthermore, using the corrected IBCs (14), this convergence is exact.Therefore, we reached our goals of solving the dispersion equation in a finite domain divided in two subdomains.
Moreover, the results of the optimization tests are very satisfying regarding a more general application of our method.Firstly, for fixed spatial and temporal discretizations, we obtained optimal coefficients for the method independently of the initial solution and the size of the subdomains (i.e., independently of the initial instant and the position of the interface).Secondly, we obtained good regression expressions for the optimal coefficient as function of ∆t and ∆x, which could allow the application of the model, with fast convergence, in other computational frameworks.

Conclusion and outlook
We presented and implemented in this paper an additive Schwarz method for the resolution of an one dimensional dispersive evolution equation, using as interface conditions between the subdomains some operators constructed based on the exact transparent boundary conditions for this equation.Although not as accurate (in the role of TBCs) as the ones proposed in the works we are based on (providing better TBCs was not our objective here), these approximate conditions stand out for its simple form and implementation and the fast convergence that they provide for the Schwarz method.Moreover, we also proposed small corrections to them, which insure that the solution of the DDM problem converges exactly to the solution of the monodomain problem.Finally, we verified that the speed of convergence depends on the time step, the mesh size and the (only) coefficient for constructing the approximate interface conditions; thus, via an optimization process, we obtained and validated regression expressions that provide the optimal coefficient (i.e., the one that provides the fastest convergence) in function of ∆t and ∆x.
Natural continuations of the work presented here would be the study of the method using more complex operators as IBCs, using for example higher-orders Padé approximations for λ 2 /s and considering different approximations for left and right boundary conditions.Moreover, we can extend this study for other problems, for instance the linearized KdV equation, which adds an advective term on the equation solved here, as well as other models of wave propagation.Finally, we can seek the development of global in time Schwarz methods, using optimized Scwharz waveform relaxation methods.

Figure 1 :
Figure 1: Scheme indicating the discretization imposed to each point in the monodomain and the DDM problems

Figure 2 :
Figure 2: Number of iterations until the convergence as function of the coefficient of the TBC, in the case of positive coefficients

Figure 3
Figure3shows the evolution of the error, as function of the iterations, for the five positive coefficients c that gave the fastest convergences, for a fixed initial instant and a fixed position of the interface.For other values of t 0 and α this graph is similar, concerning the number of iterations and the fact that the convergence is more regular for the coefficients closest to zero, compared to the other optimal coefficients.

Figure 3 :
Figure 3: Error evolution with the iterations for the fastest results (a) Fixed ∆x = 12 250 (b) Fixed ∆t = 0.02

Figure 4 :
Figure 4: Number of iterations until the convergence as function of the coefficient of the TBC (for positive coefficients)

Figure 5 :
Figure 5: Optimal coefficients as function of the time step and the space step

Figure
Figure5suggests a dependence of the optimal coefficient on (∆t) ν and (∆x) η , with 0 ≤ ν ≤ 1 and η < 0. In fact, performing some regressions with ∆t or ∆x fixed, we conclude that ν = 2 3 and η = −1 provide really well-fitted regression curves (with the coefficients of determination R 2 bigger than 0.99), both for the negative and the positive coefficients (although each one of these cases correspond to different curves).Therefore, we seek to model a function