{'fraction_non_alphanumeric': 0.03285433368628044, 'fraction_numerical': 0.009015277567690213, 'mean_word_length': 4.3732119635890765, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An almost symmetric Strang splitting scheme for nonlinear evolution equations

In this paper we consider splitting methods for the time integration of parabolic and certain classes of hyperbolic partial differential equations, where one partial flow can not be computed exactly. Instead, we use a numerical approximation based on the linearization of the vector field. This is of interest in applications as it allows us to apply splitting methods to a wider class of problems from the sciences. However, in the situation described the classic Strang splitting scheme, while still a method of second order, is not longer symmetric. This, in turn, implies that the construction of higher order methods by composition is limited to order three only. To remedy this situation, based on previous work in the context of ordinary differential equations, we construct a class of Strang splitting schemes that are symmetric up to a desired order. We show rigorously that, under suitable assumptions on the nonlinearity, these methods are of second order and can then be used to construct higher order methods by composition. In addition, we illustrate the theoretical results by conducting numerical experiments for the Brusselator system and the KdV equation.


Introduction
In this paper we consider the time discretization of parabolic and certain classes of hyperbolic partial differential equations. More specifically, we assume that the considered problem can be written as the following abstract Cauchy problem u = Au + B(u), where A is a linear (but possibly unbounded) operator. In this context splitting methods can be applied if the partial flows generated by A and B have an analytical representation, or if an efficient algorithm for finding their exact solution is known. For a review of splitting methods we refer the reader to [17]. However, the assumption that the partial flow generated by the nonlinear operator B can be solved exactly is usually a very strong requirement. Most partial differential equations which are drawn from the sciences do neither admit an analytical solution nor can they be solved exactly in an efficient manner by some algorithm from the literature.
To remedy this deficiency of classic splitting methods, we propose and analyze splitting schemes which approximate the partial flow generated by B, i.e.
by that of an inhomogeneous linear differential equation. That is, instead of equation (2), we consider the linearized problem given by where for consistency reasons we have to assume that That is, once a value u is substituted, the flow corresponding to (3) can be computed efficiently. That such a linearized problem can be solved exactly is a much less stringent condition and thus splitting methods can potentially be applied to a larger class of problems.
For the remainder of the paper let us denote the flow corresponding to equation (3) by ϕ b(u ) t , which can be written explicitly, by employing the exponential and φ 1 function, as follows ϕ b(u ) t (u(0)) = e tb(u ) u(0) + tφ 1 (tb(u )) d, where the φ 1 function is given by In this context we can formulate the (classic) Strang splitting scheme as follows Such a splitting scheme has been employed, for example, to numerically solve the Vlasov-Poisson equations (see e.g. [4] and [5]) or the Davey-Stewartson equation (see e.g. [14]). Note, however, that in both of these cases the underlying structure of the equations guarantees that the Strang splitting scheme, as stated above, is still symmetric. Therefore, the usual construction of composition methods (as described e.g. in [7]) can be carried out without modification.
However, for many problems arising in the sciences, such as nonlinear reactiondiffusion equations or the KdV equation (both of which are discussed in section 6), such a simplification can not be assumed. Then the Strang splitting scheme above is, unfortunately, no longer symmetric. The lacking symmetry of the method does not severely affect performance; however, if composition is used as a means to construct higher order methods, symmetry is usually a necessary condition.
To remedy this situation a family of Strang splitting schemes were introduced in the context of ordinary differential equation (see [6]). The construction of these schemes is based on the fact that the Lie splitting algorithm composed with its adjoint method, which we denote by L * τ 2 , is a method of second order and symmetric. It is given by However, for our purpose it it is more convenient to represent u 1 = S τ u 0 as the solution of the following implicit equation Note that to solve this equation is computationally not attractive in practice. Therefore, it has been suggested to employ i ∈ N fixed-point iterations to compute an approximation to u 1 . Let us denote the resulting scheme by S (i) τ . The corresponding starting value (for the fixed-point iteration) is given by The one-step methods S (i) τ are not symmetric but they are symmetric up to order i; that is, they satisfy the following relation S Therefore, it is possible to construct composition methods of arbitrary (even) order, where i, and thus the computational cost of S (i) τ , increases linearly with the desired order p. For a more detailed description we refer the reader to [6].
In this paper, we will provide a rigorous convergence analysis which shows that for parabolic and certain classes of hyperbolic partial differential equations the method S (i) τ is of second order for i ≥ 2. This analysis conducted in the context of partial differential equations is significantly more involved than the analysis conducted in [6]. The assumptions we make in order to prove the main result will be that the classic Strang splitting scheme given in equation (6) is convergent of order two (i.e. that the classical Strang splitting scheme can be successfully applied to the problem under consideration) and that terms of the form are bounded in a suitable norm; hereũ is, in a sense to be made precise, close to the exact solution u(t). Note that in the linear case this assumption is equivalent to the fact that b 3 u(t) can be bounded; i.e. we have to bound three applications of the operator b to the exact solution. It is useful to keep that in mind in order to compare the assumptions made here to, for example, [11], where the analysis of the exact Strang splitting scheme in the linear case is provided.
To that end we will discuss the fixed-point iteration in the stiff case (section 2), and show second order convergence of S (i) τ for the hyperbolic (section 3) and parabolic (section 4) case. In section 5 we discuss the construction of higher order composition methods. Finally, we will show in section 6 that in the case of the KdV equation as well as the Brusselator system the predicted behavior agrees with the numerical experiments conducted.

Fixed-point iteration in the stiff case
The purpose of this section is to discuss the fixed-point iteration necessary to compute S (i) τ , as outlined in section 1, in more detail. In the non-stiff case a sufficiently small τ can always be found such that the iteration converges. However, in the present discussion it is necessary to employ a condition on the, in general, unbounded nonlinearity b. We will assume that a Lipschitz condition is satisfied in a Banach space (V, · ) that includes, in a sense to me made precise in section 3, functions close to the exact solution.
Then there exists a step size τ 0 > 0, which depends only on the Lipschitz contant L and on bounds of the partial flows, such that for all τ ∈ [0, τ 0 ) the fixed-point iteration convergences. In addition, it holds that Proof. We have to consider the fixed-point iteration given by By employing the variation-of-constants formula, we can write and by using the Lipschitz continuity of b, we get Thus, our fixed-point iteration convergences, if τ 0 , as given in the statement of the theorem, satisfies CLτ 0 ≤ 1.
It should be duly noted that although the step size is limited by the above theorem this is not a CFL condition. In our case the step size does only depend on properties of b(u), the exact solution as well as the scheme under consideration. However, the bound holds independently of the specific space discretization under consideration. In fact, the above result is what one would expect for a stable explicit scheme; namely that the step size is limited by the Lipschitz constant of the right-hand side.
In the next section we will employ heavily the following result; it states that a time step conducted by applying S τ to some u 0 results in a value that is different from u 0 by at most an order of τ . For the convenience of the reader we state and proof this well-known result in the next theorem.
Lemma 2. If the Lie splitting algorithm (7) is consistent of order 1 and τ is chosen such that the fixed-point iteration (9) converges, then where the constant C is independent of τ .
Proof. From equation (8) we know that S (1) τ is the composition of two Lie steps. Therefore, it holds that By employing the uniqueness and existence result for the fixed-point iteration, we get for any i ≥ 2. This bound gives the desired result, since under the assumptions of convergence of the fixed-point iteration we know that C 1 Lτ < 1 (see Theorem 1).

Convergence for hyperbolic differential equations
First, we discuss the condition on the initial value that is necessary to achieve a given regularity in time. This will give us an idea of the required regularity assumptions on the initial value, which is needed for the schemes described here to be convergent of second order. To accomplish this, in section 3.1, a suitable toy model is analyzed. In section 3.2 we will derive, for the hyperbolic case, the conditions under which the splitting scheme S (i) τ , for i ≥ 2, is convergent of order two. We will formulate a condition that depends on the regularity of the exact solution. More formally, we have to bound a number of applications of b(u ) and its derivative to the exact solution and assume that the (classic) Strang splitting scheme M τ is convergent of order two. For hyperbolic systems which do not change the regularity of the initial value this is equivalent to assuming some regularity for the initial value.

Regularity in time
An interesting question is how much regularity on the initial value is required such that a scheme can be convergent of second order in time. To elaborate on this question let us consider the following (nonlinear) advection equation The solution of this toy model can be written as where the advection speed depends on the solution as follows Now, let us investigate the time regularity of the solution u(t, x). To that end let us assume that for the initial value u 0 it is true that u 0 ∈ C 1,α for some α ∈ (0, 1). From (10) we deduce that u(t, x) is at least once differentiable in time. The resulting function ∂ t u(t, ·) lies in C 0,α ; this can be easily deduced from equation (11). Our goal is to show that is not bounded for any β > α as s tends to t. Rewriting the above expression we get The first term is Lipschitz continuous in time and therefore can be bounded. However, for the second term we can write where by defining z = x − v(t) and w = x − v(s) and using the estimate If we now choose a u 0 ∈ C 1,α that nowhere can be improved to a better Hölder exponent β > α and is bounded away from zero (which then also holds true for u(t, ·)) the term on the right hand side diverges as z tends to w and thus we have shown that the regularity of u(t, x) in time is at most C 1,α . That such an initial value can be found is shown in [13].
As a numerical approximation to the solution of (10), we consider the exact solution of on the interval [0, τ ]. This is not a practical numerical method as we have, in general, no algorithm that allows us to compute v(τ /2) exactly. However, it can be seen as the limit of a scheme, where we approximate v by approximating u(t, 0) in some suitable manner. In this case the numerical solution, i.e. the exact solution of equation (12), is given by the midpoint rule, i.e.
, which has to be compared to the exact solution, given in (11). By using the method of characteristics we can estimate the error as follows where L is the Lipschitz constant of u 0 . Therefore, the numerical scheme can be seen as the midpoint rule approximating an integral. It is well-known that the midpoint rule yields an approximation with a local error of size Cτ min (3,k+α) , if and only if the integrand lies in C k,α .

Convergence of the iterated Strang splitting scheme
The discussion in the previous section shows that we can expect the Strang splitting scheme to be second order accurate in time only if we require that u 0 ∈ C 2,1 . In the framework discussed in section 1 this means that we have to bound at least three applications of the operator b(u ) to the exact solution.
In fact our assumption, which we will state in this section, will turn out to be slightly more complicated due to the fact that we have to bound certain derivatives of b(·) as well.
Before stating the assumption that is used to prove the consistency result, let us remind the reader that for a function f : V × W → Z which is Fréchet differentiable with respect to the first component it holds that where it is understood that d ∈ V . Thus, we can separate the condition on the function f from the direction of the derivative. We will see below that this fact simplifies the condition needed to establish the correct order of the scheme significantly. Note that, clearly, this also holds true if the function f is linear in the second component and we employ the notation outlined in section 1.
Now, let us assume the following bound where ∈ {0, 1, 2}, i, j ∈ N. The functionũ are generic, i.e. different occurrences of such functions do not need to be indistinguishable, linear combinations of the exact solution u(s) as well as L τ 2 (u(s)) and S τ (u(s)) for some s ∈ [0, T ].
The following theorem shows consistency under the additional assumption that the (classic) Strang splitting scheme M τ is of the expected order.
Theorem 3. (Consistency) Let us consider the following abstract initial value problem and suppose that assumption (13) is satisfied. If, in addition, the Strang splitting scheme M τ is consistent of order two, we can conclude that S (i) τ , for i ≥ 2, is consistent of order two, i.e.
Proof. First, let us consider a single time step of the implicit scheme S τ , i.e. u 1 = S τ u 0 , can be written as (note that u 0 is a not to be confused with the initial value; it is the exact solution at a point in time to which a splitting operator is applied) . Now let us proceed by comparing our method, i.e. S τ , to the Strang splitting scheme M τ . The difference can be estimated as follows can be written as Therefore we have to compare to the following expression Let us first compare the homogeneous parts. By expanding this part in equation (15) To expand equation (14) let us first define g(τ ) = e which gives the desired expansion for appropriately chosen constants c i . Now we have to compare the corresponding terms of equal order.
Terms of first order. By employing the fundamental theorem of calculus we can show that where v 1/2 = u 1/2 − u 0 . For θ = 1 this gives the following operator identity The same expansion can be carried out for b(u 1 ). Finally, we have to bound where for the sake of brevity the remainder terms are denoted by R b(u 1/2 ) 2 and R b(u 1 ) 2 , respectively. All the terms can be bounded by assumption (13).
Terms of second order. We have to estimate A bound for the first and second term can easily be found. The third term can be rewritten as (13) is sufficient to show the desired bound.

Lemma 2 as well as assumption
Terms of third order. We have to bound R 3 as well as b 3 (u 1/2 )φ 3 (τ b(u 1/2 )). The first bound is immediate, for the second we note that since the commutator b 3 (u 1/2 ), φ 3 (τ b(u 1/2 )) vanishes, the desired result follows easily.
Finally, let us compare the two inhomogeneous terms which are given by The term of first order vanishes and the term of second order is easily shown to yield an additional order. In addition, we can easily bound the remainder terms.
Therefore, we have shown consistency for S τ . The extension of this result to S (i) τ , for i ≥ 2, follows immediately from Theorem 1.
In the previous theorem we have established consistency and are thus in a position to show convergence. Before stating this result let us note that since we assume convergence of the (classic) Strang splitting scheme M τ , it is not necessary to solve the significantly more involved problem of nonlinear stability. Much effort has been devoted to solve this problem for a variety of partial differential equations found in the literature (see e.g. [5] for the Vlasov-Poisson equations, [16] for the Schrödinger-Poisson and cubic nonlinear Schrödinger equation, and [10] for the KdV equation). Therefore, we show in the following theorem that if convergence has already been established for the (classic) Strang splitting scheme this is also true for the iterated scheme we consider in this paper. Proof. We once again proceed by comparing S (i) τ to M τ . First, let us write where we have denoted the numerical approximation of the (classic) Strang splitting scheme M τ at time kτ by u k . This immediately gives In the equation above we have used e k to denote the difference to the (classic) Strang splitting scheme at step k, i.e.
Since by assumption M τ is a stable numerical scheme we know that Now we know that e 0 = 0 and thus it can be easily shown that the solution of the recursion inequality given in equation (16) satisfies e k+1 ≤ Cτ 2 , which is the desired result.

Convergence for parabolic differential equations
Let us note that Theorem 3 and Theorem 4 have an immediate extension to the parabolic case. However, due to the parabolic smoothing property we can weaken assumption 13 in the following sense: there exists a constant C such that for all holds for 0 ≤ α ≤ 1. Againũ represents a generic function that is constructed in the manner described in section 3, but with the additional restriction that s ≥ t holds true (see the paragraph following assumption (13)).
The following theorem states that to achieve second order in time (for α < 1) it is sufficient to assume less regularity in the initial value than is required for the same convergence order in the case hyperbolic differential equations. (17) is satisfied and that M τ is convergent of order two. Then it holds that

Theorem 5. (Convergence in the parabolic case). Suppose that assumptions
If assumption (17) holds for α < 1, we get a scheme of second order, i.e.
Proof. We follow the proof presented in [11]. By using assumption (17) we can deduce from Theorem 3 that where C depends on the norm of the initial value. The above recursion can be estimated to give is the harmonic series. Since we know that e 0 = 0 we get the bound e k+1 ≤ Cτ 2 log k, as desired.

High order composition methods
In the previous sections we have exclusively considered the iterated Strang splitting scheme S (i) τ . Now, we will describe how this (almost) symmetric Strang splitting scheme can be used to construct splitting methods of higher order. To that end we will first discuss composition for a symmetric one-step method.
Suppose that a symmetric one-step method Φ τ is of even order p. Then a method of order p + 2 can be constructed by composition (see e.g. [7, p. 43]). More specifically, under suitable conditions on γ 1 , γ 2 , γ 3 , the method Φ γ 3 τ •Φ γ 2 τ •Φ γ 1 τ is of order p+2. Therefore, we are able to construct methods of arbitrarily high even order. The cost, in terms of a single evaluation of the corresponding second order method, is given by 3 p/2−1 . For p = 4, for example, the corresponding method is the well-known triple jump scheme.
That such an approach can be extended to the case of the almost symmetric Strang splitting method S (i) τ has been shown, in the case of non-stiff ordinary differential equations, in [6]. In this section we will extend the results presented in the before mentioned paper to the time integration of partial differential equations. Recall that if the following system of equations is satisfied the composition results in a scheme of order at least p + 1. By a symmetry argument we can then deduce that the order is indeed p + 2. This symmetry argument requires that Φ * τ −Φ τ = O (τ p+3 ), where we have used Φ * τ to denote the adjoint method of Φ τ , and γ 1 = γ 3 (see [6]).

The single real solution that simultaneously satisfies (18) and (19) is given by
Note that in order to perform this scheme we have to compute at least one negative time step. In fact, this is shown to be necessary if the desired order of the composition method is strictly larger than two (see e.g. [1]). While it is possible to take negative time steps in the case of hyperbolic equations, in parabolic problems this almost certainly leads to numerical instabilities as, for example, roundoff errors are exponentially amplified in each backward step. To remedy this situation in [9] and [3] it has been pointed out that the system of equations given by (18) and (19) admits complex solutions with positive real part. We will discuss the implementational ramifications of using complex values in the entire computation in section 6.1.
Now, let us turn to the main result of this section. That is, we will indicate how the analysis conducted above can be extended from the non-stiff ordinary differential equations considered in [6] to the partial differential equations considered in this paper. For such an analysis appropriate regularity assumptions are indispensable. In particular we have to control the remainder terms in the expansions. Now, under the assumption that the symmetric splitting scheme is consistent of order p, it holds that the local error of Φ γ i τ satisfy where R i (·) denotes the (bounded) remainder term. Note that this, in fact, requires one more application of the b operator as computed in the previous section. However, it should be clear that if the solution is sufficiently regular, such an expansion can be accomplished by the straightforward, but tedious, extension of sections 3 and 4. Then it is possible to generalize the proof given above by employing Lemma 2. Let E i denote the full error, i.e.
Then we have with a bounded remainder termR i . The total error of the composition method is then given by where Q is a bounded remainder term. Thus, we get To recover the result for the non-stiff case we have to expand C(u 1 ) and C(u 2 ). This can be done, for example, if the coefficient C(·) is differentiable with respect to the (local) initial value. All the estimates derived in sections 3 and 4 have this property.
To conclude this section, let us note that the argument provided establishes that, under suitable assumptions on the regularity of the solution, the composition method is of order p + 1. Due to the fact that S (p+2) τ is symmetric up to order p + 2, the method is in fact of order p + 2.

Numerical simulations
Having established the theoretical convergence rates for the Strang splitting scheme S (i) τ and its composition in the previous sections, we now turn our attention to illustrate these results by conducting numerical simulations. To that end we have chosen the well-known Brusselator as a parabolic system. The Brusselator consists of two coupled advection-reaction equations and provides a valuable prototype for many similar, but often considerably more complicated, systems that model the interaction of chemical species. As an example for a hyperbolic system we have chosen the KdV equation in a single dimension. This equation is then split in such a manner that the Burgers type nonlinearity can be solved as a position-dependent but linear advection. For both of these examples a C++ program has been written that employs the fftw library and is parallelized using OpenMP.

A parabolic system: the Brusselator
As a example of a parabolic system we consider the Brusselator (see e.g. [8, p. 152]). This system is given by where u(t, x, y) and v(t, x, y) are the two unknowns that usually represent the concentration of a chemical species. We have chosen the reaction parameters as β = 4.4, γ = 3.4, δ = 1, and employ a relatively weak diffusion with α = 10 −2 . The equations above have to be supplemented by suitable initial conditions. For the purpose of the numerical simulations conducted in this section the following initial values are used and all computations are carried out on the domain [0, 1] 2 with periodic boundary conditions.
In the context of the splitting scheme outlined in section 1, let us define The corresponding flow can be computed efficiently by employing two discrete Fourier transforms. In addition, the nonlinearity is represented by (note that this choice is not unique) Thus, we have to compute the solution of a linear system in two variables which can be done analytically. Note that in this case we could, in fact, find, by a rather tedious calculation, an analytical solution for the complete nonlinearity. However, this is no longer possible if one considers either additional variables, i.e. additional chemical species, or alternatively a more complicated nonlinearity. However, in such a case the splitting scheme outlined here would not suffer any further difficulty or even a loss of efficiency as the nonlinearity in question still reduces to a linear system that is easily solved by standard methods.
In the numerical simulation we will compare three splitting schemes. The (classic) Strang splitting scheme M τ which is expected to be of order two. The naive triple jump scheme which is constructed by composition from M τ . Since M τ is not symmetric we expect this to be a third order scheme. Finally, we will consider the triple jump scheme constructed by composition from S τ . Since S (4) τ is symmetric up to order four we expect that this approach results in a fourth order scheme. In fact, this is the main motivation for our approach as it enables us to construct high order (for example, fourth or sixth order methods) by employing the well known composition rules. The numerical results shown in Figure 1 confirm the expected order for all numerical schemes considered.  From Figure 1 we can clearly see that if similar step sizes are taken the error made by the the fourth order iterated triple jump scheme is significantly less than the error of both the third order triple jump scheme as well as the second order (classic) Strang splitting scheme. However, to discuss whether the high order schemes constructed here can also provide a significant gain in efficiency we have to plot the simulation time as a function of the error. This is done in Figure 2.
It is shown that for high precision requirement (or equivalently long integration times) the use of the fourth order iterated triple jump scheme (with 4 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 error  iterations, i.e. based on S (4) τ ) results in a significant increase in efficiency. Also note that our fourth order scheme is superior to the naive triple jump scheme for almost any precision requirements. However, due to the overhead involved in the use of complex arithmetics for small precision requirements and short integration intervals the Strang splitting scheme is clearly the preferred choice.

A hyperbolic system: the KdV equation
As an example of a hyperbolic system we consider the KdV (Korteweg-de Vries) equation in a single dimension. It is given by It is shown in [12] and [2] that for sufficiently regular initial values the regularity of the solution to the KdV equation is not diminished as it is evolved in time (although this does not rule out the appearance of high frequency oscillations). For the purpose of this section we will consider two initial values with extremely different dynamical behavior.
First, let us consider the following initial value The exact solution is a soliton that travels to the right with speed 4, i.e., the solution is given by (see e.g. [15]) For our numerical studies we consider the domain [−20, 20] in order to limit artifacts which originate from the fact that periodic boundary conditions are imposed on the domain of finite length that is used in the numerical simulations. We integrate up to a final time T = 0.4.
Second, we consider an initial value where oscillations appear. The so called Schwartzian initial value is given by (see e.g. [14]) In this instance we integrate only up to T = 0.05 and employ the domain [−4π, 4π] with periodic boundary conditions (following the same argument as given above). Now let us turn our attention to the splitting scheme employed. In the framework of section 1 the linear problem is defined as x and can be solved efficiently by employing two discrete Fourier transforms. Instead of the full nonlinearity we numerically compute the evolution corre- This is significantly less involved than solving the full Burgers type nonlinearity. In fact, [14] states that solving the full nonlinearity invalidates splitting as a viable approach to numerically solve the KP (Kadomtsev-Petviashvili) equation (note that, compared to the KdV equation, the two dimensional KP equation includes an additional diffusive term in the additional variable). However, the linear problem given by equation (22) can be computed numerically to high precision, for example, by using the (exact) exponential of a finite difference stencil (as it is often done in the context of exponential integrators, for example). We intended to employ the Expokit package to compute the exponential of the a seven-point stencil that is used to approximate the first derivative. However, in case of the soliton solution the results were not satisfactory except for very small step sizes. Therefore, we used the SUNDIALS CVODE solver, with a prescribed tolerance of 10 −12 , for this example 1 .
The numerical results are shown in Figure 3 (for the soliton initial value given in equation (20)) and in Figure 4 (for the Schwartzian initial value given in equation (21)). Note that the numerical simulation matches the predicted order of the schemes studied very well. Thus, we conclude that the data obtained are clearly consistent with the analysis conducted in sections 3 and 5.

Conclusion and Outlook
We have rigorously shown that the proposed iterated Strang splitting scheme is of second order in time and due to its symmetry properties can be used to construct methods of arbitrary (even) order by composition. The main assumption we have made is that the classic Strang splitting method is convergent of order two (i.e., that we deal with a problem for which applying a splitting scheme of second order is sensible). Further, a technical assumption on the nonlinearity is made. This assumption reduces, in the linear limit, to the statement that we have to bound an appropriate number of application  of the operator in question to the exact solution (assumptions of that type have been used in much of the literature to show convergence of splitting methods for linear partial differential equations).
In addition, we have provided an argument demonstrating that our iterated Strang splitting can be used, similar to the case of ordinary differential equations, to construct higher order methods by composition if sufficient regularity of the exact solution can be assumed. This has been verified up to order four for both the Brusselator system (a parabolic problem) and the KdV equation (a hyperbolic partial differential equation). In both instances the iterated fourth order method is shown to provide superior performance, in case of medium to high precision requirements (or equivalently long integration times), compared to both the (classic) Strang splitting scheme as well as the (classic) triple jump scheme (which is a method of order three). We also conclude that the necessity of using complex precision arithmetics,  is computed by comparing the numerical solution for a given τ with a reference solution for which a sufficiently small time step has been chosen. In addition, for each scheme a line with slope equal to the expected order is displayed as well.
in the case of parabolic problems, does not negate the performance gain we expect from higher order methods.
In this paper we have only considered methods of order up to four. However, composition can be used to construct methods of arbitrary (even) order. To provide a clear picture of the efficiency gain expected for such high order methods, especially in the context of the semi-Lagrangian approach discussed in this paper, longer integration times as well as a space discretization with significantly more grid points has to be considered. Since we have not considered parallelization and other computing aspects in this paper, we will consider such an implementation as a subject of further research.