Singularly Perturbed Population Models with Reducible Migration Matrix: 2. Asymptotic Analysis and Numerical Simulations

The paper is concerned with asymptotic analysis of a singularly perturbed system of McKendrick equations of population with age and geographical structure. It is assumed that the migration between geographical patches occurs on a much faster time scale than the demographic processes and is described by a reducible Kolmogorov matrix. We apply a novel regularizing technique which makes the error estimates easier than that in previous papers and provide a numerical illustration of theoretical results.


Introduction
The present paper is an extension of [3][4][5][6] to reducible migration matrices and it can be treated as a (self-consistent) companion paper to [8]. All three papers are concerned with asymptotic analysis of McKendrick type equations with age and spatial structure, modelling fast migrations between the spatial patches and showing that for large migration rates the solutions of the original problem can be approximated by solutions to a much simpler aggregated model. In [8] and here we consider a more general mechanism of migrations than in [3][4][5][6], which is described in detail in Sect. 2

. However, in
The results of this paper have been obtained with partial support of the Research Funds of the Universities of KwaZulu-Natal and Zululand and the Grant N N201 605640 from the National Scientific Centre of Poland. 534 J. Banasiak et al. MJOM [8] we focused on the functional analytic aspects of analysis by showing how the model fits into the general framework of Trotter-Kato and Sova-Kurtz theorems, see [9], and how the practical implementation of these theorems can be facilitated by a careful use of classical asymptotic expansions. However, we stopped short of performing the full analysis of the layers which is necessary for the convergence to occur uniformly in time and for all initial conditions. In this paper we complete the analysis by considering all necessary layers; that is, the initial, the boundary and the corner layers provide biological interpretation of the obtained results and illustrate them on a worked example. Using this example we also numerically demonstrate the accuracy of the approximations obtained by the asymptotic analysis. We emphasize that the asymptotic analysis used in this paper yields a more detailed error analysis than that yielded by [8] but at a much higher computational cost.
To describe the results in more detail, we have to recall the model.

The Model and Some Heuristic Considerations
In recent years there has been an interest in structured population models describing populations with various levels of organization such as individual, population, community or ecosystem levels. In many cases we observe different time scales at which each level evolves. For instance, considering demography of a population together with migratory processes between various geographical patches its individuals occupy, we often find that the latter, consisting, e.g. of moving between feeding and resting locations occurs at a much faster pace than the former (which includes e.g. breeding and deaths).
To account for this, in a series of papers, [1,6,8,16], the authors consider the following version of the McKendrick system: Here, as in [6,8], u (a, t) = (u 1, (a, t), . . . , u n, (a, t)), where u i, (a, t) is the population density at time t of individuals residing in patch i and being of age a. Further, the symbol ∂ · denotes the differentiation with respect to the variable appearing as the subscript and, accordingly, Su := −∂ a u is the operator describing the ageing, M(a) = {μ ij (a)} 1≤i,j≤n is the matrix describing mortality as well as possible slow migrations. The matrix C(a) = {c ij (a)} 1≤i,j≤n is a Kolmogorov matrix (whose columns sum to 0) which describes the migration of individuals between patches and 1/ , with small > 0, accounts for the fact that the rate of migration is much higher than the demographic rates. In the boundary condition, γ denotes the operator of taking the trace at a = 0 and B(a) = {β ij (a)} 1≤i,j≤n is the maternity matrix. The problem (1)-(3) is considered as an abstract Cauchy problem in the space X n := L 1 (R + , R n ) = (L 1 (R + , R)) n . Accordingly, the domain of S is the subspace of the Sobolev space Y n := (W 1 1 (R + )) n .
Vol. 11 (2014) Singularly Perturbed Population Models 535 If the rate of migration is much faster than the demographical processes, and the number of patches large, then it is of interest to provide an 'aggregated' description of the system which 'averages' over the geographical structure. Assuming for simplicity that M and B are diagonal matrices, we can use the Kolmogorov structure of C to remove from (1) by adding the equations to get where u = n i=1 u i is the density of the total population across all patches. However, Eq. (4) is not closed unless μ 1 = · · · = μ n . To approximately close it, in applications typically it is assumed that the ratios u i /u stabilize in time and thus can be replaced by the components of a constant vector k = (k 1 , . . . , k n ) with k i ≈ u i /u for i = 1, . . . , n, called the stable patch distribution. This leads to the approximating aggregated equation where μ * = n i=1 μ i k i is the aggregated mortality. The aggregated maternity rate β * is defined in the analogous way.
The results of [1,6,16] confirm this biological heuristics but only if the migration matrix C is irreducible, see, e.g. [20]. This ensures that λ = 0 is the dominant simple eigenvalue of C with corresponding positive right eigenvector k, and a left eigenvector 1 = (1, 1, . . . , 1). These vectors are normalized to satisfy 1·k = 1 so that k indeed gives the stable patch distribution mentioned above. Then, introducing the spectral decomposition of R n corresponding to λ = 0, by R n = Span{k, w}, we can write u (t) = u (t)k + w (t) and We denote byū the solution of the full aggregated model corresponding to (5), Actually, the analysis can be done for general M and B. Then, see [6], uniformly in t on finite intervals of [0, ∞). Since λ = 0 is the strictly dominant eigenvalue, we recognize e t C(·) • w (·) as the initial layer contribution to the approximation, whose influence rapidly decreases for t > 0 and small .
It is to be emphasized that the above result only is valid for irreducible migration matrices. This assumption often proves to be too restrictive. For instance, in population theory it amounts to requiring that any geographical patch is accessible from any other, see e.g. [10]. In many cases, however, one has to consider patches that are either isolated or only admitting emigration (or only immigration). A typical example is a population with a breeding refuge which typically is inaccessible to males living in other patches but MJOM supplies them with offspring, or a pool above a waterfall which forms a barrier for fish living downstream. Such models result in systems of the form (1), but with a reducible matrix C. In such cases the total population u cannot be approximated by the solution of a scalar equation, such as Eq. (6), but instead one must use a system of equations for the population density in each of the so-called ergodic components of the system. The appropriate aggregated system, that replaces (6) if C is reducible, was constructed in [8,12], where we also used the Sova-Kurtz theory to prove the convergence of u to the solutions of this system as → 0, but only for the initial conditions from the ergodic states. The main aim of this paper was to provide a detailed description of the decomposition of the state space into ergodic and the complementary, transient, states and, by a careful analysis of the appropriate layers, to extend the convergence results to the solutions emanating from arbitrary initial conditions. The paper is concluded by a worked example and numerical simulations illustrating the accuracy of the constructed approximations.

The Migration Matrix
As we know from [6,8], the key to asymptotic analysis of (1)-(3) is understanding the large time behaviour of the solutions to the problem where u = (u 1 , . . . , u n ) ∈ R n and u 0 is an initial condition, which describes the so-called fast dynamics of (1). Here the age variable a is a parameter and thus it will be suppressed in this section. The general theory explicitly describing behaviour of systems with reducible transition matrices seems to be hard to achieve due to a complicated internal structure they may have, see, e.g. [7]. However, in the structured population models considered here, we only consider matrices C that describe migration and, as mentioned before, the principle of conservation of the total population yields that they are Kolmogorov matrices. A Kolmogorov matrix is a particular case of an ML matrix, see [20, Section 2.3], or a Metzler matrix, [17, Section 6.5], in which columns sum up to zero. In other words, a Kolmogorov matrix is the (transpose) of a Q matrix describing a continuous time Markov process with finitely many states. The theory of the large time behaviour of such processes, even in the reducible case, is relatively well developed in the probabilistic context, especially in the discrete time, see, e.g. [18,Section 8.4] or [2,Chapter 5]. Here we shall look at it from the dynamical systems point of view which seems to be less available.
Before we discuss the relevant properties of Kolmogorov matrices, let us introduce basic notation. The spectrum of C will be denoted by σ(C). The spectral bound of C is defined as s(C) := max{ λ; λ ∈ σ(C)}. It is easy to see that 1 := (1, . . . , 1) ∈ R n is a left eigenvector of any Kolmogorov matrix C, corresponding to the eigenvalue 0 which is dominant, that is, all other eigenvalues λ ∈ σ(C) satisfy λ < 0, see, e.g. [17, Section 6.5, Theorem 1].
Vol. 11 (2014) Singularly Perturbed Population Models 537 Let us recall, [11, p. 90] (note, however, that due to the particular application we have in mind, we use the notation of [11] in the transposed form), that a reducible Kolmogorov matrix C can be written, by simultaneous permutation of rows and columns, in the following normal form: where C l = [c ij ] n l ≤i,j≤N l for l ≤ m and T l = [c ij ] n l ≤i,j≤N l for m + 1 ≤ l ≤ r, where n 1 = 1, n l+1 = N l +1 and N r = n, are square irreducible matrices (with some possibly being the scalar 0). Furthermore, if r > m, then in each column of (9) at least one of the off-diagonal matrices is not equal to O. Clearly, the matrices The structure of (9) emphasizes the division of the space R n of states into two sets, namely the closed, or ergodic, states from which no transition outside is possible (1 ≤ j ≤ N m ) and transient states from which a transition to other states is possible (n m+1 ≤ j ≤ N r ).
Using the standard equivalence of positive off-diagonal matrices and nonnegative matrices, see, e.g.
. . , N l and x 1 + · · · + x m = 1. 6. There is a strictly positive right eigenvector of C corresponding to λ = 0 if and only if r = m, that is, when there are only ergodic states.
We call a left (resp. right) eigenvector corresponding to the dominant eigenvalue λ = 0, a left (resp. right) Perron eigenvector. MJOM Remark 2.1. Statements 2 and 3 of the theorem are equivalent to the existence of a strictly positive left Perron eigenvector which, in the present case, is 1. Considering, however, only Kolmogorov matrices does not compromise generality of the considerations since any matrix A with strictly positive left Perron eigenvector and s(A) = 0 is similar to a Kolmogorov matrix. Indeed, let p = (p 1 , . . . , p n ) be a strictly positive left Perron eigenvector. Denoting A discussion of Theorem 2.1 and proofs of its less known statements is given in [8]. Here we recall the formulae which are necessary for the asymptotic analysis carried out in this paper and provide an interpretation of them in the context of population theory.
Let us denote by ν i the dimension of the matrices where hereafter we adopt the convention that any summation over an empty set of indices is 0. Only m first entries in (10) are nonzero (even strictly positive) and the eigenvectors e i , i = 1, . . . , m, of Theorem 2.1, point 4., can be obtained by normalizing them and extending by 0 to n dimensional vectors. For a left Perron eigenvector y = (y 1 , . . . , y r ) of C [notation as in (10)] we clearly get This formula recursively defines y i , m+1 ≤ i ≤ r, since the first term is given explicitly by The explicit formula for y l , l ≥ m + 1, is related to the characterization of absorbing states in reducible Markov chains (see, e.g. [18,Section 8.4] in the discrete case). As the continuous case seems to be less known, we provide an Vol. 11 (2014) Singularly Perturbed Population Models 539 elementary algebraic discussion of it below. For this, we have to introduce some notation. For given m + 1 ≤ j, l ≤ r, V j,l will denote the set of all sequences of indices s := {j, s 1 , . . . , s k , l} with j < s 1 < · · · < s k < l. In other words, V j,l is the set of all potential pathes from j to l. Let n s denote the length of s. Then, identifying s 0 = j, s ns = l, we write s = {s i } 0≤i≤ns .
Here, we adopt convention that any summation over an empty set of indices is 0 and the product over such set of indices is the identity matrix. Further, in the products of matrices the indices increase from left to right.
Proof. To shorten notation we denote To prove (13), we use induction. The formula is true for l = m + 1. Indeed, then V m+1,m+1 = {m + 1}, the set over which we take the product is empty and the formula reduces to (12). Assuming now that (13) is valid for all i = m + 1, . . . , l, from (11) we have This ends the proof of (13).
Defining y we obtain the basis of the left null space of C as the span of . .
where ) and y We can strengthen, [8], the statement 5 of Theorem 2.1 to Furthermore, since T l , m + 1 ≤ l ≤ r, are positive off diagonal matrices with s(T l ) < 0, they are resolvent positive so that Using the constructed eigenvectors, we obtain the spectral decomposition of the state space as with the spectral projections onto the eigenvectors defined by Thus, the projection of u onto the null space V of C is given by We also will use the m-dimensional vector of coordinates The complementary space W defined by with the complementary projection However, in contrast to [8], there is no need to construct the full spectral basis in R n as we use another approach. Hence, we leave elements of W in coordinate free form.
Remark 2.2. In the considered meta-population it is clear that the states j = n m+1 , . . . , n, eventually will become depleted with the population moving towards states j = 1, . . . , N m . However, according to (17), the population in any final state will contain a contribution from the transient states that communicate with that state (but not from any other ergodic state). It is interesting to quantify these contributions. First, let v l be a solution of Using the fact that s( Vol. 11 (2014) Singularly Perturbed Population Models 541 Hence, (−T l ) −1 v l,0 is the vector of all individuals that originated from the initial state v l,0 and existed at any time at any state n l ≤ j ≤ N l . This interpretation is clear if T l = t l is one-dimensional. Then 1/t l is the average lifespan of an individual and thus v l,0 /t l is the average number of individuals from the initial cohort v l,0 that will be in the system over the whole time.
Recall that Now, since T j l is the matrix covering all possible transitions from the block of states l to the block j with m + 1 ≤ j, l ≤ r, and A i j is the matrix of rates, with respect to the chosen unit of time, of transfers from states m+1 ≤ j ≤ r to the states

Model with General M(a), B(a) and Reducible C(a)
In the previous section we suppressed the age variable a in the notation since in the analysis of (8) it has played the role of a parameter. Hereafter we return to the analysis of the full problem (1)-(3). Let, for j ∈ N, X j := L 1 (R + , R j ) and Y j := W 1 1 (R + , R j ). We recall that X n is our state space. We denote by | · | j the norm in the l 1 j (R j equipped with l 1 norm). By · j the denote the norm of a function u : X → l 1 j , We drop the index j in the notation of the space and the respective norm if the dimension of the space is clear from the context. For any other norm in a Banach space Z we shall write · Z . By B(Z 1 , Z 2 ) we denote the space of linear bounded operators from Z 1 to Z 2 , equipped with the standard operator norm ||| · ||| B(Z1,Z2) (with short notation ||| · ||| = ||| · ||| j = ||| · ||| B(Xj ,Xj ) , if no ambiguity arises). We consider the problem (1)-(3) with general a dependent matrices M(a) and B(a) ≥ 0 discussed in Introduction. We assume that their entries are bounded. The matrix −M(a) is assumed to be sub-Kolmogorov; that is, −M(a) is positive off-diagonal and satisfies − n i=1 μ ij (a) ≤ 0 for any 1 ≤ j ≤ n and a.e. a ∈ R + . Furthermore, we assume that a → M(a) ∈ C 1 b (R + , R n 2 ) and a → C(a) ∈ C 2 b (R + , R n 2 ) (bounded functions with bounded derivatives), that the structure (9) In what follows, slightly abusing the notation, we denote the operators in X, generated by pointwise multiplication of elements of X by matrices M(a), C(a), P(a), etc; that is, e.g. X f (·) → C(·)f (·) ∈ X, with the same letters. Due to the adopted assumptions, these operators are bounded in X. We further denote Bf = ∞ 0 B(a)f (a) da and, finally, A is the operator given by the differential expression where ω is a constant depending on B and M but independent of , see [8,15,21].
For the error estimates one needs to deal with inhomogeneous boundary conditions in problems of the form where f is a function differentiable with respect to t. If f = 0, the problem cannot be approached by only using semigroup theory but some version of the lifting technique must be used. By linearity, we can consider the case with where One of the problems with classical asymptotic approach to (1)-(3) is that the solution to it as well as the solution to the aggregated system (6) are discontinuous unless the initial value satisfies compatibility conditions described in the definition of the generator of the semigroup. Thus, in general, using the differential form of the equation is in both cases incorrect. One can deal with this difficulty as in [6] and write the problem in the integral form. Such an approach gives optimal error estimates in , but it requires cumbersome calculations and depends on the fact that the differentiation operator S is an n-fold copy of −∂ a , which restricts its applicability. Here we suggest an alternative approach which allows for working with the original differential formulation and thus also can be applied to systems with different growth rates, such as occur in size structured populations. The method depends on the following lemma: where C is a constant independent of δ and ψ.

The Asymptotic Expansion and Its Convergence
Decomposition (16) was constructed for a fixed a ∈ R + but it induces a corresponding decomposition of the space X into V and W, where V is the subspace V = N (C) = {u ∈ X; C(a)u(a) = 0, a ∈ R + } and W is defined analogously using W . Then, according to the convention, the operators P j , defined by (17), are spectral projections onto V and thus P j Cu = 0 for 1 ≤ j ≤ m.
The following calculations are similar to those carried out in [8, Section 5]. The only difference is that in the approach of this paper there is no need to write the complementary, kinetic, part w of the solution in the spectral coordinates. Thus, we will skip details of some calculations.
We have, for 1 ≤ j ≤ m, where u k is defined (17) and μ * jk := x j ·Me k . We define M * := {μ * jk } 1≤j,k≤m . It follows, [8] Writing the solution as u = v + w = m j=1 u j e j + w, we obtain Properties of e j yield m j=1 u j ∂ a e j ∈ W and, since where w = (w 1 , . . . , w n ); that is, w i , i = 1, . . . , n are coefficients of the expansion of w in the canonical basis in R n . Hence Now, applying projections P j , 1 ≤ j ≤ m, and Q to (1) we obtain the projected system in a compact form as where E v and υ = (u 1 , . . . , u m ) were defined in (19), F, G, H, J are linear bounded operators on respective subspaces of X, related to the projections of the operator M and derivatives of the spectral projections, whose particular form does not have any bearing on the calculations, is a solution to (1)-(3) and, therefore, is a solution to (29). Note that the above formula can be used to define projections (υ (t), w (t)) of T A (t) • u with arbitrary • u ∈ X; that is, for mild solutions to (1)-(3), but such (υ (t), w (t)) no longer are (classical) solutions to (29).

Deriving the Aggregated System
We obtain the aggregated system corresponding to (6) by performing an asymptotic expansion with respect to the small parameter of the solution and approximately closing the obtained hierarchy at the zeroth level. Since our aim is to obtain an approximate equation for v (or υ ), the best suited for this is the Chapman-Enskog procedure in the version introduced in [19]. According to this approach, we leave v 1, , . . . , v m, unexpanded so that where v k, = u k e k , k = 1, . . . , m. Recall also that v = E v υ , see (19). Since the formal procedure is analogous to that in [6,8], we skip the details. Inserting (31) into (29) and retaining only terms of zeroth and first order in , we obtain The invertibility of C W on W yieldsw 0 = 0 and this allows for the approximate closure of (32) by only retaining the zeroth order terms. The solution of the obtained system will be denoted byῡ; that is, The first order corrector in W, will be useful in the error estimates. System (33) has the same structure as (1)-(3). Thus, similarly, we define by A the realization of the differential expression where ω * is a constant depending on B * and M * . Furthermore, if This estimate follows since then υ(t) is a classical solution of (33) and hence for some constant C T (remember that here · = · Xm ). Then (35) follows from (33) since the coefficients of M * are bounded. Now, instead of using the solutionῡ to (33) with the initial condition • υ, as constructed in Lemma 3.2. Accordingly, the -order corrector of w is given as We note that in the final estimates we shall set = δ here; however, we use the parameter δ to indicate approximate terms which can be controlled independently of the original small parameter .

Construction of the Layers and the Error Equation
The construction of the initial, boundary and corner layers is carried out as for irreducible C, [6], and thus we just list the final results.
Boundary Layer. The boundary layer is constructed is a similar way but by rescaling the age variable a as α := a/ . Accordingly, we defineû(a, t) := (υ(α, t),ŵ(α, t)) and, expanding the components as before, at the 0 order we get where the subscript δ indicates that the kinetic part of the boundary layer is constructed for the approximation of Corner Layer. The corner layer is constructed by simultaneous rescaling of t and a: τ := t/ and α := a/ . Proceeding as before, we get with the following side conditions:

Vol. 11 (2014)
Singularly Perturbed Population Models 547 If we assume for a moment that • u ∈ D(A ) (which will be justified in the proof of Theorem 4.1) then, by construction, all terms of the error are differentiable (a.e.) and we find that the error satisfies where the inhomogeneity F in the equation is given by and the inhomogeneity H in the boundary conditions is given by Proof. Let us fix an arbitrary T > 0. First, we note that the assumption • u ∈ W 1 1 (R + , R n ) is not sufficient for u to be differentiable so that the error equation (39) cannot be directly used. Thus, we begin by modifying the solution u to (1)-(3) by the first part of Lemma 3.2 (note that here the estimates of the derivatives of the solution are not needed). Hence, let for sufficiently small η > 0. If u η is the solution to (1)-(3) with the initial condition u η (0) = • u η then, by (24), we have MJOM for t ∈ [0, T ], where C T is a constant independent of and η. Next, letῡ η denote the (mild) solution to (33) with the initial condition υ η (0) = Since the projections are continuous, we obtain on [0, T ], for some, possibly different, constants C T , C; note that the second estimate is independent of t on account of the negative type of e τ C on QX n . Since, in general, , we use again Lemma 3.2 to approximate it by and consequently with constants depending neither on δ nor η (for sufficiently small η and δ). Now, we can use the error equation (39) for u η and the (classical) solutionῡ η δ of (33) with the initial condition We note that to be completely correct we should also regularize the corner layer since the compatibility conditions are not necessarily satisfied. Since, however, the error equations do not contain derivatives of the corner layer, we skip this step for simplicity.
By linearity, we can split the error into the part E 1 , coming from F and the initial condition E 0 (with H = 0), and E 2 which is due to the boundary inhomogeneity H. The error estimates of E 2 use the representation (27) and thus they are practically identical to those for the irreducible matrix C carried out in [6]. Hence, they are only briefly sketched. On the other hand, the estimates for E 1 use the differential formulation (39) of the error equation which is simpler to handle than the integrated version used in [6] but requires some additional estimates described below.
Let us recall that the semigroup (T A (t)) t≥0 is equibounded in , see (24). Setting H = 0, E 1 can be written using the Duhamel formula and consequently, given T < ∞, for any t ∈ [0, T ] we have In what follows, we use a generic constant c only depending on the coefficients of the problem and T, but not on the initial data or the parameters η or δ.
To estimate E 0 we see that for the first row of E 0 we have Vol. 11 (2014) Singularly Perturbed Population Models 549 where we used the first equation of (45). For the term w 1,δ (a, 0) of the second row, we use (34), the assumptions on M and the first equation of (48) to get as above For the last term in E 0 , we observe thatῡ η δ is a classical solution; hence it is continuous up to the boundaries and thus γῡ η δ | t=0 = • υ η δ (0). Next, due to (23) and the assumptions on B, for some 0 < σ < −ζ, we have where we used the fact that • υ η δ ∈ Y m and thus the trace operator at 0 can be estimated by the Y m -norm. We also used second inequality in (48).
Next, let us considerF. The term ∂ aw η 0 is well defined due to the assumption that  (50) To get the contribution ofF to the error, we begin with the term We note that this calculation is only possible sincew η 1,δ is differentiable with respect to both a and t separately which is obtained by regularization of the initial conditions. Without it, in [6] we had to resort to the cumbersome integral formulation to achieve essentially the same result. Now, using regularity of C −1 W and M, we get as above All the other estimates are the same as in [6] (compare the error terms [ and However, using the estimates (46), (47) and (49), we immediately obtain (43) and (44) for arbitrary W 1 1 initial conditions. Remark 4.1. Let us summarize and introduce a more compact notation. In the construction of the approximation of u , we use the following terms: the bulk approximationū + w 1 = m j=1ū j e j + w 1 (whereῡ = (ū 1 , . . . ,ū m )), the initial layerũ =w 0 , the boundary layerû =ŵ 0 , and the corner layeȓ u = Eυ 0 +w 0 . However, it turns out that the X n = (L 1 (R + )) n norms of w 1 ,û andȗ are of order of themselves so that, finally, uniformly on finite time intervals, provided • u ∈ (W 1 1 (R + )) n . We emphasize here that there is no need for • u ∈ D(A ) or P • u ∈ D(A); these regularity requirements were necessary only for interim calculations and are not necessary for validity of (54).

A Worked Example
In this section, we will illustrate the asymptotic procedure as applied to the McKendrick models with 5 patches. To make calculations simple, we consider age-independent matrices M, B and C with It is a continuation of Example 3.1 of [8], so the initial section is only briefly sketched. The matrix C is reducible, with σ(C) = {0, −2}. The eigenvalue λ = 0 is dominant, semisimple with multiplicity 2, while λ = −2 is of multiplicity 3. Accordingly, the null space of C is also two dimensional. The nonnegative basis of (normalized) vectors for the null space of C, described in Theorem 2.1, is given by e 1 := 1 2 , 1 2 , 0, 0, 0 and e 2 := 0, 0, where {β * ij } 1≤i≤2,1≤j≤2 is given by formulae analogous to {μ * ij } 1≤i≤2,1≤j≤2 . According to Theorem 4.1, to write down the approximation we need the initial layer; that is, the solution to ∂ τw0 = C Ww0 . Using the coordinate form (58), we find the solution as where, by (56), Thus, as in (54), we can write the combined estimates (43) and (44)

Numerical Results
For numerical simulations it is convenient to write system (1)-(3) in the integral equation form. As here it will not cause any misunderstanding, we shall drop the subscript from the notation. Using the fundamental solution matrix L(a 0 , a), as in (26), we find that, similarly to (27) The numerical approximation to (61)-(62) is obtained in two steps: first, we solve Volterra integral equation (62) for u(t, 0) and then we recover u(t, a) by integrating linear ODEs along the characteristic lines using (61).
To solve (62)  where τ k = t k − t k−1 and a j , 0 ≤ j ≤ 4 are coefficients of the BDF formula.
The starting values y(·, t k ), 0 ≤ k ≤ 3, are obtained with the 2-stage Radau IIA method, see [14]. Note that the algorithm requires one evaluation of F (t) per integration step. In our simulations this is done by means of the threepoints, composite Gauss-Lobatto quadrature rule.   It can be shown that the above algorithm converges with order four to u(t, 0) for sufficiently smooth M, C and B (i.e. the global error is O(max k τ 4 k )) in any finite interval [0, T ]; moreover, the convergence is uniform for all ε > 0. Vol. 11 (2014) Singularly Perturbed Population Models 555 Figure 2. The bulk approximation error (left) and the effect of initial layer correction (right), ε = 10 −3 The described method is rather general and does not utilize the linear convolution structure of the Volterra integral operator. The side effect is that its computational cost grows quadratically with the number of grid points. That is a serious drawback. One obvious way to overcome the difficulty is to