1 Introduction
The aim of this paper is to adapt the Magnus integrator for nonautonomous homogeneous problems to a wide class of nonautonomous delay problems. We are interested in problems of the form
(QDE) 
where is a Banach space, denotes the history of the solution at time (i.e., the initial condition could also be written as ), where is an unbounded operator on and is bounded for all , and actually only depends on the restriction to for some fixed . The exact assumptions on the operators and functions involved will be detailed later.
Note that the case corresponds to a quasilinear delay equation with bounded operators, while leads to the unbounded case.
In [15], Magnus set out to solve the nonautonomous homogeneous problem (), , where and are bounded linear operators on appropriate spaces. In the case when all of the operators commute, the solution takes the simple form . In the general, noncommutative case, however, one has to correct the exponent, and the exact solution is given by the Magnus series expansion, involving integrals of commutators. The first term of this expansion corresponds to the commutative solution, and is a good approximant leading to a secondorder numerical method using the midpoint rule to approximate the integral in the exponent:
(1) 
where is an arbitrary timestep and denotes the numerical approximation to with initial value .
Convergence of the classical Magnus integrators, such as (1), has been widely studied in the literature for nonautonomous problems without delay. For finite dimensional spaces, Blanes et al. analysed the Magnus expansion, being the basis of the Magnus integrators, in [3] and [4]. Moan and Niesen gave a condition on the Magnus expansion’s convergence in [16]. Casas and Iserles investigated the expansion in case of nonlinear equations in [5]. Csomós derived a Magnustype integrator (see (3) below) for delay equations and showed its secondorder convergence and positivity preserving property in [6].
For infinite dimensional undelayed problems with inhomogeneity, González et al. derived a numerical method by using the Magnus expansion, and showed its secondorder convergence in [9] for sufficiently smooth solutions, in addition to (lower order) error bounds measured in the domain. Bátkai and Sikolya proved the firstorder (operator) norm convergence of the same method under very mild conditions for homogeneous problems in [2].
Delay equations describe processes where the time evolution of the unknown function not only depends on the actual state of the system but also on its past values. Delay differential equations find more and more frequent application in the modelling of scientific, financial, or even social phenomena, since there a delay term often naturally appears among the processes. One may here think of the role of the latent period when modelling the spread of epidemic diseases, the pregnancy period in population models, or the reaction time in any social or financial model. In our problem (QDE), the time parameters delimit the time window into the past that is influencing the present dynamics in the problem.
Since the exact solution of problem (QDE) is difficult or even impossible to compute analytically, one needs to find a way to approximate it. To this end, we first reformulate the delay problem as a nonautonomous problem, and then present a novel Magnustype integrator based on Magnus integrators introduced for nonautonomous problems.
For the reformulation as a nonautonomous problem, we first define the function as
and the operators
Then the function satisfies the problem
(2) 
Despite the operators depending on the unknown function via , the minimum delay within the term stemming from allows problem (2) to be treated as a nonautonomous Cauchy problem a posteriori. Indeed, for , is actually determined as it only depends on the known restriction . Hence solving the problem iteratively on the subintervals , , yields an explicit nonautonomous Cauchy problem for each timesegment considered, provided wellposedness is maintained along the way. Hence, we will consider problem (2) essentially as a formally nonautonomous problem on the whole time interval. The solution can then be approximated by the method derived by Magnus in [15] for such problems. This method however uses the values , which have to be themselves approximated, involving an appropriate secondorder approximation of based off of a grid that is compatible with the one used for .
To this end, one defines an arbitrary , takes the time step , and denotes the approximate value of by for all . The Magnustype integrator, which we will derive in detail in Section 3, takes then the form
(3)  
for , where possible negative indices refer to the corresponding values of the history function, i.e., for , the halfindexed terms are auxiliary values (essentially approximating for an application of the midpoint rule), and the expressions stem from the approximation of .
Under appropriate smoothness and uniform exponential bounds on the generators involved, we will prove the secondorder convergence of the Magnustype integrator (3) when applied to the abstract delay equation (QDE) on Banach spaces. Moreover, we show that the method inherits the invariance properties of the generators (e.g. positivity).
The paper is organised as follows. In Section 2 we introduce the abstract setting of evolution equations like (2). In Section 3 we recall the original Magnus integrators, summarise the main results from the literature regarding its application to nonautonomous Cauchy problems (2), and then describe how we adapt this method to our case where we have no a priori knowledge of the operators on the time interval .
Section 4 contains our main result, Theorem 23 on the Magnustype integrator’s secondorder convergence when applied to the quasilinear delay equation (QDE). Some of the assumptions needed for this convergence (or even the existence of the solution to the delay equation) – especially a uniform exponential bound for the semigroups generated by the operators – are typically not naturally achievable for all . However, many problems admit invariants and have qualitative preservation features (e.g., positivity of the solutions), and we shall show that our Magnustype integrator naturally exploits such invariants, allowing us to restrict our assumptions to some smaller, invariant closed subset of that we assume our initial history runs within. This will allow us to prove in parallel that the Magnustype integrator and the exact solution both exist for a positive amount of time (namely, ), with the former never leaving the invariant set. That in term will imply the secondorder convergence on this time interval, implying that the solution itself stays in the invariant set as well. Iterating these arguments, we obtain our results for any compact time interval.
In Section 5 we use a spacedependent epidemic model to illustrate the power of invariants (in this case both total population size and positivity) in ensuring secondorder convergence.
2 Autonomous and nonautonomous evolution equations
In this section we introduce the notions necessary to understand the Magnus integrator and our approach to the error bounds. Throughout we shall assume that is a Banach space. Our main references are Engel and Nagel [8] and Nickel [17].
Definition 1.
A family of bounded linear operators on is said to be a strongly continuous semigroup generated by the linear, closed, and densely defined operator if the following holds:

, the identity operator in ,

for all ,

the function is continuous for all ,

there exists for all .
This strongly continuous semigroup describes the solution to the Abstract Cauchy Problem
in the sense . It is known (see e.g. [8, Prop. I.5.5]) that such semigroups are exponentially bounded, i.e., there exist and such that for all .
In many processes, however, one cannot assume the generator in to be constant in time, leading to socalled nonautonomous Cauchy problems, or for short. Let be a linear operator on for every . Furthermore, let and be given. Then we consider the following nonautonomous problem for the differentiable unknown function :
The following definitions are based on [8, Section VI.9].
Definition 2.
A continuous function is called a (classical) solution of if , for all , , and for all .
Definition 3.
For a family of linear operators on the Banach space , the nonautonomous Cauchy problem is called wellposed with regularity subspaces if the following holds.

Existence: For all the subspace
is dense in .

Uniqueness: For every the solution is unique.

Continuous dependence: The solution depends continuously on and , i. e., if , with then we have uniformly for in compact subsets of , where
If, in addition, there exist constants and such that
for all and , then is called wellposed with exponentially bounded solutions.
Definition 4.
A family of linear, bounded operators on a Banach space is called an (exponentially bounded) evolution family if

for all ,

the map is strongly continuous,

for some , and all .
An evolution family is said to be contractive if we can choose and , and quasicontractive if we can choose for some .
Definition 5.
An evolution family is called an evolution family solving , if for every the regularity subspace
is dense in .
We have the following result connecting the wellposedness of to the existence of a unique evolution family solving it.
Theorem 6 ([17, Prop. 2.5]).
Let be a Banach space, a family of linear operators on . Then the nonautonomous Cauchy problem is wellposed if and only if there exists a unique evolution family solving .
As mentioned in the Introduction, our goal is to rephrase (QDE) as a nonautonomous Cauchy problem on the time interval . Of great help in establishing wellposedness is the following consequence of a result by Kato ([12, Thm. 4]).
Theorem 7.
Let be a closed interval, and a family of generators of contraction semigroups on the Banach space with common domain satisfying for all . Then the nonautonomous Cauchy problem is wellposed with regularity subspaces () and admits a contractive evolution family . In particular, for any and the initial condition leads to a unique (classical) solution with for all .
Remark 8.
The contractivity of the evolution family is actually part of the earlier Theorem 2 in the same paper.
Also, note that Kato’s result talks of evolution families and wellposedness on a possibly bounded closed interval instead of , and it is not immediately clear how the two can be connected. Let with . Since each () is a generator, it makes sense to extend to by setting when and when . Then the evolution family corresponding to the extended interval can obviously be restricted to with all the properties preserved. But also the evolution family corresponding to the restricted problem has a natural extension to as follows.
Denote by the contraction semigroup generated by , and by the contraction semigroup generated by .
If , then let and if , then let . This defines the evolution family for all in a way compatible with Definition 4 ( may have to be replaced with ).
An easy rescaling argument yields that the same holds if instead of generators of contraction semigroups we consider a family such that the generated semigroups are uniformly quasicontractive.
Corollary 9.
Let be a closed interval, and a family of generators of uniformly quasicontractive semigroups on the Banach space with common domain satisfying for all . Then the nonautonomous Cauchy problem is wellposed with regularity subspaces () and admits a quasicontractive evolution family . In particular, for any and the initial condition leads to a unique (classical) solution with for all .
3 Magnustype integrator
We saw in the Introduction that the delay equation (QDE) could formally be written as the nonautonomous abstract Cauchy problem with the solutiondependent operator . Whilst wellposed autonomous abstract Cauchy problems have their solutions given through of a oneparameter strongly continuous semigroup, problem – if wellposed and is in the regularity subspace – has its solution given through a unique twoparameter evolution family (cf. Theorem 6). Since the exact form of the evolution family is usually unknown, we need to approximate the solution in (4). To this end we define an arbitrary and take the time step . Then the approximate value of at time levels , , is denoted by .
In case of a finite dimensional space , the evolution family is the exponential of a bounded operator for . Magnus showed in [15] that the bounded operator could be expressed by the integral of an infinite series, called Magnus series. Casas and Iserles showed in [5] that the appropriate truncations of the series led to convergent approximations. The further approximation of the integral terms yields the Magnustype integrators. More precisely, by the property of the evolution family we have
for all . By approximating by as in [5], and then by with the midpoint quadrature rule, we arrive at the formula of the simplest Magnus integrator
(5) 
for all with .
In case of an infinite dimensional Banach space , we consider formally the same formula (5) where the exponential refers to the strongly continuous semigroup generated by the corresponding operator (cf. Definition 1).
Since the Magnus integrator (5) gives only an approximation to the exact solution at time for all , it is necessary to show that the approximate value converges to the exact value as the time step tends to zero, or equivalently, tends to infinity. As is usual, we will want to show that our numerical scheme yields a good approximation on any compact time interval .
Definition 10.
Let denote the exact solution to problem (QDE). Its approximation (or, equivalently, the corresponding numerical method) is called convergent of order , if there exists such that holds for all and with , where the constant is independent of and but may depend on .
We note that we defined the convergence for the sequence of the time steps in order to simplify our proofs, which could be done for all sequences and with by introducing a more complicated formalism. Since one has an initial history function (or a set of data) on the time interval , it is natural to choose a time step that is compatible with it.
In what follows we present two results from the literature about the Magnus integrator (5) for nonautonomous problems, since we will use them in our analysis.
Theorem 11 ([9, Thm. 2]).
Let and be Banach spaces with densely embedded in . We suppose that the closed linear operator is uniformly sectorial for . Moreover, we assume that the graph norm of and the norm in are equivalent. We also assume that , and in particular there then exists a constant such that
holds for all . Moreover, we introduce the notations
(6)  
and corresponding notations will also be used with instead of . Then the Magnus integrator (5) applied to problem is convergent of second order, that is, there exists a constant , being independent of and
, such that the following estimate holds for the global error:
(7) 
for all , provided that the quantities on the righthand side are welldefined.
The following theorem is essentially the quasicontractive, continuously differentiable special case of the consistency result from the proof of [2, Thm. 3.2], more specifically inequality (5) therein. It hinges on the fact that Corollary 9 implies that the wellposedness, stability and local Hölder continuity conditions of that theorem are automatically satisfied.
Theorem 12.
We consider the problem on the Banach space and suppose the following.

There exists a such that for all and . Further, , where is the generator of a strongly continuous semigroup, and for all .

The map is continuously differentiable as a map from to for some .
Then for all with we have the following estimate:
(8) 
where is the Lipschitz constant of on .
To illustrate the process of obtaining our Magnustype integrator, we shall first present a special case that already exhibits some of the new ideas involved, and then we indicate what further changes are needed to accommodate for the general case.
Example 13.
In this example we shall focus on the special case , i.e., when the delay is concentrated on a specific point of the history function. Our equation (QDE) then takes the form
(QDE’) 
In this particular case, we have in (2). Recall that for , however, its value is unknown for . Thus, we need to approximate in the Magnus integrator (5) to obtain a working method. Now the natural idea would be to use the appropriate , however falls exactly between two points of our time grid, and has to be obtained via further approximation. We therefore introduce the corresponding term as an auxiliary value, obtained via another Magnus step with half time step. In full detail, we have the following.
Definition.
Let be an arbitrary integer and . Then the Magnustype integrator for pointdelay, which yields an approximation to the solution of the pointdelay equation (QDE’) at time levels is given as follows. We introduce the notation for . Then for , we have the recursion
(9)  
Note that since they are essentially auxiliary values, the way we indexed the terms is not indicative of the time layer they correspond to (which would actually be ). Rather, the indices reflect the natural order in which one would execute the algorithm, i.e.,
despite being able to compute already at earlier stages.
In the general case of the quasilinear delay evolution equation (QDE), we have in (2). If we want to apply the Magnus integrator (5), we need to be able to approximate . Even for , all we have available is an approximation to some of its values, corresponding to the time levels in our grid. Thus we first need a discretisation of itself using only discrete values of , and then use the approximation of those values in the final form of our method. The simplest approach is to make the discretisation of compatible with the original time grid, and use the same auxiliary Magnus step as in the previous Example 13.
To this end we introduce the approximation of in the following form
(10) 
for appropriate elements , weights , and functions having properties to be detailed in Section 4. We rewrite now (QDE) as a nonautonomous problem (2) with , apply the Magnus integrator (5), and approximate as in (10) to obtain for all :
(11)  
where we used the definition of the history function in the last step. We approximate next the intermediate values using another Magnus step with time step but now by taking the leftrectangle rule when approximating the integral in the exponent, cf. [5]:
(12)  
for all with . Observe that formula (12) can be reindexed by using instead of . However, since is not defined for negative times, whenever , the values at the corresponding intermediate time levels should be obtained from the initial history function instead. By combining (11) and (12), we thus obtain the following definition.
Definition 14.
Let be an arbitrary integer and . Then the Magnustype integrator, which yields an approximation to the solution of the delay equation (QDE) at time levels is given as follows. As before, we use the notation for . Then for , we have the recursion
(13)  
where the exact properties of the weights and of the functions will be presented in Section 4.
Remark 15.
To fit Example 13 in this general formulation, set and for , with the identity for all .
Example 16.
Let us consider a delay term where a fixed delay time interval uniformly governs the dynamics. More specifically, we shall take a closer look at the delay when
i.e, we have
in (QDE).
The first thing to note is that whenever is divisible by , the endpoint
of the integral also falls on the discretisation grid, significantly simplifying things, and the odd
’s have to be treated slightly differently to fit the general framework. Alternatively, we could simply adapt the general scheme to this special situation by only considering even values for , but we shall detail the odd case nevertheless.In contrast to the pointinteraction delay in Example 13 where we could directly substitute the computed values into the delay term, we here have an integral that itself has to be numerically approximated using some appropriate quadrature. Taking into consideration the order of the error that we need to achieve for the recursive inequalities to result in a secondorder Magnustype integrator, the error of quadrature has to be of the magnitude of .
So for even , we use the composite trapezoidal rule with nodes for as
The Magnustype integrator thus has the form for :
It has the form (13) with weights and for where each is the identity. We remark that the weights sum up to (cf. (14) later on).
For odd values of , the point is not part of the grid, so we have to use the truncated approximation
instead, again with error of order , i.e., and for .
It is reasonable to wonder what happens if the delays above are not given in the convenient form where we are looking at a convex combination/normalised integral, for instance, if we were dealing with a delay of the form . Actually, this is not really an issue, as we may then define and thereby revert to the convex combination case detailed above. This is only a question of convenience, related to the interpretation of the invariance of the set . After rescaling should stay invariant under the semigroups generated by the for all , whereas without rescaling the invariance is under the ones generated by with .
4 Convergence
This section contains our main result regarding the secondorder convergence of Magnustype integrator (13) applied to the quasilinear delay equation (QDE). We have already seen how problem (QDE) can formally be written as the nonautonomous problem . Hence, Definition 10 of the convergent approximation remains valid also for the Magnustype integrator (13) applied to the delay problem (QDE).
From now on we use the notations and for the left and right derivatives of a function at zero, respectively.
We will need the following list of assumptions.
Assumptions.
Let be a Banach space, a closed subset and a dense subspace.

We have for all , where , generates a contraction semigroup on , , and the operators () are all dissipative. We shall also assume , which is no real added restriction as we may simply replace by and by for some .

The function is continuously differentiable on the set .

The function is twice continuously differentiable on the set .

For any the operator leaves invariant, and is bounded and dissipative on , where . In addition, is continuously differentiable on .

The operator is a sectorial operator generating an (analytic) contraction semigroup with sector for some , and the operators () are all such that for any the operator is dissipative.

with whenever .

with whenever .

The functions and weights () are such that there exists a constant independent of such that
(14) and the function
(15) satisfies
(16) for some constant independent of . In addition, there exists a constant such that and are Lipschitz with constant for all and .

The set satisfies
and
(17) for all .

The initial history function is in and satisfies for all and the boundary conditions
(18) (19)
Remark 17.
The condition (17) is automatically satisfied when is convex, and and for all and .
The following results will show that under appropriate smoothness assumptions on and the initial history function , the solution itself will exhibit similar smoothness properties on and . Also, we shall show that Theorem 11 is applicable to our setting.
The next two results show that our Assumptions imply that the semigroups generated by the operators () are uniformly quasicontractive, and uniformly sectorial as well.
Lemma 18.
Under Assumption (a) each operator () with domain is the generator of a contraction semigroup on . In particular, the semigroups generated by are uniformly quasicontractive, that is, for all and .
Proof.
The operators are by assumption dissipative, and with bound 0 (as they are bounded operators). Thus the claims follow directly from [8, Thm. III.2.7] and the usual rescaling argument. ∎
Lemma 19.
Proof.
The assumptions imply that for any the operator generates a contraction semigroup and is dissipative with bound 0, hence each generates a contraction semigroup as well by [8, Thm. III.2.7]. To see that the semigroup generated by is analytic on , we use the semigroup property and the fact that [8, Thm. III.2.14] implies that generates an analytic semigroup on some sector