1 Introduction
In this paper we derive a Lietype splitting integrator for abstract semilinear boundary coupled systems, and prove first order error estimates for the time integrator by extending the results of [CsEF] from the linear case. The main idea of our algorithm is to decouple the two nonlinear problems appearing in the original coupled system, while maintaining stability of the boundary coupling. More precisely, we combine the splitting scheme presented in [CsEF] with the appropriate handling of the nonlinear terms. We use techniques from operator semigroup theory to prove the firstorder convergence in the following abstract setting.
We consider the abstract semilinear boundary coupled systems of the form:
(1) 
where are linear operators on the Banach spaces and , respectively, , are suitable functions, and the two unknown functions and are related via the linear coupling operator acting between (subspaces of) and . A typical setting would be that is a tracetype operator between the space (for the bulk dynamics) and the boundary space (for the surface dynamics). The precise setting and assumptions for (1) will be described below.
This abstract framework simultaneously includes problems which have been analysed on their own as well. For instance, abstract boundary feedback systems, see [DeschLasieckaSchappacher], [DeschMilotaSchappacher], [CENN] and the references therein, fit into the above abstract framework where the equations in and
representing the bulk and boundary equations. Such examples arise, for instance, for the boundary control of partial differential equation systems, see
[LasieckaTriggianiI, LasieckaTriggianiII], and [KumpfNickel], [Engel_etal_2010, Section 3], and [Adler_etal_2017, Section 3]. These problems usually involve a bounded feedback operator acting on , which can be easily incorporated into the nonlinear term above. We further note, that semilinear parabolic equations with dynamic boundary conditions, see [Wen59, engel2005analyticity, GalG08, ColF, vazquez2011heat, Lie13, Racke2003cahn, Gal07, dynbc], etc., and diffusion processes on networkswith boundary conditions satisfying ordinary differential equations in the vertices, see
[Mugnolo, MugnoloDiss, Sikolya_flows, MugnoloRomanelli, Mugnolo_book], etc., both formally fit into this setting. In both cases, however, the feedback operator is unbounded.In this paper we propose, as a first step into this direction, a Lie splitting scheme for abstract semilinear boundary coupled systems, where the semilinear term is locally Lipschitz (and might include feedback). An important feature of our splitting method is that it separates the flows on and , i.e. separates the bulk and surface dynamics. This could prove to be a considerable computational advantage if the bulk and surface dynamics are fundamentally different (e.g. fast and slow reactions, linear–nonlinear coupling, etc.). In general, splitting methods simplify (or even make possible) the numerical treatment of complex systems. If the operator on the righthand side of the initial value problem can be written as a sum of at least two suboperators, the numerical solution is obtained from a sequence of simpler subproblems corresponding to the suboperators. We will use the Lie splitting, introduced in [BagGod57], which, from the functional analytic viewpoint, corresponds to the Lie–Trotter product formula, see [Trotter], [EngelNagel, Corollary III.5.8]. Splitting methods have been widely used in practice and analysed in the literature, see for instance the survey article [McLachlanQuispel_acta], and see also, e.g., [Strang68, JahnkeLubich, Thalhammer2008, HanO09], etc. In particular, for semilinear partial differential equations (PDEs) with dynamic boundary conditions, two bulk–surface splitting methods were proposed in [dynbc]. The numerical experiments of Section 6.3 therein illustrate that both of the proposed splitting schemes suffer from order reduction. Recently, in [dynbc_PDAE_splitting], a firstorder convergent bulk–surface Lie splitting scheme was proposed and analysed.
In the present work we start by the variation of constants formula and apply the Lie splitting to approximate the appearing linear operator semigroups. More precisely, we will identify three linear suboperators: two describing the dynamics in the bulk and on the surface, respectively, and one corresponding to the coupling. Then, either the solutions to the linear subproblems are known explicitly, or can be efficiently obtained numerically. We will show that the proposed method is firstorder convergent for boundary coupled semilinear problems. The proposed method does not suffer from order reduction, and is therefore suitable for PDEs with dynamic boundary conditions, cf. [dynbc], see the experiment in Section 5.2. However, due to the unbounded boundary feedback operator, our present results do not apply to this case directly. Nevertheless, we strongly believe that the developed techniques presented in this work provide further insight into the behaviour of operator splitting schemes of such problems. This is strengthened by our numerical experiments.
The convergence result is based on studying stability and consistency, using the procedure called Lady Windermere’s fan from [HairerWannerI, Section II.3], however, these two issues cannot be separated as in most convergence proofs, since this would lead to suboptimal error estimates. Instead, the error is rewritten using recursion formula which, using the parabolic smoothing property (see, e.g., [EngelNagel, Theorem 4.6 (c)]), leads to an induction process to ensure that the numerical solution stays within a strip around the exact solution. A particular difficulty lies in the fact that the numerical method for the linear subproblems needs to approximate a convolution term in the exact flow [CsEF], therefore the stability of these approximations cannot be merely established based on semigroup properties. Estimates from [CsEF] together with new technical results yield an abstract firstorder error estimate for semilinear problems (with a logarithmic factor in the time step), under suitable (local Lipschitztype) conditions on the nonlinearities. By this analysis within the abstract setting we gain a deep operator theoretical understanding of these methods, which are applicable for all specific models (e.g. mentioned above) fitting into the framework of (1). Numerical experiments illustrate the proved error estimates, and an experiment for dynamic boundary conditions complement our theoretical results.
The paper is organised as follows.
In Section 2 we introduce the used functional analytic framework, and derive the proposed numerical method. We also state our main result, namely, the firstorder convergence, the proof of which along with error estimates takes up Sections 3 and 4.
Section 5 presents numerical experiments illustrating and complementing our theoretical results.
2 Setting and the numerical method
We consider two Banach spaces and , sometimes referred to as the bulk and boundary space, respectively, over the complex field . The product space is endowed with the sum norm, or any other equivalent norm, rendering it a Banach space and the coordinate projections bounded. Elements in the product space will be denoted by boldface letters, e.g. for and . We first discuss a convenient framework established in [CENN] to treat linear boundary coupled problems. Then we treat the nonlinearities, derive the numerical method, and present the main result of the paper.
General framework
We will now define the abstract setting for linear boundary coupled systems, established in [CENN], i.e. for (1) with and . We will also list all our assumptions on the linear operators in (1).
The following general conditions—collected using Roman numerals—will be assumed throughout the paper:

The operator is linear.

The linear operator is surjective and bounded with respect to the graph norm of on .

The restriction of to generates a strongly continuous semigroup on .

The operator generates a strongly continuous semigroup on .

The operator matrix is closed.
We recall from [CENN, Lemma 2.2] that is invertible, and its inverse, often called the Dirichlet operator, given by
(2) 
is bounded, and that
Let us briefly recall the following example from [CsEF] (see Examples 2.7 and 2.8 therein), which is also one of the main motivating examples of [CENN]; we refer also to [GeMit11, GeMit08, BeGeMi20] for facts concerning Lipschitz domains.
Example 1 (Bounded Lipschitz domains)
Let be a bounded domain with Lipschitz boundary , and .

Consider the following operators: with domain , and the Dirichlet trace of on (see, e.g., [McL00, pp. 89–106]). Then is surjective and actually has a bounded rightinverse , which is the harmonic extension operator, i.e. for any the function solves (uniquely) the Poisson problem with inhomogeneous Dirichlet boundary condition . The operator is strictly positive and selfadjoint operator generating the Dirichletheat semigroup on .

One can also consider the Laplace–Beltrami operator on , which (with an appropriate domain) is also a strictly positive, selfadjoint operator, see [GeMit11, Theorem 2.5] or [GeMiMiMi] for details.
In summary, we see that the abstract framework of [CENN], hence of this paper, covers interesting cases of boundary coupled problems on bounded Lipschitz domains.
We now turn our attention towards the semigroup, and its generator, corresponding to the linear problem. Consider the linear operator
(3) 
For and define the convolution
(4) 
For all we also define , and using integration by parts, see [CENN], we immediately write
(5) 
We see that and are both linear operators on and bounded when is endowed with the graph norm.
The next result, recalled from [CENN], characterizes the generator property of , which in turn is in relation with the wellposedness of (1), see Section 1.1 in [MugnoloDiss].
Theorem 1 ([Cenn, Theorem 2.7])
In other words, if the conclusion of Theorem 1 holds, then the linear problem is wellposed and the solution with initial value is given by the semigroup as .
We further add to the list of general conditions (i)–(v) by further assuming:

The operators and are invertible.

The operators and generate bounded analytic semigroups.
Remark 1

By Corollary 2.8 in [CENN] the assumption in (vii) implies that is the generator of an analytic semigroup on .

The invertibility of or is merely a technical assumption which slightly simplifies the proofs and assumptions, avoiding a shifting argument.

In principle, one can drop the assumption of being the generator of an analytic semigroup. In this case minor additional assumptions on the nonlinearity are needed, and the error bound for the numerical method will look slightly differently. We will comment on this in Remark 2 below, after the proof of the main theorem.

The fact that generates a bounded analytic semigroup implies the bound , see, e.g., [EngelNagel, Theorem 4.6 (c)].
For further details on analytic semigroups we refer to the monographs [Pazy, Lunardi, EngelNagel, Haase].
The abstract semilinear problem
We now turn our attention to semilinear boundary coupled problems (1). In particular we will give our precise assumptions related to the solutions of the semilinear problem, and to the nonlinearity .
The function , , is a mild solution of the problem (1), written on as
(8) 
i.e. it satisfies the variation of constant formula:
(9) 
We further assume that the exact solution has the following properties:

The function is Lipschitz continuous on the strip
around the exact solution with constant .

The second component is Lipschitz continuous on , with constant .

For each , and .

The second component along the solution satisfies for each , and

Furthermore, is differentiable and .
The numerical method
We are now in the position to derive the numerical method. For a time step , for all , we define the numerical approximation to via the following steps.

We approximate the integral in (9) by an appropriate quadrature rule.
In what follows we describe the numerical method by using firstorder approximations in Steps 1–2, and show its firstorder convergence. We note here that the application of a correctly chosen exponential integrator could be inserted as a preliminary step, see [HochbruckOstermann]. Since it eliminates the integral’s dependence on , the quadrature rule simplifies in Step 1. This approach, however, leads to the same numerical method as Steps 1–2.
Before proceeding as proposed, for all , we rewrite formula (9) at as
(10) 
Now, according to Step 1, we approximate the integral by the left rectangle rule leading to
for any .
In Step 2, we apply the Lie splitting, which, according to [CsEF], results in the approximation of the convolution operator by an appropriate (to be specified later). Altogether, we approximate the semigroup operators by
(11) 
We remark that holds with the notations introduced in [CsEF]:
This leads to the numerical method approximating at time :
(12) 
with .
The actual form of operator depends on the underlying splitting method. Here, we will use the Lie splitting of the operator , proposed in [CsEF, Section 3]. Namely, we split up the operator with
and , , . It was shown in [CsEF, Prop. 3.2.] that the operator parts , and generate the strongly continuous semigroups
respectively, on . Then the application of the Lie splitting as leads to the formula (11) with
(13) 
Thus, the Lie splitting transfers the coupled linear problem into the sequence of simpler ones. First we solve the equation on by using the original initial condition , then we propagate the solution by , which serves as an initial condition to the homogeneous problem on . To get an approximation at , the semilinear expressions and the terms coming from the “diagonalisation” should be treated. Then the whole process needs to be cyclically performed times.
We note that the approximation can also be obtained by using an appropriate convolution quadrature, i.e. by approximating from the left (at ) and from the right (at ).
Upon plugging in the splitting approximation (13) into the convolution , and by introducing the intermediate values
the method (12) reads componentwise as
(14)  
This formulation only requires two applications of the Dirichlet operator per time step. We point out that the two terms with the Dirichlet operator can be viewed as correction terms which correct the boundary values of the bulksubflow along the splitting method.
The main result
We are now in the position to state the main result of this paper, which asserts first order (up to a logarithmic factor) error estimates for the approximations obtained by the splitting integrator (12) (with (13)) separating the bulk and surface dynamics in and .
Theorem 2
In the above setting, let be the solution of (1) subject to the conditions in Assumptions 2 and consider the approximations at time determined by the splitting method (12) (with (13)). Then there exists a and such that for any time step we have at time the error estimate
(15) 
The constant is independent of and , but depends on , on constants related to the semigroups and , as well as on the exact solution .
The proof of this result will be given in Section 4 below. In the next section we state and prove some preparatory and technical results needed for the error estimates.
Recall that the splitting method (12), written componentwise (14), decouples the bulk and surface flows, which can be extremely advantageous if the two subsystems behave in a substantially different manner. We remind that, when applied to PDEs with dynamic boundary conditions, naive splitting schemes suffer from order reduction, see [dynbc, Section 6], and a correction in [dynbc_PDAE_splitting].
We make the following remark about the logarithmic factor in the above error estimate. Inequality (15) implies that for any we have with another constant . This amounts to saying that the proposed method has convergence order arbitrarily close to , and in fact this is also what the numerical experiments show. Indeed, numerical experiments in Section 5 illustrate the firstorder error estimates of Theorem 2, including an example with dynamic boundary conditions, Section 2, without any order reductions.
3 Preparatory results
In this section we collect some general technical results which will be used later on in the convergence proof. After a short calculation, or by using the results in Section 3 of [CsEF], we obtain
(16)  
see [CsEF, equation (3.9)]. Now we are in the position to prove exponential bounds for the powers of .
Lemma 1
Proof
From the sum norm on the product space , we have
The exponential boundedness of the semigroups and , and the boundedness of directly yield
It remains to bound the term . We obtain
which completes the proof of the first statement.
If is a bounded analytic semigroup, then we improve the last estimate to
By putting the estimates together, the assertions follows.
We recall the following lemma from [CsEF].
Lemma 2 ([CsEF, Lemma 4.4])
There is a such that for every , for any , and for every we have
Using the above quadrature estimate we prove the following approximation lemma.
Lemma 3
For and we have
Proof
Lemma 4
For we have
Proof
Resorting to the Taylor expansion we have for that
which readily implies , and hence the assertion.
Lemma 5
Let be Lipschitz continuous and consider
Then for all we have
Proof
For , we have
Let and denote the projection onto the first and second coordinate in .
Lemma 6
For there is a such that for every , we have
Proof
We have
and the second asserted inequality follows at once.
On the other hand, for the first component
and we obtain
and the first inequality is also proved.
Lemma 7
For there is a such that for every , , , we have
Proof
From (16) we obtain
where denote the four terms in the order of appearance. By Lemma 4
so we need to estimate . Since has the appropriate convolution form, Lemma 5 implies
Altogether we obtain
For and we have
To estimate we recall from the proof of Lemma 1 that (for ), so that
Finally, the estimates for together yield the assertion.
4 Proof of Theorem 2
The proof of our main result is based on a recursive expression for the global error, which involves the local error and some nonlinear error terms. The recursive formula is obtained using a procedure which is sometimes called Lady Windermere’s fan [HairerWannerI, Section II.3]; our approach is inspired by [OPiazzolaWalach], [diss_Walach, Chapter 3]. The local errors are weighted by , therefore a careful accumulation estimate—heavily relying on the parabolic smoothing property—is required. In order to estimate the locally Lipschitz nonlinear terms we have to ensure that the numerical solution remains in the strip (see Assumptions 2). This will be shown using an induction process, which is outlined as follows:

We shall find and a constant such that for any if belong to the strip and , then

Since is a constant independent of and , we can take sufficiently small such that for each we have , the width of the strip , therefore by the previous step we have