Electronic Journal of Qualitative Theory of Differential Equations

We show, by way of an example, that numerical bifurcation tools for ODE yield reliable bifurcation diagrams when applied to the pseudospectral approximation of a one-parameter family of nonlinear renewal equations. The example resembles logisticand Ricker-type population equations and exhibits transcritical, Hopf and period doubling bifurcations. The reliability is demonstrated by comparing the results to those obtained by a reduction to a Hamiltonian Kaplan–Yorke system and to those obtained by direct application of collocation methods (the latter also yield estimates for positive Lyapunov exponents in the chaotic regime). We conclude that the methodology described here works well for a class of delay equations for which currently no tailor-made tools exist (and for which it is doubtful that these will ever be constructed).


Introduction
In the field of population dynamics renewal equations (i.e., integral equations of Volterra type) abound [1,29,32,35,37]. In fact they form the core of the theory of structured population models, see [20,22] and also [21,33,43,48]. The mathematical theory for such equations is extensive as well as comprehensive [16,19,31]. And yet there are relatively few simple examples for which a detailed analysis reveals interesting nonlinear dynamics.

The equation and the numerical results
Let h : R → R be a C 1 function with h(0) = 0 and h(z) ≥ 0 for z ≥ 0. We are interested in the renewal equation where A : [0, +∞) → R takes nonnegative values. We want A to have compact support and to be reflection symmetric with respect to the midpoint of the support. In order to facilitate reference to Kaplan and Yorke [34], we choose the midpoint to be τ = 2. More precisely, from now on we assume that A is a bounded, nonnegative and measurable function satisfying A(τ) = 0 for τ > 4, In particular, for the numerical simulations we combine the step function A(τ) = 1 2 , 1 ≤ τ ≤ 3, 0, 0 ≤ τ < 1 or τ > 3, (2.4) with either the quadratic nonlinearity h(z) = γz(1 − z) (2.5) or the exponential Ricker-type nonlinearity h(z) = γze −z , (2.6) where γ > 0 is a free parameter, which we will consider as the bifurcation parameter. The main advantage of the pseudospectral method considered in [2] is that a generic nonlinear delay equation is reduced to a low-dimensional system of ODE, which can be written in a convenient way using the right-hand side of the original equation and an additional block of linear equations that depends only on the discretization mesh and on the delay (hence, it is independent of the specific problem under a suitable scaling of time). The dynamical properties of the approximating system can be easily analyzed through well-established software for numerical bifurcation of ODE (e.g., MATCONT or AUTO, to name a couple). Hence the approach yields an interesting resource for users without a solid numerical background who are interested in gaining insights in the properties of delay equations. More in this regard will be said in Section 5.
The results in this paper are obtained with MATCONT [15,42], a MATLAB continuation and bifurcation package. The outcome of the bifurcation analysis with respect to γ of the pseudospectral approximation of (2.1) with the quadratic nonlinearity (2.5) is plotted in Figure 2.1a: we observe the nontrivial steady state bifurcating into a periodic orbit through a Hopf bifurcation, followed by the starting of a period doubling cascade. Some periodic trajectories for parameter values in the regime of the period doubling cascade are plotted in Figure 2.1b (panels A-C, with reference to the corresponding values in Figure 2.1a).
An important advantage of the approximation through a system of ODE is that, in principle, the trajectories of the corresponding initial value problems can be simulated by using ordinary numerical solvers (we used ode45 and ode23s, available in MATLAB). The simulations produced for parameter values in the range of the period doubling cascade (panels A-C in Figure 2.1b) agree with those approximated by MATCONT. In addition, the simulation of the initial value problem allows us to investigate the dynamical behavior beyond the range of the period doubling cascade. In particular, the system seems to show chaotic behavior for larger values of γ, see, e.g., panel D of Figure 2.1b for γ = 5, as well as Figure 2.3. The numerical bifurcation analysis of the exponential nonlinearity (2.6) has similar features as the quadratic case. The bifurcation diagram with respect to log γ is plotted in Figure 2.2a. The nontrivial equilibrium undergoes a Hopf bifurcation, followed by the starting of a period doubling cascade. Some trajectories are plotted in Figure 2.2b for parameter values in the range of the period doubling cascade (panels A-C, with reference to the corresponding values in Figure 2.2a). Again, the simulation of the initial value problem for larger values of γ shows the appearance of chaotic behavior, visible, e.g., in panel D of Figure 2.2b for log (γ) = 4.3.
There are, unfortunately, too few values of γ at period doubling points to check whether they display the Feigenbaum universality [28].
The results described above are, for the quadratic nonlinearity, visualized in the movie quadratic.mp4 * available as supplementary material. Combining the pseudospectral discretization with a different method described in Section 5 and the methodology of [9] we were able to numerically compute Lyapunov exponents for the quadratic nonlinearity. The results are depicted in Figure 2.3, and they clearly reveal the signature of chaotic behavior for values of γ beyond 4.55. Summarizing, the numerical results for both the quadratic and exponential nonlinearities support the conjecture, formulated 30 years ago by O. D. but never verified, of a period dou-  bling bifurcation which breaks the original symmetry of the periodic solutions arising from the Hopf bifurcation (see Figure 2.1b or Figure 2.2b, as well as Sections 3-4). For larger parameter values, chaos is observed. This behavior is in line with the results obtained in [49] through numerical simulations of initial value problems.
The results presented in this section demonstrate how the pseudospectral approximation allows to gain valuable information about the dynamical properties of nonlinear renewal equa-tions. The convergence of the method when increasing the number of approximating ODE is proved in [2] for the steady states and their stability. There, the authors also conjecture the same convergence properties for periodic solutions and more complicated dynamical objects together with their stability and relevant bifurcations. The proof of convergence for periodic solutions and their stability is currently ongoing work. For this reason, in Section 5 we shall validate the numerical outcome presented in this section in two ways: first, by exploiting the analytical results guaranteed by the symmetry property (2.2); second, by comparing it with the numerical results obtained through a more sophisticated and specific approach, which is based on the principle of linearized stability.

The renewal equation: some analytical aspects
In the following, we shall study equation (2.1) under the assumptions (2.2)-(2.3), and we shall pay particular attention to solutions of period 4. Note that "solutions of period 4" includes solutions of minimal period 4 k for some integer k ≥ 1. When we mean minimal period 4, we shall say so explicitly. Moreover, we shall pay particular attention to solutions with some symmetry. Our first result establishes a connection between symmetry and period 4. Proof. From (2.1) it follows that proving the claim since z is even.
Note that the renewal equation (2.1) is translation invariant: if z is a solution, so is z(· + t 0 ) for every choice of t 0 . To make a solution even, one has to choose the translate appropriately, if possible at all.
Our second result gives the connection between period 4 and symmetry in translationinvariant form for the specific kernel (2.4).
Lemma 3.2. Assume (2.4) and let z be a solution of (2.1) of period 4. Then there exists θ ∈ R such that Proof. Hence, i.e., z(t) + z(t + 2) is constant. If we call 2θ this constant value, we arrive at (3.1).
The next, rather important, result indicates the relation between the symmetric periodic solutions of (2.1) and a particular delay differential equation.
In particular, z is a solution of (2.1) if and only if the additional condition is satisfied.
In the next section we shall take (3.4) and (3.7) as the key equations for the construction of 4-periodic solutions of (2.1) when (2.4) holds. In the remainder of this section we derive Hopf bifurcation conditions for the general class of kernels A satisfying (2.2)-(2.3).
The linearized version of (2.1) in z is where p = h (z) is considered as a parameter. The corresponding characteristic equation is obtained by substituting ξ(t) = e λt and reads for an integer k ≥ 1.
For positive p, the rightmost root of (3.9) is real and hence roots on the imaginary axis are not crucial for stability. Therefore we are especially interested in negative values of p.
If (2.4) holds, p 2l = +∞ and showing that l needs to be even in order to obtain negative p. When l = 2m we can narrow (3.10) and (3.11) down to with corresponding period showing that 4 is always a period, but it is the minimal period only for the primary branch m = 0. Consider the particular nonlinear functions analyzed in Section 2. For the quadratic case (2.5), we first of all observe that the nontrivial steady state z = 1 − 1 γ is positive if and only if γ > 1 and it is stable for γ slightly larger than 1 (see [16] for the dynamical system perspective and the principle of linearized stability). Linearization in this steady state yields (3.8) with p = 2 − γ, which, combined with (3.13), allows to recover (3.14) Observe now that γ 0 = 2 + π 2 is approximated correctly by the numerical bifurcation analysis of Section 2, see Figure 2.1a. In the exponential case (2.6) the nontrivial steady state is z = log γ and the corresponding linearization is (3.8) with p = 1 − log γ, so that Hopf bifurcations occur at parameter values Observe again that log (γ 0 ) = 1 + π 2 is approximated correctly by the numerical bifurcation analysis of Section 2, see Figure 2.2a. The book [18] provides a detailed Hopf bifurcation analysis that applies also to (2.1), see [16]. Alternatively, see [17]. Also see [43, Section VI 3.2] for a formula giving the direction of the Hopf bifurcation.
We claim that the period remains 4 (possibly divided by an integer of the form 4m + 1) along the m-th branch. We sketch a proof.
The key idea is to consider, independently of the Hopf bifurcation approach described above, equation (2.1) as a fixed point problem for an integral operator H defined by on a space of even 4-periodic functions, e.g., continuous functions provided with the supremum norm. To verify that H maps even 4-periodic functions to even 4-periodic functions, define and verify that H(R(z)) = R(H(z)) if (2.2)-(2.3) hold and z has period 4. It follows that In this setting one can apply the Crandall-Rabinowitz theory for local bifurcation at a simple eigenvalue and the Krasnosel'skii theory for global bifurcation, see [36]. The result of uniqueness modulo translation, derived from the dynamical systems approach, implies that a suitably chosen translate equals the even fixed point of period 4. It follows that the period does not vary along the branch.
Note that, for A given by (2.4), the approach of the next section yields an alternative proof that the period is fixed along the branch.

Periodic solutions of delay differential equations generated by Hamiltonian systems (d'après Kaplan & Yorke)
In this section, in view of (3.4), our aim is to find 4-periodic solutions of the delay differential where the continuous function f : R → R is assumed to be odd, i.e., In fact, we shall guarantee that x has period 4 by requiring the symmetry property (3.3), i.e., x(t + 2) = −x(t). Then, (x, y) is a 4-periodic solution of the planar Hamiltonian systeṁ Conversely, if (x, y) is a periodic solution of (4.5) such that (4.4) holds, then both (4.1) and (4.3) are satisfied, provided that x and y have mean zero.

Remark 4.2.
Given the symmetry properties explained below, the condition that the mean of x equals zero can alternatively be formulated as: the orbit corresponding to (x, y) encloses the origin, see Figure 4.1b.
Because of (4.2), the Hamiltonian system (4.5) is equivariant with respect to the dihedral group of order 8 corresponding to the symmetry of a square. Let be the counterclockwise rotation over π 2 and let s = 0 1 1 0 be the reflection in the line y = x. Then r 4 = e = s 2 , for e the identity matrix, and the 8 elements of the group are e, r, r 2 , r 3 , s, rs, r 2 s, r 3 s. Assume that the graph of f has the features depicted in Figure 4.1a, i.e., assume that f ∈ C 1 and f (0) > 0, Then, system (4.5) has, in addition to (0, 0), the eight steady states (±x, 0), (0, ±x), (±x, ±x), that form a group orbit.
The qualitative phase portrait for (4.5) with f satisfying (4.6) is depicted in Figure 4.1b. In Theorem 4.1 we focused on closed orbits of period 4. This is, of course, a severe restriction. However, when the right-hand side of the delay equation involves a multiplicative parameter γ (recall (2.5) and (2.6)), we can eliminate that parameter in the Hamiltonian system by a scaling of time. It follows that any periodic solution (with the right symmetry) of the Hamiltonian system without parameter yields a periodic solution of the delay equation of exactly period 4 by choosing appropriately the value of the parameter γ. In other words, we can parametrize branches of 4-periodic solutions of the delay differential equation with the parameter specifying the members of a one-parameter family of closed orbits of a Hamiltonian system that does not involve γ.
Proof. For better readability, in the proof we do not write explicitly the dependence on ξ 0 and m when they are clear. The key point is that the symmetry and the fact that T(ξ 0 ) is the minimal period imply that We first verify (4.3): Next we verify (4.7):ẋ (t) = γξ(γt) Note that, when the period is a monotone function of ξ 0 , we can invert ξ 0 → γ m (ξ 0 ) and hence parametrize the branch with γ. But in general it is rather difficult to establish monotonicity of the period, see [13,30] and the references therein. At this point we want to mention that for (4.7) much is known about secondary bifurcations along the primary branch of symmetric 4-periodic solutions, see the survey paper [52] and [23,53]. Now recall Theorem 3.3. We want to apply Theorem 4.3 in the context of Theorem 3.3. This means that we have an additional parameter, viz. θ, and an additional equation, viz. (3.7). We can view (3.7) as determining θ as a function of ξ 0 , thus retaining the parametrization of a branch of 4-periodic solutions by ξ 0 .
Therefore, we now consider the renewal equation with parameter γ explicitly introduced and not incorporated in h (thus deviating from the earlier convention (2.1) with (2.5) or (2.6)). We define f θ by (3.5) and assume that for a certain relevant range of values of θ the Hamiltonian systeṁ has a one-parameter family of periodic orbits with minimal period T(ξ 0 , θ) intersecting the positive ξ-axis at ξ 0 . We rewrite (3.7) as Defining we can write (4.10) as and (try to) solve the latter for θ as a function of ξ 0 . Since and If h(θ) = θg(θ), we can divide out the trivial root θ = 0 and arrive at Because of (4.11), we may also write (4.13) as which expresses that, for ξ 0 = 0, θ should be a nontrivial steady state of (4.8).
Note that (4.14) For the linearized equation (3.8), we have p = γh (θ) = γ g(θ) + θg (θ) , allowing us to rewrite (4.14) as which reduces to exactly as in (3.13). Thus we see that ξ 0 = 0 corresponds, as expected, to Hopf bifurcations. The discussion above, together with Theorem 3.3 and Theorem 4.1, allows us to characterize the 4-periodic solutions of the renewal equation (2.1) arising from the Hopf bifurcation as solutions of the Hamiltonian system (4.9) and (4.11)-(4.12). In the next section we shall exploit this relation to compare the numerical periodic solutions obtained in Section 2 for (2.6) with the periodic solutions obtained by simulating the Hamiltonian system. Although the same reasoning applies to (2.5), we stress that in this case the periodic solutions are explicitly available (since the Hamiltonian system is the harmonic oscillator), see Appendix A. Hence, the quadratic case allows to compare the numerical solutions to the exact ones.

On reliability and accuracy
The pseudospectral discretization [2] starts from the reformulation of the initial value problem for a generic nonlinear delay equation as an equivalent abstract Cauchy problem, i.e., an infinite-dimensional ODE. The state of the corresponding dynamical system, which belongs to a suitable space of functions, is approximated with a polynomial of a chosen degree M, obtained by collocating the infinite-dimensional ODE on a mesh of points. In principle, the higher M, the smaller the error. In the sequel, we refer to M as discretization index.
As already mentioned, there are two main reasons that make this approach particularly suitable for, say, pragmatic users: the straightforward formulation of the approximating ODE system in terms of the original right-hand side, and the possibility of studying the bifurcation properties by means of well-established software for ODE. For these reasons we shall also informally refer to it as "pragmatic-pseudospectral" method (PP). However, a rigorous proof of convergence of the approximation to periodic solutions for increasing M is still ongoing research (included in the PhD program of F. S.). In this section, we provide numerical support for the validity of the results of Section 2 by way of two different reference settings: on the one hand, we exploit the Hamiltonian equivalence and the analytical results described in Sections 3-4 and Appendix A; on the other hand, we consider a more sophisticated numerical approach.
This second approach relies on the principle of linearized stability. As such, it requires, first, a method to compute or approximate the solution of interest and, second, a method to recover the stability information about the system linearized around this solution. Usually, the latter is obtained from the spectrum of suitable operators. Concretely, we use this approach for steady states and periodic solutions, and their stability properties (in Section 2, we briefly pushed the approach even further to demonstrate chaotic dynamics by way of positive Lyapunov exponents). Concerning the renewal equations treated here, steady states can be found analytically, so that one only needs to approximate the associated characteristic roots. These can be obtained as eigenvalues of the corresponding infinitesimal generator, via the pseudospectral techniques in [3,4]. Both reduce the generator to a finite-dimensional matrix, whose eigenvalues are taken as approximations. As for periodic solutions, instead, their analytical expression is unattainable in general (the quadratic case (2.5) represents an exception, the branch of periodic solutions arising from the Hopf bifurcation being explicitly known, see Appendix A). To remedy, we extend to renewal equations the ideas of the collocation approach proposed in [25] for delay differential equations. See also [39][40][41] for a modern and abstract treatment. With a numerical approximation of the periodic solution at hand, again through pseudospectral collocation, we construct a finite-dimensional version of the monodromy operator associated to the linearized problem, and we use its eigenvalues as approximations of the characteristic multipliers. The latter method is inspired by [6], and it implicitly assumes that Floquet theory can be extended to renewal equations, despite the current lack of a detailed elaboration. The above techniques used to compute periodic solutions and their stability indicators (as well as their extensions to chaotic trajectories) are part of ongoing research (included in the PhD program of D. L.), targeted to coupled renewal and delay differential equations in the spirit of [4].
The above description suggests that this alternative, more sophisticated technique is less suitable for users who are not confident, or less expert, in numerical analysis, since it requires different specific tools for the approximation of the solution first and of the spectrum of the linearized system next. For this reason we shall informally refer to it as "expert-pseudospectral" method (EP). At the same time, since this approach is specifically directed to delay equations, it allows to obtain considerably more accurate results and shorter computational times compared with PP. Therefore it makes sense to use EP to validate the outcome of PP.
Notice that parts of EP are very recent and still under development, and their details are indeed unpublished. Here we do not go into the theoretical justification. Neither do we describe the implementation or the details of the observed numerical convergence. We intend to present a convergence proof as well as the results of systematic testing on various problems in future work. For now we consider the matching of the results of PP and EP as convincing evidence that both are correct, at least when applied to the present problem.
In the spirit described above, we present in Figure 5.1 the periodic trajectories approximated by EP and PP compared with the reference trajectory. The latter is taken to be the analytical solution in the quadratic case (see Appendix A) and the solution computed via simulation of the Hamiltonian system (4.9) with (4.11)-(4.12) in the exponential case. Notice the remarkable overlap of the trajectories, obtained after suitable translation as necessitated by the different phase conditions imposed by the two approaches.
To make the comparison more quantitative, we study the convergence of both the trajectories and their periods in terms of increasing M. The latter is the degree of the collocation polynomial for both PP and EP. Recall from Section 3 that the concerned trajectories have exact period 4. The errors are plotted in Figures 5.2 and 5.3 for, respectively, the quadratic and the exponential case. We can clearly observe similar features: the approximations are converging faster than polynomially in 1/M. This constitutes the main strength of the pseudospectral approximation, also known as spectral accuracy (see, e.g., [50]): it guarantees high accuracies with low discretization indices. Moreover, it is evident how EP allows to obtain more accurate results than PP with lower computational effort, reaching, e.g., the machine precision for already approximately M = 20 in the quadratic case. For discretization indices up to roughly 10, the order of accuracy of the two methods is comparable, but, when increasing M further, PP is affected by the limits imposed by the continuation tolerance we set in MATCONT (TOL = 10 −8 for the quadratic case and TOL = 10 −6 for the exponential case).
The last test we present aims at confirming the accuracy of EP in the approximation of the periodic solutions and of their multipliers. We restrict the analysis to the quadratic nonlinearity. Along the branch arising from the Hopf bifurcation (panel A in Figure 2.1b), the reference solution is given by the analytical expression obtained in Appendix A. In the period doubling range (panels B-C in Figure 2.1b), the reference solutions are computed with M = 80, so that they can be considered sufficiently accurate. Moreover, we test the numerical approximation of the exact multiplier 1, which is always present due to translation invariance (see, e.g., [11]). The results reported in Figure 5.4 confirm the spectral accuracy of EP. They also show another important property of the pseudospectral approximation: the larger the discretization interval (in this case, the period), the slower the convergence. Indeed, we can observe that at least twice the nodes are necessary to reach the same accuracy after every period doubling. This is reasonably supported by standard collocation and interpolation results [10,51]: the error depends on both the length of the interval (in our case, the period) and bounds on higher derivatives (in our case, related to the number of oscillations per period). This, in general, makes the continuation of the period doubling cascade rather difficult beyond the first few doubling points.

Discussion and outlook
The paper [2] announced new prospects for the numerical bifurcation analysis of delay equations. As dreams do not always come true, it is crucial to address test problems in order to ascertain the merit for the practice of applied mathematics. The current paper serves as a first test case. We are happy to conclude that the methodology worked even better than expected: we were able to accurately compute (with only moderate computational effort and standard tools) a nontrivial one-parameter bifurcation diagram for a scalar renewal equation. This is all the more gratifying since at present, as far as we know, no tailor-made tools for numerical bifurcation analysis of renewal equations exist. So the approach reported here seems to be the one and only option. And, once again, it worked well! Subgroups of the authors, reinforced by others, plan to elaborate both the PP and the EP approach (as described in the preceding section) and to test them on various other problems, including equations with infinite delay. A more detailed exposition of the computation of Lyapunov exponents is also on the "to do" list.
Concerning the special class of renewal equations studied here, we plan to perform a two-parameter bifurcation analysis of h(z(t − τ)) dτ (6.1) with special attention to the limit → 0. Note that, formally, (6.1) reduces to the map z(t) = γh(z(t − 2)) with continuous, rather than discrete, time t. Here we expect to derive inspiration from [12] and the work of J. Mallet-Paret and R. D. Nussbaum, see, e.g., [38] and the references given there.

Acknowledgements
The

A Explicit solutions for a quadratic nonlinearity
Let us consider the equation with parameter γ > 0. We make the ansatz that is, we aim to derive conditions on θ, A and k such that (A.2) yields a solution of (A.1). Since showing that odd l correspond to negative values of γ and are therefore of no relevance for us. So we require that l = 2m and find Hopf bifurcations at (3.14). For more results in this spirit we refer to [24].