Nyström discretization of integrodifference equations: numerical continuation of periodic solutions and Floquet multipliers

Abstract

Integrodifference equations are discrete-time counterparts to reaction-diffusion equations and have various applications in, e.g., theoretical ecology. Their behavior is often illustrated using numerical simulations, which require a spatial discretization. In this paper, we establish that periodic solutions to time-periodic integrodifference equations, their stability and their Floquet spectrum persist under discretization of Nyström-type, which replaces integrals by quadrature or cubature rules. Moreover, it is shown that the convergence rates of the particular integration rules are preserved. By means of a typical model from theoretical ecology, these results are demonstrated in terms of a numerical continuation for periodic solutions and their Floquet multipliers.

1 Introduction

Integrodifference equations are popular and successful models in fields such as theoretical ecology to describe the temporal and spatial evolution of species evolving in non-overlapping generations or which are metered over discrete intervals. Their right-hand sides are typically Hammerstein integral operators (see [14, 15, 19]) and therefore fit into the more general setting of Urysohn difference equations

$$\begin{aligned} u_{t+1}(x)=\int _\Omega f_t(x,y,u_t(y))\,{\mathrm d}y\quad \text {for all }x\in \Omega \end{aligned}$$

(1.1)

with a compact set $\Omega \subset {\mathbb R}^\kappa$ studied in this paper. The values $u_t(x)$ are interpreted as the population density of species at a given point $x\in \Omega$ from the habitat $\Omega$ in the t-th generation. Note that $u_t$ can also be vector-valued when dealing with several interacting species or with stage-structured models [27]. In order to simulate the asymptotic behavior of these recursions (1.1) numerically, it is purposeful to replace the integrals by a quadrature or cubature rule yielding a corresponding Nyström discretization

$$\begin{aligned} u_{t+1}(x)=\sum _{\eta \in \Omega _n}w_\eta f_t(x,\eta ,u_t(\eta )) \end{aligned}$$

(1.2)

involving nonnegative weights $w_\eta$ and nodes $\eta \in \Omega _n$ from a finite grid $\Omega _n\subset \Omega$ (for this see [9] or [19, pp. 112ff]). Here the parameter $n\in {\mathbb N}$ indicates the quality of approximation and increasing n is supposed to yield smaller discretization errors. In contrast to alternative approaches from the numerical analysis of integral equations [9] such as projection and degenerate kernel methods, Nyström discretizations have the advantage to allow an immediate implementation. Given this, the paper at hand is concerned with

• (1) Feasible numerical simulations of integrodifference (1.1) and robust computation and continuation of their periodic solutions including the associate Floquet spectrum providing stability information,

• (2) A profound mathematical justification of these aspects.

We restrict ourselves to time-periodic Urysohn difference (1.1), which is well-motivated from applications in order to describe seasonal influences.

Concerning (1), we provide explicit algorithms (and their implementation) to simulate the forward behavior of solutions to Urysohn difference (1.1), locate and continue their periodic solutions, and to determine their stability.

Concerning (2), it is a natural question to ask how the dynamics of (1.1) and of its Nyström discretization (1.2) are related. When aiming for an answer, one is confronted with the following problem: Even for convergent integration rules the linearizations of (1.2) converge pointwise but not uniformly to linearizations of (1.1) in the limit $n\rightarrow \infty$ (cf. [13, pp. 130–131, Lemma 4.7.6]). Two strategies enable us to circumvent this problem:

• Classically, collective compactness [5] provides a suitable perturbation framework establishing persistence and convergence of fixed points [7, 30] for nonlinear integral operators defined on the continuous functions over $\Omega$.

• More recently, retreating to the state space of Hölder continuous functions over $\Omega$ guaranteed convergence results in the operator norm [22] for which, as illustrated in [21], a more conventional perturbation theory applies. Nonetheless, [21] requires rather involved assumptions to ensure that the difference (1.1) and (1.2) are well-defined when acting between spaces of Hölder functions.

In this paper, we avoid such complications and immediately work with the natural state space of continuous functions. Our basic Section  establishes that the right-hand sides of (1.1) and (1.2) are well-defined, have appropriate properties, and converge pointwise to each other as $n\rightarrow \infty$, provided the particular integration rule is convergent. The stability properties of periodic solutions are based on the linearizations of (1.1) resp. (1.2) and their Floquet spectrum. Relying on preparations in [6, 8], we therefore provide a corresponding perturbation theory for the Floquet multipliers in Section , which preserves the convergence order of the integration rule used for (1.2). In Section  we prove that generic periodic solutions to (1.1) persist under Nyström discretization preserving both stability properties when passing to (1.2), as well as the convergence order of the integration method. Here, our arguments are based on tools from [7, 30]. The actual convergence rates that can be realized in this process are derived in Section . Our concluding Section  contains some numerical tests, but also an illustration of our techniques by means of a concrete ecological model (a generalized spatial Beverton–Holt equation equipped with a Laplace kernel depending on one parameter). Using pseudo-arclength continuation, we compute a (partial) bifurcation diagram illustrating the number of periodic solutions depending on the parameter value. For the convenience of the reader, a short appendix closes our elaborations.

Notation

Given a finite set N, we write $\#N$ for its cardinality. The positive integers are denoted by ${\mathbb N}$, while ${\mathbb N}_0$ are the nonnegative integers. For the nonnegative reals, we write ${\mathbb R}_+:=[0,\infty )$ and the complex unit circle is ${\mathbb S}^1:=\left\{ z\in \mathbb {C}:\,|z|=1\right\}$. Norms on finite-dimensional spaces are denoted by $\left| \cdot \right|$. The Cartesian product $X\times Y$ of normed spaces X, Y is equipped with the product norm

$$\begin{aligned} \left\| (x,y)\right\| :=\max \bigl \{\left\| x\right\| _X,\left\| y\right\| _Y\bigr \} \end{aligned}$$

(1.3)

and we proceed accordingly on products of more than two spaces. For the open and closed balls in X with radius $r>0$ and center $x\in X$, we respectively write

$$\begin{aligned} B_r(x)&:=\left\{ y\in X:\,\left\| y-x\right\| <r\right\} ,&\bar{B}_r(x)&:=\left\{ y\in X:\,\left\| y-x\right\| \le r\right\} . \end{aligned}$$

The set of all bounded k-linear maps from the k-fold Cartesian product $X^k$ to Y is denoted by $L_k(X,Y)$; in particular $L(X,Y):=L_1(X,Y)$, $L(X):=L(X,X)$, $L_0(X,Y):=Y$ and GL(X, Y) are the invertible elements in L(X, Y). In this context, $\text {id}_X$ is the identity mapping on X, N(T) stands for the kernel, and R(T) for the range of $T\in L(X,Y)$. Furthermore, $\sigma (S)$ abbreviates the spectrum and $\sigma _p(S)$ the point spectrum of (the complexification to) $S\in L(X)$.

Given a nonempty set $A\subseteq X$, ${{\,\textrm{dist}\,}}_A(x):=\inf _{a\in A}\left\| x-a\right\|$ is the distance of a point $x\in X$ from A and $B_r(A):=\left\{ x\in X:\,{{\,\textrm{dist}\,}}_A(x)<r\right\}$ stands for the open r-neighborhood of A. We denote any subset $\mathcal {A}\subseteq {\mathbb Z}\times X$ as nonautonomous set having the fibers $\mathcal {A}(t):=\left\{ x\in X:\,(t,x)\in \mathcal {A}\right\}$, $t\in {\mathbb Z}$. In particular, the nonautonomous set

$$ \mathcal {B}_r(\phi ):=\left\{ (t,u)\in {\mathbb Z}\times X:\,\left\| u-\phi _t\right\| <r\right\} $$

is the open r-neighborhood of a sequence $\phi =(\phi _t)_{t\in {\mathbb Z}}$ in X.

If $X_1,\ldots ,X_n$ are Banach spaces and a function $f:A\subseteq X_1\times \ldots \times X_n\rightarrow Y$ depends differentiably on the ith variable, $1\le i\le n$, then the corresponding partial derivative is denoted by $D_if:A\rightarrow L(X_i,Y)$; for higher order partial derivatives, we write $D_i^kf$, $k\in {\mathbb N}_0$.

A family ${\mathscr {F}}$ of maps $F:X\rightarrow Y$ is said to be collectively compact, if the union $\bigcup _{F\in {\mathscr {F}}}F(B)\subseteq Y$ is relatively compact for each bounded $B\subset X$ and ${\mathscr {T}}\subset L(X,Y)$ is called collectively compact, if the set $\bigcup _{T\in {\mathscr {T}}}TB_1(0)\subseteq Y$ is relatively compact.

Throughout the paper, we assume that the habitat $\Omega \subset {\mathbb R}^\kappa$ is nonempty and compact. The set $C(\Omega ,{\mathbb R}^d)$ of all continuous functions $u:\Omega \rightarrow {\mathbb R}^d$ is a Banach space w.r.t. the norm $\left\| u\right\| _\infty :=\sup _{x\in \Omega }\left| u(x)\right|$ and will be abbreviated as ${C_{d}}:=C(\Omega ,{\mathbb R}^d)$. Finally, we say the function u is Hölder continuous with exponent $\alpha \in (0,1]$, if

$$ \left[ u\right] _\alpha :=\sup _{x,\bar{x}\in \Omega }\frac{\left| u(x)-u(\bar{x})\right| }{\left| x-\bar{x}\right| ^\alpha }<\infty $$

and write $C^\alpha (\Omega ,{\mathbb R}^d)$ for the set of all such functions.

2 Integrodifference equations and Nyström discretizations

Let us first introduce some basics on Nyström discretizations of integrodifference equations over a compact habitat $\Omega \subset {\mathbb R}^\kappa$. A quadrature ($\kappa =1$) or cubature rule ($\kappa >1$) is a sequence of mappings

determined by grids $\Omega _n\subset \Omega$ consisting of finitely many nodes $\eta \in \Omega _n$ and weights $w_\eta >0$; note that the dependence of $w_\eta$ on $n\in {\mathbb N}$ is suppressed here. For concrete examples, we refer to [11, pp. 361ff, Chap. 15–16] and Appendix B. One denotes $Q^n$ as

• stable, provided the weights satisfy

  $$\begin{aligned} W&:=\sup _{n\in {\mathbb N}}W_n<\infty ,&W_n&:=\sum _{\eta \in \Omega _n}w_\eta , \end{aligned}$$

  (2.1)

• convergent, if $\lim _{n\rightarrow \infty }Q^nu=\int _\Omega u(y)\,{\mathrm d}y$ holds for all $u\in {C_{d}}$. Moreover, $Q^n$ has the convergence order $\alpha >0$, if $\left| \int _\Omega u(y)\,{\mathrm d}y-Q^nu\right| \le \tfrac{C}{n^\alpha }$ for all $n\in {\mathbb N}$ with a constant $C\ge 0$ typically depending on the function u or its derivatives.

Because of [13, p. 20, Thm. 1.4.17], convergence implies stability.

A natural way to simulate an Urysohn difference (1.1) is to substitute the integral by a convergent integration rule $Q^n$ yielding a Nyström discretization (1.2). Aiming for a unified notation capturing both (1.1) and (1.2) simultaneously, we consider families of difference equations

with the right-hand sides

$$\begin{aligned} {\mathscr {F}}_t^n:{C_{d}}\rightarrow {C_{d}}, \qquad {\mathscr {F}}_t^n(u) := {\left\{ \begin{array}{ll} \int _\Omega f_t(\cdot ,y,u(y))\,{\mathrm d}y,\qquad \qquad \, n=0,\\ \sum _{\eta \in \Omega _n}w_\eta f_t(\cdot ,\eta ,u(\eta )),\qquad n\in {\mathbb N}. \end{array}\right. } \end{aligned}$$

(2.2)

In particular, $\left\{ {\mathscr {F}}_t^n\right\} _{n\in {\mathbb N}}$ is called a Nyström method corresponding to the integration rule $Q^n$, and the functions ${\mathscr {F}}_t^n(u)$, $n\in {\mathbb N}$, are the Nyström interpolants. For the sake of well-defined and smooth mappings ${\mathscr {F}}_t^n$, we require several standing assumptions on the kernel functions $f_t:\Omega \times \Omega \times {\mathbb R}^d\rightarrow {\mathbb R}^d$ in (2.2):

Hypothesis 2.1

Let $m\in {\mathbb N}_0$. Suppose there exists a $\theta _0\in {\mathbb N}$ such that

$$\begin{aligned} f_t=f_{t+\theta _0}:\Omega \times \Omega \times {\mathbb R}^d\rightarrow {\mathbb R}^d\quad \text {for all }t\in {\mathbb Z} \end{aligned}$$

(2.3)

and that the partial derivatives $D_3^kf_t:\Omega \times \Omega \times {\mathbb R}^d\rightarrow L_k({\mathbb R}^d,{\mathbb R}^d)$ exist as continuous functions for all $0\le t<\theta _0$ and $0\le k\le m$.

Remark 2.2

Let $(\phi _t^0)_{t\in {\mathbb Z}}$ be a sequence in ${C_{d}}$. Then, the above assumptions guarantee that there exist functions $b_t^k:{\mathbb R}_+\times \Omega \rightarrow {\mathbb R}_+$ such that

$$\begin{aligned} b_t^k(r,y) := \sup _{x\in \Omega }\sup _{z\in \bar{B}_r(\phi _t^0(y))} \left| D_3^kf_t(x,y,z)\right| <\infty \quad \text {for all }y\in \Omega . \end{aligned}$$

(2.4)

Proposition 2.3

(properties of ${\mathscr {F}}_t^0$) If Hypothesis  holds with $m\in {\mathbb N}_0$, then the Urysohn operators ${\mathscr {F}}_t^0:{C_{d}}\rightarrow {C_{d}}$ are well-defined, completely continuous and satisfy ${\mathscr {F}}_t^0={\mathscr {F}}_{t+\theta _0}^0$ for all $t\in {\mathbb Z}$.

Proof

This results from [23, proof of Prop. 2.6 and Cor. 2.7], where (2.2) and (2.3) yield the periodicity ${\mathscr {F}}_t^0={\mathscr {F}}_{t+\theta _0}^0$ for all $t\in {\mathbb Z}$.

A forward solution to $\Delta ^n$, $n\in {\mathbb N}_0$, is a sequence $(\phi _t)_{\tau \le t}$ satisfying the solution identity $\phi _{t+1}={\mathscr {F}}_t^n(\phi _t)$ for all $\tau \le t$ with some initial time $\tau \in {\mathbb Z}$. An entire solution is a sequence $(\phi _t)_{t\in {\mathbb Z}}$ in ${C_{d}}$ satisfying $\phi _{t+1}\equiv {\mathscr {F}}_t^n(\phi _t)$ on ${\mathbb Z}$ and $\theta$-periodic solutions additionally satisfy $\phi _t\equiv \phi _{t+\theta }$ on ${\mathbb Z}$ with some $\theta \in {\mathbb N}$.

Proposition 2.4

(smoothness of ${\mathscr {F}}_t^0$) If Hypothesis  holds with $m\in {\mathbb N}$, then the Urysohn operators ${\mathscr {F}}_t^0:{C_{d}}\rightarrow {C_{d}}$ are m-times continuously differentiable with Fréchet derivatives

$$\begin{aligned} D^k{\mathscr {F}}_t^0(u)v_1\cdots v_k = \int _\Omega D_3^kf_t(\cdot ,y,u(y))v_1(y)\cdots v_k(y)\,{\mathrm d}y \end{aligned}$$

for all $t\in {\mathbb Z}$, $u,v_1,\ldots ,v_k\in {C_{d}}$ and $1\le k\le m$. In particular, the Fréchet derivative $D{\mathscr {F}}_t^0(u)\in L({C_{d}})$ is compact.

Proof

Let $t\in {\mathbb Z}$. The claimed differentiability is shown in [23, proof of Prop. 2.15]. Because the right-hand side ${\mathscr {F}}_t^0$ is completely continuous, due to Proposition , the compactness of the operator $D{\mathscr {F}}_t^0(u)$ results from [20, p. 89, Prop. 6.5].

Next, let us establish that analytical properties of the right-hand sides ${\mathscr {F}}_t^0$ in (1.1) extend to their Nyström discretizations ${\mathscr {F}}_t^n$, $n\in {\mathbb N}$:

Proposition 2.5

(properties of ${\mathscr {F}}_t^n$) If Hypothesis  holds with $m\in {\mathbb N}_0$, then the Nyström discretizations ${\mathscr {F}}_t^n:{C_{d}}\rightarrow {C_{d}}$ are well-defined and completely continuous for each $n\in {\mathbb N}$ and satisfy ${\mathscr {F}}_t^n={\mathscr {F}}_{t+\theta _0}^n:{C_{d}}\rightarrow {C_{d}}$ for all $t\in {\mathbb Z}$. If the integration rule $Q^n$ is stable, then

• (a) $\left\{ {\mathscr {F}}_t^n(u)\right\} _{n\in {\mathbb N}}$ is equicontinuous for all $u\in {C_{d}}$,

• (b) $\left\{ {\mathscr {F}}_t^n\right\} _{n\in {\mathbb N}}$ is collectively compact.

Proof

Let $t\in {\mathbb Z}$. The grids $\Omega _n$ are a family of compact, discrete subsets of $\Omega$ and we equip them with the measure $\mu (\Omega _n):=\sum _{\eta \in \Omega _n}w_\eta$. Then, due to [23, Exams. 2.1 and 2.2], the abstract measure-theoretical integral considered in [23] becomes

$$\begin{aligned} \int _{\Omega _n}f_t(x,y,u(y))\,{\mathrm d}\mu (y) = \sum _{\eta \in \Omega _n}w_\eta f_t(x,\eta ,u(\eta ))\quad \text {for all }x\in \Omega \end{aligned}$$

(2.5)

and yields the discrete integral operators ${\mathscr {F}}_t^n$ given by (2.2) for $n\in {\mathbb N}$. Consequently, well-definedness and complete continuity result as in the proof of Proposition . Finally, the relations (2.2) and (2.3) imply that ${\mathscr {F}}_t^n$ are $\theta _0$-periodic in t.

(a) Let $\varepsilon >0$ and $u\in {C_{d}}$ with $r:=\left\| u\right\| _\infty$. Since the continuous kernel function $f_t:\Omega \times \Omega \times {\mathbb R}^d\rightarrow {\mathbb R}^d$ is uniformly continuous on every compact set $\Omega \times \Omega \times \bar{B}_r(0)$, there exists a $\delta >0$ such that $\left| x-\bar{x}\right| <\delta$ implies

$$ \left| f_t(x,y,z)-f_t(\bar{x},y,z)\right| \le \tfrac{\varepsilon }{2W} \quad \text {for all }x,\bar{x},y\in \Omega ,\,z\in \bar{B}_r(0). $$

For each $n\in {\mathbb N}$ this readily leads to

$$\begin{aligned} \left| {\mathscr {F}}_t^n(u)(x)-{\mathscr {F}}_t^n(u)(\bar{x})\right|&\overset{(2.2)}{\le } \sum _{\eta \in \Omega _n}w_\eta \left| f_t(x,\eta ,u(\eta ))-f_t(\bar{x},\eta ,u(\eta ))\right| \overset{(2.1)}{\le } \tfrac{\varepsilon }{2} \end{aligned}$$

and passing to the least upper bound over $x\in \Omega$ yields $\left\| {\mathscr {F}}_t^n(u)-{\mathscr {F}}_t^n(\bar{u})\right\| _\infty <\varepsilon$, that is, the family $\left\{ {\mathscr {F}}_t^n(u)\right\} _{n\in {\mathbb N}}\subset {C_{d}}$ is equicontinuous.

(b) Given a bounded subset $B\subset {C_{d}}$, we show that $\bigcup _{n\in {\mathbb N}}{\mathscr {F}}_t^n(B)$ is relatively compact by means of the Arzelà-Ascoli theorem [18, p. 58, Cor. 3.2]:

First, since B is bounded, there exists a $r>0$ such that $\left| u(y)\right| \le \left\| u\right\| _\infty \le r$ for all $y\in \Omega$ and $u\in B$. Moreover, the continuous function $f_t:\Omega \times \Omega \times {\mathbb R}^d\rightarrow {\mathbb R}^d$ is bounded on the compact product $\Omega \times \Omega \times \bar{B}_r(0)$ by some $K\ge 0$. This implies

$$ \left| {\mathscr {F}}_t^n(u)(x)\right| \overset{(2.2)}{\le } \sum _{\eta \in \Omega _n}w_\eta \left| f_t(x,y,u(y))\right| \le K\sum _{\eta \in \Omega _n}w_\eta \overset{(2.1)}{\le } KW \quad \text {for all }x\in \Omega $$

and passing to the supremum over $x\in \Omega$ yields $\left\| {\mathscr {F}}_t^n(u)\right\| _\infty \le KW$ for all $n\in {\mathbb N}$ and $u\in B$. This means that $\bigcup _{n\in {\mathbb N}}{\mathscr {F}}_t^n(B)$ is bounded.

Second, let $\varepsilon >0$ and choose $u\in \bigcup _{n\in {\mathbb N}}{\mathscr {F}}_t^n(B)$. Hence, there exist $n\in {\mathbb N}$, $v\in B$ with $u={\mathscr {F}}_t^n(v)$ and as in (a), we stablish that $\left| x-\bar{x}\right| <\delta$ implies

$$ \left| u(x)-u(\bar{x})\right| = \left| {\mathscr {F}}_t^n(v)(x)-{\mathscr {F}}_t^n(v)(\bar{x})\right| \le \tfrac{\varepsilon }{2} \quad \text {for all }x,\bar{x}\in \Omega , $$

which shows that $\bigcup _{n\in {\mathbb N}}{\mathscr {F}}_t^n(B)$ is equicontinuous.

Let us next address differentiability properties. A Nyström method $\left\{ {\mathscr {F}}_t^n\right\} _{n\in {\mathbb N}}$ is called equidifferentiable in a point $u\in {C_{d}}$, if each ${\mathscr {F}}_t^n$ is continuously differentiable in u with derivative $D{\mathscr {F}}_t^n(u)\in L({C_{d}})$ and

$$ \lim _{h\rightarrow 0}\tfrac{1}{\left\| h\right\| _\infty }\left\| {\mathscr {F}}_t^n(u+h)-{\mathscr {F}}_t^n(u)-D{\mathscr {F}}_t^n(u)h\right\| _\infty = 0\quad \text {for all }t\in {\mathbb Z}$$

holds uniformly in $n\in {\mathbb N}$.

Proposition 2.6

(smoothness of ${\mathscr {F}}_t^n$) If Hypothesis  holds with $m\in {\mathbb N}$, then the Nyström discretizations ${\mathscr {F}}_t^n:{C_{d}}\rightarrow {C_{d}}$ are m-times continuously differentiable for each $n\in {\mathbb N}$ with Fréchet derivatives

$$\begin{aligned} D^k{\mathscr {F}}_t^n(u)v_1\cdots v_k = \sum _{\eta \in \Omega _n}w_\eta D_3^kf_t(\cdot ,\eta ,u(\eta )) v_1(\eta )\cdots v_k(\eta ) \end{aligned}$$

(2.6)

for all $t\in {\mathbb Z}$, $u,v_1,\ldots ,v_k\in {C_{d}}$ and $1\le k\le m$. If the integration rule $Q^n$ is stable, then

• (a) the Nyström method $\left\{ {\mathscr {F}}_t^n\right\} _{n\in {\mathbb N}}$ is equidifferentiable on ${C_{d}}$ and

• (b) $\left\{ D{\mathscr {F}}_t^n(u)\right\} _{n\in {\mathbb N}}$ is equicontinuous in each $u\in {C_{d}}$ and every $D{\mathscr {F}}_t^n:{C_{d}}\rightarrow L({C_{d}})$ is uniformly continuous on bounded sets uniformly in $n\in {\mathbb N}$.

Proof

Let $t\in {\mathbb Z}$ and $u\in {C_{d}}$. The smoothness assertions on ${\mathscr {F}}_t^n$ are shown by applying [23, proof of Prop. 2.15] to (2.5). Assume that $Q^n$ is stable from now on.

(a) In the present situation the remainder term [23, proof of Prop. 2.15] becomes

$$\begin{aligned} r_0(h)&=\sup _{\vartheta \in [0,1]}\bigg \Vert \sum _{\eta \in \Omega _n}w_\eta \left[ D_3f_t(\cdot ,\eta ,(u+\vartheta h)(\eta ))-D_3f_t(\cdot ,\eta ,u(\eta ))\right] \bigg \Vert _\infty \end{aligned}$$

for functions $h\in {C_{d}}$. It follows using (2.1) that $\lim _{h\rightarrow \infty }r_0(h)=0$ holds uniformly in $n\in {\mathbb N}$. This yields the claimed equidifferentiability.

(b) Mimicking the proof of [23, Lemma 2.14], where $(u_l)_{l\in {\mathbb N}}$ is a sequence converging to $u\in {C_{d}}$ and $v\in \bar{B}_1(0)$, as in inequality [23, (2.14)] one obtains that there exists a $L_1\in {\mathbb N}$ such that $l\ge L_1$ implies

$$ \left\| D{\mathscr {F}}_t^n(u_l)v-D{\mathscr {F}}_t^n(u)v\right\| _\infty \le \sum _{\eta \in \Omega _n}w_\eta \frac{\varepsilon }{2W} < \varepsilon \quad \text {for all }n\in {\mathbb N}. $$

Passing to the supremum over $v\in \bar{B}_1(0)$ in the above estimate establishes the claimed equicontinuity of $\left\{ D{\mathscr {F}}_t^n(u)\right\} _{n\in {\mathbb N}}$. Moreover, the claimed uniform continuity on bounded subsets holds.

We finally obtain at least pointwise convergence in the local discretization error:

Proposition 2.7

(convergence of ${\mathscr {F}}_t^n$) If Hypothesis  holds with $m\in {\mathbb N}_0$ and the integration rule $Q^n$ is convergent, then

$$\begin{aligned} \lim _{n\rightarrow \infty }\left\| {\mathscr {F}}_t^n(u)-{\mathscr {F}}_t^0(u)\right\| _\infty =0\quad \text {for all }t\in {\mathbb Z},\,u\in {C_{d}}. \end{aligned}$$

(2.7)

In Section , we provide explicit convergence rates in the limit relation (2.7).

Proof

Let $t\in {\mathbb Z}$, $u\in {C_{d}}$ be fixed and define the functions $g_n:={\mathscr {F}}_t^n(u)$ for $n\in {\mathbb N}_0$. In order to apply Lemma , we observe that $f_t(x,\cdot ,u(\cdot )):\Omega \rightarrow {\mathbb R}^d$ is continuous and the convergence of $Q^n$ yields pointwise convergence

$$ \lim _{n\rightarrow \infty }g_n(x)=g_0(x)\quad \text {for all }x\in \Omega . $$

According to Proposition (a) the family $\left\{ g_n\right\} _{n\in {\mathbb N}}$ is equicontinuous. In combination, Lemma  guarantees uniform convergence, i.e., (2.7) holds.

The above Proposition  suggests that the Nyström discretizations $\Delta ^n$ capture the solution behavior of integrodifference equations $(\Delta ^0)$ at least over finite time-intervals, provided that $n\in {\mathbb N}$ is sufficiently large. Now, in order to simulate forward solutions

$$ \varphi ^0(t;\tau ,u_\tau ):= {\left\{ \begin{array}{ll} u_\tau ,&{}t=\tau ,\\ {\mathscr {F}}_{t-1}^0\circ \cdots \circ {\mathscr {F}}_\tau ^0(u_\tau ),&{}\tau <t \end{array}\right. } $$

to an Urysohn difference equation $(\Delta ^0)$ starting at time $\tau \in {\mathbb Z}$ in $u_\tau \in {C_{d}}$ by means of the discretizations $\Delta ^n$, $n\in {\mathbb N}$, first choose an integration rule $Q^n$ with weights $w_\eta >0$ and nodes $\eta \in \Omega _n$ from a grid $\Omega _n=\left\{ \eta _1,\ldots ,\eta _{q_n}\right\} \subseteq \Omega$ (cf. Appendix B for details). Given this, numerical implementations require us to represent continuous functions $u:\Omega \rightarrow {\mathbb R}^d$ as finite-dimensional vectors in ${\mathbb R}^{dq_n}$. While one can proceed accordingly in arbitrary spatial dimensions $\kappa >1$, on habitats $\Omega \subset {\mathbb R}$, we provide a detailed description: For $\kappa =1$, we represent numerical approximations of u as matrix

$$ U:=(u(\eta _1),\ldots ,u(\eta _{q_n})) = \begin{pmatrix} u_1(\eta _1) &{} \ldots &{} u_1(\eta _{q_n})\\ \vdots &{} &{} \vdots \\ u_d(\eta _1) &{} \ldots &{} u_d(\eta _{q_n}) \end{pmatrix} \in {\mathbb R}^{d\times q_n}\cong {\mathbb R}^{dq_n} $$

by means of the isomorphism $T:{\mathbb R}^{d\times q_n}\rightarrow {\mathbb R}^{dq_n}$,

$$ TU:=(u_1(\eta _1),\ldots ,u_d(\eta _1),u_1(\eta _2),\ldots ,u_d(\eta _2),\ldots ,u_1(\eta _{q_n}),\ldots ,u_d(\eta _{q_n}))^T. $$

One obtains right-hand sides ${\mathscr {F}}_t^n(u)=\sum _{j=1}^{q_n}w_{\eta _j}f_t(\cdot ,\eta _j,\upsilon ^j)$, where

$$\begin{aligned} \upsilon&= \begin{pmatrix} \upsilon ^1\\ \vdots \\ \upsilon ^{q_n} \end{pmatrix}\in {\mathbb R}^{dq_n},&\upsilon ^j&=u(\eta _j)\in {\mathbb R}^d\quad \text {for all }1\le j\le q_n \end{aligned}$$

denotes the vector representation of u, then the mapping

$$\begin{aligned} F_t^n:{\mathbb R}^{dq_n}&\rightarrow {\mathbb R}^{dq_n},&F_t^n(\upsilon )&:= \sum _{j=1}^{q_n} w_{\eta _j} \begin{pmatrix} f_t(\eta _1,\eta _j,\upsilon ^j)\\ \vdots \\ f_t(\eta _{q_n},\eta _j,\upsilon ^j) \end{pmatrix} \end{aligned}$$

(2.8)

reduces our problem to the $dq_n$-dimensional difference equation

$$\begin{aligned} \upsilon _{t+1}=F_t^n(\upsilon _t) \end{aligned}$$

(2.9)

as follows: With the initial vector $\upsilon _\tau \in {\mathbb R}^{dq_n}$ obtained from an initial function $u_\tau$, the forward solution $\upsilon _{t-1}$ to $\Delta ^n$ consequently results in the Nyström interpolant

$$ \varphi ^n(t;\tau ,u_\tau ):= {\left\{ \begin{array}{ll} u_\tau ,&{}t=\tau ,\\ \sum _{j=1}^{q_n}w_{\eta _j}f_{t-1}(\cdot ,\eta _j,\upsilon _{t-1}),&{}\tau <t \end{array}\right. } $$

serving as approximation of the actual forward solution $\varphi ^0(t;\tau ,u_\tau )\in {C_{d}}$. A pseudo-code of this procedure can be found in Algorithm 1:

Algorithm 1

Simulation of forward solutions to Urysohn difference (1.1)

Remark 2.8

(Total population) When using integrodifference equations $(\Delta ^0)$ as models in theoretical ecology [14, 15, 19], the iterates $u_t=\varphi ^0(t;\tau ,u_\tau )$ describe population densities over a habitat $\Omega$ and

$$ \overline{u_t}:=\int _\Omega u_t(y)\,{\mathrm d}y\in {\mathbb R}^d $$

can be interpreted as the total population of the habitat in generation t. Consequently,

$$ \overline{\varphi ^n(t;\tau ,u_\tau )}=\sum _{j=1}^{q_n}w_{\eta _j}\varphi ^n(t;\tau ,u_\tau )(\eta _j)\in {\mathbb R}^d $$

serve as numerical approximations for the total population sizes.

Remark 2.9

(FFT method) Rather than relying on integration formulas, another popular approach to approximately evaluate the right-hand sides in $(\Delta ^0)$ is the Fast Fourier Transformation (abbreviated as FFT, [25] and [19, pp. 106ff, Sect. 8.2]). In contrast to Nyström methods, such approximations are restricted to Hammerstein integrodifference equations of convolution type $u_{t+1}(x)=\int _\Omega k_t(x-y)g_t(u_t(y))\,{\mathrm d}y$ for all $x\in \Omega$. As pointed out in [19, pp. 105ff], in order to arrive at a similar accuracy, Nyström methods require less nodes but tend to be slower than FFT methods.

3 Linear integrodifference equations and Nyström discretization

This section focusses on linear integrodifference equations

$$ u_{t+1}(x)=\int _\Omega k_t(x,y)u_t(y)\,{\mathrm d}y\quad \text {for all }x\in \Omega , $$

where the kernels $k_t:\Omega \times \Omega \rightarrow L({\mathbb R}^d)$ are assumed to be continuous and to satisfy the periodicity assumption $k_{t+\theta }=k_t$ for all $t\in {\mathbb Z}$ and some period $\theta \in {\mathbb N}$.

After an integration rule $Q^n$ is chosen, we once again aim for a unified notation capturing both linear integral operators, as well as their Nyström discretizations and consider families of linear difference equations

with coefficient operators

$$\begin{aligned} {\mathscr {K}}_t^n&\in L({C_{d}}),&{\mathscr {K}}_t^nv&:= {\left\{ \begin{array}{ll} \int _\Omega k_t(\cdot ,y)v(y)\,{\mathrm d}y,&{}n=0,\\ \sum _{\eta \in \Omega _n}w_\eta k_t(\cdot ,\eta )v(\eta ),&{}n\in {\mathbb N}, \end{array}\right. } \end{aligned}$$

having the periodicity property ${\mathscr {K}}_t^n={\mathscr {K}}_{t+\theta }^n$ for all $t\in {\mathbb Z}$. Note that these linear operators are compact (cf. [13, p. 45, Thm. 3.2.5]).

The transition operators $\Phi ^n(t,s)\in L({C_{d}})$ of linear difference $L^n$ are defined as

$$\begin{aligned} \Phi ^n(t,s) := {\left\{ \begin{array}{ll} \text {id}_{{C_{d}}},&{}t=s,\\ {\mathscr {K}}_{t-1}^n\cdots {\mathscr {K}}_s^n,&{}s<t \end{array}\right. } \quad \text {for all }n\in {\mathbb N}_0. \end{aligned}$$

(3.1)

Given this, we introduce the compact period operators

$$ \Xi _\theta ^n:=\Phi ^n(\theta ,0)\in L({C_{d}}) $$

associated with $L^n$. Hence, the nonzero spectral points of $\Xi _0^n$ are eigenvalues. We denote them as Floquet multipliers and $\sigma _p(\Xi _\theta ^n)$ as Floquet spectrum of $L^n$. The minimal $\nu (\lambda )\in {\mathbb N}$ satisfying $N(\lambda \text {id}_{{C_{d}}}-\Xi _\theta ^n)^{\nu (\lambda )}=N(\lambda \text {id}_{{C_{d}}}-\Xi _\theta ^n)^{\nu (\lambda )+1}$ is the index and the dimension of the generalized eigenspace

$$ {{\,\textrm{Eig}\,}}_\lambda \Xi _\theta ^n:=N(\lambda \text {id}_{{C_{d}}}-\Xi _\theta ^n)^{\nu (\lambda )} $$

the algebraic multiplicity $a(\lambda ):=\dim {{\,\textrm{Eig}\,}}_\lambda \Xi _\theta ^n$ of an eigenvalue $\lambda \in \sigma _p(\Xi _\theta ^n)$.

We denote a linear difference $L^n$ as weakly hyperbolic, if $1\not \in \sigma (\Xi _\theta ^n)$, and as hyperbolic, if ${\mathbb S}^1\cap \sigma (\Xi _\theta ^n)=\emptyset$ holds. In the hyperbolic situation, the Floquet spectrum can be represented as $\sigma (\Xi _\tau ^n)=\sigma _+\dot{\cup }\sigma _-$ with spectral sets

$$\begin{aligned} \sigma _+&\subseteq B_1(0),&\sigma _-&\subseteq \mathbb {C}\setminus \bar{B}_1(0). \end{aligned}$$

The unstable spectral set $\sigma _-$ consists of finitely many Floquet multipliers with finite algebraic multiplicity. In particular, there exists a $\rho >0$ such that the neighborhoods $\bar{B}_\rho (\lambda )$, $\lambda \in \sigma _-$, are pairwise disjoint. This allows us to define the spectral projections

$$ P_\lambda ^n:=\frac{1}{2\pi i}\int _{\partial B_\rho (\lambda )}(z\text {id}_{{C_{d}}}-\Xi _\theta ^n)^{-1}\,{\mathrm d}z \quad \text {for all }\lambda \in \sigma _- $$

and to obtain $R(P_\lambda ^n)={{\,\textrm{Eig}\,}}_\lambda \Xi _\theta ^n$. Based on this, we introduce the fibers

$$\begin{aligned} \mathcal {V}_-^n(t):=\Phi ^n(t,0)\bigoplus _{\lambda \in \sigma _-}R(P_\lambda ^n), \end{aligned}$$

first for integers $t\ge 0$ and then by $\theta$-periodic continuation for $t\in {\mathbb Z}$. This yields invariant $\theta$-periodic nonautonomous sets $\mathcal {V}_-^n\subseteq {\mathbb Z}\times {C_{d}}$ (the unstable vector bundle) of $L^n$, i.e., $\mathcal {V}_-^n(t)=\mathcal {V}_-^n(t+\theta )$ for all $t\in {\mathbb Z}$. The fibers $\mathcal {V}_-^n(t)$ have constant finite dimension, which we denote as Morse index of $L^n$ and equals the sum $\sum _{\lambda \in \sigma _-}a(\lambda )$.

Using this terminology, let us study the behavior of the Floquet spectrum $\sigma _p(\Xi _\theta ^0)$ under Nyström discretization of $(L^0)$. For a Floquet multiplier $\lambda \ne 0$ of $(L^0)$, choose $\rho >0$ sufficiently small such that $B_\rho (\lambda )\cap \sigma _p(\Xi _\theta ^0)=\left\{ \lambda \right\}$. Given this, we define the set

$$ \Lambda _\rho ^n(\lambda ):=B_\rho (\lambda )\cap \sigma _p(\Xi _\theta ^n) $$

of Floquet multipliers to the Nyström discretizations $L^n$ being $\rho$-close to $\lambda$.

Proposition 3.1

(Floquet multipliers of $L^n$ under Nyström discretization) Let $Q^n$ be convergent. If $\lambda \ne 0$ is a Floquet multiplier of $(L^0)$ with index $\nu (\lambda )$ and $\left\{ e_1,\ldots ,e_{a(\lambda )}\right\}$ denotes a basis of ${{\,\textrm{Eig}\,}}_{\lambda }\Xi _\theta ^0$, then there exist $C\ge 0$ and an integer $N\in {\mathbb N}$ such that the following holds true for all $n\ge N$:

• (a) $\max _{\mu \in \Lambda _\rho ^n(\lambda )}\left| \mu -\lambda \right| \le C\max _{i=1}^{a(\lambda )}\root \nu (\lambda ) \of {\left\| [\Xi _\theta ^n-\Xi _\theta ^0]e_i\right\| _\infty }$,

• (b) if $(L^0)$ is (weakly) hyperbolic, then also $L^n$ is (weakly) hyperbolic,

• (c) $\dim \mathcal {V}_-^n=\dim \mathcal {V}_-^0$.

Proof

As preparation, let us embed the operators $\Xi _\theta ^n\in L({C_{d}})$, $n\in {\mathbb N}_0$, into the perturbation setting of [5, 6, 8] needed here. Applying Proposition (b) with the kernel function $f_t(x,y,z):=k_t(x,y)z$ shows that the family $\left\{ {\mathscr {K}}_t^n\right\} _{n\in {\mathbb N}}$ is collectively compact. Then, the set $\left\{ {\mathscr {K}}_t^nv\right\} _{n\in {\mathbb N}}$ is relatively compact, thus bounded for each $v\in {C_{d}}$, and the uniform boundedness principle [18, p. 395, Thm. 3.2] yields $c_t:=\sup _{n\in {\mathbb N}}\left\| {\mathscr {K}}_t^n\right\| <\infty$ for all $t\in {\mathbb Z}$. The limit relation

$$\begin{aligned} \lim _{n\rightarrow \infty }\Vert \Xi _\theta ^n v-\Xi _\theta ^0 v\Vert _\infty =0\quad \text {for all }v\in C_d \end{aligned}$$

(3.2)

is shown as in [21, step (I) in the proof of Thm. 2.1]. Moreover, we establish by mathematical induction over $\theta \in {\mathbb N}$ that the family $\{\Xi _\theta ^n\}_{n\in {\mathbb N}}$ is collectively compact. The argument for $\theta =1$ was given above. As induction step $\theta \rightarrow \theta +1$, we obtain $\Xi _{\theta +1}^n={\mathscr {K}}_\theta ^n\Xi _\theta ^n$ from (3.1) with the collectively compact family $\left\{ \Xi _\theta ^n\right\} _{n\in {\mathbb N}}$ and since $c_\theta$ is finite, [5, p. 59, Prop. 4.2.2.] yields the claim.

(a) Results from [8, Theorem].

(b) The related arguments to prove [21, Thm. 2.1(a) and (b)] are based on the upper semicontinuity of the spectrum requiring $\lim _{n\rightarrow \infty }\Vert \Xi _\theta ^n-\Xi _\theta ^0\Vert =0$. A corresponding substitute under the weaker assumption (3.2) provides [5, p. 65, Thm. 4.8].

(c) Is a consequence of [5, p. 69, Thm. 4.13].

Proposition  has an immediate impact on the persistence of stability:

Corollary 3.2

(stability of $L^n$ under Nyström discretization)

• (a) If $\sigma (\Xi _\theta ^0)\subset B_1(0)$ holds, then there exists an integer $N\in {\mathbb N}$ such that $L^n$ is uniformly exponentially stable for all $n\in \left\{ 0,N,N+1,\ldots \right\}$.

• (b) If $\sigma (\Xi _\theta ^0)\cap (\mathbb {C}\setminus \bar{B}_1(0))\ne \emptyset$ holds, then there exists an integer $N\in {\mathbb N}$ such that $L^n$ is unstable for all $n\in \left\{ 0,N,N+1,\ldots \right\}$.

Proof

Let $\rho _n:=\rho (\Xi _\theta ^n)$ abbreviate the spectral radius of $\Xi _\theta ^n$, $n\in {\mathbb N}_0$.

(a) By assumption, we have $\rho _0<1$. In the proof of Proposition , we have established that [5, p. 65, Thm. 4.8] applies and $\rho _n<\tfrac{1+\rho _0}{2}<1$ follows for all $n\ge N$. Hence, also $L^n$ is uniformly exponentially stable.

(b) By assumption $(L^0)$ has Morse index $\dim \mathcal {V}_-^0>0$ and Proposition (c) guarantees $\dim \mathcal {V}_-^n>0$ for all $n\ge N$. In both cases, there exists a Floquet multiplier of modulus $>1$ yielding an unstable bundle and hence instability.

For the remainder of this section, we keep $n\in {\mathbb N}_0$ fixed and assume that $\phi ^n$ is a solution to $\Delta ^n$. Under Hypothesis  for $m\in {\mathbb N}$, we conclude from Propositions 2.4 and 2.6 that ${\mathscr {F}}_t^n$ is continuously differentiable and introduce the variation equation

along $\phi ^n$. The previous terminology for general linear integrodifference $L^n$ now applies to the variation $V^n$ with coefficients ${\mathscr {K}}_t^n=D{\mathscr {F}}_t^n(\phi _t^n)$ and the solution $\phi ^n$. In particular, Proposition  and Lemma  guarantee compact coefficients $D{\mathscr {F}}_t^n(\phi _t^n)\in L({C_{d}})$ and therefore compact period operators $\Xi _\theta ^n$.

Lemma 3.3

Let $r>0$ and $n\in {\mathbb N}$. If Hypothesis  holds with $m\in {\mathbb N}$, then the variation $V^n$ have compact coefficient operators $D{\mathscr {F}}_t^n(\phi _t^n)\in L({C_{d}})$. Moreover,

$$\begin{aligned} \left\| D^k{\mathscr {F}}_t^n(u_t)v_1\cdots v_k\right\| _\infty \le W_n \sup _{\eta \in \Omega _n}b_t^k(r,\eta )\left\| v_1\right\| _\infty \cdots \left\| v_k\right\| _\infty \end{aligned}$$

(3.3)

holds and if the integration rule $Q^n$ is stable, then we even have

$$\begin{aligned} \left\| D^k{\mathscr {F}}_t^n(u_t)v_1\cdots v_k\right\| _\infty \le W\sup _{\eta \in \Omega }b_t^k(r,\eta )\left\| v_1\right\| _\infty \cdots \left\| v_k\right\| _\infty \end{aligned}$$

(3.4)

for all $t\in {\mathbb Z}$, $u\in \bar{\mathcal {B}}_r(\phi ^0)$, $v_1,\ldots ,v_k\in {C_{d}}$ and $1\le k\le m$.

Remark 3.4

For stable quadrature rules $Q^n$ it is a consequence of the mean value theorem [18, p. 341, Thm. 4.2] and (3.3) applied for $k=2$ that the Nyström method $\left\{ {\mathscr {F}}_t^n\right\} _{n\in {\mathbb N}}$ is equidifferentiable in $\phi _t^n$ (cf. [5, p. 110, Prop. 6.9]).

Proof

Let $t\in {\mathbb Z}$ and $n\in {\mathbb N}$. First, we conclude from Proposition  that ${\mathscr {F}}_t^n:{C_{d}}\rightarrow {C_{d}}$ is of class $C^m$. Second, due to Proposition  the mapping ${\mathscr {F}}_t^n$ is completely continuous and it results from [20, p. 89, Prop. 6.5] that $D{\mathscr {F}}_t^n(\phi _t)\in L({C_{d}})$ is compact.

We choose $u_s\in \bar{B}_r(\phi _s^0)$, i.e., $\left| u_s(y)-\phi _s^0(y)\right| \le r$ for all $y\in \Omega$ and obtain that

$$\begin{aligned} \left| [D^k{\mathscr {F}}_s^n(u_s)v_1\cdots v_k](x)\right|\overset{(2.6)}{\le } & {} \sum _{\eta \in \Omega _n}w_\eta \left| D_3^kf_s(x,\eta ,u_s(\eta ))\right| \left| v_1(\eta )\right| \cdots \left| v_k(\eta )\right| \\\le & {} \sum _{\eta \in \Omega _n}w_\eta \left| D_3^kf_s(x,\eta ,u_s(\eta ))\right| \left\| v_1\right\| _\infty \cdots \left\| v_k\right\| _\infty \\\overset{(2.4)}{\le } & {} \sum _{\eta \in \Omega _n}w_\eta b_s^k(r,\eta )\left\| v_1\right\| _\infty \cdots \left\| v_k\right\| _\infty \\\overset{(2.1)}{\le } & {} W_n\sup _{\eta \in \Omega _n}b_s^k(r,\eta )\left\| v_1\right\| _\infty \cdots \left\| v_k\right\| _\infty \quad \text {for all }x\in \Omega \end{aligned}$$

and $v_1,\ldots ,v_k\in {C_{d}}$. Passing to the supremum over $x\in \Omega$ yields (3.3) and passing to the least upper bound over $n\in {\mathbb N}$ implies (3.4).

4 Numerical persistence of periodic solutions

This section studies the behavior of $\theta _1$-periodic solutions $\phi ^0$ to $\theta _0$-periodic integrodifference equations $(\Delta ^0)$ under Nyström discretization. For convergent rules $Q^n$ and weakly hyperbolic $\phi ^0$, we establish that the Nyström discretizations $\Delta ^n$, $n\in {\mathbb N}$, possess a periodic solution for sufficiently large n, which converges to $\phi ^0$ as $n\rightarrow \infty$.

For the least common multiple $\theta :={{\,\textrm{lcm}\,}}\left\{ \theta _0,\theta _1\right\}$ and sequences $u=(u_t)_{t\in {\mathbb Z}}$ in ${C_{d}}$, we define

$$\begin{aligned} \hat{u}&:=(u_0,\ldots ,u_{\theta -1})\in {C_{d}}^\theta ,&(\overline{u_0,\ldots ,u_{\theta -1}})&:=(u_{t\text { mod }\theta })_{t\in {\mathbb Z}} \end{aligned}$$

and characterize periodic solutions to $\Delta ^n$ using the cyclic operators

$$\begin{aligned} \hat{\mathscr {F}}^n:{C_{d}}^\theta&\rightarrow {C_{d}}^\theta ,&\hat{\mathscr {F}}^n(\hat{u})&:= \begin{pmatrix} {\mathscr {F}}_{\theta -1}^n(u_{\theta -1})\\ {\mathscr {F}}_0^n(u_0)\\ {\mathscr {F}}_1^n(u_1)\\ \vdots \\ {\mathscr {F}}_{\theta -2}^n(u_{\theta -2}) \end{pmatrix} \quad \text {for all }n\in {\mathbb N}_0. \end{aligned}$$

(4.1)

Here, the Cartesian product ${C_{d}}^\theta$ is equipped with a product norm induced via (1.3). The next result is immediate yet basic (and we skip the proof):

Lemma 4.1

Let $n\in {\mathbb N}_0$ and $\Delta ^n$ be $\theta _0$-periodic:

• (a) If $(\phi _t)_{t\in {\mathbb Z}}$ is a $\theta$-periodic solution to $\Delta ^n$, then $\hat{\phi }\in {C_{d}}^\theta$ is a fixed point of $\hat{\mathscr {F}}^n$.

• (b) Conversely, if $\hat{\phi }\in {C_{d}}^\theta$ is a fixed point of $\hat{\mathscr {F}}^n$, then $(\overline{\phi _0,\ldots ,\phi _{\theta -1}})$ is a $\theta$-periodic solution to $\Delta ^n$.

The subsequent result is crucial for a numerical computation of the Floquet spectra of variation $V^n$ used to determine stability properties of periodic solutions. From now on let $\Xi _\theta ^n\in L(C_d)$ denote the period operator of $V^n$.

Lemma 4.2

Let $\phi ^n$, $n\in {\mathbb N}_0$, be a $\theta$-periodic solution to $\Delta ^n$. If Hypothesis  holds with $m\in {\mathbb N}$, then $\hat{\mathscr {F}}^n:{C_{d}}^\theta \rightarrow {C_{d}}^\theta$ is of class $C^m$ and the following holds:

• (a) $D\hat{\mathscr {F}}^n(\hat{\phi }^n)\in L({C_{d}}^\theta )$ is compact and of cyclic structure

  $$\begin{aligned} \hspace{-5mm} D\hat{\mathscr {F}}^n(\hat{\phi }^n) = \begin{pmatrix} 0 &{} 0 &{} \cdots &{} \cdots &{} D{\mathscr {F}}_{\theta -1}^n(\phi _{\theta -1}^n)\\ D{\mathscr {F}}_0^n(\phi _0^n) &{} 0 &{} \cdots &{} \cdots &{} 0\\ 0 &{} D{\mathscr {F}}_1^n(\phi _1^n) &{} 0 &{} \cdots &{} 0\\ \vdots &{} \ddots &{} \ddots &{} \ddots &{} \vdots \\ 0 &{} \cdots &{} 0 &{} D{\mathscr {F}}_{\theta -2}^n(\phi _{\theta -2}^n) &{} 0 \end{pmatrix}, \end{aligned}$$

  (4.2)

• (b) $\sigma (\Xi _\theta ^n)\setminus \left\{ 0\right\} =\sigma (D\hat{\mathscr {F}}^n(\hat{\phi }^n))^\theta \setminus \left\{ 0\right\}$.

Proof

Because of Propositions 2.4 and 2.6 this results from [18, p. 338].

(a) By [20, p. 89, Prop. 6.5] the derivatives $D{\mathscr {F}}_t^n(\phi _t^n)$ of the completely continuous mappings ${\mathscr {F}}_t^n$ are compact. The $D\hat{\mathscr {F}}^n(\hat{\phi }^n)$-image of the unit ball $B_1(0)^\theta$ is

$$ D{\mathscr {F}}_{\theta -1}^n(\phi _{\theta -1}^n)B_1(0)\times D{\mathscr {F}}_0^n(\phi _0^n)B_1(0)\times \ldots \times D{\mathscr {F}}_{\theta -2}^n(\phi _{\theta -2}^0)B_1(0), $$

where each component set is relatively compact since the derivatives $D{\mathscr {F}}_s^n(u_s)$ are compact. Therefore, by Tychonoff’s theorem [18, p. 37, Thm. 3.12], also the image $D\hat{\mathscr {F}}^n(\hat{\phi }^n)B_1(0)^\theta$ is relatively compact yielding the compactness assertion.

(b) See [21, Lemma 2.3].

Proposition 4.3

(unique periodic solutions) Let Hypothesis  hold with $m\in {\mathbb N}$. If $\phi ^0$ is a weakly hyperbolic, $\theta _1$-periodic solution of $(\Delta ^0)$, then there exists a $\rho _0>0$ so that $\phi ^0$ is the unique $\theta$-periodic solution of $(\Delta ^0)$ in $\mathcal {B}_{\rho _0}(\phi ^0)$.

Proof

Assume by contradiction that $1\in \sigma (D\hat{\mathscr {F}}^0(\hat{\phi }^0))$, i.e., $1\in \sigma _p(D\hat{\mathscr {F}}^0(\hat{\phi }^0))$. Then, Lemma (b) implies $1\in \sigma (D\hat{\mathscr {F}}^0(\hat{\phi }^0))^\theta =\sigma (\Xi _\theta ^0)$, contradicting our weak hyperbolicity assumption. Thus, we can conclude from [7, pp. 801–802, Property P4] that $\hat{\phi }^0$ is an isolated fixed point of $\hat{\mathscr {F}}^0$. By Lemma (b) this implies the claim.

This leads us to our main result:

Theorem 4.4

(periodic solutions under Nyström discretization) Let Hypothesis  hold with $m\ge 2$ and suppose that the integration rule $Q^n$ is convergent. If $\phi ^0$ is a weakly hyperbolic, $\theta _1$-periodic solution of $(\Delta ^0)$, then there exist $\rho \in (0,\rho _0]$ (with $\rho _0>0$ from Proposition ), $C>0$ satisfying

$$\begin{aligned} C\rho \max _{t=0}^{\theta -1}C_t(r)&<1,&C_t(r)&:=\sup _{x\in \Omega }\sup _{y\in \Omega }\sup _{z\in \bar{B}_r(\phi _t^0(y))}\left| D_3^2f_t(x,y,z)\right| \end{aligned}$$

for some $r>0$, and an integer $N_0\in {\mathbb N}$ such that the following hold for all $n\ge N_0$:

• (a) There exists a unique $\theta$-periodic solution $\phi ^n$ to $\Delta ^n$ in $\mathcal {B}_\rho (\phi ^0)$,

• (b) $\lim _{n\rightarrow \infty }\left\| \phi _t^n-\phi _t^0\right\| _\infty =0$ for all $t\in {\mathbb Z}$ and in particular

  $$\begin{aligned} \hspace{-5mm} \sup _{t\in {\mathbb Z}}\left\| \phi _t^n-\phi _t^0\right\| _\infty \le \frac{2C}{2-C\rho \max _{s=0}^{\theta -1}C_s(r)}\max _{s=0}^{\theta -1}\left\| {\mathscr {F}}_s^n(\phi _s^0)-{\mathscr {F}}_s^0(\phi _s^0)\right\| _\infty . \end{aligned}$$

  (4.3)

The actual convergence rates achieved in (4.3) depend on the quadrature rule $Q^n$, but also on smoothness properties of the kernel functions $f_t$ in $(\Delta ^0)$ (cf. Section ).

Proof

Throughout the proof keep $s\in {\mathbb Z}$ fixed. We first verify that the perturbation result [30, Thm. 1] applies to the operator $\hat{\mathscr {F}}:{C_{d}}^\theta \rightarrow {C_{d}}^\theta$ defined in (4.1).

(I) Claim: $\{\hat{\mathscr {F}}^n\}_{n\in {\mathbb N}}$ is collectively compact.

The collective compactness of $\{\hat{\mathscr {F}}^n\}_{n\in {\mathbb N}}$ means that for every bounded $\hat{B}\subset {C_{d}}^\theta$ the union $\bigcup _{n\in {\mathbb N}}\hat{\mathscr {F}}^n(\hat{B})$ is relatively compact. For this, due to Tychonoff’s theorem [18, p. 37, Thm. 3.12] and (4.1) it suffices to show that for every bounded set $B\subset {C_{d}}$ the union $\bigcup _{n\in {\mathbb N}}{\mathscr {F}}_s^n(B)\subset {C_{d}}$ is relatively compact. Since $Q^n$ is convergent, it is stable and Proposition (b) guarantees that $\bigcup _{n\in {\mathbb N}}{\mathscr {F}}_s^n(B)$ is relatively compact.

(II) Claim: $\lim _{n\rightarrow \infty }\Vert \hat{\mathscr {F}}^n(\hat{u})-\hat{\mathscr {F}}(\hat{u})\Vert =0$ for all $\hat{u}\in {C_{d}}^\theta$.

Because of (4.1) it suffices to show (2.7), which was established in Proposition .

(III) Claim: $\hat{\mathscr {F}}^n:{C_{d}}^\theta \rightarrow {C_{d}}^\theta$ is of class $C^2$ and satisfies

$$\begin{aligned} \left\| D^2\hat{\mathscr {F}}^n(\hat{u})\right\| \le W\max _{s=0}^{\theta -1}C_s(r) \quad \text {for all }n\in {\mathbb N},\,\hat{u}\in B_r(0)\subseteq {C_{d}}^\theta . \end{aligned}$$

(4.4)

Due to Proposition  the mappings ${\mathscr {F}}_t^n$, and thus $\hat{\mathscr {F}}^n$, are of class $C^2$. Let $u_s\in \bar{B}_r(\phi _s^0)$, i.e., $\left| u_s(y)-\phi _s^0(y)\right| \le r$ for $y\in \Omega$, and $v_i:=(v_i^{\theta -1},v_i^0,\ldots v_i^{\theta -2})\in {C_{d}}^\theta$ for $i=1,2$. Then, Lemma  implies

$$\begin{aligned} \left\| D^2\hat{\mathscr {F}}^n(\hat{u})v_1v_2\right\|&\overset{(1.3)}{=} \max _{s=0}^{\theta -1}\left\| D^2{\mathscr {F}}_s^n(u_s)v_1^sv_2^s\right\| _\infty \\&\overset{(3.3)}{\le } W_n\max _{s=0}^{\theta -1}\sup _{\eta \in \Omega _n}b_s^2(r,\eta )\left\| v_1^s\right\| _\infty \left\| v_2^s\right\| _\infty \end{aligned}$$

and since convergence of $Q^n$ implies stability, (2.1) yields the desired inequality $\bigl \Vert D^2\hat{\mathscr {F}}^n(\hat{u})\bigr \Vert \le W\max _{s=0}^{\theta -1}\sup _{\eta \in \Omega }b_s^2(r,\eta )$ for all $n\in {\mathbb N}$. Hence, (4.4) holds.

(IV) Claim: $\text {id}_{{C_{d}}^\theta }-D\hat{\mathscr {F}}(\hat{\phi }^0)$ is invertible.

As in the proof of Proposition , we show that $1\not \in \sigma (D\hat{\mathscr {F}}^0(\hat{\phi }^0))$, i.e., $\text {id}_{{C_{d}}^\theta }-D\hat{\mathscr {F}}(\hat{\phi }^0)$ has a bounded inverse. We now obtain the assertions of Theorem  as follows:

(a) Having established (I–III), [30, Thm. 1] guarantees that there exist $\rho \in (0,\rho _0]$ and $N_0\in {\mathbb N}$ such that $\hat{\mathscr {F}}^n$ has a unique fixed point $\hat{\phi }^n\in B_\rho (\hat{\phi }^0)$ for all $n\ge N_0$. Due to Lemma (b) this implies that the corresponding $\theta$-periodic solutions $\phi ^n$ of $\Delta ^n$ are uniquely determined in $\mathcal {B}_\rho (\phi ^0)$.

(b) It follows directly from [30, Thm. 1] that

$$\begin{aligned} \left\| \phi _t^n-\phi _t^0\right\| _\infty \overset{(1.3)}{\le } \left\| \hat{\phi }^n-\hat{\phi }^0\right\| \xrightarrow [n\rightarrow \infty ]{}0 \quad \text {for all }t\in {\mathbb Z}, \end{aligned}$$

where the proof of [7, Thm. 4, p. 552] actually contains the estimate

$$ \left\| \phi _t^n-\phi _t^0\right\| _\infty \le \left\| \hat{\phi }^n-\hat{\phi }^0\right\| \le \frac{C}{1-\tfrac{1}{2}C\rho \max _{s=0}^{\theta -1}C_s(r)} \left\| \hat{\mathscr {F}}^n(\hat{\phi }^0)-\hat{\mathscr {F}}(\hat{\phi }^0)\right\| , $$

which implies the assertion (4.3).

Due to the error analysis provided by Theorem , we can now approximate periodic solutions of $(\Delta ^0)$ by solutions of Nyström discretizations $\Delta ^n$, $n\in {\mathbb N}$. However, rather than working with the Nyström interpolants ${\mathscr {F}}_t^n:{C_{d}}\rightarrow {C_{d}}$ our computations rely on the mappings $F_t^n:{\mathbb R}^{dq_n}\rightarrow {\mathbb R}^{dq_n}$ introduced in (2.8). Since typical applications feature smooth kernel functions $f_t$ also the mappings $F_t^n$ become differentiable. Hence, instead of working with the fixed point problem $\hat{\phi }=\hat{\mathscr {F}}(\hat{\phi })$ in ${C_{d}}^\theta$, we determine zeros to the equivalent problem

$$\begin{aligned} \hat{G}^n(\hat{\upsilon })&=0,&\hat{G}^n(\hat{\upsilon })&:= \begin{pmatrix} \upsilon _0-F_{\theta -1}^n(\upsilon _{\theta -1})\\ \upsilon _1-F_0^n(\upsilon _0)\\ \upsilon _2-F_1^n(\upsilon _1)\\ \vdots \\ \upsilon _{\theta -1}-F_{\theta -2}^n(\upsilon _{\theta -2}) \end{pmatrix} \end{aligned}$$

(4.5)

in ${\mathbb R}^{\theta dq_n}$ using on a suitable (exact or inexact) Newton method.

Algorithm 2

Computation of periodic solutions to Urysohn difference (1.1)

The pseudo-code from Algorithm 2 describes an exact Newton method, but can easily be adapted to inexact methods (or further Newton-like schemes) as follows: Before the loop (Newton iteration) one defines $\hat{\upsilon }^0\leftarrow \hat{\upsilon }$ and performs Newton iteration based on the linear systems $D\hat{G}^n(\hat{\upsilon }^0)\delta =-\hat{G}^n(\hat{\upsilon })$ with fixed coefficient matrix. In both cases, as convergence criterion to interrupt the Newton iteration, we use the $\infty$-norm of $\hat{G}^n(\hat{\upsilon })$ to be sufficiently small.

Corollary 4.5

(Floquet multipliers of $(V^0)$ under Nyström discretization) If $\lambda \ne 0$ is a Floquet multiplier of $(V^0)$ with index $\nu (\lambda )\in {\mathbb N}$ and $\left\{ e_1,\ldots ,e_{a(\lambda )}\right\}$ denotes a basis of ${{\,\textrm{Eig}\,}}_\lambda \Xi _\theta ^0$, then there exist $C\ge 0$ and an integer $N_1\ge N_0$ such that the following holds for all $n\ge N_1$:

• (a) $\max _{\mu \in \Lambda _\rho ^n(\lambda )}\left| \mu -\lambda \right| \le C\max _{i=1}^{a(\lambda )}\root \nu (\lambda ) \of {\left\| [\Xi _\theta ^n-\Xi _\theta ^0]e_i\right\| _\infty }$,

• (b) $\phi ^n$ is weakly hyperbolic,

• (c) with $\phi ^0$ also $\phi ^n$ is hyperbolic having the same Morse index as $\phi ^0$.

Proof

Let $t\in {\mathbb Z}$ and $v\in {C_{d}}$. We define

$$\begin{aligned} {\mathscr {K}}_tv&:=D{\mathscr {F}}_t^0(\phi _t^0)v=\int _\Omega D_3f_t(\cdot ,y,\phi _t^0(y))v(y)\,{\mathrm d}y,\\ {\mathscr {K}}_t^nv&:=D{\mathscr {F}}_t^n(\phi _t^n)v=\sum _{\eta \in \Omega _n}w_\eta D_3f_t(\cdot ,\eta ,\phi _t^n(\eta ))v(\eta ) \quad \text {for all }n\ge N_0, \end{aligned}$$

and begin with two preparations:

(I) Claim: $\lim _{n\rightarrow \infty }\left\| {\mathscr {K}}_t^nv-{\mathscr {K}}_tv\right\| _\infty =0$

As in Proposition , we obtain $\lim _{n\rightarrow \infty }\left\| D{\mathscr {F}}_t^n(\phi _t^0)v-D{\mathscr {F}}_t^0(\phi _t^0)v\right\| _\infty =0$ and according to Proposition (b), we have $\lim _{n\rightarrow \infty }\left\| D{\mathscr {F}}_t^n(\phi _t^n)v-D{\mathscr {F}}_t^n(\phi _t^0)v\right\| _\infty =0$. These limit relations combined with the triangle inequality yields

$$ 0 \le \left\| {\mathscr {K}}_t^nv-{\mathscr {K}}_tv\right\| _\infty \le \left\| D{\mathscr {F}}_t^n(\phi _t^n)v-D{\mathscr {F}}_t^n(\phi _t^0)v\right\| _\infty + \left\| D{\mathscr {F}}_t^n(\phi _t^0)v-D{\mathscr {F}}_t^0(\phi _t^0)v\right\| _\infty $$

and therefore the claim.

(II) Claim: $\left\{ {\mathscr {K}}_t^n\right\} _{n\in {\mathbb N}}$ is collectively compact.

We apply the Arzelà-Ascoli theorem [18, p. 58, Cor. 3.2] in order to show that the union $\bigcup _{n\in {\mathbb N}}{\mathscr {K}}_t^nB_1(0)$ is relatively compact. First, it is not difficult to check that the union is bounded. Second, in order to prove equicontinuity, we choose $\varepsilon >0$ and a function $u\in \bigcup _{n\in {\mathbb N}}{\mathscr {K}}_t^nB_1(0)$, i.e., there exist $n\in {\mathbb N}$, $v_0\in B_1(0)$ such that $u=D{\mathscr {F}}_t^n(\phi _t^0)v_0$. A uniform continuity argument on the compact product $\Omega \times \Omega$ guarantees that there exists a $\delta >0$ such that $\left| x-\bar{x}\right| <\delta$ and $\left| z-\bar{z}\right| <\delta$ imply

$$ \left| D_3f_t(x,y,z)-D_3f_t(\bar{x},y,\bar{z})\right| < \tfrac{\varepsilon }{2W}\quad \text {for all }x_,\bar{x},y\in \Omega ,\,z,\bar{z}\in K $$

on each compact $K\subseteq {\mathbb R}^d$. Given this, Theorem  allows us to choose $N_1\ge N_0$ such that $\left| \phi _t^n(y)-\phi _t^0(y)\right| \le \left\| \phi _t^n-\phi _t^0\right\| _\infty <\delta$ for all $n\ge N_1$ and $y\in \Omega$. Consequently,

$$\begin{aligned} \left| u(x)-u(\bar{x})\right|\overset{(2.6)}{\le } & {} \sum _{\eta \in \Omega _n}w_\eta \left| \left[ D_3f_t(x,\eta ,\phi _t^n(\eta ))-D_3f_t(\bar{x},\eta ,\phi _t^0(\eta ))\right] v_0(\eta )\right| \\\le & {} \sum _{\eta \in \Omega _n}w_\eta \left| \left[ D_3f_t(x,\eta ,\phi _t^n(\eta ))-D_3f_t(\bar{x},\eta ,\phi _t^0(\eta ))\right] \right| \overset{(2.1)}{<} \tfrac{\varepsilon }{2} \end{aligned}$$

for all $x,\bar{x}\in \Omega$. Since $\delta >0$ is independent of $n\in {\mathbb N}$ and $v_0\in B_1(0)$, this shows equicontinuity.

Based on steps (I–II) the assertions follow as in the proof of Proposition .

Corollary 4.6

(stability of $\phi ^0$ under Nyström discretization)

• (a) If $\sigma (\Xi _\theta ^0)\subset B_1(0)$ holds, then there exists an integer $N\ge N_1$ such that $\phi ^n$ is uniformly exponentially stable for all $n\in \left\{ 0,N,N+1,\ldots \right\}$.

• (b) If $\sigma (\Xi _\theta ^0)\cap (\mathbb {C}\setminus \bar{B}_1(0))\ne \emptyset$ holds, then there exists an integer $N\ge N_1$ such that $\phi ^n$ is unstable for all $n\in \left\{ 0,N,N+1,\ldots \right\}$.

Proof

Note that the dichotomy spectrum $\Sigma (\phi ^n)$ (cf. [24]) and the Floquet spectrum of a $\theta$-periodic solution $\phi ^n$ are related by $\Sigma (\phi ^n)=\root \theta \of {\left| \sigma (\Xi _\theta ^n)^\theta \right| }\setminus \left\{ 0\right\}$ for $n=0$ or all $n\ge N_0$. This allows us to employ the stability criterion [24, Thm. 2.1]. If $n=0$, then the assertions result. If $n\ge N_0$ then Corollary  guarantees $\sigma (\Xi _\theta ^n)\subset B_1(0)$ (in case (a)) resp. $\sigma (\Xi _\theta ^n)\cap (\mathbb {C}\setminus \bar{B}_1(0))\ne \emptyset$ (in case (b)) for sufficiently large n, which guarantees the claims.

After having seen how to approximate periodic solutions $\phi ^0$ to Urysohn difference (1.1) numerically in Algorithm 2, we now determine their Floquet spectrum and hence their stability properties. This first of all requires an approximation $\phi ^n$ as $\theta$-periodic solution to the Nyström discretization (1.2) computed using Algorithm 2. Thus, the variation $V^n$ can be evaluated and error estimates relating the Floquet multipliers of $\phi ^n$ to those of $\phi ^0$ are provided by Corollary . Nevertheless, as pointed out in, e.g., [12], even in finite dimension the numerical computation of the Floquet spectrum via the eigenvalues of the period operator $\Xi _\theta ^n$ is not backward stable. Hence, we approximate the Floquet multipliers of $(V^0)$ using the cyclic operator (4.2) for a sufficiently large $n\in {\mathbb N}$. Indeed, thanks to Lemma (b) the $\theta$th powers of the eigenvalues of $D\hat{\mathscr {F}}^n(\hat{\phi }^n)$ are of the Floquet multipliers for $V^n$.

In detail, projected to finite-dimensions the periodic solution $\phi ^n$ of (1.2) corresponds to a $\theta$-periodic solution $\upsilon ^n$ of $\Delta ^n$ and coefficients $D{\mathscr {F}}_t^n(\phi _t^n)\in L(C_d)$ in the variation $V^n$ become $dq_n\times dq_n$-block matrices

$$ DF_t^n(\upsilon _t^n) = \begin{pmatrix} w_1D_3f_t(\eta _1,\eta _1,\upsilon _{t,1}^n) &{} \ldots &{} w_{q_n}D_3f_t(\eta _1,\eta _{q_n},\upsilon _{t,q_n}^n)\\ \vdots &{} &{} \vdots \\ w_1D_3f_t(\eta _{q_n},\eta _1,\upsilon _{t,1}^n) &{} \ldots &{} w_{q_n}D_3f_t(\eta _{q_n},\eta _{q_n},\upsilon _{t,q_n}^n) \end{pmatrix} $$

consisting of the blocks $w_iD_3f_t(\eta _i,\eta _j,\upsilon _{t,j}^n)\in {\mathbb R}^{d\times d}$. This, in turn, yields the cyclic $\theta dq_n\times \theta dq_n$-block matrices

$$\begin{aligned} D\hat{F}^n(\hat{\upsilon }^n) = \begin{pmatrix} 0 &{} 0 &{} \cdots &{} \cdots &{} DF_{\theta -1}^n(\upsilon _{\theta -1}^n)\\ DF_0^n(\upsilon _0^n) &{} 0 &{} \cdots &{} \cdots &{} 0\\ 0 &{} DF_1^n(\upsilon _1^n) &{} 0 &{} \cdots &{} 0\\ \vdots &{} \ddots &{} \ddots &{} \ddots &{} \vdots \\ 0 &{} \cdots &{} 0 &{} DF_{\theta -2}^n(\upsilon _{\theta -2}^n) &{} 0 \end{pmatrix} \end{aligned}$$

(4.6)

from which we approximate the Floquet multipliers of the $\theta$-periodic solution $\phi ^0$. Keeping Lemma (b) in mind we first determine the spectrum $\sigma (D\hat{F}^n(\hat{\phi }^n))$. Then, in order to respect potential numerical inaccuracies, we locate clusters of neighboring elements of $\sigma (D\hat{F}^n(\hat{\phi }^n))^\theta$. A cluster $C\subseteq \mathbb {C}$ is understood as a nonempty set of (up to) $\theta$ elements of $\sigma (D\hat{F}^n(\hat{\phi }^n))^\theta$ all contained in a sufficiently small ball. There exist up to $dq_n$ such clusters. Then, the arithmetic mean $\lambda :=\frac{1}{\#C}\sum _{\gamma \in C}\gamma$ approximates a Floquet multiplier of $V^n$, and in turn a Floquet multiplier of the initial $\theta$-periodic solution $\phi ^0$ to (1.1). We point to Algorithm 3 for an algorithmic description of this procedure.

Algorithm 3

Computation of the Floquet spectrum for a periodic solution to an Urysohn difference (1.1)

Remark 4.7

(1) Note that the functions $\phi _\tau ^n,\ldots ,\phi _{\tau +\theta -1}^n$ required in Algorithm 3 could also coincide with an exact $\theta$-periodic solution $\phi ^0$ to an Urysohn difference (1.1), provided it is known. In general, the $\phi _\tau ^n,\ldots ,\phi _{\tau +\theta -1}^n$ are computed using Algorithm 2.

(2) In applications it is often unnecessary to determine the entire spectrum of the block-cyclic matrix $D\hat{F}^n(\hat{\upsilon })$ (or all clusters), which could contain up to $\theta dq_n$ elements. In fact, we only need Floquet multipliers of large magnitude to infer stability properties. Furthermore, being eigenvalues of a compact operator, the Floquet multipliers typically accumulate near zero which complicates to identify clusters of eigenvalues having small magnitude.

5 Convergence rates

The crucial error estimate (4.3) in Theorem (b) relating the periodic solutions $\phi ^0$ of an integrodifference equation $(\Delta ^0)$ and $\phi ^n$ of the Nyström discretizations $\Delta ^n$ requires information on the local discretization error ${\mathscr {F}}_t^n(\phi _t^0)-{\mathscr {F}}_t^0(\phi _t^0)$. This section establishes that the convergence order of an integration rule $Q^n$ extends to these differences for sufficiently smooth kernel functions $f_t$. Furthermore, a similar statement also holds for the error estimates provided by Corollary (a) concerning the Floquet multipliers under Nyström discretization.

5.1 Quadrature rules for Hölder-continuous kernel functions

Assume that $\Omega =[a,b]$ is an interval with reals $a<b$. In [21, Exam. 3.1] it is shown that the composite midpoint, trapezoidal and Simpson rules $Q^n$ have convergence order $\alpha \in (0,1]$, i.e., there exists a $c_0\ge 0$ such that

$$\begin{aligned} \left| \int _a^bu(y)\,{\mathrm d}y-Q^nu\right| \le \frac{c_0}{n^\alpha }[u]_\alpha \quad \text {for all }u\in C^\alpha [a,b]. \end{aligned}$$

(5.1)

In addition to the standard assumptions from Hypothesis , we suppose that $f_t$ fulfills the following Hölder condition: There exist $L_t\ge 0$ and $\alpha \in (0,1]$ such that

$$ \left| f_t(x,y,z)-f_t(x,\bar{y},\bar{z})\right| \le L_t\max \left\{ \left| y-\bar{y}\right| ^\alpha ,\left| z-\bar{z}\right| \right\} $$

for all $x,y,\bar{y}\in [a,b]$ and $z,\bar{z}\in \phi _t^0([a,b])$. In case $\phi _t^0\in C^\alpha [a,b]$ this readily implies

$$\begin{aligned} \left| f_t(x,y,\phi _t^0(y))-f_t(x,\bar{y},\phi _t^0(\bar{y}))\right|&\le L_t\max \left\{ \left| y-\bar{y}\right| ^\alpha ,\left| \phi _t^0(y)-\phi _t^0(\bar{y})\right| \right\} \\&\le L_t\max \left\{ 1,[\phi _t^0]_\alpha \right\} \left| y-\bar{y}\right| ^\alpha \quad \text {for all }y,\bar{y}\in [a,b], \end{aligned}$$

i.e., the function $f_t(x,\cdot ,\phi _t^0(\cdot )):[a,b]\rightarrow {\mathbb R}^d$ satisfies a Hölder condition with exponent $\alpha$ uniformly in $x\in [a,b]$. Consequently, using the consistency condition (5.1) this yields the error estimate

$$ \left\| {\mathscr {F}}_t^n(\phi _t^0)-{\mathscr {F}}_t^0(\phi _t^0)\right\| _\infty \le \frac{c_0}{n^\alpha }L_t\max \left\{ 1,[\phi _t^0]_\alpha \right\} \quad \text {for all }t\in {\mathbb Z}$$

and one obtains a convergence order $\alpha \in (0,1]$.

5.2 Cubature and differentiable kernel functions

For $u\in C^{d_p}(\Omega )$ assume that $Q^n$ allows a quadrature error of the form

$$\begin{aligned} \left| \int _a^b u(y)\,{\mathrm d}y-Q^n u\right| \le \frac{c_p}{n^{d_p}}\sup _{x\in \Omega }\left| D^{d_p}u(x)\right| . \end{aligned}$$

(5.2)

The values of the real $c_p$ and $d_p\in {\mathbb N}$ for common quadrature rules are summarized in Appendix B. It is convenient to define the mappings

$$\begin{aligned} F_t:\Omega \times \Omega&\rightarrow {\mathbb R}^d,&F_t(x,y)&:=f_t(x,y,\phi _t^0(y))\quad \text {for all }t\in {\mathbb Z},\\ \Phi ^t:\Omega&\rightarrow \Omega \times {\mathbb R}^d,&\Phi ^t(y)&:=(y,\phi _t^0(y)). \end{aligned}$$

If the derivatives of $f_t(x,\cdot ):\Omega \times {\mathbb R}^d\rightarrow {\mathbb R}^d$ exist up to order $d_p$ as continuous functions and $\phi _t^0\in C^{d_p}(\Omega )$, then the higher order chain rule [29] implies

$$\begin{aligned} D_2^lF_t(x,y)h_1\cdots h_j\\ = \sum _{j=1}^l\sum _{(N_1,\ldots ,N_j)\in P_j^<(l)} D_{(2,3)}^jf_t(x,y,\phi _t^0(y))D^{\#N_1}\Phi ^t(y)h_{N_1}\ldots D^{\#N_j}\Phi ^t(y)h_{N_j} \end{aligned}$$

for all $0\le l\le d_p$, $x,y\in \Omega$. Here, with given $j,l\in {\mathbb N}$, we define the ordered partitions of $\left\{ 1,\ldots ,l\right\}$ having length j by

$$\begin{aligned} P_j^<(l):=\left\{ (N_1,\ldots ,N_j)\left| \, \begin{array}{l} N_i\subseteq \left\{ 1,\ldots ,l\right\} \text { and } N_i\ne \emptyset \text { for }i\in \left\{ 1,\ldots ,j\right\} ,\\ N_1\cup \ldots \cup N_j=\left\{ 1,\ldots ,l\right\} ,\\ N_i\cap N_k=\emptyset \text { for }i\ne k,\,i,k\in \left\{ 1,\ldots ,j\right\} ,\\ \max N_i<\max N_{i+1}\text { for }i\in \left\{ 1,\ldots ,j-1\right\} \end{array}\!\! \right. \right\} . \end{aligned}$$

If $N=\left\{ n_1,\ldots ,n_k\right\} \subseteq \left\{ 1,\ldots ,l\right\}$ with $k\in {\mathbb N}$, $k\le l$, then we briefly write $D^j\Phi ^t(y)h_N:=D^j\Phi ^t(y)h_{n_1}\cdots h_{n_j}$ for $y\in \Omega$ and $h_1,\ldots ,h_l\in {\mathbb R}^\kappa$.

In conclusion, we arrive at the error estimate

$$ \left\| {\mathscr {F}}_t^n(\phi _t^0)-{\mathscr {F}}_t^0(\phi _t^0)\right\| _\infty \le \frac{c_p}{n^{d_p}} \sup _{x,y\in \Omega } \left| D_2^{d_p}F_t(x,y)\right| \quad \text {for all }t\in {\mathbb Z}$$

guaranteeing a convergence order $d_p$.

6 Numerical simulations of periodic integrodifference equations

In applications, problems will typically depend on parameters. Rather than $\Delta ^n$, we therefore deal with $\theta _0$-periodic difference equations

involving a parameter $\alpha$, which in ecological applications might be the carrying capacity of a habitat, or a dispersal or growth rate of a population.

As a consequence, also the finite-dimensional difference (2.9) depend on $\alpha$ and so do the problems

$$\begin{aligned} \hat{G}^n(\hat{\upsilon },\alpha )&=0,&\hat{G}^n(\hat{\upsilon },\alpha )&:= \begin{pmatrix} \upsilon _0-F_{\theta -1}^n(\upsilon _{\theta -1},\alpha )\\ \upsilon _1-F_0^n(\upsilon _0,\alpha )\\ \upsilon _2-F_1^n(\upsilon _1,\alpha )\\ \vdots \\ \upsilon _{\theta -1}-F_{\theta -2}^n(\upsilon _{\theta -2},\alpha ) \end{pmatrix}. \end{aligned}$$

(6.1)

Whence, the zeros of $\hat{G}^n$ characterizing the periodic solutions of (2.9) are functions of $\alpha$. Let us now outline how to determine such branches of periodic solutions using pseudo-arclength continuation (cf., for instance, [1]). For parameters $\alpha \in {\mathbb R}$ this leads to functions $\hat{G}^n:{\mathbb R}^{\theta dq_n}\times {\mathbb R}\rightarrow {\mathbb R}^{\theta dq_n}$.

• (0) Set $k=0$, choose $h>0$ (step size), $K\in {\mathbb N}$ (number of steps). Given an initial parameter value $\alpha _0\in {\mathbb R}$, we solve

  $$\begin{aligned} {\left\{ \begin{array}{ll} \hat{G}^n(\hat{\upsilon }^0,\alpha _0)=0,\\ D\hat{G}^n(\hat{\upsilon }^0,\alpha _0)z_0=0 \end{array}\right. } \end{aligned}$$

  (6.2)

  for $\hat{\upsilon }^0\in {\mathbb R}^{\theta dq_n}$, $z_0=(\zeta _0,\delta _0)\in {\mathbb R}^{\theta dq_n}\times {\mathbb R}$. To ensure the existence of a unique solution, we additionally impose that the second equation in (6.2) (for $z_0$) is supplemented with the additional condition $z_0^Tz_0=1$. For the final component $\delta _0$ in $z_0$, we require $\delta _0>0$ when interested in forward continuation resp. $\delta _0<0$ when interested in backward continuation. If $D\hat{G}^n(\hat{\upsilon }^0,\alpha _0)$ has full rank, then the implicit function theorem shows that there exists a smooth curve $(\hat{\upsilon },\alpha ):I\rightarrow {\mathbb R}^{\theta dq_n}\times {\mathbb R}$ defined on a neighborhood I of 0 so that $\hat{G}^n(\hat{\upsilon }(s),\alpha (s))\equiv 0$ on I and $\hat{\upsilon }(0)=\hat{\upsilon }^0$, $\alpha (0)=\alpha _0$ and $(\hat{\upsilon },\alpha )'(0)=z_0$.

• (1) Predictor step: Set $\tilde{\upsilon }^{k+1}=\hat{\upsilon }^k+h\zeta _k$ and $\tilde{\alpha }_{k+1}:=\alpha _k+h\delta _k$.

• (2) Corrector step: Iteratively solve

  $$\begin{aligned} {\left\{ \begin{array}{ll} \hat{G}^n(\hat{\upsilon }^{k+1},\alpha _{k+1})=0,\\ z_k^T(\left( {\begin{array}{c}\hat{\upsilon }^{k+1}\\ \alpha _{k+1}\end{array}}\right) -\left( {\begin{array}{c}\tilde{\upsilon }^{k+1}\\ \tilde{\alpha }_{k+1}\end{array}}\right) )=0 \end{array}\right. } \end{aligned}$$

  (6.3)

  for $\hat{\upsilon }^{k+1}\in {\mathbb R}^{\theta dq_n}$ and $\alpha _{k+1}\in {\mathbb R}$ with initial values $\tilde{\upsilon }^{k+1}$ and $\tilde{\alpha }_{k+1}$.

• (3) Determine the new tangent direction $z_{k+1}=(\zeta _{k+1},\delta _{k+1})\in {\mathbb R}^{\theta dq_n}\times {\mathbb R}$ from the linear equation

  $$\begin{aligned} {\left\{ \begin{array}{ll} D\hat{G}^n(\hat{\upsilon }^{k+1},\alpha _{k+1})z_{k+1}=0,\\ z_k^Tz_{k+1}=1. \end{array}\right. } \end{aligned}$$

  (6.4)

• (4) Increase k by 1 and go to step (1), while $k\le K$

Pseudo-code for self-starting pseudo-arclength continuation is given in Algorithm 4.

Remark 6.1

(numerical implementation) (1) Characteristic for the tasks tackled in this paper is that the mappings $\hat{G}^n$ and the linearizations have a (block-) cyclic structure (4.5) resp. (4.6). For these problems, there are tailor-made methods to determine corresponding eigenvalues in a numerically stable way. We refer to the periodic QR algorithm (cf. [16]) suitable for small and mid-sized values $dq_n<50$, i.e., dimension $\kappa =1$) and the periodic Krylov-Schur algorithm (see [17] for dimensions $\kappa >1$).

(2) We implemented Algorithms 1–4 in Python (and C++) using numerical solvers from the PETSc¹ toolbox. For the solution of nonlinear (4.5) and (6.1), (6.2), and (6.3), we rely on line search Newton methods (SNESNEWTONLS with KSPGMRES and PCNONE). The associate eigenvalue problems (4.6) are tackled using routines from the SLEPc² library. In particular, we are working with our own implementation of Algorithm 4. We are nevertheless well-aware that there exist continuation packages in Python (such as PyNCT³). They might be sufficiently flexible to interact with our routines to solve nonlinear equations and eigenvalue problems.

Algorithm 4

Pseudo-arclength continuation of periodic solutions to $(\Delta _\alpha ^0)$

6.1 Logistic equations with separable kernels

Over the habitat $\Omega =\left[ -\tfrac{L}{2},\tfrac{L}{2}\right]$, we consider logistic integrodifference equations

$$\begin{aligned} u_{t+1}(x)=a_t\int _{-\tfrac{L}{2}}^{\tfrac{L}{2}}k(x,y)u_t(y)(1-u_t(y))\,{\mathrm d}y\quad \text {for all }x\in \Omega \end{aligned}$$

(6.5)

with a periodic sequence $(a_t)_{t\in {\mathbb Z}}$ of positive reals and $L>0$. We choose the kernel

$$ k(x,y):=\frac{\pi }{4\,L} {\left\{ \begin{array}{ll} \cos \left( \tfrac{\pi }{2\,L}(x-y)\right) ,&{}\left| x-y\right| \le L,\\ 0,&{}\text {else}, \end{array}\right. } $$

which is degenerate (cf. [14]), since it allows the representation

$$\begin{aligned} k(x,y)&= \cos \left( \tfrac{\pi }{2L}x\right) \cos \left( \tfrac{\pi }{2L}y\right) +\sin \left( \tfrac{\pi }{2L}x\right) \sin \left( \tfrac{\pi }{2L}y\right) \\&= \cos \left( \tfrac{\pi }{2L}x\right) e_1(y)+\sin \left( \tfrac{\pi }{2L}x\right) e_2(y) \quad \text {for all }x,y\in \Omega ,\,|x-y|\le L \end{aligned}$$

with the functions $e_1,e_2:[-\tfrac{L}{2},\tfrac{L}{2}]\rightarrow {\mathbb R}$,

$$\begin{aligned} e_1(x)&:=\cos \left( \tfrac{\pi x}{2L}\right) ,&e_2(x)&:=\sin \left( \tfrac{\pi x}{2L}\right) . \end{aligned}$$

Note that this degenerate kernel k inherits the smoothness of the $\cos$-function on the domain of integration, i.e., it is real-analytical. The ansatz $u_t:=(v_1)_te_1+(v_2)_te_2$ in (6.5) leads to the periodic difference equation

$$\begin{aligned} v_{t+1}&=h_t(v_t),&h_t(v)&:= \frac{a_t}{24} \begin{pmatrix} (3(\pi +2)v_1-10\sqrt{2}v_1^{2}-2\sqrt{2}v_2^{2}\\ {\left( 3\pi - 4\sqrt{2}v_1 - 6\right) } v_2 \end{pmatrix} \end{aligned}$$

in ${\mathbb R}^2$ being independent of L. In particular, for the 3-periodic coefficient sequence $(a_t)_{t\in {\mathbb Z}}$ with $a_0=4$, $a_1=6$, $a_2=8$ it has three 3-periodic solutions, which can be obtained by solving the fixed-point equation $h_2\bigl (h_1(h_0(v,0))\bigr )=\left( {\begin{array}{c}v\\ 0\end{array}}\right)$ for $v\in {\mathbb R}$, among which we consider

$$ v_t^*= {\left\{ \begin{array}{ll} (0.5646522179381946,0),&{}t\mod 3=0,\\ (0.7001113170615100,0),&{}t\mod 3=1,\\ (0.9668027543129443,0),&{}t\mod 3=2 \end{array}\right. } \quad \text {for all }t\in {\mathbb Z}. $$

Consequently, the integrodifference (6.5) has the 3-periodic solution

$$ \phi _t^*= e_1 {\left\{ \begin{array}{ll} 0.5646522179381946,&{}t\mod 3=0,\\ 0.7001113170615100,&{}t\mod 3=1,\\ 0.9668027543129443,&{}t\mod 3=2 \end{array}\right. } \quad \text {for all }t\in {\mathbb Z}, $$

which serves as a reference solution. We compute 3-periodic solutions for Nyström discretizations $\Delta ^n$ of the logistic integrodifference (6.5) equipped with the above kernel based on the cyclic systems (4.5). The initial vector is set to the exact solution evaluated at the nodes of the respective quadrature rule. We display the error

$$ {{\,\textrm{err}\,}}_n:=\max _{t=0}^2\left\| \phi _{t}^{n} - \phi _t^*\right\| _\infty $$

with respect to the Nyström interpolant $\phi _t^{n}(x)$ for the different quadrature rules from Appendix B as a function of n in Fig. 1.

Fig. 1

Error ${{\,\textrm{err}\,}}_n$ in the Nyström approximations for the 3-periodic solution $\phi ^*$ of the logistic integrodifference (6.5) with $L=10$ and different quadrature rules

This confirms the theoretical results obtained in Theorem (b) combined with Section .

6.2 Generalized Beverton–Holt equation with Laplace kernel

Over a compact habitat $\Omega \subset {\mathbb R}^\kappa$ let us consider the generalized Beverton–Holt integrodifference equation

$$\begin{aligned} u_{t+1}(x)=\alpha \int _\Omega k(x,y)\frac{c_tu_t(y)}{1+u_t(y)^4}\,{\mathrm d}y\quad \text {for all }x\in \Omega \end{aligned}$$

(6.6)

involving a $\theta _0$-periodic sequence $(c_t)_{t\in {\mathbb Z}}$ of positive reals and the bifurcation parameter $\alpha >0$. The following kernel will be central for our analysis:

Example 6.1

(Laplace kernel in 1d) For $a,L>0$ and $\Omega =[-\tfrac{L}{2},\tfrac{L}{2}]$ the Laplace kernel

$$ k(x,y):=\frac{a}{2}e^{-a\left| x-y\right| }\quad \text {for all }x,y,\in {\mathbb R}$$

leads to the linear integral operator

$$ {\mathscr {K}}u(x):=\frac{a}{2}\int _{-L/2}^{L/2}e^{-a\left| x-y\right| }u(y)\,{\mathrm d}y\quad \text {for all }x\in \Omega , $$

whose eigenpairs $(\lambda ^j,\xi ^j)$, $j\in {\mathbb N}_0$, can be explicitly computed (cf. [26, Appendix 2]) as follows: If the transcendental equation

• $\tan (\tfrac{aL}{2}\nu )=\tfrac{1}{\nu }$ has the positive solutions $\nu _0<\nu _2<\ldots$,

• $\cot (\tfrac{aL}{2}\nu )=-\tfrac{1}{\nu }$ has the positive solutions $\nu _1<\nu _3<\ldots$,

then $\sigma ({\mathscr {K}})\setminus \left\{ 0\right\} =\left\{ \tfrac{1}{1+\nu _j^2}\in {\mathbb R}:\,j\in {\mathbb N}_0\right\} \subset (0,1)$ with the associate eigenfunctions

$$\begin{aligned} \xi ^j:[-\tfrac{L}{2},\tfrac{L}{2}]&\rightarrow {\mathbb R},&\xi ^j(x)&:= {\left\{ \begin{array}{ll} \cos (a\nu _jx),&{}j\text { is even},\\ \sin (a\nu _jx),&{}j\text { is odd} \end{array}\right. } \quad \text {for all }j\in {\mathbb N}_0. \end{aligned}$$

Note that $\sigma ({\mathscr {K}})$ depends only on the product aL, and not on a and L individually.

Example 6.2

(Laplace kernel in 2d) For $a_i,L_i>0$, $i=1,2$, and the rectangular habitat $\Omega =[-\tfrac{L_1}{2},\tfrac{L_1}{2}]\times [-\tfrac{L_2}{2},\tfrac{L_2}{2}]$ consider the Laplace kernel

$$ k(x,y):=\frac{a_1a_2}{4}e^{-a_1\left| x_1-y_1\right| -a_2\left| x_2-y_2\right| }\quad \text {for all }x,y,\in {\mathbb R}^2 $$

yielding the linear integral operator

$$ {\mathscr {K}}u(x):=\frac{a_1a_2}{4}\int _{-L_1/2}^{L_1/2}\int _{-L_2/2}^{L_2/2}e^{-a_1\left| x_1-y_1\right| -a_2\left| x_2-y_2\right| }u(y)\,{\mathrm d}(y_1,y_2)\quad \text {for all }x\in \Omega . $$

In this setting, we are able to compute the eigenpairs of ${\mathscr {K}}$ via the Fredholm operators ${\mathscr {K}}_iu(x):=\frac{a_i}{2}\int _{-L_i/2}^{L_i/2}e^{-a_i\left| x-y\right| }u(y)\,{\mathrm d}y$, whose eigenpairs $(\lambda _i^j,\xi _i^j)$ with $i=1,2$ and $j\in {\mathbb N}_0$ are known from Example . Thereto, we use the ansatz

$$ \xi ^{j_1,j_2}(x):=\xi _1^{j_1}(x_1)\xi _2^{j_2}(x_2)\quad \text {for all }x=(x_1,x_2)\in \Omega $$

to obtain the countably many eigenpairs $(\lambda _1^{j_1}\lambda _2^{j_2},\xi ^{j_1,j_2})$, $j_1,j_2\in {\mathbb N}_0$, of ${\mathscr {K}}$ from the eigenvalue-eigenvector relation

$$ {\mathscr {K}}\xi ^{j_1,j_2} = ({\mathscr {K}}_1\xi _1^{j_1})({\mathscr {K}}_2\xi _2^{j_2}) = (\lambda _1^{j_1}\xi _1^{j_1})(\lambda _2^{j_2}\xi _2^{j_2}) = \lambda _1^{j_1}\lambda _2^{j_2}\xi ^{j_1,j_2}. $$

Fig. 2

Feigenbaum diagrams for the generalized Beverton–Holt integrodifference (6.6) with parameters (6.7) for $\kappa =1$ (left) and $\kappa =2$ (right). They were obtained from Nyström discretizations based on the trapezoidal rule with 257 resp. 1089 nodes

In order to illustrate how the dynamics of generalized Beverton–Holt integrodifference (6.6) depends on $\alpha >0$, let us choose $\theta _0=3$ and the parameters

$$\begin{aligned} a&=1,\,L=2\text { for }\kappa =1,{} & {} {\left\{ \begin{array}{ll} a_1=a_2=1,\\ L_1=L_2=2 \end{array}\right. }\text { for }\kappa =2,&c&=(\overline{1,2,\tfrac{1}{2}}) \end{aligned}$$

(6.7)

resulting in the habitats $\Omega =[-1,1]$ resp. $\Omega =[-1,1]^2$. This yields the Feigenbaum diagrams displayed in Fig. 2 for the total populations (cf. Remark ) after a transient phase. They indicate that the trivial solution loses its exponential stability at some parameter $\alpha _1^0>0$ and bifurcates into a nontrivial branch $\phi ^1(\alpha )$ of 3-periodic solutions along which, in turn, a period-doubling scenario unfolds. The colors in Fig. 2 encode the generation modulo 3.

Table 1 Theoretical and achieved convergence order in the Nyström approximations for the 3-periodic solution $\phi ^*$ of the logistic integrodifference (6.5) with $L=10$

In the following, we substantiate these experimental observations in more detail (Table 1).

6.2.1 Behavior along the trivial solution

Linearizing the generalized $\theta _0$-periodic Beverton–Holt (6.6) along the trivial solution branch $\phi ^0(\alpha )\equiv 0$ yields the variation equation

$$ u_{t+1}(x)=\alpha \int _\Omega c_tk(x,y)u_t(y)\,{\mathrm d}y\quad \text {for all }x\in \Omega . $$

With the period operator $\Xi _{\theta _0}^0=\alpha ^{\theta _0}c_{\theta _0-1}\cdots c_0{\mathscr {K}}^{\theta _0}$ the spectral mapping theorem yields the Floquet spectrum $\sigma (\Xi _{\theta _0}^0)=\alpha ^{\theta _0}c_{\theta _0-1}\cdots c_0\sigma ({\mathscr {K}})^{\theta _0}$. Hence, there exists a strictly increasing sequence of parameter values $(\alpha _k^0)_{k\in {\mathbb N}}$ such that the Morse index of the trivial solution increases at each $\alpha _k^0$. For $\theta _0=3$ and the parameters (6.7) the dependence of the Floquet multipliers on $\alpha$ is illustrated in Fig. 3.

Fig. 3

Cubic parabolas illustrating the dependence of the six largest Floquet multipliers $\lambda _i^n(\alpha )$ for the generalized Beverton–Holt integrodifference (6.6) with parameters (6.7) along the trivial solution on $\alpha$ for $\kappa =1$ (left) and $\kappa =2$ (right)

The branches of real Floquet multipliers form the cubic parabolas displayed in Fig. 3 whose intersections with the horizontal line equal to 1 yield the critical values $\alpha _k^0$ given in Table 2. Because the spectrum of the Laplace kernels is explicitly known, we can determine those parameter values explicitly.

Table 2 Approximate critical parameter values $\alpha _k^0$ for the generalized Beverton–Holt integrodifference (6.6) with period $\theta _0=3$, parameters (6.7) along the trivial solution for $\kappa \in \left\{ 1,2\right\}$

Table 2 confirms the observation from the Feigenbaum diagrams in Fig. 2 that the loss of stability for the trivial solution occurs for smaller $\alpha$-values when $\kappa =1$.

The explicit knowledge of the Floquet spectrum for Laplace kernels allows us to test our numerical algorithms: We approximate Floquet multipliers using the cyclic operator with the quadrature resp. cubature rules from Appendix B. Choosing again the parameters (6.7) and $\alpha =1$ the resulting error diagrams relating the absolute error to the number of integration nodes are displayed in Fig. 4.

Again the logarithmic scale on both axes in these diagrams yields that the slope of the graphs is an approximation to the convergence order. Their values are given in Table 3. Interestingly, the experimental convergence rate is independent of the integration rule, namely roughly 2 for $\kappa =1$ and about 1 for $\kappa =2$ dimensions. This deficit compared to the theoretically possible values (cf. Proposition ) stated in Appendix B is due to the lacking smoothness of the Laplace kernel, which is merely Lipschitz.

Remark 6.2

(integration over the Laplace kernel) Better convergence rates than those in Table 3 might be obtained by subdividing the integration interval $\Omega :=[-\tfrac{L}{2},\tfrac{L}{2}]$ from Example  into subintervals

$$ \Omega =[-\tfrac{L}{2},x]\cup [x,\tfrac{L}{2}]\quad \text {for all }x\in \Omega _n $$

respectively the integration rectangle $\Omega :=[-\tfrac{L_1}{2},\tfrac{L_1}{2}]\times [-\tfrac{L_2}{2},\tfrac{L_2}{2}]$ from Example  into subrectangles

$$\begin{aligned} \Omega =&\left\{ y\in \Omega :\,x_1\le y_1,x_2\le y_2\right\} \cup \left\{ y\in \Omega :\,x_1\le y_1,x_2\ge y_2\right\} \\&\cup \left\{ y\in \Omega :\,x_1\ge y_1,x_2\le y_2\right\} \cup \left\{ y\in \Omega :\,x_1\ge y_1,x_2\ge y_2\right\} \quad \text {for all }x\in \Omega _n, \end{aligned}$$

over which the Laplace kernels are sufficiently smooth. Of course, these subdivision techniques generalize to further kernels lacking smoothness (cf. [19, p. 30, Table 3.1]).

As a consequence, since higher-order integration rules do not yield an improvement in convergence, we proceed with a composite trapezoidal rule with 257 resp. $1089=(32+1)^2$ nodes for our Nyström discretization.

Fig. 4

Error in the dominant Floquet multiplier for the generalized Beverton–Holt integrodifference (6.6) with parameters (6.7) and $\alpha =1$ along the trivial solution for $\kappa =1$ (left) and $\kappa =2$ (right) with different integration rules

Table 3 Convergence rate of the dominant Floquet multiplier for the generalized Beverton–Holt integrodifference (6.6) with parameters (6.7) and $\alpha =1$ along the trivial solution for $\kappa \in \left\{ 1,2\right\}$

Fig. 5

Solution branches for the generalized Beverton–Holt (6.6) with parameters (6.7) and $\kappa =1$, where green color indicates exponential stability and red indicates unstable solutions: Left: Nontrivial branch $\phi ^1(\alpha )$ of 3-periodic solutions bifurcating from the trivial branch $\phi ^0(\alpha )$ at the critical parameter $\alpha =\alpha _1^0$: $\phi _{3t}^1(\alpha )$ (top), $\phi _{3t+1}^1(\alpha )$ (center) and $\phi _{3t+2}^1(\alpha )$ (bottom). Right: Branches of 6-periodic solutions $\phi ^2(\alpha )$ bifurcating from the 3-periodic solution $\phi ^1(\alpha )$ at the critical parameter $\alpha =\alpha _1^1$: $\phi _{3t}^{2}(\alpha )$ (top), $\phi _{3t+1}^{2}(\alpha )$ (center) and $\phi _{3t+2}^{2}(\alpha )$ (bottom)

Fig. 6

Dependence of the six largest Floquet multipliers $\lambda _i^n(\alpha )$ for the generalized Beverton–Holt integrodifference (6.6) with parameters (6.7) along the nontrivial branch of 3-periodic solutions $\phi ^1(\alpha )$ on $\alpha$ for $\kappa =1$ (left) and $\kappa =2$ (right)

6.2.2 Behavior along the first nontrivial branch

At the critical value $\alpha _1^0$ the trivial branch of the generalized Beverton–Holt integrodifference (6.6) loses its exponential stability and a nontrivial branch $\phi ^1(\alpha )$ of $\theta _0$-periodic solutions bifurcates from $\phi ^0(\alpha )$. We illustrate $\phi ^1(\alpha )$ in Fig. 5 (left).

Linearizing (6.6) along $\phi ^1(\alpha )$ leads to the $\theta _0$-periodic variation equation

$$ u_{t+1}(x)=\alpha \int _\Omega c_tk(x,y)\frac{1-3\phi _t^1(\alpha )(y)^4}{(1+\phi _t^1(\alpha )(y)^4)^2}u_t(y)\,{\mathrm d}y\quad \text {for all }x\in \Omega . $$

It may have the period operator $\Xi _{\theta _0}^0$ with now unknown spectrum. We thus approximate it numerically based on Algorithm 3 and illustrate its dependence on $\alpha$ in Fig. 6.

Based on this diagram, one obtains the critical value $\alpha _1^1$ listed in Table 4, where the eigenvalue branches cross the horizontal $-1$ indicating a period-doubling along the branch $\phi ^1(\alpha )$. While the critical values in Table 2 could be obtained explicitly based on the eigenvalues of the Laplace kernel, we now compute them using bisection while following the solution branch $\phi ^1(\alpha )$ and monitoring the Floquet multipliers.

The corresponding branches of the 6-periodic solutions $\phi ^2(\alpha )$ bifurcating supercritically at $\alpha _1^1$ from $\phi ^1(\alpha )$ for $\kappa =1$ are illustrated in Fig. 5 (right).

Finally, we complement the Feigenbaum diagram in Fig. 2 with the (incomplete) bifurcation diagram for the 3-periodic generalized Beverton–Holt integrodifference (6.6) given in Fig. 7.

Table 4 Approximate critical parameter value $\alpha _1^1$ for the generalized Beverton–Holt integrodifference (6.6) with period $\theta _0=3$, parameters (6.7) along the solution branch $\phi ^1(\alpha )$ for $\kappa \in \left\{ 1,2\right\}$

7 Perspectives

For simplicity reasons and for the sake of immediate implementations, we restricted our analysis to Nyström methods. More advanced techniques yielding superconvergence, in particular when approximating periodic solutions of integrodifference equations, were developed in [2, 3]. Thereto, one first discretizes the state space $C_d$ using discontinuous piecewise polynomials (of degree at most $\gamma$) and then evaluates the remaining integrals using Nyström methods (of degree $\gamma$). This results in a convergence order $2\gamma$ (or higher when the kernel functions are sufficiently smooth).

Fig. 7

Bifurcation diagrams for the generalized Beverton–Holt integrodifference (6.6) with parameters (6.7) for $\kappa =1$ (left) and $\kappa =2$ (right). Exponentially stable solutions are in green, while unstable branches are marked in red

For various problems the kernel function $f_t$ in (1.1) is singular or at least non-smooth along the diagonal $x=y$ (like the Laplace kernels from Examples 6.1 and 6.2). As illustrated in Table 3 this dooms most conventional quadrature schemes to slow convergence. Nevertheless, the kernel typically has a singularity of known type. This allows to develop specific quadrature formulas based on corrections of the trapezoidal rule [4, 28]. For an alternative approach, we also refer to [10].

Appendix

1.1 A Miscellaneous tools

Suppose that $(\Omega ,d)$ is a compact metric space and that (X, d) is a metric space (denoting the metric on both spaces by d should cause no confusion). A family $\mathcal {G}$ of mapping $g:\Omega \rightarrow X$ is called equicontinuous, if for every $\varepsilon >0$ there exists a $\delta >0$, such that for every $g\in \mathcal {G}$ one has the implication

$$ d(x,\bar{x})<\delta \quad \Rightarrow \quad d(g(x),g(\bar{x}))<\varepsilon \quad \text {for all }x,\bar{x}\in \Omega . $$

Lemma A.1

Let a sequence $(g_n)_{n\in {\mathbb N}}$ of continuous functions $g_n:\Omega \rightarrow X$ converge pointwise to a function $g:\Omega \rightarrow X$. If the family $\left\{ g_n\right\} _{n\in {\mathbb N}}$ is equicontinuous, then $(g_n)_{n\in {\mathbb N}}$ converges uniformly to g.

Proof

Let $\varepsilon >0$. Since $\left\{ g_n\right\} _{n\in {\mathbb N}}$ is equicontinuous, there exists a $\delta _0>0$ such that $d(x,\bar{x})<\delta _0$ implies $d(g_n(x),g_n(\bar{x}))<\tfrac{\varepsilon }{3}$ for all $n\in {\mathbb N}$, $x,\bar{x}\in \Omega$.

(I) We show that $g:\Omega \rightarrow X$ is uniformly continuous. For $x,\bar{x}\in \Omega$, the pointwise convergence implies that there exists a $N\in {\mathbb N}$ such that $d(g(x),g_n(x))<\tfrac{\varepsilon }{3}$ and $d(g(\bar{x}),g_n(\bar{x}))<\tfrac{\varepsilon }{3}$ hold. This gives $d(g(x),g(\bar{x}))\le d(g(x),g_n(x))+d(g_n(x),g_n(\bar{x}))+d(g_n(\bar{x}),g(\bar{x}))<\varepsilon$, provided that $d(x,\bar{x})<\delta _0$.

(II) We establish uniform convergence. Since g is uniformly continuous, there is a $\delta \in (0,\delta _0)$ such that $d(x,\bar{x})<\delta$ leads to the inequality $d(g(x),g(\bar{x}))<\tfrac{\varepsilon }{3}$ for $x,\bar{x}\in \Omega$. Let $\mathcal {U}$ be an open cover of $\Omega$ consisting of $\tfrac{\delta }{2}$-neighborhoods. Because $\Omega$ is compact, there exist finitely many points $x_1,\ldots ,x_p\in \Omega$ with $\Omega \subseteq \bigcup _{j=1}^pB_{\delta /2}(x_j)$. The pointwise convergence implies that there exists a $N\in {\mathbb N}$ such that $d(g_n(x_j),g(x_j))<\tfrac{\varepsilon }{3}$ holds for $n\ge N$, $1\le j\le p$. Given $x\in \Omega$, we choose j with $x\in B_{\delta /2}(x_j)$ and consequently the estimate $d(g_n(x),g(x))\le d(g_n(x),g_n(x_j))+d(g_n(x_j),g(x_j))+d(g(x_j),g(x))<\varepsilon$ yields the claim.

1.2 B Quadrature and product cubature rules

Quadrature rules

Let $a<b$ be reals, $n\in {\mathbb N}$, and $h:=\tfrac{b-a}{n}$. The quadrature rules employed in this paper are supposed to satisfy (cf. (5.2))

$$ \left| \int _a^b u-\sum _{j=1}^{q_n}w_ju(\eta _j)\right| \le \frac{c_p}{n^{d_p}}\sup _{x\in [a,b]}\left| u^{(d_p)}(x)\right| \quad \text {for all }u\in C^{d_p}[a,b]. $$

Among them, our implementations rely on the following rules:

• Composite midpoint rule (MID)

  $$ \int _a^bu=h\sum _{j=0}^{n-1}u(a+\tfrac{2j+1}{2}h)+\tfrac{(b-a)^3}{24n^2}u''(\eta ) $$

• Composite trapezoidal rule (TRP)

  $$ \int _a^bu = \tfrac{h}{2}\left( u(a)+\sum _{j=1}^{n-1}u(a+jh)+u(b)\right) -\tfrac{(b-a)^3}{12n^2}u''(\eta ) $$

• Composite Simpson rule (SMP)

  $$ \int _a^bu = \tfrac{h}{6}\left( u(a)+4\sum _{j=0}^{n-1}u(a+\tfrac{2j+1}{2}h)+2\sum _{j=1}^{n-1}u(a+jh)+u(b)\right) - \tfrac{(b-a)^5}{2880n^4}u^{(4)}(\eta ) $$

• Composite Clenshaw-Curtis rule order 6 (CC6): There exists a real $C\ge 0$ such that

  $$\begin{aligned} \int _a^bu = \tfrac{h}{30} \Biggl ( \sum _{j=0}^{n-1}(u(a + \tfrac{2j+1}{2}h \pm \tfrac{h}{2}) + 8(u(a + \tfrac{2j+1}{2}h + \cos (\tfrac{\pi }{4})\tfrac{h}{2}) \\ + u(a + \tfrac{2j+1}{2}h + \cos (\tfrac{3\pi }{4})\tfrac{h}{2})) + 12u(a + \tfrac{2j+1}{2}h)) \Biggr ) + C\tfrac{(b - a)^{7}}{n^{6}}u^{(6)}(\eta ) \end{aligned}$$

• Composite Gauß-Legendre rule order 8 (GL8)

  $$\begin{aligned} \int _a^bu = \tfrac{h}{72} \Biggl (\sum _{j=0}^{n-1}((18 - \sqrt{30})u(a+\tfrac{2j+1}{2}h \pm \sqrt{\tfrac{3}{7} + \tfrac{2}{7}\sqrt{\tfrac{6}{5}}}\tfrac{h}{2}) \\ + (18 + \sqrt{30})u(a+\tfrac{2j+1}{2}h \pm \sqrt{\tfrac{3}{7} - \tfrac{2}{7}\sqrt{\tfrac{6}{5}}}\tfrac{h}{2})) \Biggr ) + \tfrac{(b - a)^{9}}{6945750n^{8}}u^{(8)}(\eta ) \end{aligned}$$

Besides these composite rules, for smooth integrands u also nodes and weights for Gauß-Legendre quadrature of arbitrary order (and its variants) perform well. For their computation we refer to, e.g., [11, pp. 389ff].

Cubature rules

Let $a_1<b_1$, $a_2<b_2$ be reals, $n_1,n_2\in {\mathbb N}$ and $h_i:=\tfrac{b_i-a_i}{n_i}$ for $i=1,2$. For double integrals we use product cubature rules derived from the above quadrature rules, i.e.,

$$ \int _{a_1}^{b_1}\int _{a_2}^{b_2}u(y_1,y_2)\,{\mathrm d}y_2\,{\mathrm d}y_2 \approx \sum _{j_1=1}^{q_{n_1}}\sum _{j_2=1}^{q_{n_2}}w_{j_1}^1w_{j_2}^2u(\eta _{j_1},\eta _{j_2}). $$