Integrating Semilinear Wave Problems with Time-Dependent Boundary Values Using Arbitrarily High-Order Splitting Methods

: The initial boundary-value problem associated to a semilinear wave equation with time-dependent boundary values was approximated by using the method of lines. Time integration is achieved by means of an explicit time method obtained from an arbitrarily high-order splitting scheme. We propose a technique to incorporate the boundary values that is more accurate than the one obtained in the standard way, which is clearly seen in the numerical experiments. We prove the consistency and convergence, with the same order of the splitting method, of the full discretization carried out with this technique. Although we performed mathematical analysis under the hypothesis that the source term was Lipschitz-continuous, numerical experiments show that this technique works in more general cases.


Introduction
We consider full discretizations by means of the method of lines of a semilinear second order in time-evolutionary problems with time-dependent boundary values. For this, we first discretize in space by means of finite differences, obtaining a system of ordinary differential equations.
For time integration, we rewrite the semidiscrete system as a first-order system in time and apply an arbitrary splitting scheme. Useful descriptions of splitting methods can be found in review articles [1][2][3]. Splitting schemes are especially useful in the field of geometric integration. In fact, splitting integrators preserve the structural properties of the original problem flow for as long as the flow of intermediate problems does. The good performance of geometric integrators in the long-term integration of systems of Hamiltonian ODEs is well-demonstrated in [4,5].
In this paper, we thus obtain a time integrator that is explicit and has the advantage of being cheap to implement, but with the disadvantage that its stability interval is finite. However, for these second-order in-time evolutionary problems, the stability condition is acceptable, and the step size in time and space may be taken to be of a similar size. It is also possible to use implicit methods, (see, for example, [6]), where Gautschi methods are studied avoiding the order-reduction phenomenon that appears with these methods.
The way in which a splitting method works requires three steps [3]: first, by choosing how to split the problem into several simpler intermediate problems, integrating each intermediate problem either exactly or approximately, and lastly composing the solution of the intermediate problems to obtain an approximation of a certain order of the original problem.
Denoting by h the step size in space and by k the step size in time, and separating the problem, we integrate into two intermediate problems with exact flows given by Φ [1] h,t and Φ [2] h,t . We consider a general splitting integrator with m stages and coefficients a j , b j , j = 1, . . . , m, Ψ h,k = Φ [2] h,b m k Φ [1] h,a m k . . . Φ [1] h,a 2 k Φ [2] h,b 1 k Φ [1] h,a 1 k . (1) In the case of the standard method of lines, each must be integrated in time (either exactly or numerically using a sufficiently accurate quadrature formula); however, for this, the term due to the space discretization of the nonvanishing boundary values must be addressed in the same way as the source term is. Since each stage of the splitting method applied to the spatial discretization is exactly integrated, we deduced that the optimal order was achieved, at least for a fixed thickness of the spatial discretization.
We propose in this paper a technique to cheaply and effectively incorporate the boundary values to the time integration carried out by the splitting method. This technique is consistent with these values being used to approximate the spatial differential operator on the boundary. We prove consistency and convergence with optimal order of the full discretization obtained with our technique under the hypothesis that the nonlinear term of the original second order in time problem was Lipschitz-continuous; however, numerical experiments showed that this hypothesis is not necessary in practice. Moreover, numerical experiments clarified the superiority of this technique compared to the use of the standard line method, with minor errors refining the discretization in both space and time.
Throughout this paper, several constants that are independent of the time step of the time integration could be likewise denoted (usually with the letter C, and possibly with some subscript).
The paper is organized as follows. The studied problem and spatial discretization are introduced in Section 2. Section 3 is devoted to the standard method of lines, time discretization being performed with a splitting method. In Section 4, we explain the alternative method that we propose to incorporate boundary values when implementing the splitting method. Numerical experiments that clearly show the better accuracy of the proposed method versus the standard method of lines are carried out in Section 5. Mathematical analysis of the convergence is developed in Section 6, where consistency is proved, and in Section 7, where convergence is stated along a brief review of the needed linear stability.

Partial Differential Equation
Let X and Y be Hilbert spaces, D(A) ⊂ X a dense subspace, and A : D(A) ⊂ X → X, B : D(A) ⊂ X → Y two closed linear operators. We consider the abstract second-order in time semilinear equation given by where source term f : [0, T] × D(A) → X is a smooth function that is generally nonlinear.
In practice, data f and g, solution u, and operators A and B are defined on a domain Ω ⊂ R n , and they could depend on spatial variables. We did not make this dependence explicit in the abstract formulation (2) in order to simplify the notation.
We make the following hypotheses on operators A and B: possesses a unique solution denoted by x = K(0)v, and there exists a constant C, such that linear operator (A4) Solution u in (2) satisfies u(t) ∈ D(A) for t ∈ [0, T] and is smooth enough in time.
(A5) Source term f (t, u) is a Lipschitz-continuous function with respect to variable u.

Spatial Discretization
Our first step to discretize (2) by means of the method of lines is spatial discretization. Let h ∈ (0, h 0 ] be a parameter that is used to measure the thickness of spatial discretization. We assumed that X h is a family of finite-dimensional spaces that approximate X. The discrete norm in X h is denoted by · h . Moreover, there is a subspace X h,0 ⊂ X h where the elements of D(A 0 ) are well-approximated by using P h : D(A) ⊂ X → X h,0 , that is, we assumed that P h u is the best approximation when u ∈ D(A 0 ). The boundary values are discretized by means of the linear operator Q h : Operator A is approximated by using operators A h : In practice, if we look for solution u ∈ D(A) of steady-state problem where F ∈ X and g ∈ Y, we cannot obtain P h u. Instead, we can compute R h u satisfying In order to discretize the source term, we suppose that function f can be defined in space X h,0 . That is, we can consider f : [0, T] × X h,0 → X h,0 , and for each u ∈ X and t ∈ [0, T].
With this spatial discretization, we obtain semidiscrete ordinary differential system We make the following hypotheses: (H1) There exists a constant C independent of h, such that, for u ∈ D(A) and small enough h, (H2) Operator A h,0 is symmetric and negative definite. Let S h,0 be the symmetric and positive definite operator, such that S 2 h,0 = −A h,0 . We also assumed that A h,0 and S h,0 were invertible and their inverses were uniformly bounded on h. (H3) There exists a subspace Z ⊂ D(A) with norm · Z , such that With these hypotheses, we can prove that the solution of (5) is a good approximation of the one of (2). Theorem 1. We assumed Hypotheses (A1-A5) and (H1-H4), that where C only depends on T, u, u and Lipschitz constant L.

Proof.
Applying P h to (2), considering (4), and making the difference with (5), Rewriting (6) as a first-order in-time problem, Using that cos(tS h,0 ) , the variation of the constant formula and the vanishing initial conditions, we deduce that Now, using (H2-H4), and that u ∈ C([0, T], Z), we have Applying the Gronwall lemma, we obtain for t ∈ [0, T], Therefore, we can obtain a good approximation of the solution of (2) by considering the solution of (5) with a small enough value of h. Next, we use a time integrator to achieve full discretization. For this, we consider a splitting scheme. In the next two sections, we study two different ways of incorporating the boundary values with the full discretization.

Full Discretization: Standard Method of Lines
First, we rewrite (5) as a first-order differential problem. We denote u h (t) = [u 1,h (t), u 2,h (t)] = [u h (t), u h (t)] and obtain system of which the exact flow is denoted by u h (t) = Φ h,t u h (0). Second, we apply to (7) a splitting scheme in the usual way. We then choose a split of (7) in two intermediate problems. The first is of which the exact flow is denoted by z h (t) = Φ I h h,t z h (0), and the second is of which the exact flow is denoted by Second, we consider z 1,j,h,n (0) = z 1,j,h,n (a j k), z 2,j,h,n (0) = z 2,j,h,n (a j k), and we advance a step b j k in time, z 1,j,h,n (b j k) = z 1,j,h,n (a j k), z 2,j,h,n (b j k) = z 2,j,h,n (a j k) + b j kA h,0 z 1,j,h,n (a j k) With this full discretization, we obtain order p in time since this is the order of the splitting method that we use, and the flows of the intermediate problems are exactly calculated when both integrals must be exactly calculated. In any case, we can always use a quadrature rule with the same accuracy as that of the splitting method. Convergence must be obtained in the discrete energy norm, and a suitable stability hypothesis is needed, similarly to the case of the discretization studied in the next section.

Full Discretization: An Alternative Way to Incorporate Boundary Values
The standard method of lines studied in Section 3 seems to be optimal, especially when integrals in (11) are exactly calculated. However, the two integrals had very different origins, since one of them came from the source term and the other one from the boundary value. Furthermore, the integral that came from the discretization of the boundary values in (11) can be arbitrarily large when spatial discretization is refined. This is because operator A h arises from the approximation at the boundary of differential operator A, which is unbounded.
We now introduce full discretization that more suitably incorporates the boundary values, as we show in the numerical experiments in Section 5. For this purpose, we need a more consistent notation with the fact that we discretize an initial boundary value problem (cf. [11]).
Suppose that In this way, semidiscrete Problem (5) can be rewritten as which is more similar to the original problem.
With this notation, we rewrite Problem (12) as first-order differential system of which the exact flow is given by u h (t) = Φ h,t u h (0), as in Section 3.
For time discretization, we consider the same splitting scheme as that in Section 3, and the challenge is to obtain a suitable way with which to incorporate the boundary values. For this, we split Problem (13) of which the exact flow is given by , and the second is We supposed that we computed . As in Section 3, we use a general splitting method with m stages, with coefficients a j , b j , j = 1, . . . , m. Therefore, To more explicitly write the j-stage in (16), we first consider initial-value problem and we assign boundary value Then, where using (18), we can calculate and we deduce that its solution at s = b j k is

Remark 2.
If we compare Formulas (11) and (19), the only difference is the treatment of the boundary values. The two ways of dealing with the boundary are for the standard method of lines and for the one that we propose, respectively.
Obviously, the second option is much simpler since the first option may even require to numerically evaluate the integral. We see in the numerical experiments in the next section that the second option also allows for obtaining much more precise results.
We prove in Sections 6 and 7 that, if the full discretization described in this section is used to approximate Problem (12), then the method is convergent, and the optimal order p of the splitting method is achieved.

Numerical Experiments
For the numerical experiments in this section, we consider splitting integrators of m stages, with coefficients a j , b j , j = 1, . . . , m, and b m = 0, which are particular cases of the symmetric ones [3]. More specifically, we use the Strang method with order p = 2, and coefficients a 1 = a 2 = 0.5 and b 1 = 1 and two other methods that are particular cases of symmetric-splitting methods, whose coefficients are given in the following way: Let l ∈ N be an even number; then, we consider a method with m = l + 1 stages satisfying For l = 6, we consider method Ψ S4 with order p = 4 obtained from and, for l = 10, method Ψ S6 with order p = 6, given by coefficients

Numerical Experiment: Test 1
We consider test problem where This problem can be written in abstract format (2) by taking We consider grid We take X h = R J+2 , and elements u h ∈ X h are denoted by u h = [u 0 , u 1 , . . . , u J , u J+1 ] T . In this way, X h,0 = {u h ∈ X h such that u 0 = u J+1 = 0} but, for the sake of simplicity, we use X h,0 = R J and their elements are denoted by Operator P h : D(A) → X h,0 is given by P h u = [u(x 1 ), . . . , u(x J )] T , and operator Q h : Y → X h,b is given by Q h (u 0 , u J+1 ) = (u 0 , 0, . . . , 0, u J+1 ). Operator A h arises from the approximation of the second-order derivative in space by using central finite differences, that is, Now, by using the previous notation, we can write which is a symmetric and definite negative operator on X h,0 . Hypothesis (H3) can be verified in the standard way by using Taylor series for the local truncation error, with ε h = O(h 2 ) and Z = C 2 (0, 1). Taking into account that the exact solution of (20) is a second-order polynomial in variable x, there is no spatial error, and we can focus on the error due to time discretization.
In this problem, for For the three symmetric-splitting methods, we compare the errors in the energy norm for the choice of boundary B(g) and for the EX(g) option in Remark 2, for values 1, 3, and 5 of µ. As Figure 1 shows, in all cases, option B(g) (solid line) obtained smaller errors than EX(g) did (dashed line). Moreover, the difference was more noticeable for Ψ S4 and Ψ S6 . Dependence on the size of boundary function g was also observed; for example, the difference between errors of B(g) and EX(g) was larger for µ = 1 than that for values µ = 3 and µ = 5, where the values taken by g were smaller.
In addition, as shown in Figure 1 the slopes of the lines were 2 (Strang), 4 (Ψ S4 ) and 6 (Ψ S6 ), which coincides with the expected optimal order of the three methods.

Numerical Experiment: Test 2
We consider test problem where x ∈ [0, 1], t ∈ [0, T], of which the exact solution is u(x, t) = e −t 2 (x 2 + 1). In this problem, a primitive of function f (t, u) = − sin(u) + e −t 2 ((−2 + 4t 2 )(x 2 + 1) − 2) + sin(e −t 2 (x 2 + 1)) cannot be exactly expressed using elementary functions, so we used a quadrature formula of appropriate order. In the calculations to obtain the data in Figure 2 we used the 3-point Gaussian quadrature of order 6, denoted by G3(g). Errors are compared for Strang's method (black), Ψ S4 (red), and Ψ S6 (blue), and for options B(g) (solid line) and G3(g) (dashed line). The slopes of the lines indicate order 2 for Strang's method, order 4 for Ψ S4 , and order 6 for Ψ S6 . Errors are always smaller for option B(g) than those for option G3(g). Moreover, this difference was more pronounced as the order of the method increased.

Numerical Experiment: Test 3
Now, we study the behavior of the error of both methods when spatial discretization is refined. The source term of semidiscrete Problem (5) grows when h → 0 due to the boundary . However, it is expected that this growth has no influence when it is treated as part of the discretization of operator A at the boundary, as in the method we propose in Section 4.
For this, we consider test problem where x ∈ [0, 1], t ∈ [0, T], of which the exact solution is u(x, t) = e −t+x . In this experiment, spatial error dominates, and we expected to observe order 2 of spatial discretization. Figure 3 shows errors for the splitting using B(g) in continuous line and EX(g) in dashed line, in red for Ψ S4 and in blue Ψ S6 . The values of h = k = 1/N were used for N = 25, 50, 100, 200. In this way, we remained in the stability interval of both methods. In the two methods with B(g), order 2 of spatial discretization was observed; when using the EX(g) option, the errors were larger, and order 2 was lost when h decreased.

Numerical Experiment: Test 4
Although theoretical results were shown for the case where f was Lipschitz-continuous, we see in a couple of examples that it also works if f is locally Lipschitz-continuous. We now consider test problem where x ∈ [0, 1], t ∈ [0, T], of which the exact solution is u(x, t) = e −t (x 2 + 1). Table 1 shows the time errors for the three symmetric-splitting methods with h = 1/100 for final time T = 1. The evolution of log 2 (error(k)/error(k/2)) is displayed in Table 2. Orders 2, 4, and 6 are shown for Strang's method, Ψ S4 , and Ψ S6 , respectively. The little loss of order in the lower-right corner was due to the influence of rounding errors.
Operator P h : D(A) → X h,0 is given by P h u = [u(x 1 , y 1 ), u(x 1 , y 2 ), . . . , u(x J , y J−1 ), u(x J , y J )] T , and A h arises from the approximation of the Laplacian operator by using central finite differences in each spatial direction, that is, considering second-order spatial discretization where I J is the identity matrix, and Taking into account that the exact solution is a second-order polynomial in variables x, y, there is no spatial error; then, we can again focus on the error due to time discretization. Table 3 shows the time errors for the three symmetric-splitting methods, with h = 1/100 and for final time T = 1. The missing value for the Strang method and k = 1/100 was due to the instability of the numerical solution in this case; see Section 7.1. The evolution of log 2 (error(k)/error(k/2)) is displayed in Table 4. Orders 2, 4, and 6 are shown for Strang's method, Ψ S4 and Ψ S6 , respectively.

Consistency Correctly Incorporating Boundary Values
Here, we deduce consistency in the energy norm of the implementation of a splitting method with the boundary values that we chose in Section 4.
As a first step in the proof of consistency, we introduce ordinary differential system We have u 3,h (t) = Q h g (t), u 3,h (0) = Q h g(0). and, therefore, Then, for the first two components of (25), we deduce that which is the same problem as (7) (and as (12) with the notation of Section 4). We also deduce that being u h (t) the solution of (5). We split (25) into two intermediate problems that are similar to the ones used in Section 4, and we applied to it the same splitting method. The solution of (25) was approximated with order p of the splitting method. This particularly is true for the two first components that match those of the solution of (12). Therefore, to prove consistency, it suffices to see that the obtained approximations for the first two components are the same as those described in Section 4 with the choice of boundary values made in (18).
We choose the following split of Problem (25). The first intermediate problem is of which the exact flow is denoted as v h (t) = Φ [1] h,t v h (0), and the second is of which the exact flow is denoted as w h (t) = Φ [2] h,t w h (0). Assuming that approximation u h,n = [u 1,h,n , u 2,h,n , u 3,h,n ] T ≈ [u 1,h (t n ), u 2,h (t n ), u 3,h (t n )] T is already calculated and, as in previous sections, we apply a general splitting of order p, h,a j k w j−1,h,n , j = 1, . . . , m, where w j,h,n = [w 1,j,h,n , w 2,j,h,n , w 3,j,h,n ] T . The performance of a time step k > 0 is as follows. For each j = 1, . . . , m, the first problem to be solved is Since the third component of (31) provides the boundary value (18), the first and second components of (31) are the same as (17). Then, the second problem is w 1,j,h,n (s) = 0, s ∈ [0, b j k], w 2,j,h,n (s) = A h,0 w 1,j,h,n + A h w 3,j,h,n + f (t n + k j−1 ∑ l=1 b l + s, w 1,j,h,n (s)), w 3,j,h,n (s) = 0, w 1,j,h,n (0) = v 1,j,h,n (a j k), w 2,j,h,n (0) = v 2,j,h,n (a j k), w 3,j,h,n (0) = v 3,j,h,n (a j k), whose solution at s = b j k is The first two components of (32) are the same as those of (19).

Remark 3.
Although the splitting approximation (30) provides the same approximation as (16) through its first two components, it is not convenient to use it in practice. It is more useful to use the implementation of Section 4, which can be carried out with minimal modifications to the standard method of lines. Now, we consider u h (t) = [u 1,h (t), u 2,h (t), u 3,h (t)] T , the solution of (25). We can easily deduce that is the solution of ordinary differential system We deduce that u 3,h (t) = Q h g(t), and that Therefore, the solution of Problem (33) is the appropriate one to calculate the energy norm of the solution of (25). Now, we see that approximating the solution of Problem (33) by means of a splitting method is equivalent to applying the same change of variables to (30), as is stated in (37).

The first intermediate problem is
of which the exact flow is denoted as v h (t) = Φ [1] t v h (0) and the second intermediate problem of which the exact flow is denoted as w h (t) = Φ [2] t w h (0). As in the previous sections, we use a general splitting method with m stages and order p, with coefficients a j , b j , j = 1, . . . , m. Therefore, h,a j k w j−1,h,n , j = 1, . . . , m, where w j,h,n = [ w 1,j,h,n , w 2,j,h,n , w 3,j,h,n ] T . We now prove that where u h,n = [u 1,h,n , u 2,h,n , u 3,h,n ] T is the solution of the splitting given by (30). For j = 1, . . . , m, the following problems are solved.

Theorem 2.
We assume that the time discretization of (12) is obtained by means of a splitting method of order p, applied as described in Formulas (16)- (17), with the choice of intermediate boundary values (18). Let u h (t n ) = [u 1,h (t n ), u 2,h (t n )] T be the value at t n of the solution of (13), and let u h,n = [u 1,h,n , u 2,h,n ] T be its approximation obtained with the splitting method (16) by taking a step of size k starting from the exact value u h (t n−1 ) = [u 1,h (t n−1 ), u 2,h (t n−1 )] T . Then, local error ρ h,n = u h (t n ) − u h,n satisfies where u h,n = [ u 1,h,n , u 2,h,n , u 3,h,n ] T is the approximation obtained at t = t n with the splitting method (36) by taking a step of size k starting from the exact value at t n−1 of the solution of (33) u h (t n−1 ) = [ u 1,h (t n−1 ), u 2,h (t n−1 ), u 3,h (t n−1 )] T . Therefore, the result is a consequence that p is the order of the splitting method.

Stability
To achieve convergence in energy norm, we need time discretization to be stable. In our case, it is sufficient to have linear stability. That is, it is enough that time discretization with the splitting method is stable for fully homogeneous linear problem corresponding to the space discretization of (2) with vanishing boundary values and source term.
To test linear stability, we first apply the splitting method to the harmonic oscillator y + λ 2 y = 0, λ > 0. We denote [p, q] T = [λy, y ] T and we consider the standard splitting If we now apply a splitting method with a time step k > 0, we obtain numerical method where ω = kλ > 0. Matrix R, of which the elements are polynomials in the ω variable, is called the stability matrix and is given by an it can be computed as where the product of matrices must be calculated in the correct order. To obtain stability for the harmonic oscillator, the boundedness of the powers of stability matrix (43) is required. For this, the following definition of "stability interval" is very useful (cf. [12,13]).
In this way, we have linear stability for test Problem (42) when kλ < ω * . The larger the value of ω * is, the less restrictive the stability condition is.
Regarding Problem (41), linear stability in the energy norm states that the powers of matrix are bounded in the matrix norm induced by the discrete norm in X h,0 , that is, if T > 0 is fixed, where C is a constant independent of h, n and k when nk ≤ T. Therefore, we need that k|λ h,0 | < ω * for all λ h,0 ∈ ρ(S h,0 ). Taking into account that S h,0 is symmetric and positive definite, any of its eigenvalues are positive; therefore, it is enough that is satisfied, where λ * h,0 is the largest eigenvalue of S h,0 .

Remark 5.
We can now deduce the ratio between parameters k and h that must be satisfied to have stability in the energy norm for the numerical experiments in Section 5. The eigenvalues of the tridiagonal matrix diag(1, −2, 1), which appears in the matrix (21) used in tests from 1 to 4, are given by −2 + 2 cos(jπ/(J + 1)), j = 1, 2, ..., J, and they all belong to the interval (−4, 0) (see for example [16,17]). We conclude that the largest eigenvalue of S h,0 satisfies λ * h,0 < 2 h and, to achieve stability, it suffices that Similarly, for the two-dimensional problem in Test 5, the eigenvalues of the block tridiagonal matrix appearing in (24) are included in interval (−8, 0), which means that, in this case, the largest and, to achieve stability, it suffices that Stability conditions (46) and (47) for the three splitting methods considered in the numerical experiments of Section 5 are given in Table 5. The missing value for the Strang method for h = 1/100 and k = 1/100 in Table 3 is due to the instability of the numerical solution because the stability condition was not fulfilled.

Convergence
We now study the convergence of the full discretization of Section 4.
Theorem 3. Assume that the time discretization of (12) is obtained by means of a splitting method of order p, applied as described in Formulas (16)- (17), with the choice of intermediate boundary values (18). Let u h (t n ) = [u 1,h (t n ), u 2,h (t n )] T be the value at t n of solution of (13), and let u h,n = [u 1,h,n , u 2,h,n ] T be its approximation obtained with the splitting method (16). Assume also that (45) and the linear stability condition (44) are satisfied. Then, global error T be the first two components of (33). We showed in On the other hand, let u h,n = [ u 1,h,n , u 2,h,n ] T be the first two components of (36). We also showed in Section 6 that u h,n = [S h,0 u 1,h,n , u 2,h,n ] T . Therefore, we can obtain u h,n from Section 4, using Equations (17) By using recursive reasoning, We now define w 0,h,n = u h (t n ), and To use an inductive reasoning, we first consider where, if (45) is satisfied and taking into account Hypothesis (H2), constant C 1,1 is independent on h, k, and n.
The global error is given by e h,n = e 1,h,n e 2,h,n = u 1,h (t n ) u 2,h (t n ) − u 1,h,n u 2,h,n = u 1,h (t n ) u 2,h (t n ) − u 1,h,n u 2,h,n + u 1,h,n u 2,h,n − u 1,h,n u 2,h,n = ae h,n + R(kS h,0 ) e h,n−1 + F n where ρ h,n is the semidiscrete local error at t = t n . Therefore, we deduce that e h,n = n ∑ j=1 R n−j (kS h,0 ) ae h,j + n ∑ j=1 R n−j (kS h,0 )F j .
The proof of convergence is achieved by using the discrete Gronwall lemma. Funding: This research was funded by the Ministerio de Ciencia y Educación, grant number PGC2018-101443-B-I00, and the first author by Consejería de Educación, Junta de Castilla y León y Leónand Feder funds, grant number VA193P20.

Conflicts of Interest:
The authors declare no conflict of interest.