Electronic Journal of Qualitative Theory of Differential Equations

Abstract. We present an application of a recently developed algorithm for rigorous integration forward in time of delay differential equations (DDEs) to a computer assisted proof of the existence of several periodic orbits in a DDE obtained by a singular perturbation limit method from the classical logistic map. The proofs are done near the parameter value for which multistability was numerically observed.


Introduction
In this work we are dealing with the delay differential equation (DDE) of the form: x ∈ R. (1.1) Equation (1.1) is a special form of the delayed feedback equation: which is extensively studied in the literature under various assumptions on the feedback function f , both analytically and numerically.Numerical studies show that complex dynamics is common in nonlinear systems of the form (1.2).A famous example was proposed by Mackey and Glass, for which there is an evidence of a series of period doubling bifurcations that apparently lead to a chaotic attractor [17].Losson, Mackey and Longtin numerically studied multistability for Equation (1.1), which was named there as so-called "logistic" DDE, as it is obtained by a singular perturbation limit procedure on the famous logistic map [16].The authors of [16] found that the basins of attraction for different stable periodic orbits may have a fractal structure.
Analytic investigation of dynamics in Equation (1.2) usually proves extremely difficult for arbitrary feedback function f , so many analytical results, despite being very impressive, are proved only for f close to a piecewise constant function.For such systems we can find several important results.Mallet-Paret and Sell proposed to use discrete Lyapunov functionals to prove a Poincaré-Bendixson type theorem in the case of monotone systems close to a step function [18].The conditions for existence and the shape of a global attractor were studied by Krisztin, Walther and Wu for systems with positive monotone feedback i.e. for monotone functions such that x • f (x) > 0 for x = 0, see [8] and references therein.Under the same assumptions, Krisztin and Vas constructed equations for which large-amplitude, slowly oscillatory periodic solutions (LSOPs) exist, such that they revolve around more than one equilibrium [9].A global picture of the attractor was also obtained.Vas showed that a positive monotone feedback function f may be chosen such that the global attractor contains any number of unstable LSOPs.[27].Earlier, Vas also showed that systems with negative monotone feedback (i.e.x • f (x) < 0 for x = 0) can exhibit multistability and that the number of stable periodic orbits can be infinite [26].Lani-Wayda and Walther were able to construct systems of the form ẋ = f (x(t − 1)) for which they proved the existence of dynamics conjugated to a symbol shift (Smale's horseshoe) [10].Srzednicki and Lani-Wayda proved the existence of multiple periodic orbits and the existence of chaos for some tooth-shaped (piecewise affine) periodic function f by the use of the generalized Lefshetz fixed point theorem [11].There are some results for equations with unimodal feedback function (i.e.smooth f that has a unique extremum), such as in Equation (1.1) see [14,15] and references therein.
Most of the mentioned works concentrate on the construction of special kinds of feedback functions for which it is possible to analytically compute the action of a semiflow associated with (1.2).Then, the system may be studied by an application of arguments from the geometric theory of dynamical systems, namely, by construction of suitable Poincaré maps.A possible approach to treat arbitrary equations such as (1.1) is to use rigorous (validated) numerical methods.In recent years, there were many computer assisted proofs of various properties of dynamical systems ranging from ordinary differential equations to (dissipative) partial differential equations, see for example [2,19,21,25,31] and references therein.A computer assisted proof is a computer program which rigorously checks assumptions of abstract theorems.Such results, in the context of continuous dynamical systems, usually require some kind of an algorithm for rigorous integration of the flow forward in time.
In this paper we present an application of the recently developed [22,23] forward in time integration algorithm for rigorous solving of Initial Value Problems (IVPs) of the form: By the rigorous integration we understand a computer procedure which, given some a priori bounds on ψ, produces bounds for the solution x(t) to (1.3) with x| [−1,0] ≡ ψ (we assume the r.h.s.f to be such that this solution is unique).Using the integrator, we rigorously construct images of sets by Poincaré maps, then we show, under some assumptions, that those maps are compact.Finally, we apply Schauder fixed point theorem to obtain existence of periodic orbits.We do not prove that the orbits are attracting, as this would require some C 1 -computations (see [30]) which are for now beyond the scope of our methods.
There are several papers that deal with computer assisted proofs of periodic solutions to DDEs [12,13,28], however their approach is very different from our method.These works transform the question of the existence of periodic orbits into a boundary value problem (BVP).The BVP is then solved by using either Newton-Kantorovich theorem [12,13] or the local Brouwer degree [28].It is clear that rigorous integration in the phase space may be used to obtain a more diverse spectrum of results such as connecting orbits and the existence of chaos.
Recently Bánhelyi et al. made progress towards proving the long-standing Wright's conjecture, by combining analytical and computer-assisted arguments [3].In their paper, the global attractivity of the stationary solution to Wright's equation was shown to be equivalent to the non-existence of slowly oscillating periodic (SOP) solutions with amplitudes greater than a certain constant ε > 0. Next, a procedure to find a set of bounding piecewise constant functions for any SOP solution to the equation was presented.Finally, the bounding functions were used to show that the SOP solution, if it exists, has an amplitude less than ε.The idea of obtaining bounds for a solution in a form suitable for rigorous computations is similar in spirit to our methods, however we use more general representation of the phase space.Moreover, the authors of [3] use the validated numerical integration to obtain the bounds of the solution only over a fixed interval and it is not clear if their method may be adopted to construct Poincaré maps.
This paper is organized as follows.In Section 2, we summarize the algorithm for rigorous construction of Poincaré maps for DDEs.In Section 3, we apply our method to prove the existence of several attracting periodic orbits in (1.1) for various values of parameter λ near 5.81, for which multistability occurs.In Section 4, we apply our method to find numerical approximations of apparently unstable orbits in the system by direct application of Newton algorithm to a nonrigorous version of our integrator.Finally in Section 5, we discuss possible future developments of the rigorous integration methods.

Notation
In this work we use the following notation.For a function f : R → R, by f (k) we denote the k-th derivative of f .By f [k] we denote the term 1  k! • f (k) .In the context of piecewise smooth maps by f (k) (t − ) and f (k) (t + ) we denote left and right one sided derivatives of f w.r.t.t, respectively.
For a given set A by cl(A) and int A we denote the closure and interior of A, respectively (in a given topology e.g.defined by the norm in the considered Banach space).
We will abuse notation and we will write C r = C r ([−τ, 0], R) for the space of all functions of class C r over [−τ, 0] equipped with the supremum norm: For v ∈ R n by π i v for i ∈ {1, 2, . . ., n} we denote the projection of v onto the i-th coordinate.For vectors u, v ∈ R n by u • v we denote the standard scalar product: We will sometimes use the usual notation v i = π i v, but only when the context is absolutely clear, as not to confuse it with the notation of We call A an interval set (a product of closed intervals in R n ).For any A ⊂ R n we denote by hull(A) a minimal interval set, such that A ⊂ hull(A).If A ⊂ R is bounded then hull(A) = [inf(A), sup(A)].For any set X by m(X) we denote the midpoint of interval set hull(X) and by diam(X) the diameter of hull(X).

Set arithmetic and rigorous numerics
By its very nature, the numerical computations done in computers are finite, so only a small subset of real numbers can be stored in an exact manner.We call them representable numbers.Other real values must be approximated to the nearest representable number.The commonly used standard for those operations is IEEE 754 [7].The idea of rigorous interval numerics is to replace operators and elementary functions (+, −, ÷, ×, sin, cos, exp, etc.) with their setvalued versions such that they operate on intervals.The interval counterpart [ ] of a given operator gurantees that the result contains all possible outcomes of applying the operator on any two possible arguments, i.e.
Therefore, by using interval arithmetic in our computations, we are sure that they contain all possible results -they are validated.To be able to implement such interval operations on computers, we need closedness, that is, we require [a 1 , b 1 ][ ][a 2 , b 2 ] to be again an interval with ends being representable numbers.Fortunately, this can be simply done under the IEEE 754 standard [6], and many software packages already exist (see [1] and references therein).We use [4] for this purpose.
To stress the fact that we are using set arithmetic we will often put a variable name in square brackets, e.g., we will write [r] to denote a set in R m .Usually it will happen in formulas used in algorithms.If both variables r and [r] are used simultaneously then usually r represents a value in [r], however this is not implied by default and it will be always stated explicitly.Please note that the notion [r] has no additional meaning and [r] may be simply regarded as a name of another variable.
Further, we will put explicit numerical values into square brackets to stress the fact that in their place we are using an interval with ends being representable numbers that contains the given value.For example [0.3333] should be understand as a representable interval obtained by the rigorous numerical computation of the division operator: [3333, 3333] ÷ [10000, 10000].If the explicit value in a numerical formula is put without brackets then it is in fact a closest representable IEEE floating point number in a given precision.This happen in formulation of some computer-assisted proofs in Sections 3 and 4.

Rigorous forward in time integration and construction of Poincaré maps for DDEs
For the rest of this paper we assume that τ = 1 in (1.3).All computations can be easily redone for any delay τ.
The semiflow ϕ associated with Equation (1.3) is defined as: where x ψ : [−1, a ψ ) → R denotes the solution to the Cauchy problem (1.3)-(1.4),such that the solution exists for all t < a ψ .Under some mild assumptions on the r.h.s.f of (1.3) x ψ is unique for all initial function ψ and the (local) semiflow ϕ is continuous on C 0 ([−1, 0], R) [5].Also, from the classical method of steps, one can obtain the following lemma, which we state without proof.

Lemma 2.1 (Smoothing property). Assume f is of class
In our computer assisted proofs, we follow a standard algorithm to construct Poincaré maps for flows: 1. Select a finite representation of the initial condition x 0 ∈ C 0 , i.e. a set X 0 ⊂ C 0 such that x 0 ∈ X 0 and X 0 can be described with a finite number of parameters.
2. Define a section S in the phase space such that x 0 ∈ S and the semiflow ϕ is transversal to S in the neighbourhood of x 0 .
3. Compute estimates on the transition time to section t S ∈ [t S ] for t S > 0 such that ϕ(t S , x 0 ) ∈ S.
4. Return a (finite representation of) set X t S such that ϕ(t S , X 0 ) ⊂ X t S .
Our rigorous integration scheme will be given by a family of operators Step 4 can be implemented as: where Notice that, the value of q can be obtained by counting the number of iterations of the operator I h before the algorithm detects a crossing of the section.Then [ε 1 , ε 2 ] may be found by application of any bisection algorithm.
In [22,23], a rigorous algorithm for construction of I h for the constant step size h = 1 p for some p ∈ N + .was proposed.Also, under some assumptions, it is possible to obtain a rigorous version of I ε for 0 < ε < h.The original algorithm extensively relies on the smoothing property of DDEs (Lemma 2.1) so from now on we assume that r.h.s.f of Equation (1.3) is "sufficiently smooth" for various expressions to make sense.In fact, we can assume that f ∈ C ∞ as this is usually true in the context of the classical equations defined by a finite composition of elementary C ∞ functions.The "logistic" (Eq.(1.1)) and Mackey-Glass equations are both good examples.
Let us recall basic definitions and results from [23].We assume that integers n ≥ 0 and p > 0 are given and we set h = 1 p .
A (p, n)-representation of g is a pair ḡ = (a, B) such that The graphical presentation of the idea of a (p, n)-representation of some exemplary function is given in Figure 2.1.The following notation is used for any (p, n)-representation ḡ: The term ḡi, [k] is called the (i, k)-th coefficient of the representation and ḡi,[n+1] the i-th remainder of the representation.The interval set B is called the remainder of the representation.We will call the constant M = p • (n + 2) + 1 the size of the (p, n)-representation.The parameters n and p are omitted when known from the context.
The two concepts are strongly linked but, formally, a (p, n)-representation is an object in a finite dimensional space, where a (p, n)-f-set contains elements of an infinite dimensional functional space.
The next lemma follows simply from the definitions of (p, n)-representations and (p, n)-fsets.[23]).Let G, F be (p, n)-f-sets with representations Ḡ and F, respectively.The following statements hold.

•
x [1] (t) =0.125 = x [1] and 1 2 • x (2) = x [2] , respectively.The remainder terms B are presented graphically as blue boxes, while the values of a are presented as blue dots.In e) the coefficients a and B are organized into a matrix that resembles their position w.r.t.indices i and k, e.g.element xi,[k] lies in i-th column and k-th row.For the sake of the presentation, numbers are rounded up to three decimal places.The reconstruction of the (p, n)-f-set represented by x is shown as blue region in d).We see that x is contained in fset( x).An animated version of this figure is available in PDF format in [24].
where x0 , xh are (p, n)-representations, such that ϕ(h, fset( x0 )) ⊂ fset( xh ) and h = 1 p .Also, under assumption that the iteration time is "long enough", they proposed a family of algorithms I ε for 0 < ε < h for which the following is true: We skip the details of the construction here.We would like to note that the condition for a "long enough" iteration is n • p ≤ q.If we denote the transition time q • h + ε in (2.6) by t, then the condition transforms into t ≥ n.By Lemma 2.1 it guarantees that ϕ(t, x) is of class at least C n (we remind that we assumed τ = 1).
Finally, we recall the theorem which guarantee that Poincaré map is well defined and it is a compact operator.Definition 2.7 (Definition 6 in [23]).Let n ∈ N, s : C n → R be a bounded linear operator and a ∈ R. We define a global C n -section as a hyperplane: (2.7) Any convex, bounded subset S a ⊂ S a is called a local C n -section (or simply a section).Let W ⊂ C n be a convex open set such that A section S a is said to be transversal on W if We will refer to (2.8) as the transversality condition.
Definition 2.8 (Definition 7 and Theorem 10 in [23]).Let ω ≥ n + 1, S be a section, V ⊂ C 0 and assume there exist t 1 , t 2 such that: • S is transversal on W, Then, the transition map to the section S after (at least) ω is well defied and given by: where t S (x) is a unique time such that ϕ(t S (x), x) ∈ S (see Theorem 10 in [23]).If V ⊂ S, then the map P ≥ω will be called the Poincaré return map on the section S after ω.
Note that the condition on the transition time given by ω is n + 1 instead of n (as compared to Theorem 2.6).This extra iteration time guarantees that the derivative of P(x) w.r.t.t is well defined and of class at least C n .
The main result that allows us to use the rigorous bounds in the actual proofs is as follows.
Theorem 2.9 (Theorem 11 in [23]).Let S be a C n -section and let V ⊂ S be bounded in C n .Assume

Computer assisted proofs of attracting periodic orbits in the logistic equation
We start this section with recall of Schauder's fixed point theorem [20,29].We will use it to prove existence of several periodic orbits to (1.1) for various values of parameter λ.
Theorem 3.1 (Schauder's fixed point theorem).Let X be a Banach space, let V ⊂ X be a nonempty, convex, bounded set and let P : V → X be a continuous mapping such that P(V) ⊂ K ⊂ V and K is compact.Then the map P has a fixed point in V.
We will be working in the space X = C n .We will use routines described in the previous Section to obtain bounds on the set K in Theorem 3.1 for some a priori estimated set V ∈ C n p ∩ C n .Theorem 2.9 will guarantee compactness of operator P, thus compactness of the set K.
Since we are using (p, n)-representations with n = 4 and p = 32 to describe sets in the phase space, it would not be feasible to present data from the proofs directly, as the number of coefficients is equal to the size of the (p, n)-representation: M = 193.Also, it would be impossible to provide reasonable initial data by hand.Therefore, we discuss shortly how we automatically obtain bounds on the set of initial functions V on a suitably chosen section S a .

(p, n)-sections
We follow notation of [23].Since we are using (p, n)-representations to describe functions in C n , it is advisable to define sections in the phase space in such a way that it would be easy to rigorously check if x ∈ S for all functions represented by a given (p, n)-f-set (or its (p, n)-representation).A straightforward way is to require s in Definition 2.7 to depend only on representation coefficients xi, [k] for k ≤ n.This motivates the following definition.Definition 3.2.Let li,k ∈ R for (i, k) ∈ C = {1, . . ., p} × {0, . . ., n} ∪ {(0, 0)}.We assume that at least one li,k is not equal to zero.We define a bounded linear operator s(x) by: A section S a = {x ∈ C n , s(x) = a} (see Definition 2.7) is called a (p, n)-section.A subset S a ⊂ S a is called a local (p, n)-section.A vector l ∈ R M−p+1 (where M is the size of the (p, n)-representation) such that π I(i,k) l = li,k for (i, k) ∈ C is called the representation of a (p, n)-section.
Note that s(x) ≡ l • x for the (p, n)-representation of a function x.

Selection of the initial conditions
To find good initial conditions, we will use an approximate flow φ defined by: where Φ h , Φ ε denotes "truncated" versions of the rigorous integration algorithms I h and I ε , respectively.Algorithms Φ h and Φ ε are called truncated, as they do not take into consideration the effect of remainder terms B of (p, n)-representations (see Section 2.2 of [23] for details).Despite this, we hope that a periodic solution x to (3.2) would be close enough to the true solution of (1.3), so that a small neighbourhood [ x] of x can be selected as V in Theorem 3.1.
We start the selection of the initial function by computing φ(T 0 , x) with x = A • sin(B • t) + C for some A, B, C ∈ R (discussed later in Section 4) and T 0 ∈ R + large enough so that the solution xT 0 stabilizes as discussed in [16].Then, we construct a Poincaré map P with a (p, n)-section defined by s(x) = π I(0,0) x and a = π I(0,0) xT 0 .For the map P we apply Newton algorithm to find x0 ∈ R m such that P( x0 ) − x0 ≤ 10 −13 .The vector x0 will be the middle point of the (p, n)-representation of the initial set V.
In the next step we choose a (p, n)-section such that it minimizes the difference between transition time t S for different x ∈ V.A heuristic algorithm to find such section exists (see Section 3.4 in [23]).
In order to demonstrate the idea behind the heuristic algorithm, let us consider for a while an ODE of the form: Let x 0 be a periodic orbit of period T of the flow ϕ induced by (3.3).Then f (x 0 ) is the right eigenvector of the monodromy matrix A = ∂ϕ ∂x (T, x 0 ) with eigenvalue λ = 1.Let l be the left eigenvector of A corresponding to λ = 1, i.e. l T • A = l T .If x 0 is hyperbolic then the left and right eigenvectors of A are uniquely defined up to a multiplier and we can choose l such that Let l ⊥ be the tangent space to l, i.e. l ⊥ = {x ∈ R n : l • x = 0}).Then, by direct computation, we get: Therefore, in the ODE setting, choosing left eigenvector of the monodromy matrix A to define (linear) section guarantees that it is transversal to the flow and the return time is constant in the first order approximation.We use this fact to obtain a good candidate for a section by computing the left eigenvector corresponding to eigenvalue 1 of the matrix ∂ φ ∂x (T, x 0 ).Let l ∈ R M−p+1 be this vector.We define the matrix This is a non-singular matrix.Let C be a matrix obtained from the procedure of orthonormalisation of Ĉ by the Gram-Schmidt algorithm, so that C is an orthonormal base for the phase space of R m .We choose the initial set V to be: for some r 0 ∈ R m , π 1 r 0 = 0. Notice, that the set V is contained on the (p, n)-section defined by l as a hyperplane in the phase space.Coefficients of r 0 are chosen experimentally to follow an exponential law in k, i.e. r i,[k] 0 ≈ d • q k for some q > 1 and d ∈ (0, 1).They are further refined during the computations to obtain better initial bounds.A method of selecting the middle point x0 for various values of the parameter λ in (1.1) is discussed in details in Section 4.

Computer assisted proofs
Before we formulate theorems, we once more recall the problem of exact representation of numbers in computer programs.Below we state theorems with floating-point (real) numbers presented up to four digits after the decimal point, so they are not exactly the ones which the software for computer-assisted proof produces.We have rounded them in such a way that they are valid bounds (upper or lower where appropriate) for the values produced in rigorous numerical procedures.Please also remember, that if we write [5.8], then we mean by that a minimal interval with end points being representable numbers such that it contains 5.8.
The computer program and initial data used to obtain estimates in theorems may be found online [24].There, one can find also a detailed instruction on compiling and running the programs.The compilation requires CAPD library [4], as we use it in underlying rigorous interval numerical procedures.Source codes were tested on a standard laptop PC with Intel ® Core ™ I7-2860QM CPU (2.50 GHz), 16 GB RAM under 64-bit Linux operating system (Ubuntu 12.04 LTS) and C/C++ compiler gcc version 4.6.3.In the sequel, we will refer to this configuration as a reference machine.On this configuration the computation time for each of the proofs is approximately one minute.
Before we state the theorems, we note that it is straightforward to check that ξ 0 = 0 and ξ 1 = λ−1 λ are the stationary solutions to (1.1) for λ = 0.Moreover, if x is a bounded solution to (1.1) with initial condition x 0 and x 0 ≡ x T for some T > 0, then x may be extended to the whole R where it will satisfy (1.1) for any t.By x(R) we denote the image of the extended function over R.
Now we can proceed to state the theorems, which are proved with computer assistance.
Proof.We choose the (32,4)-representation of the phase space and we set λ = [5.5] in (1.1).We choose the (p, n)-section and we construct representation of the (p, n)-f-set of initial conditions V as described in Section 3.2.For those initial conditions we perform the rigorous construction of the image of Poincaré map P ≥ω as described in Section 2. Let W be the (p, n)f-set enclosing the image of P ≥ω obtained from the integration procedure, i.e.P(V) • 209 = q > (n + 1) so by the Theorem 2.9 the Poincaré map is well defined and compact on V.
Next, the program checks inclusion W ⊂ V (see output of the program for details).By Schauder theorem we get the existence of a fixed point in the (p, n)-f-set V which corresponds to a periodic orbit (as ξ 0 , ξ 1 ∈ V).
Since x is periodic, it suffices to compute lower and upper estimates on the value of x(t) over an interval containing basic period to obtain estimates of inf t∈R x(t) and sup t∈R x(t).The program verifies that x(t), The data for this proof can be found in a .zipfile at [24] under logistic/5.5/0_B.initial.The estimates obtained on the reference machine can be found under logistic/5.5/0_C.proof.
The approximate orbit which defines the middle point of the initial set V in Theorem 3.3 is presented in Figure 3 Proof.We proceed in the same manner as in the proof of Theorem 3.3.We present those numbers that are needed to verify assumptions: • s( ẋ) > 0.2167 for x ∈ W ∩ C n , which guarantees transversality, • 0.0163 < inf t∈R x(t) < 0.0195, 1.1449 < sup t∈R x(t) < 1.1461, ξ 1 ∈ [0.8260, 0.8261] The data for this proof can be found in a .zipfile at [24] under logistic/5.75/1_B.initial.The estimates obtained on the reference machine can be found under logistic/5.75/1_C.proof.
The approximate orbit which defines the middle point of the initial set V in Theorem 3.4 is presented in Figure 3 Proof.We proceed in the same manner as in the proof of Theorem 3.3.We present only numbers needed to verify assumptions: • s( ẋ) > 0.4290 for x ∈ W ∩ C n , which guarantees transversality, The data for this proof can be found in a .zipfile at [24] under logistic/5.8/1_B.initial.
The estimates obtained on the reference machine can be found under logistic/5.8/1_C.proof.
The approximate orbit which defines the middle point of the initial set V in Theorem 3.5 is presented in Figure 3.3.

Further numerical investigations of the logistic equation
As we have stated in Section 2, we have generated initial conditions by a very long, nonrigorous integration of (1.1) for initial functions of the form A • sin(B • t) + C.  We call the solutions a-type and b-type.The a-type solutions tend to the periodic orbit with for larger values of parameter λ.For each value of the parameter λ ∈ {5.5 + 0.1 • i : i ∈ 0, . . ., 3} ∪ {5.75 + 0.01 • i : i ∈ 0, . . ., 10} we computed 100-th iterates of a 0 and b 0 : a 100 := Phi 100 h (a 0 ), b 100 := Phi 100 h (b 0 ).Then, Newton algorithm was applied to a 100 and b 100 to obtain final candidates a and b such that P(x) − x < 10 −13 for x ∈ {a, b}.
We ran an automatic procedure to verify (non-rigorously) if the orbit is attracting or repelling by computing (approximate) eigenvalues of the Jacobian of P at x ∈ {a, b}.The resulting orbits are shown in Tables 4.1, 4.2 and 4.3.With green box we have marked solutions for which we were able to get a computer assisted proof of existence.An a-type periodic solution for λ = 5.5 for which we have managed to get the computer assisted proof is shown in Figure 3.1.Solutions presented in Tables 4.1 and 4.2 were classified as attracting, where in Table 4.3 we have gathered initial conditions obtained from Newton algorithm that are apparently unstable (numerical computations show that Jacobian of P at x has one unstable eigenvalue).
We were not able to obtain a computer assisted proof in the case of a-type periodic solution when λ > 5.7.The reason for this is that the approximate monodromy matrix of ϕ at x have not contained an eigenvalue 1 in its spectrum.Thus, we were not able to select a good representation of the (p, n)-section.The observed dominant approximate eigenvalue was always complex with the magnitude ≈ 1.We have tried to repeat computations with larger (p, n)-representations (a (128,4)-representation and a (32,10)-representation), but without success.However, for the (128,4)-representation, the magnitude of the imaginary part of the dominant eigenvalue was much less than in the case of (32,4)-representation.This suggest that the lack of eigenvalue 1 may be only an artefact of numerical computations.We plan to carry out the computations for even larger representations of the phase space, but due to the time constraints, we need first to optimise the procedure of computation of the monodromy matrix as suggested in [22].With the current implementation we estimate that computations on (256, 4)-representation would take more than a year.

Summary
In this work we have applied a recently developed rigorous procedure for forward in time integration of DDEs to prove existence of several periodic orbits to the "logistic" equation (1.1), for which it is difficult to obtain such results by analytical means.
There is much work left to be done.First of all, we were not able to prove the existence of two coexisting attracting periodic orbits, observed numerically in [16].A possible approach Table 4.2: Numerically attracting periodic orbits in (1.1) for higher values of the λ parameter.Orbits were obtained by integration of the initial conditions (4.1) and (4.2) (upper and lower row respectively) forward in time up to t = 100 (i.e. 100 • τ).Then, the pictures were obtained by plotting the numerical solution over time t ∈ [0, 20].It seems that for λ > 5.83 a chaotic attractor appeared.However, using Newton algorithm, we were able to find apparently unstable periodic orbit as seen in Table 4. Table 4.3: Some apparently unstable periodic orbits (with one unstable direction) in (1.1) for various λ.
would be to increase the size of the (p, n)-representation, but for this we need to optimise our algorithms.We have already made some suggestions in [22].
Further, we do not have tools to prove that the existing orbits are in fact attracting.The tools for this are already present for ODEs [30] and in the future work we will try to translate them to the DDE setting.
An additional future goal is to establish existence of unstable periodic orbits.We have shown that our algorithms may be used to find good numerical approximations to such orbits (examples in Table 4.3).The tools needed in the context of infinite dimensional flows to prove existence of unstable hyperbolic periodic solutions are available [31] and it should be easy to incorporate them in our algorithms.

Figure 3 . 1 :Theorem 3 . 4 .
Figure 3.1: Plot of the solution from Theorem 3.3.The dots mark the shift by a next full delay, where the associated number is the corresponding time.

Figure 3 . 2 :Theorem 3 . 5 .
Figure 3.2: Plot of the solution from Theorem 3.4.The dots mark the shift by a next full delay, where the associated number is the corresponding time.

Figure 3 . 3 :
Figure 3.3: Plot of the solution from Theorem 3.5.The dots mark the shift by a next full delay, where the associated number is the corresponding time.

Table 4 .
1: Numerical investigations of logistic equation (1.1) with various λ.With a green box we have marked periodic solutions we were able to prove rigorously.We see that the double oscillation in the first orbit (upper row) shrinks.Probably, it is the reason why we were not able to prove their existence, as they get closer and closer to the unstable periodic orbit, as seen in Table4.3.