## Abstract

Obtaining information about pest-insect population size is an important problem of pest monitoring and control. Usually, this problem has to be solved based on scarce spatial data about the population density. The problem of monitoring can thus be linked to a more general mathematical problem of numerical integration on a coarse grid. Numerical integration on coarse grids has rarely been considered in literature as it is usually assumed that the grid can be refined. However, this is not the case in ecological monitoring where fine grids are not available. In this paper, we introduce a method of numerical integration that allows one to accurately evaluate an integral on a coarse grid. The method is tested on several functions with different properties to show its effectiveness. We then use the method to obtain an estimate of the population size for different population distributions and show that an ecologically reasonable accuracy can be achieved on a very coarse grid consisting of just a few points. Finally, we summarize our mathematical findings as a protocol of ecological monitoring, thus sending a clear and practically important message to ecologists and pest-control specialists.

## 1. Introduction

Monitoring of pest insects is an essential part of the integrated pest control and management program (Metcalf 1982; Burn 1987). Its purpose is to provide information about the population size of the pest species in a given agricultural field or area. An increase in the population size above a certain threshold is regarded as a signal to apply pesticides (Stern 1973; Ester & van Rozen 2005). Delays in pesticide treatment may result in pest population outbreak and hence in irretrievable damage to the agricultural product. On the other hand, application of pesticides is costly and brings considerable damage to the environment (Jepson & Thacker 1990). Therefore, there is a need for reliable methods to estimate the pest population size in order to avoid unjustified pesticide application and yet to prevent pest outbreaks. The demand for such methods is getting stronger nowadays as an observed increase in the average annual temperature is expected to result in more pest outbreaks (UK Department of Environment, Food and Rural Affairs 2009).

Depending on the pest biology, the information about pest population size is usually collected through sampling (Mayor & Davies 1976; Blackshaw 1983) or trapping (Holland *et al.* 1999; Alexander *et al.* 2005). In particular, in the case of trapping, a typical procedure is as follows. A number of traps are installed in the field. Generally speaking, the position of the traps can be arbitrary (Mayor & Davies 1976; Alexander *et al.* 2005), but usually they are placed at the nodes of a rectangular grid (Holland *et al.* 1999; Ferguson *et al.* 2003). The traps are exposed for a certain time (e.g. for a week), and after that they are emptied and the caught insects are counted. The trap counts obtained in this way provide information about the pest population density at the position of the traps (Byers *et al.* 1989; Raworth & Choi 2001).

Once the samples are collected, the next step is to estimate the population size, i.e. the total number of the insects in the field. There are a variety of statistical methods that are used for this purpose based on the sample counts (Cochran 1977; Seber 1982; Sutherland 1996). However, they all have their limitations that make them ineffective in routine pest-monitoring procedures. A relatively accurate method known as ‘mark–release–recapture’ requires considerable additional effort (such as insect marking and recapture) that works well in scientific studies, but can hardly be afforded in nation-wide monitoring programs. Simpler but less accurate approaches are usually heavily based on an assumption about the species spatial distribution, e.g. assuming it to be approximately homogeneous throughout the monitored area (Sutherland 1996, p. 86).

In this paper, we endeavour to develop an alternative approach to obtain the population size. In order to derive an estimate of the population size, the space-discrete function (the population density) obtained in trap counts should be integrated over the whole monitored area by a chosen method of numerical integration. Indeed, from a mathematical viewpoint, installing traps in a domain where sampled data are collected means that the discrete integrand function is defined on a (uniform) computational grid. This is a conventional problem of numerical integration, as the integration of sampled data problem frequently arises in experimental work as well as in computational applications. However, the integration of a discrete function usually implies that the grid over the domain where the integration is to be taken can be made sufficiently fine to provide required accuracy. Moreover, if/where necessary, integration can be taken over a sequence of refined grids, which gives thorough information about the integration error, convergency rate, etc.

The situation in the pest-monitoring problem is basically different. The issue is that, in routine monitoring, the number of traps installed in a field cannot be made large. Although this is not a principal limitation of the method, its practical reasons are important and cannot be neglected. First, installment of many traps per unit agricultural area would, by itself, create considerable damage to the agricultural product and hence would make the whole procedure rather senseless. The second important factor that must be taken into account is that financial and labour resources available for monitoring are always limited. Thus, the number of traps installed over a typical British agricultural field (with a characteristic size on the order of a few hundred metres) and/or the number of samples taken very rarely exceed a few dozen (Blackshaw 1983; Holland *et al.* 1999; Ferguson *et al.* 2003). Correspondingly, when traps are placed at the nodes of a rectangular grid, the number of grid nodes in each direction is usually taken to be between 5 and 10 (Holland *et al.* 1999; Alexander *et al.* 2005). Finally, another factor that makes a difference is that repeated trapping with an increased number of traps is not available in ecological applications because of the impossibility of reproducing the initial conditions. Obviously, the integrand function would be essentially different if a ‘refined grid’ of traps would be generated on the same agricultural field but later in time.

The above arguments show that the problem of pest monitoring requires numerical integration of a discrete function obtained on a very coarse grid. The primary goal of this paper is therefore to address the problem of pest monitoring using the tools and methods of computational mathematics, in particular, using the method of numerical integration recently developed in Petrovskaya & Venturino (2010). To the best of our knowledge, although many accurate and efficient methods for spatial data reconstruction are well known and documented in the literature (e.g. Davis 1975; Piessens *et al.* 1983), the problem of numerical integration on a coarse grid with a fixed number of grid points has not been properly addressed yet. Here, we are especially interested in the following aspects of the problem:

— how to obtain the integral evaluate with the maximum accuracy possible on a given coarse grid, especially in cases when the integrand may have a non-trivial spatial structure, such as in regions with a large solution gradient, irregular spatial oscillations, etc. and

— what can be a criterion to distinguish between ‘coarse’ and ‘fine’ grids in the problem, i.e. whether it is possible to make any conclusion about the reliability of numerical integration on a grid with a fixed number of grid nodes where asymptotic error estimates are not available in the problem.

Since the purpose of this paper is to address the problem of coarse-grid integration conceptually rather than to develop a complete algorithm ready to use in ecological practice, we restrict our analysis to one-dimensional problems. It was shown in Alexander *et al.* (2005) that, in some ‘real-life’ ecological cases, placing traps along a line can be a more efficient sampling strategy than placing them on a two-dimensional grid; thus, our analysis is not entirely abstract.

The paper is organized as follows. In §2, we introduce a method to perform accurate numerical integration on a coarse grid. Several integrands with different properties are considered in §3 in order to validate the method’s capability. In §4, we introduce a generic population-dynamics model to generate different spatial population distributions and to show that implementation of our method allows one to achieve reasonable accuracy on a coarse grid. A further discussion of the method properties and limitations, as well as the message for its use in ecological applications, will be provided in §5.

## 2. Numerical integration on a coarse grid: a method

The idea of the Gauss–Legendre rule with interpolation (GLI) that we use for numerical integration on a coarse grid is to employ the Gauss–Legendre rule at each subinterval of a uniform grid in order to increase the accuracy of computations. It is well known (e.g. Davis & Rabinowitz 1975) that the general *M*-point Gauss–Legendre rule is exact for polynomial functions of degree ≤2*M*−1. The data at the roots of the Legendre polynomials are not available to us by the problem statement. However, we can interpolate those data from a uniform grid and further use them for numerical integration. A detailed discussion of this integration technique is provided in Petrovskaya & Venturino (2010). Below we briefly describe the results of work (Petrovskaya & Venturino 2010) that will further be used in this paper.

Consider a set of points *x*_{n}, *n*=0,1,…,*N*, in the domain [*a*,*b*], where we require that *x*_{0}=*a*, *x*_{n+1}=*x*_{n}+*h*, *i*=0,…,*N*−1, and refer to set {*x*_{n}} as a uniform computational grid with nodes *x*_{n}. The grid step size *h* is defined as *h*=(*b*−*a*)/*N*. Consider a function *f*(*x*) and let data *f*_{n}=*f*(*x*_{n}) be known at each point *x*_{i}. A general problem of numerical integration is that the integral
2.1
should be approximated by the sum to meet the condition
2.2
where *ϵ* is the given tolerance.

Consider the integral in the domain [*X*_{1},*X*_{2}], where . The standard approach in numerical integration is to expand the integrand *f*(*x*) as
where *ϕ*_{k}(*x*) are polynomial basis functions, *f*(*x*_{k}) are function values at interpolation nodes and *r*(*x*) is the remainder term. Substituting the above expansion into the integral, we arrive at

Thus, the integral is where weights are and the integral remainder is given by .

Let us now assume that we do not know the precise values *f*(*x*_{k}), *k*=0,1,…,*K*. Instead, the function *f*(*x*) is defined with the error *ϵ*_{k} at each point *x*_{k},
We then have
2.3
where the new remainder *R** is

Consider a numerical integration method on a uniform grid that has the integration error *R*_{u}. Now let us interpolate the function *f*(*x*) from a uniform grid at points *x*_{k} of a non-uniform grid with interpolation error *ϵ*_{k} and evaluate integral *I* as equation (2.3). Obviously, the Gauss–Legendre formula of integration can then be applied, if the points *x*_{k} are chosen to be the roots of the Legendre polynomials. The basic requirement for such an approach is that the interpolation error for the Gauss–Legendre formula should not be larger than the integration error for a numerical integration method (e.g. Simpson’s rule) employed on a uniform grid. Since we want the integral value (2.3) made on the non-uniform grid to be more accurate than the numerical integration on the original uniform grid, we require that
2.4

Let the remainder of each method depend on grid step size *h* as |*R*|≤*Ch*^{r1}, |*R*_{u}|≤*C*_{u}*h*^{r2}, where *r*_{1},*r*_{2}>1 and *C*_{u} and *C* are constants. Then, the condition (2.4) holds, if we have another constant *C*_{ϵ} and |*R*_{ϵ}|≤*C*_{ϵ}*h*^{r}, where we require that when .

Let us note here that in the problem of ecological monitoring, the number of grid points is fixed and *h* is usually not small. However, for a consistent analysis of our method, here and in the next section, we rather consider the problem in a general context of numerical integration and therefore assume that the grid can be refined.

The error *ϵ*_{k} depends on the method chosen to approximate the function *f*(*x*) at points {*x*_{k}}. We use the Lagrange interpolation in order to compute the function values on a non-uniform grid. Below, we briefly recall a general Lagrange interpolation procedure (e.g. Davis 1975).

Suppose that we want to interpolate a smooth function *g*(*x*) over an arbitrary interval [*X*_{1},*X*_{2}]. The Lagrange interpolation requires the knowledge of values of the function *g*(*x*) at some points , *q*=0,…,*p*, that belong to the interval [*X*_{1},*X*_{2}]. An interpolating polynomial of degree *p* is then defined at the interval [*X*_{1},*X*_{2}] as
2.5
where the basis functions are given by
For instance, if we interpolate a function *g*(*x*) by a polynomial of degree *p*=5 over the interval [*X*_{1},*X*_{2}], we need six points where the values are defined at points .

Let us assume that points are defined as
where *h*=(*X*_{2}−*X*_{1})/*p*. We also assume that *g*(*x*) and the derivatives of *g*(*x*) up to the order (*p*+1) are continuous and bounded on the interval [*X*_{1},*X*_{2}]. Then, the error of the Lagrange interpolation is
2.6
where *C*_{p} is constant.

Let us apply the interpolation (2.5) in our problem. By the problem statement, the function *g*(*x*) in equation (2.5) is given by the integrand *f*(*x*). Consider the subinterval *e*_{n}=[*x*_{n},*x*_{n+1}], *n*=0,…,*N*−1, of a uniform grid and let , *k*=1,…,*M*, be the abscissae of the Gauss–Legendre formula at the subinterval *e*_{n}. Since , *k*=1,…,*M*, are interior points of the subinterval *e*_{n}, the integrand *f*(*x*) is not defined at points . Hence, we interpolate the values of the function *f*(*x*) from a uniform grid to points as , where the notation reads that the function *f*(*x*) is interpolated by a Lagrange polynomial of degree *p* in order to reconstruct the function values at *k* interior points of the subinterval *e*_{n}.

The values of the function *f*(*x*) required in formula (2.5) are available to us at nodes of a uniform grid only. Therefore, we should take points , *q*=0,…,*p*, in the formula (2.5) as nodes of a uniform grid generated in the problem. For instance, consider the interpolation by a polynomial of degree *p*=5 over the subinterval *e*_{n}=[*x*_{n},*x*_{n+1}]. The six points we need for the interpolation are then selected over a uniform grid as , . Obviously, if the boundary subinterval *e*_{0}=[*x*_{0},*x*_{1}] or *e*_{N−1}=[*x*_{N−1},*x*_{N}] is considered, then the interpolation stencil should be shifted to the right or to the left, respectively.

Once the function values have been defined at the abscissae of the Gauss–Legendre formula at the subinterval *e*_{n}=[*x*_{n},*x*_{n+1}], we can apply the Gauss–Legendre integration rule at the subinterval *e*_{n} and get the evaluate . We then compute the integral over the domain [*a*,*b*] as . Since the remainder of the *M*-point Gauss–Legendre rule at each subinterval *e*_{n} is given by
(e.g. Davis & Rabinowitz 1975) the error (2.6) becomes dominant in the evaluation (2.4) for *M*>*p*/2. The error (2.6) should then be compared with the error *R*_{u} of a numerical integration method on a uniform grid. Let us, for instance, use the compound Simpson’s rule, then the error *R*_{u} is given by
2.7
where *C*_{S} is a constant (Davis & Rabinowitz 1975). It can be readily seen from equations (2.6) and (2.7) that the polynomial degree in equation (2.5) should be chosen as *p*≥4 to make the error (2.6) smaller than that in the compound Simpson’s method. This should result in more accurate numerical integration in comparison with Simpson’s rule used on the same grid.

Finally, let us note here that Simpson’s rule belongs to the well-known family of Newton–Cotes (N–C) integration formulas that replace an integrand function with an interpolating polynomial on uniform grids (Davis & Rabinowitz 1975). Simpson’s rule requires one to replace the integrand with a quadratic function, while higher order N–C formulas exploit higher order polynomials in order to evaluate the integral. However, a problem with the N–C rules that considerably limits their use in practical applications is that increasing the number *N* of grid points used for numerical integration does not necessarily lead to a more accurate integral evaluate. The example considered in Davis & Rabinowitz (1975) demonstrated that increasing the number *N* from *N*=2 up to 21 resulted in a larger integration error. Another serious drawback of higher order N–C formulas is that it is often difficult to apply a compound N–C rule with the polynomial degree *p*>2 on a uniform grid with an arbitrary number of grid subintervals, as the total number *N* of grid cells is required to be a multiple of *p*.

The main difference between our method and the N–C integration rules is that we do not directly integrate the Lagrange interpolant on a uniform grid, but use it for the evaluation of the function at the abscissae of the Gauss–Legendre formula. The above difference makes our method more flexible in comparison with higher order N–C rules as the compound method can now be applied on a grid with arbitrary number of grid nodes. The GLI method should also remain accurate when we vary the number of grid subintervals, as we use a local interpolation stencil that does not depend on the total number of grid nodes. For instance, if we have a grid of *N*=23 subintervals, the N–C rules allow us to use either a single polynomial of degree *p*=23 over the grid or a compound trapezoidal rule (the first rule in the N–C family corresponding to the polynomial degree *p*=1). The integral evaluate will not be accurate in both cases (see the discussion in Davis & Rabinowitz 1975). On the other hand, our approach makes it possible to use a compound rule with the polynomial degree *p*=5 on such a grid, the accuracy of integration being controlled by equation (2.6).

## 3. The coarse-grid problem: test cases

While the detailed validation of the GLI method has been performed in Petrovskaya & Venturino (2010), here we consider several standard test cases that illustrate our approach and discuss the accuracy of the GLI method. For the purpose of our discussion, we need to look at the convergence rate of the GLI method, as the study of the convergence rate is essential for the understanding of the coarse-grid problem. Hence, we assume that a coarse uniform grid of *N*=8 cells^{1} can be uniformly refined by halving each cell. The integral is then computed on a sequence of uniformly refined grids where the error norm ||*e*|| is defined on each grid in the sequence as
3.1
In the error definition (3.1), *I*_{p} is the exact integral and *I*_{a} is the approximate integral obtained by a chosen method (i.e. the GLI method with a chosen polynomial degree or Simpson’s rule). Obviously, the above assumption does not work in real-life applications where we have a single grid at hand, but we need it in order to compare the convergence rate for the GLI method with that for Simpson’s rule.

Note that the test cases considered in this section are not supposed to have an immediate ecological meaning, but they are important to validate our mathematical approach.

We first consider the integrand at the interval [0,*π*] shown in figure 1*a*. The convergence curve for this integrand function is shown in figure 2*a*, where the polynomial degrees *p*=1,3,5 have been considered in the Lagrange interpolant (2.5) used in the GLI method. The convergence rate for each value of *p* is shown in the figure, along with the convergence curve obtained for Simpson’s rule. The error norm ||*e*|| as a function of the number *N* of grid cells is presented in the figure on a logarithmic scale.

The above test case illustrates well the properties of the GLI method. The convergence for *p*=1 is slower than the convergence for Simpson’s rule. The convergence line for *p*=3 is parallel to that for Simpson’s rule, while the interpolation by a higher order polynomial *p*=5 results in a faster convergence rate, as predicted by the estimate (2.6).

We now consider the function
3.2
where the parameter *μ*>0. The function (3.2) is often thought of as the solution to the following ‘convection-dominated’ boundary-value problem (e.g. Carey 1997)
3.3

The solution *f*(*x*) depends on the parameter *μ* as shown in figure 1*b*, where *f*(*x*) is displayed for *μ*=1.0 and 0.005. It can be seen from the figure that a narrow region of the steep gradient (‘a boundary layer’) is present near the point *x*=1.0 for small values of *μ*. It is readily seen that the width |*r*_{b}| of the boundary layer can be estimated as |*r*_{b}|∼*O*(*μ*).

The convergence history for the integrand (3.2) is shown in figure 2*b*,*c* for *μ*=1.0 and 0.005, respectively. The convergence graphs in both cases are obtained for the GLI method with *p*=5 and for Simpson’s method. The results presented in figure 2*b* are very similar to those in figure 2*a*. According to the error estimate (2.6), the convergence of the GLI method with higher polynomial degree is better than the convergence rate for Simpson’s rule.

Meanwhile, the situation becomes different when we decrease the parameter *μ* and a sharp boundary layer appears near the right endpoint. The initial coarse grid with grid step size *h*=0.125 cannot capture the region of the steep gradient that has the width |*r*_{b}|∼0.005. As a result, the error of numerical integration is very large for both the GLI and Simpson’s method. Obviously, the value of the integral obtained on such a grid cannot be considered as reliable in practical applications and a finer grid would be required in the problem to resolve the integrand.

Let us make uniform refinement of the initial grid. Since we use approximation by polynomials of degree *p*=5 in the GLI method, we can expect a faster convergence rate in comparison with Simpson’s rule on refined grids. However, it can be seen from the figure that, after several refinement cycles, both the methods still have a slow convergence rate. The convergence of the GLI method is not getting better than the convergence of Simpson’s method until a sufficient number of nodes have been inserted in the boundary-layer region. This occurs when the grid is refined up to the number of grid cells marked as *N** in figure 2*c*, where *N**≈200 for the resolution of the boundary layer defined by *μ*=0.005. Thus, we suggest that any grid with *N*<*N** should be considered as a coarse grid that does not provide us with a reliable value of the integral. In other words, the asymptotic error estimate (2.6) does not hold on any grid with *N*<*N**. Accordingly, any grid with *N*>*N** is a fine grid where numerical integration will give us a reliable integral evaluate.

The problem is further illustrated by consideration of the following fast-oscillating function:
3.4
shown in figure 1*c*. The convergence of the numerical integration methods is presented in figure 2*d*. Again, the integration error on the initial coarse grid is very large for both the GLI *p*=5 and Simpson’s method. The transition from coarse to fine grids happens after several cycles of grid refinement; it is readily seen that a grid of *N**≈180 can be regarded as sufficiently refined. The convergence rate cannot be predicted in the domain of coarse grids, while asymptotic error estimates hold on fine grids with *N*>*N**. An important conclusion that follows directly from the discussion above is that using higher order polynomial degrees in Lagrange interpolations is not advantageous on coarse grids, as it does not reduce the integration error. The integration error of the GLI method remains approximately the same as the error of Simpson’s rule, until the grid is refined up to the number *N** of subintervals.

It is worth mentioning here that, in some cases, the number *N** of grid points when the grid becomes ‘sufficiently refined’, i.e. the convergence rate reaches its asymptotical value, can be estimated heuristically. Indeed, let Δ*x* be a characteristic width of a spatial heterogeneity described by a given integrand, e.g. the width of the boundary layer in equation (3.2) or the width of a single hump in equation (3.4). Then a uniform grid is not refined enough until at least one grid point (i.e. one detector or one trap, if we think about the problem from a practical viewpoint) falls into the heterogeneity region. Therefore,
3.5
where 1 in the numerator stands for the length of the whole domain. Here, *k*≥1 is a numerical coefficient depending on the type of the heterogeneity: while a boundary layer can be ‘resolved’ by stationing just one grid point in the region of the heterogeneity, one needs at least three grid points to resolve a hump. That can be seen from the following argument. If *f*(*x*) is a monotone function at the interval [*x*_{1},*x*_{2}], then we consider the function values at two points, e.g. *x*_{1}+*δ* and *x*_{2}−*δ*, where 0<*δ*<0.5(*x*_{1}+*x*_{2}), as the minimum ‘amount of information’ required to reconstruct *f*(*x*) over [*x*_{1},*x*_{2}], as these data are sufficient for linear polynomial approximation of *f*(*x*). Consequently, for a function with a ‘hump’ (that is an extremum point *x*_{2}∈[*x*_{1},*x*_{3}]) we need at least three points, e.g. *x*_{1}+*δ*, *x*_{2} and *x*_{3}−*δ*, as this will result in linear approximations at each subintervals [*x*_{1},*x*_{2}] and [*x*_{2},*x*_{3}] where *f*(*x*) is a monotone function.

Correspondingly, from equation (3.5) we arrive at *N**≈200 for function (3.2) (with *μ*=0.005) and at *N**≈150 for function (3.4), which is in good agreement with the results of the more accurate analysis shown in figure 2*c*,*d*.

Summarizing the above, our approach is based on the comparison between the two methods of numerical integration available on a given uniform grid. On fine grids, the GLI method with *p*=5 is proved to have better accuracy against Simpson’s method. Therefore, if the error of Simpson’s method appears to be smaller or compatible with that in the GLI method on a given grid with the fixed number of grid nodes, then we can conclude that the grid is too coarse to provide us with a reliable integral evaluate. In the next section, we extend our results to the coarse-grid problem in ecological applications.

## 4. Ecological model

In the view of the suggested ecological applications, the integrand functions considered in the previous section described somewhat abstract test cases. The problem is that their spatial structure is difficult to interpret in ecological terms, nor could it be related to the values of any ecologically meaningful parameters. Therefore, applicability of our method to the problem of ecological monitoring yet remains to be demonstrated.

Ideally, the capability of the method to deal with an integrand on a coarse grid should be checked against appropriate ecological data. However, we are not currently aware about real-world data that could be used for this purpose: in order to make a sensible comparison between our method and data, we must know not only the trap counts but the actual population size of the pest population, which can hardly be available. Thus, in this study, instead of using real field data on trap counts, we use predictions of a generic model describing spatiotemporal dynamics of a pest-insect population. Mathematically, such a system is described by a system of coupled diffusion–reaction equations (Murray 1989). Namely, we consider the following system: 4.1 and 4.2 which is also known as the Rosenzweig–MacArthur model.

The model describes a prey–predator system where the prey is the pest insect and the predator is its consumer (e.g. an avian population). The functions *U* and *V* in equations (4.1) and (4.2) are the densities of prey and predator, respectively, at time *T* and position *X*, *T*>0 and 0<*X*<*L* where *L* is the typical size of a given agricultural field. Coefficients *D*_{1} and *D*_{2} describe species diffusivity due to the movement of the individuals, parameter *K* is the prey carrying capacity, *α* is the prey per capita linear growth rate, *A* is the predator attack rate, *H* is the half-saturation prey density, *κ* is the food assimilation efficiency coefficient and *M* is the predator mortality. For the sake of simplicity, below we assume that *D*_{1}=*D*_{2}=*D*.

Let us emphasize that the model (4.1) and (4.2) is well validated; in an appropriate parameter range, it is in good agreement with experimental data (Malchow *et al.* 2008) and therefore can be used for generation of ecologically realistic data.

Note that, since the model (4.1) and (4.2) does not contain any spatial heterogeneity, the test cases considered below apply more to the population distribution on the scale of a single field rather than to a landscape scale where heterogeneity cannot be neglected.

In order to obtain the population distributions, system (4.1) and (4.2) is solved numerically. For convenience, we first introduce dimensionless variables. It is readily seen that the domain length *L* makes a convenient scale for the coordinate *X*, *x*=*X*/*L*, where the new variable *x* is dimensionless. In a similar way, the prey maximum per capita growth rate *α* gives a scale for the time, *t*=*αT*, and the carrying capacity *K* makes a convenient scale for the prey population density, *u*=*U*/*K*. Introducing dimensionless predator density as *v*=*V* *A*/(*αK*), equations (4.1) and (4.2) take the following form at the interval 0<*x*<1:
4.3
and
4.4
where *k*=*κA*/*α*, *m*=*M*/*α*, *h*=*H*/*K* and *d*=*D*/(*L*^{2}*α*) are dimensionless parameters.

Note that, since the principal goal of the pest control is to prevent pest outbreaks, we cannot neglect the reaction terms in equations (4.3) and (4.4) that describe the population growth and predation. In its turn, interaction between reaction and diffusion is known to result in pattern formation; e.g. Malchow *et al.* (2008). The properties of the patterns essentially depend on the type of kinetics of the corresponding non-spatial system, i.e. system (4.3) and (4.4) without the diffusion terms. Here, we focus on the case when the kinetics is oscillatory, which is in agreement with observations on various predator–prey systems (Turchin 2003; Sherratt & Smith 2008). It is readily seen (Malchow *et al.* 2008) that the non-spatial system has a unique positive equilibrium that becomes unstable for *h*<(*k*−*m*)/(*k*+*m*). Once it loses its stability, a limit cycle appears through Hopf bifurcation. The pattern appearing in the spatial system (4.3) and (4.4) will be oscillatory as well. Importantly, the properties of the pattern^{2} then depend on the value of dimensionless diffusivity *d*. In particular, for *d* being on the order of 1 or larger, the solution *u*(*x*,*t*) will be a monotonous function of *x* (figure 3*a*), which means that the local population oscillations are almost synchronized over the entire domain.

However, oscillations at different positions can become de-synchronized for *d*≪1 (see Petrovskii *et al.* 2003). In the latter case, the initial conditions *u*(*x*,0) and *v*(*x*,0) evolve to an ensemble of irregular humps and hollows (figure 3*d*). The smaller the *d* is, the larger the number of humps is. For an intermediate value of *d*, the pattern can consist of just one or a few peaks only; see figure 3*b*,*c*. Also, for a fixed value of *d*, the complexity of the pattern (e.g. the number of humps) tends to increase with time; compare figure 3*a* with figure 3*b* (obtained for *d*=10^{−4}), and figure 3*c* with figure 3*d* (obtained for *d*=10^{−5}). In terms of the ecologically meaningful parameters of the original system (4.1) and (4.2) (see the definition of *d* below equation (4.4)), it means that a slowly diffusing pest population is more likely to form a spatial pattern than a fast one. Similarly, patterns are more likely to be observed in a large field than in a small one.

Available data show that, for common insect pests in Britain, a typical size of spatial heterogeneity in their distribution is between 30 and 50 m (Holland *et al.* 1999; Alexander *et al.* 2005). Therefore, the distributions shown in figure 3*a*,*b* are likely to be observed in an agricultural field with a characteristic size of the order of 100 m, while figure 3*c*,*d* would correspond to a larger field with a size of several hundred metres.

Now we are going to check whether our method of numerical integration on a coarse grid can be applied to the population distributions shown in figure 3. Note that, unlike the test cases considered in the previous section, we do not know the exact value of the integral when the pest population density represented by the function *u*(*x*) is integrated. To overcome this difficulty, we first numerically solve the system (4.3) and (4.4) on a very fine uniform grid *G*_{f} of *N*_{f}=2^{15}≡32 768 grid cells, and consider the result as the ‘exact’ solution to the problem. The corresponding value of the solution integral is then considered as the exact integral computation that will be compared with the value of the integral on a coarse grid. We also need to look at the convergence rate for each of the distributions in order to conclude about the reliability of computations on the initial coarse grid. For this purpose, we generate a sequence of uniformly refined grids where the ‘exact’ solution obtained on the fine grid with *N*_{f} grid subintervals should be available on each grid in the sequence. Hence, we consider a projection of fine grid *N*_{f} onto a uniform coarse grid *G*_{c} of *N*_{c} grid cells. The number *N*_{c} is defined as *N*_{c}=*sN*_{0}, where *N*_{0}=8 is the number of grid cells on the initial grid and the scaling coefficient *s*=2^{m},*m*=0,1,2,…,12. The nodal coordinates , *k*=1,…,*N*_{c}, and , *i*=1,…,*N*_{f}, considered on the coarse and fine grid, respectively, are related as , where *i*=*sk*. The solution at nodes of a coarse grid of *N*_{c} subintervals is then readily available and the integration error
4.5
can be easily computed. In the definition above, the integral *I*_{Nf} computed on the fine grid is considered as the exact integral for a given function *u*(*x*) and *I*_{Nc} is an integral evaluate obtained on a grid with *N*_{c} subintervals.

The convergence history for the solution given by a monotone function (see figure 3*a*) is shown in figure 4*a*. At a first glance, the convergence curve presented in the figure contradicts the asymptotic estimate (2.6). The discussion of the previous section demonstrates that, for a monotone function without any humps, boundary layers or fast-oscillating regions, one can expect better convergence of the GLI method with *p*=5 starting from the initial grid of *N*_{0}=8 nodes. However, it can be seen from the figure that the convergence of the GLI method with *p*=5 is the same as the convergence of Simpson’s rule, unless a very fine grid is generated. This contradiction can be readily resolved. Namely, it appears that the monotone solution *u*(*x*) is close to a cubic polynomial; see figure 5*a*. In this case, the GLI method with *p*=5 is identical to GLI *p*=3 method, as *p*=3 is the highest polynomial degree of the interpolation. Consequently, the convergence of the GLI method is the same as the convergence of Simpson’s rule, unless the function is resolved on a very fine grid.

Meanwhile, it is worth noting here that the integration error for the monotone solution *u*(*x*) is very small even on coarse grids, as we have ||*e*||≈5×10^{−3} already on the initial grid. Such accuracy is much higher than the usual accuracy of ecological data;^{3} therefore, the coarse grid with *N*=8 appears to be sufficient for the problem of ecological monitoring.

The next test case is given by the function *u*(*x*) shown in figure 3*b*. The convergence history for this function is presented in figure 4*b*. This test case clearly demonstrates the advantage of the GLI method with a higher order polynomial used for the interpolation. The conclusion that we make from the analysis of the convergence graph is that, for the GLI method, the initial grid of *N*=8 subintervals appears to be sufficiently refined and results of integration on this grid are reliable as asymptotic error estimates can be applied.

The convergence test for the function of figure 3*c* is presented in figure 4*c*. A somewhat finer grid should now be generated to resolve the humps that are present in the solution. It follows from the consideration of the convergence graphs for the GLI *p*=5 and Simpson’s method that a grid of *N*=32 cells is sufficient to obtain a reliable integral evaluate.

Finally, we consider the multi-peak solution shown in figure 3*d*. As we have already seen in §3, several cycles of grid refinement may be required to resolve solution oscillations. Indeed, it is readily seen from the convergence graph that asymptotic error estimates only begin to work when a grid is refined to *N*=128 subintervals. It is worth mentioning, however, that the GLI method provides a very good accuracy already for *N*=32.

In §3, we have shown that if some *a priori* information is available with regard to the solution properties, a simple heuristic formula (3.5) can be used to estimate the minimum number of grid points on a sufficiently refined grid. Now we are going to validate the estimate (3.5) for more the ecologically realistic solutions shown in figure 3. Consider, for instance, a typical ecological distribution given by a multi-peak pattern (Malchow *et al.* 2008). It was shown in Petrovskii & Malchow (2001) and Petrovskii *et al.* (2003) that a characteristic hump width is
4.6
where Re(*λ*) is the real part of the eigenvalues of the non-spatial system (4.3) and (4.4) linearized in the vicinity of the (unstable) positive equilibrium and *σ* is a numerical coefficient of the order of 1. The value of Re(*λ*) depends on the parameters; however, it appears to be approximately constant in a wide parameter range (Petrovskii & Malchow 2001), so that equation (4.6) can be written as
4.7
where *ω*≈25. For the parameters of figure 3*d*, we thus obtain Δ*x*≈0.08, which is in very good agreement with the properties of the pattern in the middle of the domain. Applying equation (3.5) with *k*=3, we obtain *N*≈38, and this is indeed the number of grid cells when the solution becomes sufficiently resolved; see above.

We want to emphasize, however, that, for all the ecological test cases considered in this section, integration with a reasonable accuracy can already be done on the initial coarse grid of just *N*=8 cells. Even for the complex multi-peak spatial pattern shown in figure 3*d*, integration on the initial coarse grid gives an estimate of the population size with the accuracy *e*≈0.3, which is acceptable for ecological applications. This seems to be a surprising and rather counter-intuitive result.

## 5. Discussion and conclusions

A central problem of ecological monitoring is to obtain an estimate of the pest population size in a given area under conditions when the population density is only known at a few locations. In our paper, we have considered this problem as a problem of numerical integration on a coarse grid. A criterion of grid coarseness has been elaborated where we use the information available on a given grid to conclude about reliability of results of numerical integration obtained for that grid. Correspondingly, a coarse grid is defined as a grid where one cannot apply asymptotic error estimates to evaluate the integration error. Hence, we elaborate an implicit error estimate that involves direct comparison of two methods of integration on a given grid. It is worth noting here that, as the grid coarseness is evaluated in terms of the integration error rather than by the number of grid nodes, its definition depends strongly on an integrand function. It has been demonstrated in the paper that a grid considered as coarse for one integrand function can be a fine grid for another integrand.

The purpose of this study has been twofold. First, we aimed to develop an approach that makes it possible to decide when a given grid becomes sufficiently refined, so that the asymptotical estimates of the error of numerical integration can be applied. That has been done by considering the relative accuracy of Simpson’s method and the GLI method. In practice, asymptotical estimates usually apply to grids with the number of nodes of the order of 100 or more. Once the asymptotical estimates become valid, the error of numerical integrations normally remains within 1 per cent (often much less than that, depending on the particulars of the integrand).

However, such detailed grids are not available in ecological applications. Correspondingly, the second goal of our study was to reveal what can be the actual integration accuracy on coarse grids. Surprisingly and counter-intuitively, we have found that an ecologically acceptable accuracy on the order of 10 per cent or less can be reached on a very coarse grid, even when the integrand is a function with a complicated spatial structure (cf. figure 3*d*). We regard this empirical result as an indirect proof that an estimate of the pest population size can indeed be obtained by means of numerical integration.

We would like to emphasize again that the integration error is, of course, not available in practical applications. Instead, some information about grid coarseness should be used prior to making field measurements in order to give recommendations with regard to the number of detectors or traps to be installed. For this purpose, a heuristic estimate (3.5) can be useful, provided some *a priori* information about the solution structure is available. The latter can be obtained either from the results of previous field studies on the same pest or from mathematical modelling. The test cases considered in the paper demonstrated that the heuristic estimate (3.5) is in good agreement with a more rigorous estimate of grid coarseness based on the convergency history.

Throughout this paper, we have assumed that the traps are positioned at the nodes of a uniform grid, as is usually the case with ecological applications (Holland *et al.* 1999; Ferguson *et al.* 2003). However, our results show that the accuracy of the integration should be significantly higher if the integrand is known at the abscissae of the Gauss–Legendre formula. It therefore sends out a clear message to ecologists: the accuracy of the monitoring can be improved without any need to install additional traps, i.e. without using additional resources, simply by choosing the trap position as prescribed by the Gauss–Legendre formula.

### (a) Ecological summary

A message to ecologists and/or pest-control specialists that immediately follows from our results is that the GLI method of numerical integration can be an effective tool to estimate the pest population size for pest monitoring and control, even when only scarce information is available.

In order to enhance application of our results to real-world ecological problems, our mathematical findings can be formulated as a protocol of ecological monitoring, i.e. as a succession of simple steps that should be taken to obtain a reliable estimate of the population size of a given pest insect.

— Secure the information about the typical scale of spatial heterogeneity in the population distribution of given species. Common pest species are well studied and therefore this information is easily available from the literature (e.g. Holland

*et al.*1999; Ferguson*et al.*2003; Ester & van Rozen 2005).— Apply the estimate (3.5) in order to calculate the minimum number of traps/sampling units. If/where possible, the traps should be positioned in the field according to the Gauss–Legendre formula.

— Once the traps have been exposed and the trap catches obtained, apply a method of numerical integration (preferably the GLI method) to obtain the population size, with the accuracy of integration normally well within 30 per cent, even on a very coarse grid.

### (b) Concluding remarks

Our work leaves a number of open questions. One direction of our future research should be consideration of two-dimensional problems: indeed, the pest insects normally live in a field rather than on a line. A point of particular interest is the effect of the system’s geometry, i.e. when the field shape cannot be approximated as a rectangle. Also, the case when the trap position is arbitrary (and not in the nodes of a rectangular grid) needs to be considered.

The second issue is the effect of measurement errors and/or noise. In terms of ecological monitoring, noise would result in irregular spatial fluctuations in the population density across the domain, in addition to those produced by deterministic factors. Although we have shown in the above that, in principle, our method can be applied to a fast-oscillating solution, the question yet remains whether the situation will be the same when the perturbation is truly stochastic, i.e. spatially aperiodic.

Also, further validation of the GLI method of numerical integration will be needed, both by analytical tools and by considering other test cases, in order to better outline the limits of its effectiveness and capability. In particular, an interesting question is whether our approach remains valid when applied to functions with discontinuities or singularities, even that such functions are not ecologically meaningful.

Finally, we want to mention here that application of the GLI method of numerical integration and the criterion of grid coarseness that we developed in this paper is not restricted to ecological applications. Minimization of the number of measurements or the number of detectors/sampling units is a challenging issue imposed by current budget cuts and restrictions in many practical problems. Hence, a general problem of evaluating the value of a functional where *u* is a continuous function of space, its value being known only in a few points, arises in many real-life applications, e.g. in weather forecasting, and our results may appear to be useful across different disciplines.

## Acknowledgements

S.P. acknowledges support from the UK East Midlands Development Agency through Innovation Fellowship 6667 ‘QuantiFly’. The authors are thankful to Ms Tara Smith who read the manuscript and made several comments aimed at improving the style and grammar.

## Footnotes

↵1 Where

*N*=8 is consistent with the ecological monitoring routine.↵2 At least, for any

*t*not too small, in order to avoid the effect of the initial conditions.↵3 In many ecological studies, acceptable accuracy of estimates of the population size or density is ||

*e*||∼1; e.g. Pascual & Kareiva (1996) and Sherratt & Smith (2008).

- Received January 17, 2010.
- Accepted April 1, 2010.

- © 2010 The Royal Society