Rate effects on the growth of centres

Entropy maximising spatial interaction models have been widely exploited in a range of disciplines and applications: from trade and migration flows to the spread of riots and the understanding of spatial patterns in archaeological sites of interest. When embedded into a dynamic system and framed in the context of a retail model, the dynamics of centre growth poses an interesting mathematical problem, with bifurcations and phase changes, which may be addressed analytically. In this paper, we present some analysis of the continuous retail model and the corresponding discrete version, which yields insights into the effect of space on the evolving system, and an understanding of why certain retail centres are more successful than others. The slowly developing growths and the fast explosive growths that are of particular concern are explained in detail.


Introduction
The dynamics of urban centre growth has been the subject of a great deal of academic interest since the 1960s. Although a number of modelling methods now exist (Batty, 2010) the so-called production-constrained spatial interaction model has been, and remains, a particularly popular tool in this context thanks to its performance in empirical studies and its flexible framework O'Kelly, 2010). The model has its origins in the shopping model initially developed by Huff (1963) and Lakshmanan and Hansen (1965) and shares many similarities to the gravity models widely used in location theory and urban planning (Isard, 1975;Senior, 1979). It was later formalised by Wilson and Wilson (1971), embedded into a dynamic Lotka-Volterra type system (Harris and Wilson, 1978) and has since been adapted to a range of applications in theory and practice, including trade (Fry, 2012), migration (Dennett and Wilson, 2013), higher education (Singleton et al., 2012), policing strategy (Baudains et al., 2013;Davies et al., 2013), analysis of archaeological sites (Davies et al., 2014) and distribution of city sizes (Favaro and Pumain, 2011).

222
H. M. Fry and F. T. Smith In its simplest form, the model is capable of predicting the most likely set of flows between any origin and destination. Techniques to calibrate the model against empirical data are well established (Batty and Mackie, 1972) and have been widely exploited for a number of years, particularly in retail, with many large retail firms developing their own in-house versions of the model. Indeed, UK corporations such as Sainsbury's and Tesco saw unprecedented growth in the latter part of the 20th Century thanks in part to their adoption and development of geographic information systems and spatial interaction models (Birkin et al., 2010). As the application to a real-world setting has mostly been performed by those in industry, there are relatively few empirical studies in the academic literature (Johnston and Thrift, 1993). However, the model has been validated in the following studies: Clarke and Wilson (1983), Gould (1972), Haggett et al. (1977), Rae (2009), Guy (1991), and Openshaw (1976), while a useful review of the issues encountered in real-world application may be found in Birkin et al. (2010). As discussed in the review and elsewhere, the model is not without its limitations. It cannot take agglomeration effects into account, nor address a retailer's option of relocation over organic growth (Fotheringham and Knudsen, 1986). Nonetheless, it has historically provided a simple and useful framework with which to study urban growth.
Despite the model's popularity, relatively little analysis or similar insight into the dynamics of the model exists in the literature. Harris and Wilson (1978) discussed the stationary points and their stabilities, and others have performed similar investigations on extensions to the model (see, for example (Fotheringham and Knudsen, 1986)) but the majority of studies focus on simulation and adaptation of the original model. The present paper aims to explore the dynamics in more detail by using asymptotic analysis and matching. Our focus will be restricted to certain special cases, chosen to allow a greater understanding of the underlying mechanisms and shedding light on continuum and discrete descriptions and their interplay when many entities are initially present.
Although the potential adaptations are wide, it is perhaps simplest to frame our investigation in the context of retail as the model was originally written. Consider a discrete spatial region with residential areas labelled i and shopping centres j. If X i is the vector of total money spent by shoppers from residential area i, the aim initially is to determine Y ij which is the matrix of money spent in shopping centre j by residents of area i. This is achieved by assuming that a shopper from i will choose to spend money in j according to some utility, U ij : where Z j is the floorspace of centre j and d ij is a term representing cost to travel between i and j, often taken as Euclidean distance. Here, α and β are constants which weigh the relative importance of distance and size to a shopper. The logarithmic term in (1.1) is used to reflect a returns to scale of customer perception of retail floorspace. This particular form of utility function follows from Harris and Wilson (1978) is widely used throughout urban modelling and has been shown to perform well in empirical tests (Clarke and Wilson, 1983;Dearden et al., 2015;Tang et al., 2015). By framing the spending flows as a problem in information theory, it can be shown that the most likely outcome for Y ij may be determined via entropy maximisation.
Rate effects on the growth of centres

223
The derivation in the present context is based on a analogy with equilibrium statistical mechanics (see Senior, 1979) and the probability that a shopper from i will choose centre j takes the Boltzman form, given by so that the flows of money from i to j may be written as This form for Y ij makes good intuitive sense: the exponential form of the utility equation (1.1) appears in the numerator, while the denominator acts as a competition term, equivalent to a partition function in statistical physics, comparing the relative attractiveness of centre j to the finite set of alternatives. It is worth noting at this stage that the form of equation (1.2) is identical to that of the conditional logit discrete choice model (McFadden, 1980). In that setting, the utility function equation (1.1) includes an extreme value error term, which notionally takes into account fluctuations in individual taste.
The expression for Y ij in equation (1.3) is known as the singly constrained spatial interaction model and serves as a one-period model in many applications as previously discussed. However, (Harris and Wilson, 1978) noted that since equation (1.3) provides a mechanism with which to track money flows through a system it is also possible to calculate the money flows into a given retail centre, D j , by summing Y ij over i (1.4) By assuming that the per-unit costs, k, to run a centre of size Z j are uniform across the system of interest we may conclude that the profits of j are given simply by D j − kZ j . Finally, assuming that the per-unit growth (or decline) of Z j in a given period is linearly related to profit, we arrive at the final form of the retail model (Harris and Wilson, 1978) 1 for some constant > 0, to ensure a positive relationship between profits and growth. is typically small. The task then is to solve equation (1.5) for a given number of centres.
In the majority of applications equation (1.5) is used in a dimensional form. However, it is possible to choose a substitution which allows for full removal of all units from the equation, leaving the intrinsic properties in the form of just one dimensionless parameter which characterises the system: in this case a quantitative measure of relative expense. By measuring against this parameter we may directly compare systems with the same properties, regardless of scale, allowing analogies to be theoretically drawn between inner city and suburban centres. The substitutions required for non-dimensionalisation are as follows: where (Z * , X * , d * ) are typical dimensional values of (Z, X, d), leaving (Z,X,d) as dimensionless and O(1). These leave the form of the dynamics equation largely unchanged, where the tildes have been dropped for ease of notation and the dimensionless parameter κ is given by (1.8) The numerator of (1.8) speaks of a typical total running cost, while the denominator is a typical amount of money which a residential area has available to spend. The ratio of the two characterises how difficult it is for a retail centre within the system to grow. High values of κ are relatively expensive systems (from the point of view of the retailer) with lots of competition, high running costs, and/or populations with low disposable income. Low κ means that the population has high spending, and/or the relatively cheap system has little competition and/or low running costs. In this way, for a given utility function equation (1.1) we see that the dynamics of a system like central London -with high running costs and competition -will share similarities with a collection of rural retail centres where the money flows are small, as both are characterised by a large κ.
The non-dimensional model in equation (1.7) will form the basis of investigation for what follows. In Section 2, we discuss characteristics that the discrete form can promote and hence reduce the system into a one-dimensional continuous framework. In doing so, it becomes possible to scale out space entirely from the model and employ asymptotic techniques to offer insight into the remaining dynamics in three special cases: α 1 in Section 3; α near unity in Section 4 and α 1 in Section 5. We conclude with a discussion in Section 6.

Continuous and discrete versions
For different values of , equation (1.7) has been known to exhibit chaotic behaviour. This notion is illustrated by the examples given in Figure 1 This feature of sensitivity to the value of in the dynamics has similarities to the wellstudied logistic map. The mapping was first popularised in 1976 by Robert May (May, 1976) and demonstrates how complex behaviour can arise from a seemingly simple nonlinear dynamic equation. The similarities to the present retail model become particularly apparent on application of the following transformation: to the difference equation form of (1.7). In such a case, the retail model is equivalent to the logistic map where the coefficient is Working under the perhaps bold assumption that D j may be treated as independent of Z j enables the findings of May to be considered relevant here. This assumption is purely for simplification purposes, but does offer useful insights into the behaviour of the system • If D j < 0: the retail population will eventually die, regardless of the initial conditions, although in this setting the value of D j is bounded below by zero, since revenue may not be negative.
• If 0 < D j < 1: the retailers will quickly approach the value Z j = D j /κ regardless of the initial conditions. The rate of convergence is linearly related to .
• If 1 < D j < 2: the population will also eventually approach the same value but will fluctuate around that value for some time. The rate of convergence is linear, except for D j = 2 when it is dramatically slow, less than linear.
• If D j > 2: the retail sizes will demonstrate oscillatory behaviour for almost all initial conditions. At approximately D j = 2.57 the onset of chaotic behaviour occurs and oscillations no longer have finite period.

H. M. Fry and F. T. Smith
Simulations for a range of D j are shown in Figure 1(a) the results reflect what is found in the ranges above. In practice, the model forces D j > 0; and numerical examples are difficult to achieve once D j > 2 and the system becomes very volatile.
Chaotic behaviour in Z j may be avoided by considering the continuous-time version of equation (1.7). First by relating to δt and then letting δt → 0, equation (1.7) becomes but we may go further still by assuming that the number of homes and centres is relatively large, or formally tends to infinity. In that case, the continuous version of the model, on which we will base a good deal but not all of our analysis, is given by The integrals here are overall x, y in effect, with x, y denoting the spatial dimensions of the system. The task now is to solve the integro-partial differential equation (2.5) for Z(x, y, t) subject to initial conditions. To make some progress in our early investigations, we examine equation (2.5) as a one dimensional problem and define Z along s between 0 and L, with one residential area of size X at s = 0. In such a case, the governing equation reduces to the coupled form It is worth noting here that in the limit L → ∞, F(β, t) represents the Laplace transform of Z α , in effect resolving Z α into its moments. The coupled system of equations (2.6) and (2.7) has steady states at Z = 0 or when found by integrating the term in the large brackets of equation (2.6). In addition, by taking the ratio of any two arbitrary points along s when Z = 0, it is relatively straightforward to show that the non-trivial steady state will take the form This expression equation (2.9) may be combined with equation (2.8) to give an exact solution for Z(0) and hence the entire steady state, namely which holds only where α = 1. The form of the steady state equations (2.9) and (2.10) for the system equations (2.6) and (2.7) gives some early insights into the effects of α and β: in cases where β/(α − 1) is small, the leading order steady state of Z(s) is X/κL and all retail centres become equally sized. If β/(α − 1) 1, however, Z(s) will have its maximum value at s = 0 of Z(0) = βX/(κ(α − 1)) and decay exponentially from that point, independent of initial conditions. The interdependence of α and β limits the knowledge that can be gained from exploring equation (2.10). Ideally, we would want to decouple these two parameters to understand the mechanisms at play which evolve to reach this steady state. It would be relatively simple to convert the system to β = 1 for example, but in such a case the initial conditions will still exert their effects on their own scale laterally, limiting the insights we could gain from such a simplification.
Fortunately, it appears that with a system description presented as in equations (2.6) and (2.7) it is possible to scale out β entirely, separating out the effects of space and centre size, and providing a much simpler basis for analysis. This is achieved by transforming Z = Z exp (λs) and switching to a different time scale: t =t exp (−λs), where λ = β/ (α − 1). By also rewritings = exp (λs) ,L = exp (λL) ,X = λX,F = λF, ( 2.11) the governing equations (2.6) and (2.7) become and so formally the dynamics ofZ is effectively equivalent to the β = 0 scenario of equations (2.6) and (2.7). Even though a radial dependence on s will still exist when marching forward in transformed time, the new reduced system of equations (2.12) and (2.13) is without an explicit spatial structure: a potentially powerful result. Initial conditions will also exert their influence of course, and are particularly important here as the system is without explicit boundary conditions. Regardless, this zero-beta system of equations (2.12) and (2.13) is an important and central case which will provide much insight into the underlying dynamics. A further simplification may be made without loss of generality, by scaling out the magnitude of the single remaining centreX and the parameter κ, both of which are now constant for a given system. This is achieved by replacings = κs/X,L = κL/X and switching to a time scalet = κt. By dropping all bars for ease of notation then, the final reduced system on which we base most of the remainder of our analysis is: The steady states of this space-free formulation are Z = 0 and Z = 1/L and uniform in s,

228
H. M. Fry and F. T. Smith found using the equivalent form of equation (2.8). Moreover, it is relatively straightforward to compute the dynamics of equations (2.14) and (2.15) for small times. The disturbance from an initial distribution Z(s, 0) = Z 0 (s) say, takes the form for any α value, with the initial F 0 being known from the integral of Z α 0 over s. Then, from equations (2.14) and (2.15) the departures Z 1 , F 1 are determined by the balances Beyond this rough overview, it is possible to explore the dynamics of the beta-zero form of the retail model equations (2.14) and (2.15) further by employing an asymptotic approach to three special ranges, yielding analytic insights into the system. First, we explore the case of α 1 in Section 3, then α near one in Section 4 and finally α 1 in Section 5.

Alpha small
The range where α 1 relates to a retail system with diminishing returns for centre size, where customers have little or no preference for shops with larger floorspace. Specialist items and high-end products are traditionally thought to belong to this category, where the consumer considers other factors such as price, distance and quality as more important than the size of an outlet. Under this simplifying assumption, the space-free form of the retail model equations (2.14) and (2.15) may be reduced to a solvable form, offering insights into the dynamics of the full model. If we assume a form for Z and F in equations (2.14) and (2.15) of Z = Z 0 + αZ 1 + · · · and F = F 0 + αF 1 + · · · respectively, with coefficients being of order unity, we may exploit the expansion so that the full form for F follows neatly

2)
Rate effects on the growth of centres

229
The dominant balance in equation (2.14) may be solved exactly to give an analytic expression for the leading order term in centre size Z: for some C 0 (s) which is straightforward to determine by initial conditions, namely C 0 (s) = Z 0 (s, 0) −1 − L. If α = 0 then the system decays exponentially to a uniform state of Z 0 = 1/L where all revenue is equally shared amongst retailers, regardless of the initial state of Z, and there is no concentrated growth of any one centre. The next order balance of equation (2.14) leaves an expression for Z 1 which may be solved numerically where F 1 is the O(α) term in equation (3.2). However, we may go further. If the system begins with a uniform initial condition to leading order and Z 0 (s, 0) is a constant, then C 0 will also be constant and Z 0 = Z 0 (t) only. In this special case, F 1 = L ln Z 0 and equation (3.4) simplifies to give an exact solution for Z 1 Here, C 1 = C 1 (s) and is dictated by initial conditions. This solution for Z 1 does not actually require that Z 1 (s, 0) is constant, but Z 1 will decay to zero with time, tending to a solution of all centres equal. An example is shown in Figure 1(b) where Z 1 (s, 0) = sin(2πs) for α = 0.05, δs = 10 −4 and δt = 0.01. The asymptotic solutions from equations (3.3) and (3.5) are combined and plotted as circles and compared to the full time-marched solution from equations (2.14) and (2.15). The asymptotics perform well in this instance with little difference between the two solutions. The condition of a constant Z(s, 0) is a little restricting however, and it would be useful to explore the results if we allow ourselves to numerically calculate Z 1 from equation (3.4). An example of this scenario is given in Figure 2 in which Z 0 and Z are plotted for comparison in (a) and Z 1 , as found from equation (3.4), is shown in (b). Here, Z 0 (s, 0) = Z(s, 0) = s sin(2πs) + 1, δs = 10 −4 , δt = 0.01 and α = 0.01. The small value of α means the leading order approximation of Z 0 does well to match the full solution, but the dynamics of Z 1 , which are independent of α, demonstrate some interesting insights. Z 1 acts to prevent Z from reaching a uniform steady state by decreasing the values of the centres with smaller initial sizes and boosting those centres which began with a larger floorspace. This effect decays over time as Z 1 , like Z 0 , becomes constant throughout.
It is interesting to note how the two solutions, the asymptotic from equations (3.3) and (3.4) and time-marched from equations (2.14) and (2.15), compare across time for different values of α. Figure 3 shows the absolute difference at s = 0.8 between the two solutions for the same parameter values and initial conditions as in Figure 2. The plot demonstrates that the asymptotic solution holds well until around α = 0.3 when the two solutions begin to deviate. Beyond that, the residual is largest in all cases at around t = 3 230 H. M. Fry and F. T. Smith when ∂Z/∂t is greatest, but the two solutions converge once more as t increases, giving useful long-time results for values of α up to 0.6.

Alpha near unity
The critical value for α in the original discrete retail model equation (1.5) is known to be unity (Harris and Wilson, 1978), which indicates the point at which steady states transition from uniform centre size to a 'winner takes all' style system. In a real world context, calibration of this original model equation (1.5) tends to yield values of α close to unity (Clarke and Wilson, 1983). This critical point is found to be present in the space-free model of equations (2.14) and (2.15), presenting a regime which, as we shall demonstrate in the present section, produces some interesting linear and non-linear evolutionary properties. First, when α = 1 exactly, Z = Z 1 (s, t) and F = F 1 (t) say, the governing dynamics in equation (2.14) forms a Riccati equation. The behaviour of the system can be explored via the implied separable equation By setting the right-hand side to equal H(t) and setting H(0) = 0, we may define B(s) = Z 1 (s, 0) −1 (known from initial conditions), so that a complete solution in implicit where dot denotes a derivative with respect to time and equation (4.2) is equivalent to equation (2.15), since: Given an initial state of the system, Z 1 (s, 0), a steady state will be reached as t → ∞ anḋ H → 0. This will force and provides an explicit expression for the final form of Z 1 (s, ∞) : where the value of H(∞) will be determined by the integral contribution of: The steady state of Z 1 in equation (4.5) demonstrates the dependence on the initial conditions of the system. The results here are consistent with the existing literature on the full model (Clarke and Wilson, 1983;Harris and Wilson, 1978), agree with the small-time features discussed near the end of Section 2 and confirm that in a sense boundary conditions are not required. If, instead, α is close to unity, the dynamics may be explored by supposing that: for δ 1 and α 1 = O(1). The resulting approximations to the governing framework are sensitive to time and an asymptotic approach here produces early and later stages of evolution for centre size Z as discussed in the following sub-sections.

The early stage
While time is O(1), the solution to the coupled governing space-free system equations (2.14) and (2.15) can be explored using a perturbed modification to Z and F of the form (4.8) (4.9) with δ as in equation (4.7). The dominant balance of the dynamic equation (2.14) gives an expression identical to that seen in the α = 1 case of the previous section in equation (4.1). The solutions and steady states given there for Z 1 (s, t) and F 1 (t) in equations (4.5), (4.4) and (4.3) are also relevant here. The secondary terms in equations (4.8) and (4.9) may be determined by looking at an order δ balance of equations (2.14) and (2.15), yielding expressions for Z 2 and F 2 respectively (4.10) At large times, we expect F 1 → 1 (as in equation (4.4)) so that the first term in equation (4.10) will tend to zero. If F 2 (t) also tends to a steady state, that is F 2 (t) → F 2 (∞), the second term in equation (4.10) will become constant in time and we may define implying a large-time response from Z 2 (s, t) as: (4.13) For this solution to hold, we require the temporal component of equation ( where Z 1 (s, ∞) is given in equation (4.5).
As t → ∞ then, |Z 2 | will continue to grow linearly with time, while |Z 1 | remains of order unity. This causes a large time breakdown of the asymptotics in equation (4.8) as non-linear effects begin to re-emerge. This second stage occurs as t becomes of order 1/δ and is explored in the following section.

The middle stage
The first later stage occurs when t =t/δ, with δ as defined in equation (4.7). Scaled time runs forward from zero where a match can be shown with the results from the previous section. Here, a simplification of the governing system in equations (2.14) and (2.15) can be made by application of the following substitutions for Z and F: Z =z(s,t) + · · · , (4.16) (4.17) A leading order balance of the dynamic equation in (2.14) again requires (4.18) as seen in previous sections. Here however, a secondary balance, after suitable manipulation ofz α , presents a non-linear pde forz: without forcingF 2 to reach a steady state. A complete analytical solution of the non-linear pair equation (4.19) with equation (4.18) cannot be obtained readily, but it is relatively straightforward to perform a numerical study of the system. In a similar manner to the method of the previous section, taking the integral of equation (4.19) with respect to s from 0 to L makes the left-hand side zero, by equation (4.18). Rearranging the result then offers an expression forF 2 in terms ofz(s,t)  Figure 4 for when α 1 = 1. In each case, the initial condition Z = (1 + s sin(πs)) is used so that ∂Z/∂s is small,   range of δs and δt regimes are tested in the numerics as shown, which offers insight into some interesting behaviour. In all cases, the dynamics is slow at first, with little deviation from the initial conditions (shown as a black-dashed line) by t = 1. At t = 2, the values of Z which were initially at the site of local maxima have begun to increase, while Z elsewhere decreases to accommodate the integral condition in equation (4.18). Beyond t = 2, the main peak gathers momentum and increases rapidly, dominating the contribution to the integral in equation (4.18) and forcing Z to decrease for all other values of s. By t = 2.8, only a single point in each numerical solution has a non-zero value -which raises an interesting scenario.
The area under the curve is fixed by equation (4.18), regardless of δt and δs, but the triangular shape which is created, in effect by discretising this continuous system for the numerical treatment, has a base whose width depends on the discretisation used. In effect then, the final value of the single non-zero point in Z(s) must depend on δs, as verified in Figure 5 where the steady-state solutions increase for smaller values of δs.
This perhaps somewhat unusual result appears through the interplay between discrete and continuous versions of the model. The computational terminal values for bothF 2 and max(z) can be calculated analyticallȳ (4.21) and are indeed determined by the numerical discretisation applied. These terminal forms equation (4.21) are also shown in Figure 5 as a single dashed line for each case, but would each be infinite in a true representation of the continuous system. Returning to the continuous system in equation (4.19) with equation (4.18), the logarithmic term in equation (4.19) will balance withF 2 near a peak to leading order, so that the term in the bracket becomes O(1) and ∂z ∂t balances withz 2 . This suggests that formallyz ∼ 1/t and indicates an ending in finite time. Specifically, we expect a solution of the formz for some positive constant A and blow-up at t = t 0 . In this scenario,F 2 tends to ∞ as verified to a large extent by the computations.

H. M. Fry and F. T. Smith
At other positions of s, away from the peak,z is not large enough for the logarithmic term to play a role in equation (4.19). Instead, ∂z/∂t =z 2 ln(t −t 0 ), suggesting a terminal form as follows:z (4.24) and thusz tends to 1/C ast → ∞. This fits together well with the findings from the numerical simulations and presents a near-complete description of the system for α values near unity. The finite-time breakup above indicates that discrete aspects come back into play at positions close to the peak.

The final stage
The growth ofz close to a peak causesF 2 to grow logarithmically with t as in equation (4.23) and violate the assumptions in the series expansion of equation (4.17). Simultaneously, the growth effect ofz at the peak, s = s 0 say, is no longer unity to leading order and we begin to see exponentially large values for Z within a tiny range of values for s, close to s 0 , and in an exponentially small time interval. Specifically, where s−s 0 = O(e −n/δ ) for some n of order unity, we expect Z, F and t to take the form (4.25) where δ is as defined in equation (4.7). These substitutions make Z α−1 = e α1t * predominantly which, when applied with equation (4.25) to the governing equations of the system equation (2.14) and (2.15), show that the leading order balance occurs between the contributions within the parentheses of equation (2.14). To the first approximation then, we expect (4.26) This final stage solution fits with that of the middle stage, where close to a peak equations (4.17) and (4.23) combine to give F = 1 − δB ln t 0 −t + · · · for some constant B. Given thatt = δt = e −t * /δ +t 0 (as in equation (4.26)), F from the middle stage may be rewritten as F = 1 + Bt * which matches with the solution for F in equation (4.26) when t * tends to zero from above and B = α 1 . The exponential growth in F suggests that the −Z 2 term in equation (2.14) will dominate the dynamics outside of this small range of s, causing all other centres to decay to zero. The exponential growth in Z near the peak then serves to satisfy the integral condition of equation (2.14) and indeed the solution continues to hold to large positive effective times t * . It is worth noting that the 'exponential' growth in both F and Z is actually in terms of the exponentially defined t * in equation (4.25), meaning that the true temporal behaviour is algebraic.
The overall effect of the three stages, early, middle and final, is that of delayed but eventually explosive growth in centre size, arising about a single location.
Rate effects on the growth of centres 237

Large alpha
The final range to explore is that of α 1. This relates to a retail system where high importance is placed on the size of a retail site and large centres tend to prosper. The original production-constrained retail system, given in equation (1.5), has been shown to permit stable steady states with α > 1 for a range of small k (Clarke and Wilson, 1983) (Harris and Wilson, 1978), a result that has been also found to exist in a real-world setting. In the present case, we restrict our analysis to the space-free system given in equations (2.14) and (2.15) which instead allows large centres to grow explosively over a comparatively short time scale as we shall see. Translation between results from the reduced form used here and the original formulation in equation (1.5) is discussed in Section 6.
For a sufficiently large α there are considerable changes in the magnitude of the Z α contribution for values of Z either side of unity. To mitigate against this, it is convenient to work with the following expansion for Z and F: With α 1, we expect a fast time response and so switch to a shorter time-scale t = T α . Together, these serve to reduce the balance of equations (2.14) and (2.15) to to leading order. For points along s where z < 0, centres ultimately decay linearly at a constant rate. Wherever z > 0 however, the exponential term in the dynamics equation (5.3) dominates, which begins to imply how dramatic the rapid growth of larger centres is within the system. Further insight is gained by a solution to equation (5.3), achieved by putting  throughout s and a global maximum at around s = 0.62. The constants (δs, δt, L) were set at (10 −4 , 10 −4 , 1) respectively throughout the numerics.
Despite the relatively small value for α used in the numerical example, the asymptotic solution performs well, matching the full time-marched system from equations (2.14) and (2.15) even as the dominant peak begins to form. As expected from equation (5.3), the peak emerges at the site of the initial global maximum and grows rapidly. Away from the peak, F 0 and Z become constant, yielding linear decay in z. As z → −∞ at these points, Z → 0. This process continues until a single non-zero point remains in the numerically discretised solution, as in Section 4.

Discussion
An asymptotic analysis of the continuous, one-dimensional, single source version of the model the retail model has revealed a number of interesting insights into the dynamics that underpin the full system. The analysis, based on the coupled equations (2.14) and (2.15) is equivalent to a continuous β = 0 case of (1.5). Such a simplification presents an opportunity to explore the effects of the parameter α and its role in the dynamics of centre size Z.
If α = 0 an analytic solution to Z was found in which each centre decays exponentially to a uniform state of Z = 1/L regardless of the initial conditions. For positive α 1, the higher order terms in an asymptotic analysis serve to slow the approach to a steady state. Small α will boost the centres which were initially smaller and decrease the size of those which began with a higher value of Z. This effect itself will decay in time, so that whenever α < 1 all centres will eventually become equal.
The critical value of α = 1 has been discussed previously in the literature. In this setting, it marks the transition from a steady state of uniform centre size to a system where the 'winner takes all'. Exploring the dynamics close to unity by setting α = 1 + δα 1 demonstrates a delayed but eventually explosive growth of Z near a peak for α 1 > 0. This growth will be exponential in asymptotic time, although algebraic in true temporal terms. Away from the peak, centres will decay to a constant value.
When α 1 an analytic solution is again possible. Here, the system experiences rapid growth at the site of the original global maximum and linear decay elsewhere.
While these behaviours match the characteristic dynamics of real-world retail systems (Clarke et al., 1998), by scaling out β from the system we introduced a radial dependence on time in equations (2.12) and (2.13) that precludes a straightforward translation between the special cases outlined here and the full system of (1.5). Nonetheless, for a qualitative comparison with the results given here, some simulations of the original model are included in Figures 7(a) and (b).
The simulations use an initial condition of Z = {j : j ∈ N, 0 < j 6 10} with the centres positioned along the integers of the x axis from x = 0 to x = 9. A single source of X = 55 is positioned at (x, y) = (0, 1) and = 10 −3 is used. To allow for a direct comparison with the simplified model in (2.14) and (2.15) a value of β = 0 is applied throughout. Figure 7(a) demonstrates the evolution of the mean centre size over time for the cases of α = 0.1, 0.9, 1.1, 2.0.
When α = 0.1 the speed of the exponential decay to a steady state for each centre is dependent on their initial size. Hence, the total size of the system drops quickly in Figure  7(a) as the initially large centres shrink much more rapidly than the growth of the initially small centres. As time continues the system approaches its steady state of uniform centre size, a result which is reflected in Figure 7(b) which shows the standard deviation of Z j .
The cases of α = 0.9 and α = 1.1 have clear similarities and both show the characteristic slow-burning dynamics described in Section 4. Over a much longer time scale than the small α case, the standard deviation of centre size does continue to grow and decay for α > 1, α < 1 respectively, as predicted by the asymptotics.
Finally, when α = 2 the centre with the largest initial value of Z experiences rapid growth, while away from the peak the other centres slowly decay. This is reflected in  (c) Is the mean centre size for an identical system, except with β = 1 and the respective standard deviation is given in (d). Figure 7(b) where the standard deviation of the system slowly reaches a steady state of a single non-zero centre. It should be noted that the finite terminal value of this steady state is a result of the discrete system and would be infinite in the continuous model (see Section 4).
An extension of the current framework to take the spatial structure of the system into account (or rather, a direct translation between equations (2.14), (2.15) and (2.6), (2.7)) would likely demonstrate that a stable, multi-centre steady state is possible for α > 1, although this is beyond the scope of the current study. However, the simulations of were Figures 7(a) and (b) were repeated with β = 1 and are presented in Figures 7(c) and (d) to offer some insight into the corresponding full system. From these graphs, it becomes clear that the inclusion of space provides additional mechanisms that are worthy of further investigation.
The analytical solutions for Z, and their respective discrete simulations, also reveal the dependence of the model on initial conditions across a range of values of α, while an interplay between discrete and continuous system suggests some caution in the middle stage when α is near unity. Both these are important considerations in real-world application.