# Computer-assisted validation of the existence of periodic orbits in the
Brusselator system

Jakub Banaśkiewicz†,⋆, Piotr Kalita⋆, Piotr Zgliczyński⋆ †Faculty of
Mathematics and Computer Science, Jagiellonian University, ul. Łojasiewicza 6,
30-348, Kraków, Poland  ⋆AGH University of Science and Technology, al.
Mickiewicza 30, 30-059 Kraków, Poland<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

We investigate the Brusselator system with diffusion and Dirichlet boundary
conditions on one dimensional space interval. Our proof demonstrates that, for
certain parameter values, a periodic orbit exists. This proof is computer-
assisted and rooted in the rigorous integration of partial differential
equations. Additionally, we present the evidence of the occurrence of period-
doubling bifurcation.

The work of all three authors was supported by National Science Center (NCN)
of Poland under project No. UMO-2016/22/A/ST1/00077. The research of JB for
this publication has been supported by a grant from the Antropocen Priority
Research Area under the Strategic Programme Excellence Initiative at
Jagiellonian University. The work of PK was also partially supported by NCN of
Poland under project No. DEC-2017/25/B/ST1/00302 and by Ministerio de Ciencia
e Innovación of Kingdom of Spain under project No. PID2021-122991NB-C21.

##  1\. Introduction

In this paper, we study the dynamics of the Brusselator system with diffusion,
described by the following initial and boundary value problem consisting of
two mutually coupled partial differential equations.

$\begin{cases}u_{t}=d_{1}u_{xx}-(B+1)u+u^{2}v+A\sin(x)\;\text{for}\;(x,t)\in(0,\pi)\times(0,\infty),\\\
v_{t}=d_{2}v_{xx}+Bu-u^{2}v\;\text{for}\;(x,t)\in(0,\pi)\times(0,\infty),\\\
u(t,x)=v(t,x)=0\;\text{for}\;(x,t)\in\\{0,\pi\\}\times(0,\infty),\\\
u(0,x)=u^{0}(x),\ v(0,x)=v^{0}(x)\ \textrm{for}\ x\in(0,\pi).\end{cases}$
(1.1)

The system models autocatalytic reactions in form

$\displaystyle A\rightarrow u,$ $\displaystyle 2u+v\rightarrow 3u,$
$\displaystyle B+u\rightarrow v+D,$ $\displaystyle u\rightarrow E.$

The coefficients $A,B$ in the system are given and they correspond to
densities of the substances $A$ and $B$ in the above reactions. The unknown
functions $u(t,x)$ and $v(t,x)$ describe the densities of two substances which
we also denote by $u$ and $v.$ The system (1.1) is an extension of the planar
Brusselator ODE (1.2) where the substances $u$ and $v$ are homogeneously
spread in the domain. The terms $u_{xx}$ and $v_{xx}$ represent the diffusion
of substances. The coefficients $d_{1},d_{2}$ are the corresponding diffusion
rates.

If we drop the diffusion and the dependence on the variable $x$ in the term
$A\sin(x)$ from the system (1.1), we obtain the following planar ODE

$\begin{cases}u^{\prime}=-(B+1)u+u^{2}v+A\ \ \text{for}\ \ t\in\mathbb{R},\\\
v^{\prime}=Bu-u^{2}v\ \ \text{for}\ \ t\in\mathbb{R}.\end{cases}$ (1.2)

The planarity of the above system implies that the invariant sets consist of
fixed points, periodic orbits, and heteroclinic connections between them. In
fact, it is known that in (1.2), there can exist an attracting periodic orbit
that arises via the Hopf bifurcation [10, Theorem 3]. The analytical results
about the Brusselator system of PDEs with diffusion are limited. In the
article [23] the existence of the global attractor for the Brusselator system
on the 3-dimensional domain is proved. While this global attractor is known to
exist, the question about its structure, which pertains to the understanding
of the problem dynamics, remains unanswered. Some partial analytical results
about the dynamics are available for Neumann boundary conditions, where the
homogeneous steady state is known from the solutions of the corresponding ODE
(1.2). In this case, one can linearize the system of PDEs in its vicinity.
Such results are available for example in [5] and [20]. The case with
Dirichlet conditions, which we consider, appears to be much more challenging.
In [3], the stability analysis of the steady state was carried out for
Dirichlet non-homogeneous conditions. In such case there exists a nonzero
homogeneous steady state. However, this type of analysis is not possible for
the problem we are dealing with, as the homogeneous steady state does not
exist (it would have to be equal to zero). Based on numerical observations,
the system (1.1) possesses a periodic orbit for some range of parameters
$d_{1},d_{2},A,B$. We anticipate that this periodic orbit arises from a
mechanism similar to the one known for the planar ODE (1.2), namely through a
Hopf bifurcation.

We deal with the apparent impossibility of obtaining purely analytical results
on the periodic orbit existence by using the computer assisted techniques.

Specifically, we perform a computer-assisted proof of the following theorem.

######  Theorem 1.1.

For the parameters $d_{1}=0.2$, $d_{2}=0.02$, $A=1$, and $B=2$, the
Brusselator system has a time-periodic orbit.

Figure 1. The numerical approximation of the time-periodic orbit obtained in
Theorem 1.1. Blue and orange surfaces correspond to plots of $u(t,x)$ and
$v(t,x)$ respectively.

It is worth noting that while numerical simulations suggest that this orbit
from Theorem 1.1 appears to be attracting, it has not been rigorously proven
to be so.

The results on the existence of periodic orbit for the Brusselator were also
obtained recently, together with the proof of the Hopf bifurcation, also using
the computer assisted techniques, in the paper [1]. The approach employed
there is based on the Newton–Kantorovich method. The author demonstrates that
a specific Newton-type operator has a fixed point, enabling him to establish
the existence of a periodic orbit. This class of methods has been successfully
applied for many problems governed by PDEs (see for example [2, 4, 18, 19]).
The monograph [17] contains detailed description and up to date overview of
these methods together with numerous applications. In our proof of Theorem 1.1
we are using different method. Namely, our approach is based on the algorithm
of rigorous forward integration of the dissipative systems. This method was
developed in articles [25, 26] and applied there for the Kuramoto–Shivasinski
equation. Similarly, as in the above works, we express the solution in terms
of the Fourier series. Specifically, as the solutions $u(t,x)$ and $v(t,x)$
satisfy the Dirichlet boundary conditions on $[0,\pi]$, they are represented
as the sum of sine Fourier series

$u(t,x)=\sum_{i=1}^{\infty}u_{i}(t)\sin(ix),\quad
v(t,x)=\sum_{i=1}^{\infty}v_{i}(t)\sin(ix).$ (1.3)

The ability to work with all Fourier coefficients facilitates a
straightforward and efficient integration of the Brusselator system. It is
noteworthy, however, that in Theorem 4.2, we obtain the same periodic orbit as
in [1], cross-validating the accuracy of the two distinct methods. Other
approaches of representing the solution are also possible. In the work [12]
the solution to the Burgers equation was represented by the first order finite
element basis together with the estimates on the norms in the Sobolev spaces.

Note that the constant term $A$ from the ODE (1.2) has been replaced by the
term $A\sin(x)$ in our PDE extension of the Brusselator system. In principle,
we could consider the PDE version of the problem with the constant term $A$,
that is the following system

$\begin{cases}u_{t}=d_{1}u_{xx}-(B+1)u+u^{2}v+A\;\text{for}\;(x,t)\in(0,\pi)\times(0,\infty),\\\
v_{t}=d_{2}v_{xx}+Bu-u^{2}v\;\text{for}\;(x,t)\in(0,\pi)\times(0,\infty),\\\
u(t,x)=v(t,x)=0\;\text{for}\;(x,t)\in\\{0,\pi\\}\times(0,\infty),\\\
u(0,x)=u^{0}(x),\ v(0,x)=v^{0}(x)\ \textrm{for}\ x\in(0,\pi).\end{cases}$
(1.4)

In such a case, we observe that $u_{xx}(t,x)=-\frac{A}{d_{1}}\not=0$ for the
boundary points $x\in\\{0,\pi\\}.$ and for every $t\in\mathbb{R}^{+},$ which
means that compatibly conditions are not met. Let us represent $u$ in term of
the sine Fourier series (1.3). Assume that the series
$\sum_{i=1}^{\infty}|u_{i}(t)|i^{2}$ is convergent. Then, we can differentiate
the Fourier series twice, and we obtain

$u_{xx}(t,x)=\sum_{i=1}^{\infty}-u_{i}(t)i^{2}\sin(ix)$

Therefore $u_{xx}(t,x)=0$ for $x\in\\{0,\pi\\},$ which is not true, and hence,
by contradiction, the series $\sum_{i=1}^{\infty}|u_{i}(t)|i^{2}$ has to
diverge. Therefore, the Fourier series cannot converge to the solution fast.
This unwelcome effect does not occur when we consider term $A\sin(x)$ instead
of $A$ in the system (1.1), because the function $Asin(x)$ on $[0,\pi]$ is a
restriction of a smooth, odd and $2\pi$ periodic function.

The novelty of the present paper is the proof of the periodic orbit existence
for the Brusselator systems using the rigorous forward integration techniques.
This integration algorithm is, according to our knowledge, applied by us for
the first time for a system of PDEs: in our case two mutually coupled
nonlinear parabolic PDEs with polynomials of order $3$ in the nonlinear term.
We underline that the rigorous integration scheme which we use is the same as
in [25, 26], where it was used for the Kuramoto–Sivashinski equation with odd-
periodic boundary conditions. We show its applicability for the problem with
higher degree of nonlinearity: is our case the two equations of the system are
coupled through the cubic term, while the nonlinearity in the
Kuramoto–Sivashinsky equation is quadratic. The key concept which makes it
possible for the integration scheme to work is the same in our case as in [25,
26]. Namely, the dissipativity of the leading linear operator together with
appropriate a priori estimates for the nonlinearity, which is of lower order,
allow the linear terms to dominate over the nonlinear ones at appropriately
high modes in the Fourier expansion. This allows us to treat the tail of the
Fourier expansion uniformly, by controlling a polynomial decay of the
coefficients in every time step, cf. Lemma 2.7. Our techniques hold potential
for wider applicability. To this end, in Section 5, we establish estimates on
the convolution of sine and cosine Fourier series, which can be utilized to
calculate nonlinear terms for a general dissipative system with polynomial
nonlinearities in one spatial dimension. These estimates are a crucial
component in rigorous integration algorithms for such systems. A further
innovation is that in our theoretical results, which ensure the validity of
the algorithm, we don’t require the use of Galerkin projections of the
solution. Instead, we work directly with the solution of the PDE system on the
level of abstract theorems as can be seen in Sections 2 and 3. This simplifies
our assumptions and hence it makes the results more accessible to the
dissipative PDE community. Additionally, we examine the limitations of our
algorithm for PDEs with nonlinearities that do not meet the compatibility
condition on the boundary: this is the case even with the constant $A$ in
place of $A\sin(x)$ in the equation for $u_{t}$ in the Brusselator system
(1.1). We illustrate this issue by considering the problem governed by the
diffusive logistic equation (2.16) as an example.

In Section 4, we observe that for a sufficiently large parameter $B$, the
system exhibits slow-fast behavior, as expected for the Brusselator system and
known in the ODE case. This effect is stronger for higher Fourier modes. To
demonstrate this, we establish the existence of periodic orbits for parameters
$d_{1}=0.2,\;d_{1}=0.02,\;A=1,\;B=2+\frac{i}{10}$ for $i\in{0,\dots,11}$.
Figures 3 and 4 show that some of these orbits exhibit slow-fast behavior.

In [1] Arioli observed a period doubling bifurcation, a phenomenon that cannot
occur in the planar ODE (1.2). Thus, the dynamics of (1.1) is expected to be
more complicated than that of the planar ODE (1.2). Although we do not
rigorously prove the bifurcation, we also show that the minimal period of the
found orbits approximately doubles with a small increase in the parameter $B$,
as seen in Theorem 4.3. In this range of parameter $B$ there should also exist
one unstable orbit with the period that is not doubled. To prove the existence
of this kind of orbit we could use the rigorous solver for the Brusselator
system together with the concept of $h$-sets and covering relations (see for
example [22, Section 2.1] and [26, Section 10.2]). Such approach has been
successfully used to prove the existence of unstable periodic orbits before
[26, Theorem 45, Theorem 46]).

Other nontrivial dynamics of the Brusselator system, such as the existence of
2-dimensional attracting tori and chaos, were numerically investigated in [8].
Our numerical observations support the existence of 2-dimensional attracting
tori for small diffusion parameters, although a rigorous proof of their
existence remains an open problem. There are many avenues for further research
on this topic. Numerical simulations indicate that the periodic orbit
established in Theorem 1.1 is attracting. However, providing a rigorous
computer-assisted proof of this observation is challenging, as it requires a
rigorous $C^{1}$ calculation, i.e. the integration of the variational equation
for the Brusselator system.

The structure of the article is as follows. In Section 2 we describe the
algorithm of rigorous integration for dissipative equations. In Section 3 we
describe the computer assisted proof of Theorem 1.1. In Section 3.2, we
address the algorithm for computing the Poincaré map and prove Theorem 3.2,
which pertains to the fixed point of this map which corresponds to the
periodic orbit of the system. In the remaining part of Section 3 we describe
the validation of the assumptions of this theorem for the Brusselator system.
Section 4 contains numerical and rigorous results for various parameters of
the Brusselator system. Finally, in Section 5, we provide results on the
algebra of infinite series, which is utilized in the algorithms.

The code which was used in computer assisted proofs is published at GitHub [7]
and based on CAPD library [6, 15].

##  2\. Algorithm of integration

In this section we present our version of the technique of integration for
infinite dimensional dissipative systems proposed in [25], where it has been
used for the Kuramoto–Shivashinsky equation. This method relies of the
rigorous integration of a differential inclusion and it can be used for many
dissipative problems in mathematical physics. We discuss it in the abstract
setting but some details will be specified for the Brusselator system. The
other approach based on the automatic differentation is presented in [22]. We
summarize the content of this section. We start, in Section 2.1, with the
formulation of the abstract problem for which the algorithm can be applied,
and in Section 2.2 we discuss its realization for Brusselator PDEs. In Section
2.3 we briefly describe the goal of the algorithm and its main steps. The sets
of states and their representation are discussed in Section 2.4, and the way
to compute the nonlinearities present in the system on those sets is described
in Section 2.5. In Section 2.6, we outline the steps for determining the
enclosure and provide a justification for its correctness. The algorithm of
evolution of sets is discussed in Section 2.7.

###  2.1. Abstract Problem

Let $H$ be a real Hilbert space with the scalar product $\left<{.},{.}\right>$
and $Y$ a be a Banach space which continuously and densely embeds in $H$, that
is $\left\|{x}\right\|_{H}\leq C\left\|{x}\right\|_{Y}$ for $x\in Y.$ We
assume that $\\{e_{1},e_{2},\ldots\\}$ is an orthogonal basis of $H$ such that
$e_{i}\in Y$ for every $i\in\mathbb{N}.$

For a given $x\in H$ by $x_{i}$ we will denote the Fourier coefficient
$x_{i}=\frac{\left<{e_{i}},{x}\right>}{\left<{e_{i}},{e_{i}}\right>}.$ We
consider following problem

$\begin{cases}\frac{d}{dt}x(t)=Lx(t)+f(x(t))=F(x(t)),\\\
x(0)=x^{0},\end{cases}$ (2.1)

where $L$ is a diagonal operator such that $Le_{i}=\lambda_{i}e_{i}$ where
$\lambda\in\mathbb{R}\setminus\\{0\\}$, which generates a $C^{0}$ semigroup
$e^{tL}:Y\to Y$. We have $e^{Lt}e_{i}=e^{\lambda_{i}t}e_{i}.$ We assume that
$x^{0}\in Y$, and $f:Y\to Y$ is a continuous mapping. For a given $x\in Y$, we
will use the notation $f_{i}(x)$ and $F_{i}(x)$ for $\frac{\langle
f(x),e_{i}\rangle}{\left<{e_{i}},{e_{i}}\right>}$ and
$f_{i}(x)+\lambda_{i}x_{i}$, respectively.

The following lemma provides the criteria for local in time existence and
uniqueness of mild solutions to problem (2.1).

######  Lemma 2.1.

Assume that

  1. (A1)

For some $C_{1}>0$ there holds $\left\|{e^{tL}}\right\|_{Y}\leq e^{C_{1}t}.$

  2. (A2)

For every $R>0$ there exists $C(R)>0$ such that for every $x,y\in Y$ such that
$\left\|{x}\right\|_{Y},\left\|{y}\right\|_{Y}\leq R$ there holds
$\left\|{f(x)-f(y)}\right\|_{Y}\leq C(R)\left\|{x-y}\right\|_{Y}.$

Then for every initial data $x^{0}\in Y$ there exists the unique local in time
solution to problem (2.1) understood in the following sense

$x(t)=e^{Lt}x^{0}+\int_{0}^{t}e^{L(t-s)}f(x(s))\,ds,$ (2.2)

where the equality in $Y$ is supposed to hold for every $t\in[0,T]$, where $T$
may depend on $x_{0}$.

Instead of proving Lemma 2.1, we will show more general result which implies
it. Namely, we can replace (A1) and (A2) with more general conditions. The
following result generalises Lemma 2.1 to the case when is $f$ is a continuous
map from $Y$ to $Y^{1},$ where $Y^{1}$ is a Banach space such that $Y$ is
continuously embedded in $Y^{1},$ and $Y^{1}$ is continuously embedded in $H.$
This more general result is useful, for example, if the nonlinear term in the
problem depends not only on the value of the solution but also on the values
of its spatial derivatives, which is the case for the Burgers or
Kuramoto–Shivashinsky equations.

######  Lemma 2.2.

Assume that

  1. (B1)

For some $C_{1},C_{2}>0$ and every $t\in[0,\infty)$ there holds
$\left\|{e^{tL}}\right\|_{Y}\leq C_{1}e^{C_{2}t}.$

  2. (B2)

For every $R>0$ there exists $C(R)>0$ such that for every $x,y\in Y$ such that
$\left\|{x}\right\|_{Y},\left\|{y}\right\|_{Y}\leq R$ there holds
$\left\|{f(x)-f(y)}\right\|_{Y^{1}}\leq C(R)\left\|{x-y}\right\|_{Y}.$

  3. (B3)

The semigroup $e^{tL}$ can be extended to the $C^{0}$ semigroup on $Y^{1}$.
There exist constants $C_{3},C_{4}>0$ and $\gamma\in[0,1)$ such that for every
$z\in Y^{1}$ and $t\in(0,\infty)$ there holds
$\left\|{e^{tL}z}\right\|_{Y}\leq
C_{3}e^{C_{4}t}\frac{1}{t^{\gamma}}\left\|{z}\right\|_{Y^{1}}$.

Then, for every initial data $x^{0}\in Y$, there exists a unique time-local
solution to problem (2.1) in the following sense:

$x(t)=e^{Lt}x^{0}+\int_{0}^{t}e^{L(t-s)}f(x(s))\,ds,$

where the equality holds in $Y$ for all $t\in[0,T]$, where $T$ may depend on
$x_{0}$.

###### Proof.

For a given $x^{0}\in Y$ and $\delta>0$ consider the set

$S_{\delta}:=\\{y\in C([0,\delta];Y)\colon y(0)=x^{0}\;\text{and for every}\ \
t\in[0,\delta]\text{ we have }\left\|{y(t)}\right\|_{Y}\leq
1+C_{1}\left\|{x^{0}}\right\|_{Y}\\},$

and define the mapping $T:C([0,\delta];Y)\to C([0,\delta];Y)$ by the formula

$T(y)(t)=e^{Lt}x^{0}+\int_{0}^{t}e^{L(t-s)}f(y(s))ds.$

The space $C([0,\delta];Y)$ is equipped with the norm
$\sup_{t\in[0,\delta]}\left\|{y(t)}\right\|_{Y}$. We have

$\displaystyle\left\|{T(y)(t)}\right\|_{Y}$ $\displaystyle\leq
C_{1}e^{C_{2}t}\left\|{x^{0}}\right\|_{Y}+\int_{0}^{t}\left\|{e^{L(t-s)}f(y(s))}\right\|_{Y}ds$
$\displaystyle\leq
C_{1}e^{C_{2}t}\left\|{x^{0}}\right\|_{Y}+C_{4}\int_{0}^{t}e^{(t-s)C_{3}}\frac{1}{(t-s)^{\gamma}}\left\|{f(y(s))}\right\|_{Y_{1}}ds$
$\displaystyle\leq
C_{1}e^{C_{2}t}\left\|{x^{0}}\right\|_{Y}+C_{4}e^{tC_{3}}\int_{0}^{t}\frac{1}{(t-s)^{\gamma}}ds\sup_{s\in[0,t]}\left\|{f(y(s))}\right\|_{Y_{1}}$
$\displaystyle\leq
C_{1}e^{C_{2}t}\left\|{x^{0}}\right\|_{Y}+\frac{t^{1-\gamma}}{1-\gamma}C_{4}e^{tC_{3}}\left(C(R)\sup_{s\in[0,t]}\left\|{y(t)}\right\|_{Y}+\left\|{f(0)}\right\|_{Y_{1}}\right),$

where $R=1+C_{1}\left\|{x^{0}}\right\|_{Y}.$ If we pick $\delta>0$ such that

$e^{C_{2}\delta}\leq
1+\frac{1}{2C_{1}\left\|{x_{0}}\right\|_{Y}},\quad\text{and}\quad\delta^{1-\gamma}e^{\delta
C_{3}}\leq\frac{1-\gamma}{2C_{4}\left(C(R)R+\left\|{f(0)}\right\|_{Y_{1}}\right)}$

we have that $T(S_{\delta})\subset S_{\delta}.$ We have also

$\displaystyle\left\|{T(y_{1})(t)-T(y_{2})(t)}\right\|_{Y}\leq
C_{4}e^{tC_{3}}C(R)\frac{t^{1-\gamma}}{1-\gamma}\sup_{s\in[0,t]}\left\|{y_{1}(s)-y_{2}(s)}\right\|_{Y}.$

If we take $\delta$ such that

$\delta^{1-\gamma}e^{\delta C_{3}}\leq\frac{1-\gamma}{2C_{4}C(R)},$

the mapping $T$ is also a contraction in the set $S_{\delta}.$ From Banach
fixed point theorem we that $T$ has a unique fixed point which is a solution
to (2.1). ∎

######  Lemma 2.3.

Assume $(A1)-(A2)$ or $(B1)-(B3)$. For every $x^{0}\in Y$ there exist
$t_{max}(x^{0})\in(0,\infty]$ such that for the unique solution
$x:[0,t_{max}(x^{0}))\to Y$ of (2.1) the interval $[0,t_{max}(x^{0}))$ is
maximal interval of existence of this solution. Additionally if we consider
the set

$\Omega:=\left\\{(t,x)\in\mathbb{R}^{+}\times
Y:t\in\left[0,t_{max}(x)\right)\right\\},$

then the function $\varphi:\Omega\to Y$ given by formula

$\varphi(t,x^{0})=e^{Lt}x^{0}+\int_{0}^{t}e^{L(t-s)}f(x(s))\,ds.$ (2.3)

defines a local semigroup.

######  Remark 2.1.

Let $x:[0,\tau]\to Y$ satisfy (2.2). Assume $(A1)-(A2)$ or $(B1)-(B3)$. Then
for every $i\in\mathbb{N},$ the Fourier coefficients $x_{i}$ satisfy the
following non-autonomous ODE

$\frac{dx_{i}}{dt}(t)=\lambda_{i}x_{i}(t)+f_{i}(x(t)),\quad\text{for every
$t\in(0,\tau)$}.$

###### Proof.

If $x(t)$ satisfies (2.2) then for every $i\in\mathbb{N}$ we have

$x_{i}(t)=e^{t\lambda_{i}}x^{0}_{i}+\int_{0}^{t}e^{(t-s)\lambda_{i}}f_{i}(x(s))ds.$
(2.4)

Observe that if $(A1)-(A2)$ or $(B1)-(B3)$ hold, then for every
$i\in\mathbb{N}$ the function $f_{i}:Y\to\mathbb{R}$ is locally Lipschitz.
Hence, in the formula (2.4) the function under integral is continuous.
Therefore we can differentiate this formula and get the ODE for $i$-th Fourier
coefficient. ∎

###  2.2. The Brusselator system.

We discuss how to represent the Brusselator system in the abstract framework
presented in the previous section. We use the notation $L^{2}$ for
$L^{2}(0,\pi)$ equipped with the norm
$\left\|{u}\right\|_{L^{2}}=\sqrt{\int_{0}^{\pi}u(x)^{2}\,dx}$ and $C_{0}$ for
$\\{u\in C([0,\pi])\,:\ u(0)=u(\pi)=0\\}$ equipped with the norm
$\|u\|_{C_{0}}=\max_{x\in[0,\pi]}\\{|u(x)|\\}$, and we consider the following
two product spaces: the Hilbert space

$H=L^{2}\times L^{2}\ \ \textrm{with the norm}\ \
\left\|{(u,v)}\right\|_{H}^{2}=\left\|{u}\right\|_{L^{2}}^{2}+\left\|{v}\right\|_{L^{2}}^{2}\
\ \textrm{for}\ \ (u,v)\in H$

and the Banach space

$Y=C_{0}\times C_{0}\ \ \textrm{with the norm}\ \
\left\|{(u,v)}\right\|_{Y}=\max\\{\left\|{u}\right\|_{C_{0}},\left\|{v}\right\|_{C_{0}}\\}\
\ \textrm{for}\ \ (u,v)\in Y.$

In the space $H$ system of functions $\\{e_{k}\\}_{k=1}^{\infty}$ defined in
following way

$e_{2k-1}=(\sin(kx),0)\quad e_{2k}=(0,\sin(kx))\ \ \textrm{for}\ \
k\in\\{1,2,\ldots\\}.$

is the orthogonal basis. For $u,v\in L^{2}$ we denote by $u_{k}$ and $v_{k}$
the $k$-th coefficients in the Fourier expansion in the sine series, of $u$
and $v$ respectively. We define the operator

$L:D(L)\supset Y\to Y\ \ \textrm{as}\ \
L(u,v)=(d_{1}u_{xx}-(B+1)u,d_{2}v_{xx}),$

where $D(L)=\\{(u,v)\in H^{1}_{0}\times H^{1}_{0}\colon\ L(u,v)\in Y\\}$. The
operator $L$ defines a $C^{0}$ semigroup on $Y$, denoted by $e^{tL}$, cf. [11,
Proposition 2.6.7 and Theorem 3.1.1]. Observe that
$Le_{2k-1}=-(d_{1}k^{2}+B+1)e_{2k-1}$ and $Le_{2k}=-d_{2}k^{2}e_{2k},$ so
$\\{e_{k}\\}_{k=1}^{\infty}$ are the eigenfunctions of $L.$ We define
$f(u,v)=(u^{2}v+A\sin(x),Bu-u^{2}v).$ We can write the Brusselator system
(1.1) as the following abstract problem

$\begin{cases}\frac{d}{dt}(u(t),v(t))=L(u(t),v(t))+f(u(t),v(t)),\\\
(u(0),v(0))=(u^{0},v^{0}).\end{cases}$ (2.5)

We apply Lemma 2.1 to the above system, which gives the following result.

######  Theorem 2.4.

For every $(u^{0},v^{0})\in Y$ there exists a function
$(u,v):[0,t_{max}(u_{0},v_{0}))\to Y$ which is the unique solution to the
(2.5), satisfying the Duhamel formula

$(u(t),v(t))=(u^{0},u^{0})+\int_{0}^{t}e^{L(t-s)}f(u(s),v(s))ds.$

###### Proof.

We demonstrate that for the Brusselator problem (2.5) assumptions (A1) and
(A2) hold and we can use Lemma 2.1. Namely

$\left\|{e^{Lt}}\right\|_{Y}\leq 1,\text{ for every }t\in[0,\infty),$

which follows from the maximum principle for the heat equation. Furthermore

$\left\|{u^{2}v-\bar{u}^{2}\bar{v}}\right\|_{C_{0}}\leq\left\|{u}\right\|^{2}_{C_{0}}\left\|{v-\bar{v}}\right\|_{C_{0}}+\left\|{u+\bar{u}}\right\|_{C_{0}}\left\|{\bar{v}}\right\|_{C_{0}}\left\|{u-\bar{u}}\right\|_{C_{0}}\
\ \textrm{for}\ \ u,v\in C_{0}.$

So for every $R>0$ there exist $C(R)$ such that for every
$\left\|{(u,v)}\right\|_{Y},\left\|{(\bar{u},\bar{v})}\right\|_{Y}\leq R$

$\left\|{f(u,v)-f(\bar{u},\bar{v})}\right\|_{Y}\leq
C(R)\left\|{(u,v)-(\bar{u},\bar{v})}\right\|_{Y},$

which concludes the proof. ∎

From Remark 2.1 and the formula for expanding the expression $u^{2}v$ in terms
of sine Fourier series, we obtain the following lemma.

######  Lemma 2.5.

Let $(u(t),v(t))$ solve the system (1.1). For every $k\in\mathbb{N}$ there
holds

$\begin{cases}\frac{d}{dt}u_{k}=-u_{k}(d_{1}k^{2}+1+B)+N(u,v)_{k}+A\delta_{1k},\\\
\frac{d}{dt}v_{k}=-v_{k}d_{2}k^{2}+u_{k}B-N(u,v)_{k},\end{cases}$ (2.6)

where

$N(u,v)_{k}=\frac{1}{4}\sum_{i_{1}+i_{2}+i_{3}=k}u_{|i_{1}|}u_{|i_{2}|}v_{|i_{3}|}\text{sgn}(-i_{1}i_{2}i_{3})\text{
with $i_{1},i_{2},i_{3}\in\mathbb{Z}\setminus\\{0\\}$}.$ (2.7)

Finally, we present the result that states the fact that the space of
functions with only odd Fourier coefficients being nonzero for functions $u$
and $v$ is forward-invariant and corresponds to the space of functions that
are symmetric with respect to the point $x=\frac{\pi}{2}$.

######  Proposition 2.2.

Space $W:=\\{(u,v)\in Y:u_{2i}=v_{2i}=0,\;\text{for}\;i\in\mathbb{N}\\}$ is
invariant for system (1.1). Specifically, if $(u_{0},v_{0})\in W$ then
$u(t,\frac{\pi}{2}+x)=u(t,\frac{\pi}{2}-x)$ and
$v(t,\frac{\pi}{2}+x)=v(t,\frac{\pi}{2}-x)$ for every
$t\in[0,t_{max}(u_{0},v_{0}))$ and almost every $x\in[0,\frac{\pi}{2}].$

######  Remark 2.3.

In the proof of Theorem 2.4 we verify that assumptions (A1) and (A2) hold for
the Brusselator system. If the lower order nonlinearity depends on the
derivatives of the unknown, then we need (B1)–(B3). Indeed, consider the
Kuramoto–Sivashinsky equation

$u_{t}=-\nu u_{xxxx}-u_{xx}+(u^{2})_{x}\ \ \textrm{for}\ \
(x,t)\in(0,\pi)\times(0,\infty),$

with odd-periodic boundary conditions
$u(0,t)=u(\pi,t)=u_{xx}(0,t)=u_{xx}(\pi,t)=0$ studied in [25, 26]. We assume
that the constant $\nu$ is positive. The system
$\\{\sin(kx)\\}_{k=1}^{\infty}$ constitutes the orthogonal in $L^{2}(0,\pi)$
basis of eigenfunctions of the leading linear operator $Lu=-\nu
u_{xxxx}-u_{xx}$ with the eigenvalues $\lambda_{k}=-\nu k^{4}+k^{2}$. To
verify (B1)-(B3) we take $H=L^{2}(0,\pi)$, $Y^{1}=\\{u\in
H^{3}(0,\pi)\,:u(0)=u(\pi)=u_{xx}(0)=u_{xx}(\pi)=0\\}$, and
$Y=H^{4}(0,\pi)\cap Y^{1}$. If $u=\sum_{k=1}^{\infty}u_{k}\sin(kx)$, then

$\|u\|_{Y^{1}}^{2}=\frac{2}{\pi}\sum_{k=1}^{\infty}k^{6}|u_{k}|^{2},\ \
\|u\|_{Y}^{2}=\frac{2}{\pi}\sum_{k=1}^{\infty}k^{8}|u_{k}|^{2}.$

The operator $L$ is diagonal and the evolution of the $k$-th mode via the
linear semigroup $e^{tL}$ is given by the formula

$u_{k}(t)=u_{k}(0)e^{(-\nu k^{4}+k^{2})t}.$

It is easy to verify that the function $k\to-\nu k^{4}+k^{2}$ has its maximum
equal to $\frac{1}{4\nu}$ at $k=\frac{1}{\sqrt{2\nu}}$. This leads to the
estimates $\|e^{tL}\|_{\mathcal{L}(Y^{1};Y^{1})}\leq e^{\frac{t}{4\nu}}$ and
$\|e^{tL}\|_{\mathcal{L}(Y;Y)}\leq e^{\frac{t}{4\nu}}$. Verification of (B3)
follows the concept of [21, Lemma 3.1]. Indeed, assuming that $u^{0}\in Y^{1}$
is the initial data, we obtain

$\|e^{tL}u^{0}\|^{2}_{Y}=\frac{2}{\pi}\sum_{k=1}^{\infty}k^{6}|u^{0}_{k}|^{2}k^{2}e^{2(k^{2}-\nu
k^{4})t}.$

A straightforward computation which involves the maximization over $k\geq 0$
shows that

$ke^{(k^{2}-\nu k^{4})t}\leq\frac{C}{\sqrt[4]{\nu t}}e^{\frac{t}{\nu}},$

where $C$ is independent of $k,\nu,t$. We deduce that

$\|e^{tL}u^{0}\|_{Y}\leq\frac{C}{\sqrt[4]{\nu
t}}e^{\frac{t}{\nu}}\|u^{0}\|_{Y^{1}},$

and (B3) is proved. To get (B2) it is enough that

$\|2uu_{x}-2vv_{x}\|_{H^{3}}\leq C(R)\|u-v\|_{H^{4}},$

where $\|u\|_{H^{4}},\|v\|_{H^{4}}\leq R$, which is straightforward to verify.

###  2.3. Overview of algorithm

First of all, since the phase space $Y$ of our abstract problem (2.1) is
infinite dimensional, we need a suitable representation of sets from this
space. Once we have such representation, the key concept of the method is, for
a given set $X(0)$ of initial data and a time-step $\tau$, to effectively
construct another set $X(\tau)$, such that it is guaranteed that every
solution starting from $X(0)$ at time $t=0$ belongs to $X(\tau)$ at time
$t=\tau$. Of course sets $X(0)$ and $X(\tau)$ must be described in previously
defined representation. In other words, if $\varphi:\Omega\to\mathbb{R}$ is
the local semigroup governed by the solutions of (2.1) we need to be able to
construct the set $X(\tau)$ such that $\varphi(\tau,X(0))\subset X(\tau)$.
Furthermore, as we do not exclude the possibility of blow-up arbitrarily, the
algorithm should also ensure that the value of $\varphi(\tau,x)$ is well-
defined for every $x\in X(0)$. At the same time, the set $X(\tau)$ should be
as small as possible, as we iterate the above procedure to find sets that
contain all solutions originating from a given set of initial data at large
times. The chosen representation always results in an overestimation at each
iteration. Because we need to represent sets $X(0)$ and $X(\tau)$ in the
computer memory, which is finite, we represent those sets as finite objects.
In the abstract problem (2.1) we assume that our phase space $Y$ is embedded
in a Hilbert space $H$ with the basis $\\{e_{k}\\}_{k=1}^{\infty}.$ We
represent $H=H_{P}\oplus H_{Q}.$ where
$H_{P}=\text{span}\\{e_{1},\dots,e_{n}\\}\cong\mathbb{R}^{n}$ and $H_{Q}$ is
an orthogonal complement of $H_{P}$ in $H.$ By $P,Q$ we will denote orthogonal
projections on the spaces $H_{P}$ and $H_{Q}$ respectively. We represent the
sets $X(0),X(\tau)\subset Y$ as

$X(0)=X_{P}(0)+X_{Q}(0)\qquad X(\tau)=X_{P}(\tau)+X_{Q}(\tau),$

where $X_{P}(0),X_{P}(\tau)\subset H_{P}$ are sets in finite dimensional space
and $X_{Q}(0),X_{Q}(\tau)\subset H_{Q}$ are infinite dimensional sets, which
need some finite representation. We realize this representation by giving some
inequalities, which are uniform with respect to coefficients in the Fourier
representation. Such decomposition of the set will be called a $P,Q$
representation. Now, essentially, we divide algorithm into two parts

  1. (1)

Find the set $X([0,\tau])$, given in the same $P,Q$ representation as the set
$X(0)$, such that every the solution to the considered problem satisfies
$\varphi(t,X(0))\subset X([0,\tau])$ for every $t\in[0,\tau]$. We will
equivalently write $\varphi([0,t],X(0))\subset X([0,\tau])$ and we will call
such set an enclosure. The details of finding this enclosure are given in
Section 2.6.

  2. (2)

Use the obtained enclosure to find the set $X(\tau).$ To this end, we use the
following procedures to find separately the sets $X_{P}(\tau)$ and
$X_{Q}(\tau).$

     * •

We formulate a differential inclusion in $\mathbb{R}^{n}$ for the $P$
component of the solution. In practice this $P$ component consists of a finite
number of Fourier coefficients (with respect to the space variable) of the
solutions of the system (2.6). The influence from the omitted variables, i.e.
the ones from $Q$, is estimated from the enclosure obtained in the first step.
The differential inclusion is integrated rigorously over time interval $\tau$
and for initial values belonging to the set $X_{P}(0)$. As a result we obtain
the bounds for the coordinates of $X_{P}(\tau)$. There are two possibilities
to get these bounds:

       * –

Study the evolution of variables belonging to $P$ separately, coordinate by
coordinate, by solving the linear differential inequalities. That is, for
every $i\in\mathbb{N}$ we are estimating the evolution of the Fourier modes
from the equation

$\frac{d}{dt}x_{i}(t)=\lambda_{i}x_{i}(t)+f_{i}(x(t)).$

       * –

Solve rigorously the following vector differential inclusion, obtained by
considering the Galerkin projection of the problem (2.1) on the space $H_{P},$
and estimating the influence the omitted terms in the equation through the
multivalued expression $I$.

$\frac{d}{dt}Px\in PF(Px)+I.$

The rigorous integration algorithm for finite dimensional vector inclusions
such as the one above is described in [14].

We intersect the estimates obtained by two above techniques in order to obtain
the sharper bounds.

     * •

Use the a priori estimates coming from the dissipativity of the linear part of
the problem to obtain the representation of $X_{Q}(\tau)$. The influence of
the nonlinear terms is estimated from the enclosure found in the first step.

The details of this step of the algorithm are given in Section 2.7.

###  2.4. Representation of sets

The important role in the algorithm will be played by the sequences
$\\{V_{i}\\}_{i=1}^{\infty}$ of intervals which we will call infinite interval
vectors. For a given infinite interval vector $V$ we will denote the $i$-th
interval by $V_{i}$ and its left and right ends by $V_{i}^{-},V_{i}^{+}$,
respectively. If for every sequence $\\{v_{i}\\}_{i=1}^{\infty}$, such that
$v_{i}\in V_{i}$, the series $\sum_{i=1}^{\infty}e_{i}v_{i}$ converges in $H$,
we will call the set $\\{\sum_{i=1}^{\infty}e_{i}v_{i}\in H:v_{i}\in V_{i}\\}$
a representation of infinite interval vector and we will say that the infinite
vector is representable in the space $H.$ It is possible that infinite
interval vector $V$ does not represent subset of $H$ as it can happen that the
series $\sum_{i=1}^{\infty}e_{i}v_{i}$ where $v_{i}\in V_{i}$ does not
converge in $H$. Whenever it will not lead to confusion we will use the same
nonantion for infinite interval vectors and their representations.

We define several useful operation on the infinite interval vectors. First,
for infinite vector $V$ we define the quantities

$V^{-}=\\{[V^{-}_{i},V^{-}_{i}]\\}_{i=1}^{\infty},\quad
V^{+}=\\{[V^{+}_{i},V^{+}_{i}]\\}_{i=1}^{\infty}.$

For a given interval $I$ the we define multiplication an infinite interval
vector $V$ by the interval $I$ as

$I\,V=V\,I=\\{IV_{i}\\}_{i=1}^{\infty}.$

For two infinite intervals vectors $V$ and $W$ we define their sum and
element-wise product as

$V+W=\\{V_{i}+W_{i}\\}_{i=1}^{\infty},\quad
V*W=\\{V_{i}W_{i}\\}_{i=1}^{\infty}.$

We say that vector $V$ is a subset of $W$ and denote by $V\subset W$ if and
only if $V_{i}\subset W_{i}$ for every $i\in\mathbb{N}.$ Additionally we
define $V\subset_{\text{int}}W$ if and only if $V_{i}\subset\text{int}W_{i}$
for every $i\in\mathbb{N}.$ We define the convex hull of two infinite
intervals vectors as

$\text{conv}\\{V,W\\}=\\{\text{conv}\\{V_{i}\cup W_{i}\\}\\}_{i=1}^{\infty}.$

The intersection of two infinite vectors is defined in the following way

$V\cap W=\\{V_{i}\cap W_{i}\\}_{i=1}^{\infty}.$

Note that all above operations make sense for all infinite interval vectors,
and not only the representable ones.

In the algorithm we consider such sets $X=X_{P}+X_{Q}$, for which there exist
infinite interval vectors whose representation contains the set $X.$
Specifically, for the Brusselator system we will work with the pairs of
interval infinite vectors which are given in the form
$(U,V)=\\{(U_{i},V_{i})\\}_{i\in\mathbb{N}^{+}}$ where $U_{i}$ and $V_{i}$ are
the intervals. Such vectors can be easily re-indexed into the form described
previously. We will work with class of infinite interval vectors, with
polynomial estimates on the tail. That means that for some $n\in\mathbb{N}$
and $s\in\mathbb{R}$ and every sequence $u_{i}\in U_{i}$ and $v_{i}\in V_{i}$
we have

$u_{i}\in\frac{[C_{U}^{-},C_{U}^{+}]}{k^{s}},\quad
v_{i}\in\frac{[C_{V}^{-},C_{V}^{+}]}{k^{s}},\quad\text{for $i\geq n$},$ (2.8)

where $C_{U}^{-}\leq C_{U}^{+}$ and $C_{V}^{-}\leq C_{V}^{+}$ are constants.
In this manner, the tail of the Fourier expansion (for $i\geq n$) can be
represented by specifying the decay rate $s$ of the coefficients and four
additional constants $C_{U}^{-},C_{U}^{+},C_{V}^{-},C_{V}^{+}$. Lemmas 5.5 and
5.6 are helpful in the implementation of operations of element-wise
multiplication and addition for such class of infinite interval vectors.

Since for the Bruselator, we consider the system of two PDEs, the states are
the sets $X=X_{P}+X_{Q}$ which are the subsets of $Y$. Their elements are
pairs $(u,v).$ To represent these pairs we use the Fourier basis which
consists of the eigenfunctions of the operator $L$, i.e. $(\sin(kx),0)$ and
$(0,\sin(kx))$ for $k\in\mathbb{N}^{+}$. The finite dimensional space $H_{P},$
in which the sets $X_{P}$ are always contained, is equal to
$\text{span}\\{(\sin(x),0),(0,\sin(x)),\ldots,(\sin(kx),0),(0,\sin(kx))\\}.$

The set $X$ has to be a subset of some representable infinite interval vector
vector $(U,V)$. In the algorithm, this means that $X_{P}$ is a subset of some
cube in $H_{P}.$ The cubes are the simplest examples of possible
representations of the set $X_{P}.$ More sophisticated parallelepiped-type
objects can also used. They can reduce overestimation of the integration
results for the rigorous ODE solvers and overcome so calling wrapping effect
(for example see [15]). The element $(u,v)\in H_{Q}$ belongs to $X_{Q}$ if it
satisfies (2.8) with $s>1.$ This guarantees that $X_{Q}$ is a subset of $Y$.
In that case we can estimate the result of series multiplication by using
Lemmas 5.1 and 5.2. Finally, we note that we may restrict our space to the
space of functions $u,v$ with only odd nonzero coefficients. Then the set $X$
is a subset of $W:=\\{(u,v)\in
Y:u_{2i}=v_{2i}=0,\;\text{for}\;i\in\mathbb{N}\\}.$ The representation of set
$X$ is roughly the same, except that we have to enforce that $u_{i}=v_{i}=0$
if $i$ is even.

###  2.5. Computation of nonlinear terms.

In the course of the algorithm for a given set $X=X_{P}+X_{Q}$, we need to
compute the set $X^{1}$ such that $f(X)\subset X^{1}.$ This set, represented
as $X^{1}=X_{P}^{1}+X_{Q}^{1}$ constitutes the estimates for $f(X)$ and hence
it should be as small as possible. The set $X^{1}$ is used in further steps of
algorithm. For the Brusselator problem we have that

$f(u,v)=(u^{2}v+A\sin(x),Bu-u^{2}v).$

For functions $(u,v)\in Y$ the components of $f$ can be represented in the
following sine Fourier series with the coefficients $a_{i},b_{i}$ dependent on
$u,v$

$u^{2}v+A\sin(x)=\sum_{i=1}^{\infty}a_{i}\sin(ix),\quad
Bu-u^{2}v=\sum_{i=1}^{\infty}b_{i}\sin(ix).$

The set $X^{1}_{P}$ is represented as the cube (or parallelogram) in $H_{P}$
and $X^{1}_{Q}$ is described by the polynomial decay of Fourier coefficients.
The first step of finding $X^{1}$ is estimating the square of $u$ which is
represented in the cosine Fourier series. For this purpose we use Lemma 5.1.
Having computed the coefficients of $u^{2}$ we need to find the coefficients
of $(u^{2})v.$ To this end we use Lemma 5.2. Finally, we use Lemma 5.5 to
compute the representation of sums of particular terms which appear in the
definition of $f$. We also need to compute the image $L(X)$ but as $L$ is a
diagonal operator we only need to multiply every given coefficient by the
corresponding eigenvalue of $L.$ The result of such multiplication is given in
Lemma 5.6. Additionally in the algorithm we need the decomposition

$f(p+q)=f(p)+f_{2}(p,q),$

where $p\in H_{P}$, $q\in H_{Q}$ and $f_{2}(p,q)=f(p+q)-f(p)$. This
decomposition is required for the formulation of the differential inclusion.
For the Brusselator system we can write
$f(u_{P}+u_{Q},v_{P}+v_{Q})=f(u_{P},v_{P})+f_{2}(u_{P},v_{P},u_{Q},v_{Q})$
where

$f_{2}(u_{P},v_{P},u_{Q},v_{Q})=((2u_{P}u_{Q}+u_{Q}^{2})(v_{P}+v_{Q})+u_{P}^{2}v_{Q},Bu_{Q}-(2u_{P}u_{Q}+u_{Q}^{2})(v_{P}+v_{Q})-u_{P}^{2}v_{Q}).$

###  2.6. Computation of the enclosure.

We start this section with the definition of a enclosure.

######  Definition 2.4.

The set $X([0,\tau])$ is a enclosure of the set $X^{0}\subset Y$ for time
$\tau>0$ if $\varphi(t,X^{0})\subset X([0,\tau])$ for every $t\in[0,\tau]$.

The following Lemma can be used in order to validate if for the given set of
initial data $X^{0}$ the set $X^{0}+Z$ is an enclosure.

######  Lemma 2.6.

Assume that (A1) and (A2) hold and let $\\{X^{0}_{i}\\}_{i=1}^{\infty}$ be a
countable family of intervals $X^{0}_{i}=[x_{i}^{-},x_{i}^{+}]$ such that the
set $X^{0}:=\\{{\sum_{i=1}^{\infty}e_{i}x_{i}:x_{i}\in X_{i}^{0}}\\}$ is
bounded in $Y$. Moreover, define another set
$Z:=\\{{\sum_{i=1}^{\infty}e_{i}z_{i}:z_{i}\in Z_{i}}\\},$ also bounded in
$Y$, where $Z_{i}=[z^{-}_{i},z^{+}_{i}]$ are intervals containing zero. Let
$x^{0}\in X^{0}$ and $n\in\mathbb{N}$ . We assume that for every
$i\in\mathbb{N}$ there holds

$[g_{i}^{-},g_{i}^{+}]\cap[h_{i}^{-},h_{i}^{+}]\subset\text{int }Z_{i},$ (2.9)

where

$g^{-}_{i}=\min_{t\in[0,\tau]}\left[(e^{\lambda_{i}t}-1)\left(\frac{f^{-}_{i}}{\lambda_{i}}+x^{-}_{i}\right)\right],\quad
g^{+}_{i}=\max_{t\in[0,\tau]}\left[(e^{\lambda_{i}t}-1)\left(\frac{f^{+}_{i}}{\lambda_{i}}+x^{+}_{i}\right)\right],$
(2.10)
$h^{-}_{i}=\min_{t\in[0,\tau]}t(\lambda_{i}(x_{i}^{-}+z_{i}^{-})+f_{i}^{-}),\quad
h^{+}_{i}=\max_{t\in[0,\tau]}t(\lambda_{i}(x_{i}^{+}+z_{i}^{+})+f_{i}^{+}),$
(2.11)

and $f^{-}_{i},f^{+}_{i}\in\mathbb{R}$ satisfy

$f^{-}_{i}\leq f_{i}(X^{0}+Z)\leq f^{+}_{i}.$ (2.12)

Then for every $x^{0}\in X^{0}$ there exists a continuous function
$x:[0,\tau]\to Y$ which is a unique solution to $\eqref{eq:AbstractProblem}.$
Moreover for every $t\in[0,\tau]$ we have

  1. (1)

$x_{i}(t)\in x^{0}_{i}+[0,t](f_{i}(X^{0}+Z)+\lambda_{i}(X^{0}+Z))$ for every
$i\in\mathbb{N},$

  2. (2)

$x^{-}_{i}+g_{i}^{-}\leq x_{i}(t)\leq x^{+}_{i}+g_{i}^{+}$ for every
$i\in\mathbb{N},$

  3. (3)

$x_{i}(t)\in
e^{t\lambda_{i}}x^{0}_{i}+\frac{e^{\lambda_{i}t}-1}{\lambda_{i}}[f^{-}_{i},f^{+}_{i}]$
for every $i\in\mathbb{N}.$

###### Proof.

We define the operator $T:C([0,\tau];Y)\to C([0,\tau];Y)$ in the following way

$T(g)(t)=e^{Lt}x^{0}+\int_{0}^{t}e^{L(t-s)}f(g(s))ds.$ (2.13)

We consider the set

$S_{\tau}=\\{g\in C([0,\tau];Y):g(0)=x^{0}\;\text{and for every $t\in[0,\tau]$
we have }g(t)\in X^{0}+Z\\}.$

We prove first that for every $x^{0}\in X^{0}$ the mapping $T$ leads from
$S_{\tau}$ to itself. We observe that $y_{i}(t)=T_{i}(g)(t)$ is a solution to
the non-autonomous ODE

$\frac{d}{dt}y_{i}(t)=\lambda_{i}y_{i}(t)+f_{i}(g(t)),$ (2.14)

for every $i\in\mathbb{N.}$ We will prove that
$y_{i}(t)\in[x^{-}_{i}+z^{-}_{i},x^{+}_{i}+z^{+}_{i}],$ for every
$t\in[0,\tau].$ For the sake of contradiction assume that $\tau_{1}<\tau,$
where
$\tau_{1}:=\sup_{t\in[0,\tau]}\\{y_{i}(s)\in[x^{-}_{i}+z^{-}_{i},x^{+}_{i}+z^{+}_{i}]\
\text{for every }s\leq t\\}.$ We observe that

$y_{i}(t)=x_{0}+\int_{0}^{t}\lambda y_{i}(s)+f_{i}(g(s))ds,\quad
y_{i}(t)=e^{\lambda_{i}t}x_{0}+\int_{0}^{t}e^{\lambda(t-s)}f_{i}(g(s))ds.$

So for $t\leq\tau_{1}$ we have

$y_{i}(t)\leq
x^{0}_{i}+\int_{0}^{t}(\lambda(x^{+}_{i}+z^{+}_{i})+f_{i}^{+})ds\leq
x^{+}_{i}+h_{i}^{+}.$

Similarly for $t\leq\tau_{1}$ we have $x^{-}_{i}+h_{i}^{-}\leq y_{i}(t).$ We
observe that for $t\leq\tau_{1}$

$y_{i}(t)\leq
e^{\lambda_{i}t}x^{0}_{i}+\int_{0}^{t}e^{\lambda_{i}(t-s)}f_{i}^{+}ds\leq
e^{\lambda_{i}t}x^{+}_{i}+\frac{e^{\lambda_{i}t}-1}{\lambda_{i}}f_{i}^{+}=x^{+}_{i}+g_{i}^{+}.$

We also see that $x^{-}_{i}+g_{i}^{-}\leq y_{i}(t).$ This implies that for
every $t\in[0,\tau_{1}]$ we have
$x^{-}_{i}+z^{-}_{i}<y_{i}(t)<x^{+}_{i}+z^{+}_{i}.$ By the continuity of
$y_{i}(t)$ we can find $\delta>0$ such that
$y_{i}(t)\in[x^{-}_{i}+z^{-}_{i},x^{+}_{i}+z^{+}_{i}]$ for every
$t\leq\tau_{1}+\delta.$ This is a contradiction, so $\tau=\tau_{1}$. We deduce
that $T(S_{\tau})\subset S_{\tau}.$ To show that $T$ is a contraction we equip
the space $C([0,\tau],Y)$ with the norm
$\|y\|_{\alpha}=\sup_{t\in[0,\tau]}\left\|{y(t)}\right\|_{Y}e^{-t\alpha},$
where $\alpha>0$ is an appropriately chosen positive constant. For
$g_{1},g_{2}\in S_{\tau}$ we have

$\displaystyle\left\|\int_{0}^{t}e^{L(t-s)}(f(g_{2}(s))-f(g_{1}(s)))ds\right\|_{Y}$
$\displaystyle\leq
e^{tC_{1}}\int_{0}^{t}e^{s\alpha}e^{-s\alpha}\left\|{f(g_{1}(s))-f(g_{2}(s)))}\right\|_{Y}ds$
$\displaystyle\leq\frac{e^{tC_{1}}e^{t\alpha}}{\alpha}\sup_{s\in[0,t]}\left(e^{-s\alpha}\left\|{f(g_{1}(s))-f(g_{2}(s))}\right\|_{Y}\right)$
$\displaystyle\leq\frac{C(R)e^{tC_{1}}e^{t\alpha}}{\alpha}\sup_{s\in[0,t]}e^{-s\alpha}\left\|{g_{1}(s)-g_{2}(s)}\right\|_{Y}.$

If we pick $\alpha=\sup_{s\in[0,\tau]}\frac{1}{2C(R)e^{sC_{1}}}$ then

$\left\|T(g_{1})-T(g_{2})\right\|_{\alpha}\leq\sup_{s\in[0,\tau]}\frac{C(R)e^{sC_{1}}}{\alpha}\left\|{g_{1}-g_{2}}\right\|_{\alpha}\leq\frac{1}{2}\left\|{g_{1}-g_{2}}\right\|_{\alpha}.$

where $R>0$ is a bound in the norm $\left\|{.}\right\|_{Y}$ of elements of the
set $X_{0}+Z.$ This shows that $T$ is a contraction. From the Banach fixed
point theorem, we deduce that $T$ must have a unique fixed point in
$S_{\tau}$. Now, formulae (1)-(3) are straightforward, which completes the
proof. ∎

######  Remark 2.5.

The assertion of Lemma 2.6 holds if we replace conditions $(A1)-(A2)$ by their
more general counterparts $(B1)-(B3)$.

###### Proof.

We consider the same operator $T$ and the set $S_{\tau}$ as in proof of Lemma
2.6. The argument that $T(S_{\tau})\subset S_{\tau}$ is the same, only the
proof of contractivity changes. We will show that in this more general case we
can find $\alpha>0$ such that $T$ is a contraction with respect to the norm
$\left\|{y}\right\|_{\alpha}=\sup_{t\in[0,\tau]}e^{-t\alpha}\left\|{y(t)}\right\|_{Y}.$
For $g_{1},g_{2}\in S_{\tau}$ we compute

$\displaystyle\left\|\int_{0}^{t}e^{L(t-s)}(f(g_{2}(s))-f(g_{1}(s)))ds\right\|_{Y}$
$\displaystyle\leq
C_{3}e^{tC_{4}}\int_{0}^{t}\frac{1}{(t-s)^{\gamma}}e^{s\alpha}e^{-s\alpha}\left\|{f(g_{1}(s))-f(g_{2}(s)))}\right\|_{Y^{1}}ds$
$\displaystyle\leq C_{3}e^{\tau
C_{4}}C(R)\left\|{g_{1}-g_{2}}\right\|_{\alpha}\int_{0}^{t}\frac{1}{(t-s)^{\gamma}}e^{\alpha
s}ds$ $\displaystyle=C_{3}e^{\tau
C_{4}}C(R)e^{t\alpha}\left\|{g_{1}-g_{2}}\right\|_{\alpha}\int_{0}^{t}\frac{1}{s^{\gamma}}e^{-s\alpha}ds.$

where $R>0$ is a bound of the norm $\left\|{\cdot}\right\|_{Y}$ satisfied by
all elements of the set $X_{0}+Z.$ For every $\delta\leq t$ we have

$\int_{0}^{t}\frac{1}{s^{\gamma}}e^{-s\alpha}ds\leq\int_{0}^{\delta}\frac{1}{s^{\gamma}}ds+\frac{1}{\delta^{\gamma}}\int_{\delta}^{t}e^{-s\alpha}ds\leq\frac{\delta^{1-\gamma}}{1-\gamma}+\frac{1}{\delta^{\gamma}}\frac{1}{\alpha}$

So, if we pick $\delta,\alpha$ such that

$\delta^{1-\gamma}\leq\frac{1}{4(1-\gamma)}C_{3}e^{\tau
C_{4}}C(R),\quad\alpha\geq\frac{1}{4\delta^{\gamma}}C_{3}e^{\tau C_{4}}C(R),$

then for every $t\leq\tau$ we have

$\left\|\int_{0}^{t}e^{L(t-s)}(f(g_{2}(s))-f(g_{1}(s)))ds\right\|_{Y}\leq\frac{1}{2}e^{t\alpha}\left\|{g_{1}-g_{2}}\right\|_{\alpha}.$

Hence, $T$ is a contraction on the set $S_{\delta}$ equipped with the norm
$\left\|{\cdot}\right\|_{\alpha}$. The Banach fixed point theorem gives the
unique fixed point of the map $T$, and the proof is complete. ∎

In Lemma 2.6 the sets $X^{0},Z$ are representable infinite interval vectors.
If $X^{0}$ is not a representable infinite interval vector, we can still find
a enclosure, provided we can choose the representable interval vector $X^{1}$
such that $X^{0}\subset X^{1}.$ The following algorithm takes as an input the
infinite interval vector $X^{0},$ the infinite interval vector $Z$ such $0\in
Z_{i}$ for every $i\in\mathbb{N},$ and the step $\tau>0.$ The algorithm
validates if $X^{0}+Z$ is an enclosure for $\tau$ by checking assumption (2.9)
of Lemma 2.6.

  1. (1)

Find the infinite vectors $V^{f}$ and $V^{L}$ such that $f(X+Z)\subset V^{f}$
and $L(X+Z)\subset V^{L}.$

  2. (2)

Construct the infinite vector $V^{\text{ND}}$ such that
$[0,\tau](V^{L}+V^{f})\subset V^{\text{ND}}$

  3. (3)

Compute the infinite vectors $V^{E_{1}}$ and $V^{L^{-1}}$ which satisfy

$e^{[0,\tau]\lambda_{i}}-1\subset(V^{E_{1}})_{i},\qquad\frac{1}{\lambda_{i}}\in(V^{L^{-1}})_{i}\
\ \textrm{for every}\ \ i\in\\{1,\dots\\}.$

  4. (4)

Compute the infinite vectors $V^{\text{D1}}$ and $V^{\text{D2}}$ which satisfy

$V^{E_{1}}*((V^{f})^{+}*V^{L^{-1}}+X^{+})\subset V^{\text{D1}},\quad
V^{E_{1}}*((V^{f})^{-}*V^{L^{-1}}+X^{-})\subset V^{\text{D2}}.$

  5. (5)

Find infinite interval vectors $V^{\text{D}}$ such that

$\text{conv}\\{V^{\text{D1}},V^{\text{D2}}\\}\subset V^{\text{D}}.$

  6. (6)

Compute $Z^{1}$ such that $V^{\text{ND}}\cap V^{\text{D}}\subset Z^{1}.$

  7. (7)

Check if $Z^{1}\subset_{\text{int}}Z$ holds.

For every $i\in\mathbb{N}$ the interval $[h_{i}^{-},h_{i}^{+}]$ from (2.11) is
contained in the interval $V^{ND}_{i},$ where $V^{ND}$ is an infinite interval
vector obtained in step (2). Similarly for every $i\in\mathbb{N}$ the interval
$[g_{i}^{-},g_{i}^{+}]$ from (2.10) is contained in the interval $V^{D}_{i},$
where $V^{D}$ is an infinite interval vector from step (5). If condition in
step (7) is satisfied, we have that $V^{ND}_{i}\cap
V^{ND}_{i}\subset\text{int}Z_{i}$ for every $i\in\mathbb{N},$ which implies
that assumptions of Lemma 2.6 hold and the set $X^{0}+Z$ is an enclosure for
our initial data $X^{0}$ and the time step $\tau>0$ .

If condition $(7)$ does not hold, we can modify the set $Z$ and repeat the
validation procedure. The reasonable guess is to take $Z=[0,c]Z^{1},$ where
$c>1.$ Another possibility is to decrease the time step. The following lemma
shows that for certain class of sets it is always possible to find the
enclosure using the above algorithm with sufficiently small time step $\tau$.

######  Lemma 2.7.

Let $\\{a_{i}\\}_{i=1}^{\infty}$ be a given sequence of positive numbers.
Assume that

  * (I)

Conditions (A1)-(A2) or (B1)-(B3) hold.

  * (II)

For some $i_{0}\in\mathbb{N}$ we have $\lambda_{i}<0$ for every $i\geq i_{0}.$

  * (III)

For every bounded set $A\subset Y$ such that

$\sup_{x\in A}|x_{i}|=O(a_{i}),$

the following holds

$\sup_{x\in A}\left|\frac{f_{i}(x)}{\lambda_{i}}\right|=o(a_{i}).$

Let $\\{X^{0}_{i}\\}_{i=1}^{\infty}$ and $\\{Z_{i}\\}_{i=1}^{\infty}$ be
sequences of intervals such that

  * (IV)

The set $X^{0}:=\\{{\sum_{i=1}^{\infty}e_{i}x_{i}:x_{i}\in X_{i}^{0}}\\}$ is
bounded in $Y$. Intervals $X_{i}^{0}$ contain zero for every $i\geq i_{0},$
where $i_{0}\in\mathbb{N}$ and

$\quad\sup_{x_{i}\in X_{i}^{0}}|x_{i}|=O(a_{i}).$

  * (V)

The set $Z:=\\{{\sum_{i=1}^{\infty}e_{i}z_{i}:z_{i}\in Z_{i}}\\}$ is bounded
in $Y$. Every interval $Z_{i}$ contains zero and

$z_{i}^{+}=\Theta(a_{i})\quad\text{and}\quad z_{i}^{-}=\Theta(a_{i}).$

Under these assumptions there exists $\tau>0$ such that condition (2.9) of
Lemma 2.6 is satisfied, and, in consequence, for every $x^{0}\in X^{0}$ there
exists a continuous function $x:[0,\tau]\to Y$ which is a unique solution to
(2.1) and the following estimates hold for every $t\in[0,\tau]$

  1. (1)

$x_{i}(t)\in x^{0}_{i}+[0,t](f_{i}(X^{0}+Z)+\lambda_{i}(X^{0}+Z))$ for every
$i\in\mathbb{N},$

  2. (2)

$x^{-}_{i}+g_{i}^{-}\leq x_{i}(t)\leq x^{+}_{i}+g_{i}^{+}$ for every
$i\in\mathbb{N},$

  3. (3)

$x_{i}(t)\in
e^{t\lambda_{i}}x^{0}_{i}+\frac{e^{\lambda_{i}t}-1}{\lambda_{i}}[f^{-}_{i},f^{+}_{i}]$
for every $i\in\mathbb{N},$

where $g_{i}^{-},g_{i}^{+}$ are given by (2.10) and $f_{i}^{-},f_{i}^{+}$ are
given by (2.12).

###### Proof.

Observe that we can find $i_{1}$ such that

$[g_{i}^{-},g_{i}^{+}]\subset\text{int}\;Z_{i}=(z_{i}^{-},z_{i}^{+}),$ (2.15)

for every $i\geq i_{1}$ and $\tau>0.$ Indeed for big enough $i$ we have

$g_{i}^{+}=\max_{t\in[0,\tau]}\left[(e^{\lambda_{i}t}-1)\left(\frac{f^{+}_{i}}{\lambda_{i}}+x^{+}_{i}\right)\right]\leq\max_{t\in[0,\tau]}\left[(e^{\lambda_{i}t}-1)\frac{f^{+}_{i}}{\lambda_{i}}\right]\leq\left|\frac{f_{i}^{+}}{\lambda_{i}}\right|<z_{i}^{+}.$

The first inequality follows from the fact that $x_{i}^{+}$ is positive and
term $(e^{\lambda_{i}t}-1)$ is negative for sufficiently large $i$. The second
inequality follows from the assumption on $\lambda_{i}.$ The third one is a
consequence of the fact that $\sup_{x\in
X}\left|\frac{f_{i}(x)}{\lambda_{i}}\right|=o(a_{i})$ and
$z^{+}_{i}=\Theta(a_{i}).$ We have shown that for $i\geq i_{1}$ we have the
inclusion (2.15). For the remaining, ”low”, indexes $i$, we take $\tau>0$ such
that

$\tau<\sup_{i<i_{1}}\frac{|z_{i}^{-}|}{\left|(\lambda_{i}(x_{i}^{-}+z_{i}^{-})+f_{i}^{-})\right|+1}\quad\text{and}\quad\tau<\sup_{i<i_{1}}\frac{|z_{i}^{+}|}{\left|(\lambda_{i}(x_{i}^{+}+z_{i}^{+})+f_{i}^{+})\right|+1}.$

For such choice of $\tau$ we have

$[h_{i}^{-},h_{i}^{+}]\subset\text{int}\;Z_{i}\quad\text{for}\quad i<i_{1},$

which ends the proof. ∎

######  Remark 2.6.

The key assumption of the above lemma is (III). The sequence $a_{i}$ signifies
some decay of Fourier coefficients of the elements $x$ of a set $A$. The decay
of Fourier coefficients of $\frac{f_{i}(x)}{\lambda_{i}}$ for $x\in A$ must be
essentially faster than $a_{i}$. Assume for example that
$a_{i}=\frac{1}{i^{s}}$ and that $f$ is a polynomial. Then, the decay of the
coefficients of $f_{i}(x)$ is the same as that of $x$, that is also
$\frac{1}{i^{s}}$ (see Section 5). If the leading operator is dissipative,
i.e. $\lambda_{i}\to-\infty$, then, after the division by $\lambda_{i}$ this
decay will be essentially faster then $a_{i}$. Thus, dissipativity of $L$
guarantees that the enclosure can be always found, and a step of the rigorous
integration algorithm can be performed. The drawback is, that the length of a
time-step $\tau$ can be very small in Lemma 2.7. We expect that such situation
holds for dissipative problems with finite type blow-up, for instance for the
problems governed by the Fujita equation $u_{t}=u_{xx}+|u|^{p-1}u$, cf. [9].

Using the above lemma we deduce that for the Brusselator system we can always
find an enclousure for any set described in the Section 2.4, by choosing
sufficiently short time-step.

######  Remark 2.7.

For every $C>0,\;\varepsilon>0,\;k>1$ there exists $\tau>0$ such that for
every initial data $(u^{0},v^{0})\in Y$ satisfying

$|u^{0}_{k}|\leq\frac{C}{|k|^{s}}\quad\text{and}\quad|v^{0}_{k}|\leq\frac{C}{|k|^{s}},$

the solution $(u(t),v(t))$ to (1.1) satisfies

$|u_{k}(t)|\leq\frac{C+\varepsilon}{|k|^{s}}\quad\text{and}\quad|v_{k}(t)|\leq\frac{C+\varepsilon}{|k|^{s}}\quad\text{for}\;t\in[0,\tau].$

One can easily construct an example of a problem without this property. To
this end let us consider the logistic model of population growth with
diffusion and homogeneous Dirichlet boundary conditions. The corresponding
equation together with the initial and boundary conditions are the following

$\begin{cases}u_{t}=du_{xx}+u-u^{2}\
\text{for}\;(x,t)\in[0,\pi]\times(0,\infty),\\\
u(t,x)=0\;\text{for}\;(x,t)\in\\{0,\pi\\}\times(0,\infty),\\\
u(0,x)=u^{0}(x),\end{cases}$ (2.16)

for $d>0.$ As we impose the Dirichlet boundary condition, from the equation we
observe that $u_{xx}$ is always $0$ at the boundary. But, in general, the
fourth derivative $u_{xxxx}$ can have nonzero values at the boundary. This
implies that the coefficients in the expansion of the solution in the sine
trigonometric series cannot have the arbitrarily fast polynomial decay.

We can write the equation in above example in the form (2.1) with
$L=du_{xx}+u$ and $f(u)=u^{2}.$ For the initial condition
$u_{0}=\sqrt{\frac{\pi}{8}}\sin(x)$ we have that

$f(u_{0})=\frac{\pi}{8}(1-\cos(2x))=\sum_{i=1}^{\infty}\frac{(-1+(-1)^{i})}{-4n+n^{3}}\sin(ix).$

So using Lemma 2.6 we have a chance to find the time step $\tau>0$ and $C>0$
such that the solution satisfies

$|u_{k}(t)|\leq\frac{C}{k^{s}}\quad\text{for}\quad t\in[0,\tau],$

only for $s\in(1,5).$ Similar difficulty would occur if in the considered
Brusselator system the term $\sin(x)$ would be replaced with the original term
from the Brusselator ODE, that is the constant $1$.

###  2.7. Evolution of sets

We assume that $X([0,\tau])=X_{P}([0,\tau])+X_{Q}([0,\tau])$ is an enclosure
for the initial data $X=X_{P}+X_{Q}$ and the time step $\tau>0.$ The following
procedure is used to find the set $X(\tau)$ such that $S(\tau)X\subset
X(\tau).$

  1. (1)

Compute the infinite interval vector $V$ such that
$f_{2}(X_{P}([0,\tau]),X_{Q}([0,\tau]))\subset V$.

  2. (2)

Solve the system of differential inclusions

$\frac{d}{dt}Px\in LPx+Pf(Px)+PV,$ (2.17)

with initial condition in the set $X_{P}$. Note that if the initial data
belongs to the set $X_{P}$ which has more complicated structure then a vector
of intervals (for example it can be a parallelepiped) this step, can be
realized without enclosing the initial data in an interval vector. This will
result in sharper estimates. As the result, the solver will generate the set
$X_{P1}$ such that $Px(\tau)\in X_{P1}$ for every $x(0)\in X_{P}+X_{Q}$.

  3. (3)

Compute the infinite interval vector $V^{2}$ such that $f(X([0,\tau]))\subset
V^{2}.$

  4. (4)

Compute the infinite vectors $V^{E_{1}},V^{E_{2}}$ and $V^{L^{-1}}$ which
satisfy

$e^{\tau\lambda_{i}}-1\in(V^{E_{1}})_{i},\qquad
e^{\tau\lambda_{i}}\in(V^{E_{2}})_{i},\qquad\frac{1}{\lambda_{i}}\in(V^{L^{-1}})_{i}\
\ \textrm{for every}\ \ i\in\\{1,\dots\\}.$

  5. (5)

Compute $V^{3}$ such that $V^{E_{2}}*X+V^{E_{1}}*V^{L^{-1}}*V^{2}\subset
V_{3}$.

  6. (6)

Return the set $(PV_{3}\cap X_{P1})+QV_{3}$.

Note that the interval vectors $V^{E_{1}}$ both in the algorithm for finding
the enclosure and evolving the set are not representable. However, we still
need the data structures to represent them as infinite vectors and perform the
operations such as addition or elementwise multiplication. Particular way of
representing such sets and their operations with polynomial decay or growth of
coefficients is addressed in Section 5.

The detailed description how to rigorously solve the differential inclusion
can be found in [14]. For the step (2) we can consider differential inclusion
only on part of variables represented explicitly. This can be beneficial for
the computational time.

As we are using the infinite interval vectors (2.8) we need to determine the
decay rate $s$ for $V^{E_{2}}.$ It is possible to impose arbitrarily fast
polynomial decay in this term. On the other hand, the vector $V^{E_{1}}$ can
be represented by (2.8) with $s=0$. Finally, the representation of
$V^{L^{-1}}$ decays with some given $s_{1}>0$ determined by the decay of the
inverses of the eigenvalues of $L$ (which have to decay to zero as the
considered problem is disspative). So, $V_{3}$ can have higher $s$ then the
initial data $X$. Maximal increase of the rate $s$ in $V_{3}$ is equal to
$s_{1}$. Theoretically in every time step it is possible to increase the decay
$s$ by the value of $s_{1}.$ But it can lead to overestimates on some
variables so it is sometimes beneficial to keep the old $s.$ In the code we
are using heuristic algorithm which estimates upper bound of resulting series
to decide if it worth to increase the exponent.

##  3\. Computer assisted proof of periodic orbit existence

In this section we describe the computer assisted proof of Theorem 1.1.

###  3.1. Overview of the proof

The proof of Theorem 1.1 is based on the Schauder fixed point theorem. We
check that for some previously prepared Poincaré map $\mathcal{P}$ the image
of appropriately chosen compact initial set $X^{0}$ is contained in itself,
i.e. $\mathcal{P}(X^{0})\subset X^{0}.$ To validate the inclusion, we apply
the previously described integration algorithm and Lemma 3.1 to address the
issue of crossing the section. As the preliminary step, we construct a set of
initial data $X^{0}$ and a section $l$ from the analysis of results of
approximate numerical integration of the Galerkin projection both of the
Brusselator equation and its variational equation.

###  3.2. Poincaré map - crossing the section

The problem of finding the periodic orbit for the system (2.1) is reduced to
finding the fixed point of the Poincaré map. In this section we present the
algorithm of rigorous computatation of the Poincaré map and justify its
correctness. While the same algorithm is for ODEs can be found in [24, Section
5], and its infinite dimensional version appears in [25, Section 3], we
present it for the exposition completeness. Note that although the results of
the present section follow closely the concepts of [24, 25], the results in
[25, Lemma 6 and Theorem 8] use the Brouwer fixed point theorem, and the
argument on passing to the limit in Galerkin projections to get the fixed
point, and our Theorem 3.2 uses the Schauder fixed point theorem.

We assume as in Section 2 that we consider abstract problem (2.1). The spaces
$H,Y$ are also the same as in Section 2. By $\varphi$ we define the local
semiflow which is given by solutions of (2.1). The following lemma allows us
to check if evolution of initial set transversally intersects the section,
which is a kernel of some affine map.

######  Lemma 3.1.

Let $l:Y\to\mathbb{R}$ be a function given by the formula
$l(x)=\sum_{i=1}^{n}\alpha_{i}(x_{i}-\beta_{i})$, where
$\alpha_{1}\ldots,\alpha_{n},\beta_{1},\ldots,\beta_{n}\in\mathbb{R}.$ Let $X$
be a bounded set in $Y.$ We assume that for some $\tau>0$ the following
conditions hold:

  1. (1)

For every $x\in X$ we have $l(x)<0\;\text{and}\;l(\varphi(\tau,x))>0.$

  2. (2)

For every $x\in\varphi([0,\tau],X)$ we have
$\sum_{i=1}^{n}\alpha_{i}F_{i}(x)>0.$

Then for every $x\in X$ there exists a unique $\tau_{l}(x)\in(0,\tau)$ such
that $l(\varphi(\tau_{l}(x),x))=0.$ Moreover the function
$\tau_{l}:X\to(0,\tau)$ is continuous in the norm of the space $Y.$

###### Proof.

Let $x\in X.$ We observe that condition (1) and the continuity of the flow
imply that there exists $\tau_{l}(x)$ such that $l(\varphi(\tau_{l}(x),x))=0.$
By condition (2) we deduce that the function $g(t)=l(\varphi(t,x))$ is
increasing. Indeed we have

$\frac{dg(t)}{dt}=\sum_{i=1}^{n}\alpha_{i}\frac{d\varphi_{i}(t,x)}{dt}=\sum_{i=1}^{n}\alpha_{i}F_{i}(\varphi(t,x))>0.$

So $\tau_{l}(x)$ has to be unique zero for $g$ on the interval $[0,\tau]$. We
prove the continuity of $\tau_{l}.$ Let $\varepsilon>0.$ We define
$g(\tau_{l}(x)-\varepsilon)=A_{1}<0$ and $g(\tau_{l}(x)+\varepsilon)=A_{2}>0.$
From the continuity of the flow and the function $l,$ we can pick $\delta>0$
such that for every $y\in B_{Y}(x,\delta)\cap X$ we have

$\sup_{s\in[0,\tau]}|l(\varphi(s,x))-l(\varphi(s,y))|<\frac{\min\\{-A_{1},A_{2}\\}}{2}.$

So for $y\in B_{Y}(x,\delta)\cap X$ we obtain
$l(\varphi(\tau_{l}(x)-\varepsilon,y))<0$ and
$l(\varphi(\tau_{l}(x)+\varepsilon,y))>0.$ Consequently there exists
$\tau_{l}(y)\in\tau_{l}(x)+(-\varepsilon,\varepsilon)$ such that
$l(\varphi(\tau_{l}(y),y))=0,$ so the funtion $\tau_{l}$ is continuous, which
concludes the proof. ∎

If assumption of Lemma 3.1 are satisfied then for every $x\in X$ we have
$l(\varphi(\tau_{l}(x),x))=0.$ Hence we can define the map Poincaré map
$\mathcal{P}:X\to Y$ by the formula $\mathcal{P}(x)=\varphi(\tau_{l}(x),x).$
The map $\mathcal{P}$ is continuous in the norm of the Banach space $Y$. If
$\varphi(t,X^{0})\subset X$ then we can define the Poincaré map
$\mathcal{P}:X^{0}\to Y$ by the formula
$\mathcal{P}(x)=\varphi(\tau_{l}(\varphi(t,x)),\varphi(t,x)).$

We will briefly describe the algorithm which estimates the image of the set
$X^{0}=X^{0}_{P}+X^{0}_{Q}$ by the Poincare map. We assume that the section
$l(x)$ is given by the same formula as in Lemma 3.1 and does not depend on
values in space $H_{Q},$ that is for every $x\in Y$ we have that $l(x)=l(Px).$
The algorithm is following.

  1. (1)

Set $t_{\text{prev}}:=0$ and $t_{\text{curr}}:=0.$

  2. (2)

Check if for every $x\in\varphi(t_{\text{prev}},X^{0})$ there holds $l(x)<0$
and for every $x\in\varphi(t_{\text{curr}},X^{0})$ there holds $l(x)>0.$

    1. (a)

If no, then change $t_{\text{prev}}:=t_{\text{curr}}$ and
$t_{\text{curr}}:=t_{\text{curr}}+\tau,$ where $\tau>0$ is a given time step.
Then go back to step $(2).$

    2. (b)

If yes, then try to minimise the difference $t_{\text{curr}}-t_{\text{prev}}$
for which condition in (2) holds.

  3. (3)

Compute the enclosure $X^{E}$ such that
$\varphi([t_{\text{prev}},t_{\text{curr}}],X^{0})\subset X^{E}.$

  4. (4)

Check that for every $x\in X^{E}$ we have $\sum_{i=1}^{n}\alpha_{i}F_{i}(x)>0$
and if this condition is satisfied return $X^{E}.$

Images $\varphi(t_{\text{prev}},X^{0})$ and $\varphi(t_{\text{curr}},X^{0})$
can be estimated from the algorithm of integration described in Section 2.
Computation of evaluations of $l(\varphi(t_{\text{prev}},X^{0}))$ and
$l(\varphi(t_{\text{prev}},X^{0})))$ depends only on the $P$ part of
representations of sets $\varphi(t_{\text{prev}},X^{0})$ and
$\varphi(t_{\text{curr}},X^{0})$. We refer to [13, Algorithm 1,3] for further
details about evaluation of $l.$ Step (b) is not necessary but it
significantly improves estimation of the algorithm result. For minimization of
the crossing time we can use the rigorous version of the bisection method or
the the rigorous Newton method [13, Algorithm 5, Lemma 8].

Steps (2) and (4) of above algorithm assure that assumptions (1) and (2) of
Lemma 3.1 are satisfied. In step (4) we can additionally return estimation of
the set $[B](X^{E}_{P}-[y^{0}])$ were $[B]$ is some interval matrix and
$[y^{0}]$ is some interval vector. This allows to compute the resulting
estimates in some system of coordinates defined on the section and ignore the
component which is normal to the section. For more details of computation of
$[B](X^{E}_{P}-y^{0})$ see [13, Algorithm 6].

The following result allows us to deduce the existence of periodic orbit.
Assumptions of this theorem can be checked with the use of previously
described algorithm.

######  Theorem 3.2.

Let $X^{0}$ be a nonempty, compact, convex subset of $Y$, such that for every
$x^{0}\in X^{0}$ we have $l(x^{0})=0$. Assume that for $t>0$ and $\tau>0$ the
following holds:

  1. (1)

For every $x^{0}\in X^{0}$ we have $l(\varphi(t,x^{0}))<0$ and
$l(\varphi(t+\tau,x^{0}))>0.$

  2. (2)

For every $x\in\varphi([t,t+\tau],X^{0}),$ we have
$\sum_{i=1}^{n}\alpha_{i}F_{i}(x)>0.$

  3. (3)

For every $x\in\varphi([t,t+\tau],X^{0})$ such that $l(x)=0$ we have $x\in
X^{0}.$

Then there exist $x^{*}\in\varphi([t,t+\tau],X^{0})$ and $T\in(t,t+\tau)$ such
that $l(x^{*})=0$ and $\varphi(T,x^{*})=x^{*}$.

###### Proof.

From Lemma 3.1 we observe that for every $x\in X^{0}$ we have
$l(\varphi(\tau_{l}(\varphi(t,x)),\varphi(t,x)))=0.$ Hence we can define the
map $\mathcal{P}:X^{0}\to Y$ by the formula
$\mathcal{P}(x^{0})=\varphi(\tau_{l}(\varphi(t,x)),\varphi(t,x)).$ The map
$\mathcal{P}$ is continuous in the norm of the Banach space $Y$. From
assumption (3) we see that $\mathcal{P}(X^{0})\subset X^{0}.$ As $X^{0}$ is
compact and convex, the Schauder fixed-point theorem ensures the existence of
a fixed point $x^{*}$. For this point we have
$x^{*}\in\mathcal{P}(X^{0})\subset\varphi([t,t+\tau],X^{0})$ which concludes
the proof. ∎

######  Remark 3.1.

If in Theorem 3.2 we additionally assume that for some $t_{1}<t_{2}<t+\tau$:

  1. (1)

For every $x\in\varphi([0,t_{1}],X_{0})$ and
$x\in\varphi([t_{2},t+\tau],X_{0})$ we have that
$\sum_{i=1}^{n}\alpha_{i}F_{i}(x)>0$ .

  2. (2)

For every $t\in[t_{1},t_{2}]$ we have $\varphi(t,X^{0})\cap X^{0}=\emptyset$.

Then $T$ has to be a fundamental period for $x^{*}$.

###  3.3. Numerical approximation of periodic orbit

The first step in constructing the initial data and the section which will be
used in validating the assumptions of Theorem 3.2, is finding the
approximation of periodic using the Galerkin projection of (2.5). This means
that we need to find an initial data $x^{*}=(u^{*},v^{*})$ and a time $T^{*}$
such that the solution to the system

$\begin{cases}\frac{d}{dt}P(u(t),v(t))=LP(u(t),v(t))+Pf(P(u(t),v(t))),\\\
(u(0),v(0))=(u^{*},v^{*}).\end{cases}$ (3.1)

is close to a periodic solution and $T^{*}$ is close to its period.

(a) Projection of the numerical approximation of periodic orbit into two first
modes $u_{1},v_{1},$ for parameters $d_{1},d_{2},A,B$ from Theorem 1.1.

(b) Blue function corresponds to function $u^{*}(x)$ and orange to the
function $v^{*}(x),$ around which we build the set of initial data.

Figure 2. First two modes and initial data for the approximate periodic
solution found by nonrigorous computations.

As, for the Brusselator system, we observe that periodic orbit is numerically
attracting, it is enough to find approximation of attracting fixed point of
some Poincare map. We need to ensure that the section which defines this
Poincare map intersects periodic orbit of (3.1). Additionally we are searching
for $(u^{*},v^{*})$ for which only odd Fourier coefficients in the sine series
are nonzero. We chose to project the Brusselator system (1.1) on the subspace
of $L^{2}\times L^{2}$ spanned by the functions
$\\{(\sin(x),0),(\sin(3x),0),\ldots,(\sin(17x),0)\\}$ and
$\\{(0,\sin(x)),(0,\sin(3x)),\ldots,(0,\sin(17x))\\}.$

With this procedure we have found the following approximations of the fixed
point of Poincaré map.

$\displaystyle u^{*}(x)$
$\displaystyle=10^{-1}*6.999\sin(x)-10^{-2}*8.170\sin(3x)+-10^{-3}*5.377\sin(5x)+10^{-2}*1.325\sin(7x)+$
$\displaystyle+10^{-3}*1.050\sin(9x)-10^{-4}*2.585\sin(11x)-10^{-6}*1.764\sin(13x)+10^{-7}*5.029\sin(15x)$
$\displaystyle+10^{-8}*2.779\sin(17x),$

and

$\displaystyle v^{*}(x)$
$\displaystyle=3.869\sin(x)+1.136\sin(3x)+10^{-1}*1.017\sin(5x)-10^{-3}*9.291\sin(7x)$
$\displaystyle-10^{-3}*1.297\sin(9x)+10^{-4}*1.960\sin(11x)+10^{-5}*1.993\sin(13x)-10^{-6}*4.109\sin(15x)$
$\displaystyle-10^{-7}*3.147\sin(17x).$

The coefficients of $u^{*}$ and $v^{*}$ are written up to three decimal
places.

###  3.4. Defining a Poincaré map and constructing an initial set

We need to define a Poincaré map and an initial set that will be used to
validate Theorem 1.1. The section is a mapping $l:Y\to\mathbb{R}$ given by the
formula $l(x)=\sum_{i=1}^{n}\alpha_{i}(x_{i}-x^{*}_{i}),$ where the $x^{*}$ is
the point from the numerical approximation of the periodic orbit. The numbers
$\alpha_{i}\in\mathbb{R}$ are chosen to assure that transversality condition
(2) of Lemma 3.1 holds. The method to choose the coefficients in the optimal
way, such that the time of passing through the section is minimal is given in
[13, Theorem 18]. We use this approach. To describe the initial set
$X^{0}=X_{P}^{0}+X_{Q}^{0}$ we define separately two sets $X_{P}^{0}$ and
$X_{Q}^{0}.$ The set $X_{P}^{0}$ is defined as follows

$X_{P}^{0}=x^{*}+Ar^{0},$

where $x^{*}$ is a (noninterval) vector, $A$ is a (noninterval) square matrix
and $r^{0}$ is an interval vector. Note, that since both $r^{0}$ is an
interval vector, $x^{*}$ is a vector, and $A$ is a matrix, $X_{P}^{0}$ is a
bounded and convex set.

The first column $c_{1}$ of the matrix $A$ is equal to
$(\alpha_{1},\ldots,\alpha_{n}).$ The subsequent columns $c_{2},\ldots,c_{n}$
of this matrix constitute the coordinate system on section. They are assumed
to satisfy the following two conditions

  1. (1)

For every $i\in\\{2,\ldots,n\\}$ we have $l(c_{i}+x^{*})=0,$

  2. (2)

For the linear functional $\hat{l}(x)=l(x+x^{*})$ we have that
$\text{span}\\{c_{2},\ldots,c_{n}\\}=P(\ker\hat{l}).$

The first coefficient of $r^{0}$ can be set to interval $[0,0]$ as our initial
set should be on the previously defined section. The next coefficients
describe size of the set in the section and can be fixed for example to
intervals $[-\delta,\delta]$ where $\delta>0$. In the algorithm of the
evolution of the set, we will take as the initial data the set
$\overline{X}_{P}^{0}$ which is a superset of $X_{P}^{0}$ and is defined by
$\overline{X}_{P}^{0}=[x^{*}]+[A]r^{0}$, where $[A]$ is an interval matrix
containing $A$ and $[x^{*}]$ is an interval vector containing $x^{*}$. The
first column of $[A]$ is the interval vector containing $c_{1}$. The remaining
interval columns denoted by $\\{C_{i}\\}_{i=2}^{n}$ are constructed in such a
way that they are guaranteed to contain vectors which constitute the
coordinate system $\\{c_{i}\\}_{i=2}^{n}$ on the section, i.e. $c_{i}\in
C_{i}$ for $i\in\\{2,\ldots,n\\}$. In construction of columns
$C_{2},\ldots,C_{n}$ of the matrix $[A]$ we can use numerical approximation of
the eigenvectors of the Poincaré map given by the section $l$. Additionally
every numerical approximation of these eigenvectors is rigorously projected on
the first column in order to assure that condition $(1)$ is satisfied for a
certain $c_{i}\in C_{i}$. We also compute an interval matrix $[B]$ which is a
rigorous interval inverse of the interval matrix $[A].$ Existence of this
matrix ensures that condition (2) is satisfied. Matrix $[B]$ is also used in
further steps. The set $X^{0}_{Q}$ is defined by infinite interval vectors
with the polynomial decay of coefficients as presented in (2.8).

###  3.5. Computer assisted proof

With the previously defined section $l$ and the initial data $X^{0}$ we check
that the assumptions of Theorem 3.2 are satisfied. We use the algorithm
described in Section 3.2 to estimate the image $\mathcal{P}(X^{0})\subset
X^{1}_{P}+X^{1}_{Q}.$ This algorithm guarantees that assumptions (1) and (2)
of Theorem 3.2 are satisfied. The set $X^{1}_{P}$ is returned in the form

$X^{1}_{P}=[x^{*}]+[A]q^{0}.$

To obtain the interval vector $q^{0}$ we use the algorithm of computing a
Poincaré map with the matrix $[B],$ which is an interval inverse of the matrix
$[A],$ and $[y^{0}]$ equal to $[x^{*}]$. To validate the assumption (3) of
Theorem 3.2 it is enough to check that

  * (C1)

we have $q^{0}_{i}\subset r^{0}_{i}$ for $i\in\\{2,\ldots,n\\},$

  * (C2)

we have $X^{1}_{Q}\subset X^{0}_{Q}.$

In the computer assisted proof of Theorem 1.1 we set

$H_{P}:=\text{span}\bigcup_{\stackrel{{\scriptstyle 1\leq k\leq
59}}{{k\,\textrm{mod}\,2=1}}}\\{(\sin(kx),0),(0,\sin(kx))\\},$

We define $H_{Q}$ as the orthogonal complement of $H_{P},$ in the space

$H:=\left\\{(u,v)\in L^{2}\times
L^{2}:\,u(x)=\sum_{k=1}^{\infty}u_{2k-1}\sin((2k-1)x),\,v(x)=\sum_{k=1}^{\infty}v_{2k-1}\sin((2k-1)x)\right\\}.$

Observe that $\text{dim}(H_{P})=60.$ The set $X_{P}^{0}$ is defined using the
procedure described in Section 3.4 applied to the Brusselator system. This is
a bounded, closed and convex set in a finite dimensional space $H_{P}$. To see
that $q^{0}_{i}\subset r^{0}_{i}$ implies that $X^{0}_{P}\subset X^{1}_{P}$
note that $q^{0}\subset r^{0}$ implies that

$[B](X_{P}^{1}-[x^{*}])=q^{0}\subset r^{0}=A^{-1}(X^{0}_{P}-x^{*}).$

But, since $A^{-1}\in[B]$ and $x^{*}\in[x^{*}]$, the last inclusion implies
that

$A^{-1}(X_{P}^{1}-x^{*})\subset A^{-1}(X^{0}_{P}-x^{*}),$

which guarantees the required inclusion $X_{P}^{1}\subset X_{P}^{0}$. The part
$X_{Q}^{0}$ is given as follows

$X_{Q}^{0}:=\left\\{(u,v)\in
H^{Q}:\;u_{k}\in\frac{[-1,1]}{|k|^{5}},\;v_{k}\in\frac{[-1,1]}{|k|^{5}}\quad\text{for}\;k>59\right\\}.$

The set $X^{0}=X_{P}^{0}+X_{Q}^{0}$ must be compact in $Y=C_{0}\times C_{0}$.
As it is closed in $Y$, it is sufficient to show that it is bounded in
$H^{1}_{0}\times H^{1}_{0}$. To this end assume that $(u,v)\in X^{0}$. Then

$\displaystyle(u_{x},v_{x})=\left(\sum_{\underset{k\;\textrm{mod}\;2=1}{k=1}}^{\infty}ku_{k}\sin(kx),\sum_{\underset{k\;\textrm{mod}\;2=1}{k=1}}^{\infty}ku_{k}\sin(kx)\right),$

and

$\displaystyle\|(u,v)\|_{H^{1}_{0}\times
H^{1}_{0}}^{2}=\|(u_{x},v_{x})\|_{L^{2}\times
L^{2}}^{2}=\frac{\pi}{2}\sum_{\underset{k\;\textrm{mod}\;2=1}{k=1}}^{\infty}k^{2}(|u_{k}|^{2}+|v_{k}|^{2})$
$\displaystyle\ \ \ \ \
\leq\frac{59^{2}\pi}{2}\sum_{\underset{k\;\textrm{mod}\;2=1}{k=1}}^{59}(|u_{k}|^{2}+|v_{k}|^{2})+\frac{\pi}{2}\sum_{\underset{k\;\textrm{mod}\;2=1}{k=61}}^{\infty}\frac{2}{k^{8}}.$

Since the last quantity is bounded uniformly with respect to the choice of
$(u,v)\in X^{0}$, we get the required compactness of this set in $Y$. After
the computation of Poincaré map we get the set $X_{P}^{1}+X_{Q}^{1}$ with the
following $Q$ part

$X_{Q}^{1}:=\left\\{(u,v)\in H^{Q}:\;u_{i}\in
10^{-14}\frac{[-3.46474,3.46474]}{|i|^{6.5}},\;v_{i}\in
10^{-13}\frac{[-3.9024,3.9024]}{|i|^{6.5}}\quad\text{for}\;i>59\right\\}.$

The conditions (C1), (C2) are satisfied, so the main Theorem 1.1 is validated.
The fact that (C1) holds is demonstrated in Tab. 1 where for brevity only the
first $10$ coordinates are depicted. Due to the extra information on the
periodic orbit obtained in the computer assisted proof, we can provide the
following extended version of the main theorem.

######  Theorem 3.3.

For parameters $d_{1}=0.2,\;d_{2}=0.02,\;A=1,\;B=2$ the Brusselator system has
a periodic solution $(\bar{u}(t,x),\bar{v}(t,x)),$ with the period
$T\in[7.69666,7.69667].$ The functions $\bar{u}(t,x),\bar{v}(t,x)$ are
symmetric with respect to the point $x=\frac{\pi}{2}.$ Moreover the following
estimates are true

$\displaystyle\sup_{t\in[0,T]}\left\|{\bar{u}(t)}\right\|_{L^{2}}$
$\displaystyle\leq
1.27261,\;\sup_{t\in[0,T]}\left\|{\bar{v}(t)}\right\|_{L^{2}}\leq 5.05587,$
$\displaystyle\sup_{t\in[0,T]}\left\|{\bar{u}_{x}(t)}\right\|_{L^{2}}$
$\displaystyle\leq
1.35194,\;\sup_{t\in[0,T]}\left\|{\bar{v}_{x}(t)}\right\|_{L^{2}}\leq
6.75405,$
$\displaystyle\sup_{t\in[0,T]}\left\|{\bar{u}(t)-{u}^{*}(t)}\right\|_{L^{2}}$
$\displaystyle\leq
0.00049664,\;\sup_{t\in[0,T]}\left\|{\bar{v}(t)-{v}^{*}(t)}\right\|_{L^{2}}\leq
0.000955005,$
$\displaystyle\sup_{t\in[0,T]}\left\|{\bar{u}_{x}(t)-{u}_{x}^{*}(t)}\right\|_{L^{2}}$
$\displaystyle\leq
0.000546005,\;\sup_{t\in[0,T]}\left\|{\bar{v}_{x}(t)-{v}_{x}^{*}(t)}\right\|_{L^{2}}\leq
0.00141263,$

where $(u^{*}(t,x),v^{*}(t,x))$ is the solution to the Brusselator system with
the initial data $u^{*}(0,x)=u^{*}(x),\;v^{*}(0,x)=v^{*}(x),$ and the same
parameters as above.

|  $\displaystyle r_{0}$ |  $\displaystyle q_{0}$
---|---|---
1. |  $\displaystyle[0,0]$ |  $\displaystyle 10^{-6}[-1.10436,1.21675]$
2. |  $\displaystyle 10^{-5}[-1\ ,\ 1]$ |  $\displaystyle 10^{-6}[-6.87377,6.87371]$
3. |  $\displaystyle 10^{-5}[-1\ ,\ 1]$ |  $\displaystyle 10^{-6}[-1.1776,1.17766]$
4. |  $\displaystyle 10^{-5}[-1\ ,\ 1]$ |  $\displaystyle 10^{-11}[-2.8735,-2.36735]$
5. |  $\displaystyle 10^{-5}[-1\ ,\ 1]$ |  $\displaystyle 10^{-9}[-5.40597,5.43755]$
6. |  $\displaystyle 10^{-5}[-1\ ,\ 1]$ |  $\displaystyle 10^{-9}[-1.60821,1.59499]$
7. |  $\displaystyle 10^{-5}[-1\ ,\ 1]$ |  $\displaystyle 10^{-10}[-1.56502,1.59018]$
8. |  $\displaystyle 10^{-5}[-1\ ,\ 1]$ |  $\displaystyle 10^{-11}[-4.72352,6.65989]$
9. |  $\displaystyle 10^{-5}[-1\ ,\ 1]$ |  $\displaystyle 10^{-11}[-5.29417,-1.37926]$
10. |  $\displaystyle 10^{-5}[-1\ ,\ 1]$ |  $\displaystyle 10^{-11}[0.131355,2.59156]$
Table 1. First 10 coordinates of interval vectors $r_{0}$ and computed $q_{0}$
which correspond to the sets $X^{0}_{P}$ and $X^{1}_{P}$ in the computer
assisted proof of Theorem 1.1 for parameters $d_{1}=0.2,d_{2}=0.02,A=1,B=2.$

##  4\. Numerical and rigorous results for other parameter values

In this section we discuss some numerical observations and rigorous results
concerning the Brusselator system with various parameter values. We conduct
the computer assisted proof of existence of periodic orbit for the parameter
values from the set
$\mathcal{A}=\mathcal{A}_{1}\cup\mathcal{A}_{2}\cup\mathcal{A}_{3},$ of
parameters $(d_{1},d_{2},A,B)$, where

$\displaystyle\mathcal{A}_{1}=\left\\{\left(0.2,0.02,1,2+\frac{i}{10}\right):i\in\\{0,\ldots,11\\}\right\\},$
$\displaystyle\mathcal{A}_{2}=\left\\{\left(1,\frac{1}{64},1,2.71\right),\left(1,\frac{1}{64},1,2.83\right),\left(1,\frac{1}{64},1,2.84\right)\right\\},\quad\mathcal{A}_{3}={\\{(0.02,0.02,1,2)\\}}.$

We have rigorously proved the existence of the periodic orbits for all
parameters from the set $\mathcal{A}$. Values in the set $\mathcal{A}_{1}$
correspond to the slow-fast behavior of the system, ones in the set
$\mathcal{A}_{2}$ correspond to the period-doubling bifurcation and cross-
validation with the results of Arioli [1], and the paramters in
$\mathcal{A}_{3}$ are related with the numerical experiments where we find
attracting torus.

In the grid set $\mathcal{A}_{1}$ we fix the parameters $d_{1},d_{2},A$ and we
increase the parameter $B$ by $0.1$ starting from its value $B=2$, which
corresponds to the periodic orbit found in our main Theorem 1.1, and ending at
the value $B=3.1$. For the corresponding planar ODE (1.2), it is known [16]
that the slow-fast dynamics of the system increases when the parameter $B$
grows. We observe the same intensification of slow-fast behavior upon the
increase of parameter $B$ also for the Brusselator with diffusion. This is
depicted in Fig. 3 and 4. The slow-fast behavior is especially visible in the
higher Fourier modes, cf. Fig. 5. We also stress that as the value of $B$
increases, we require more modes to accurately represent the solution.
Therefore, to successfully carry out the computer-assisted proof, we must
increase the dimension of the inclusion as the value of $B$ increases, and,
consequently the computation time lengthens. This is shown in Tab.2.

Value of parameter $\displaystyle B$ |  Number of variables in the inclusion (2.17) |  Computation time in seconds
---|---|---
$\displaystyle 2.1$ |  $\displaystyle 14$ |  $\displaystyle 298$ $\displaystyle s$
$\displaystyle 2.3$ |  $\displaystyle 22$ |  $\displaystyle 332$ $\displaystyle s$
$\displaystyle 2.5$ |  $\displaystyle 22$ |  $\displaystyle 341$ $\displaystyle s$
$\displaystyle 2.7$ |  $\displaystyle 30$ |  $\displaystyle 477\ s$
$\displaystyle 2.9$ |  $\displaystyle 32$ |  $\displaystyle 534\ s$
$\displaystyle 3.1$ |  $\displaystyle 40$ |  $\displaystyle 771\ s$
Table 2. Selected values of parameter $B$ used in the computer assisted proof
of Theorem 4.1, the dimensions of corresponding inclusion (2.17) and the
computation times needed to rigorously validate the periodic orbit existence.
The computations were single-threaded and realized on the processor Intel(R)
Core(TM) i5-4200M CPU @ 2.50GHz.

The result on the existence of the priodic orbits for parameters in
$\mathcal{A}_{1}$ is contained in the next theorem.

######  Theorem 4.1.

The Brusselator system has a periodic orbit for
$(d_{1},d_{2},A,B)\in\mathcal{A}_{1}$.

Figure 3. Plots of projection of numerical approximation of the periodic
orbits into first modes of $u$ and $v$ (left) and the evolution in time of the
first Fourier modes (right) for $d_{1}=0.2$, $d_{2}=0.02$, $B=2.1$ (top) and
$B=2.3$ (bottom), $A=1.$

Figure 4. Plots of projection of numerical approximation of the periodic
orbits into first modes of $u$ and $v$ (left) and the evolution in time of the
first Fourier modes (right) for
$d_{1}=0.2,\;d_{2}=0.02,\;B=2.5,\;2.7,\;2.9,\;3.1$ (from top to bottom),
$A=1.$

Figure 5. The evolution of fourth mode of $u$ for $B=2.1$ (left picture) and
$B=3.1$ (right picture).

In [1, Theorem 2] Arioli used computer assisted method and proved the
existence of periodic orbit for parameters
$d_{1}=1,\,d_{2}=\frac{1}{64},\,B\in[2.6993,2.7419],\,A=1.$ We have checked
his result for $B=2.71$, and we found that both results correspond to each
other in terms of the found period of the orbit. Our resuls, which succesfully
cross-validates both approaches, is contained in the following theorem.

######  Theorem 4.2.

For parameters $d_{1}=1,\;d_{2}=\frac{1}{64},\;A=1,\;B=2.71$ the Brusselator
system has a periodic solution $(\bar{u}(t,x),\bar{v}(t,x)),$ with the period
$T\in[10.4549,10.455].$ The functions $\bar{u}(t,x),\bar{v}(t,x)$ are
symmetric with respect to the point $x=\frac{\pi}{2}.$ Moreover the following
estimates are true

$\displaystyle\sup_{t\in[0,T]}\left\|{\bar{u}(t)}\right\|_{L^{2}}$
$\displaystyle\leq
0.600569,\;\sup_{t\in[0,T]}\left\|{\bar{v}(t)}\right\|_{L^{2}}\leq 5.05587,$
$\displaystyle\sup_{t\in[0,T]}\left\|{\bar{u}_{x}(t)}\right\|_{L^{2}}$
$\displaystyle\leq
10.0529,\;\sup_{t\in[0,T]}\left\|{\bar{v}_{x}(t)}\right\|_{L^{2}}\leq
11.5792,$
$\displaystyle\sup_{t\in[0,T]}\left\|{\bar{u}(t)-{u}^{*}(t)}\right\|_{L^{2}}$
$\displaystyle\leq
10^{-5}*8.69037,\;\sup_{t\in[0,T]}\left\|{\bar{v}(t)-{v}^{*}(t)}\right\|_{L^{2}}\leq
0.000295375,$
$\displaystyle\sup_{t\in[0,T]}\left\|{\bar{u}_{x}(t)-{u}_{x}^{*}(t)}\right\|_{L^{2}}$
$\displaystyle\leq
10^{-5}*8.99234,\;\sup_{t\in[0,T]}\left\|{\bar{v}_{x}(t)-{v}_{x}^{*}(t)}\right\|_{L^{2}}\leq
0.000450742,$

where $(u^{*}(t,x),v^{*}(t,x))$ is the solution to the Brusselator system with
the initial data

$\displaystyle u^{*}(0,x)$
$\displaystyle=0.43\sin(x)-0.0231361\sin(3x)-0.00129933\sin(5x)+10^{-5}*7.48643\sin(7x)$
$\displaystyle+10^{-6}*7.466799\sin(9x)-10^{-7}*2.998759\sin(11x)-10^{-8}*3.89921\sin(13x)$
$\displaystyle+10^{-9}*1.15918\sin(15x)+10^{-10}*1.98398\sin(17x),$
$\displaystyle v^{*}(0,x)$
$\displaystyle=7.85996\sin(x)+1.59666\sin(3x)+0.091348\sin(5x)-0.0041776\sin(7x)$
$\displaystyle-10^{-4}*4.73146\sin(9x)+10^{-5}*1.66984\sin(11x)+10^{-6}*2.42261\sin(13x)$
$\displaystyle-10^{-8}*6.5405\sin(15x)-10^{-8}*1.23104\sin(17x),$

and the same parameters as above.

(a) Projection of the numerical approximation of periodic orbit on two first
modes $u_{1},v_{1},$ for parameters from Theorem 4.2.

(b) Blue plot corresponds to $u^{*}(0,x)$ and orange to the $v^{*}(0,x)$ in
Theorem 4.2.

Figure 6. Periodic orbit from Theorem 4.2, i.e. for $B=2.71$.

Figure 7. Numerically attracting periodic orbit for $B=2.83$ (left) and
$B=2.84$ (right).

Arioli observed a period doubling bifurcation, a phenomenon that cannot occur
in the planar ODE (1.2). Thus, the dynamics of (1.1) is expected to be more
complicated than that of the planar ODE (1.2). Although we do not rigorously
prove the bifurcation, we show that the minimal period of the found orbits
approximately doubles with a small increase in the parameter $B$, as seen in
Theorem 4.3. Specifically, as we increase $B$ and keep other parameters fixed
$d_{1}=1,\;d_{2}=\frac{1}{64},\;A=1$, from numerical simulations we observe
that the system goes through the period doubling bifurcation. The bifurcation
appears to occur between $B=2.83$ and $B=2.84$. We have proven the following
theorem about periodic orbits for these parameters.

######  Theorem 4.3.

For $d_{1}=1,\;d_{2}=\frac{1}{64}\;,$ $B=2.83\;A=1$ there exist a periodic
orbit with its fundamental period in the interval $[13.2128,13.2130]$. For
$d_{1}=1,\;d_{2}=\frac{1}{64}\;,$ $B=2.84,\;A=1$ there exist a periodic orbit
with its fundamental period in the interval $[27.2436,27.2439]$.

Note that we only obtain rigorous computer assisted proofs for the parameter
values $B$ before and after the expected bifurcation. We believe that
conducting computer assisted proof of the bifurcation existence would be
interesting and challenging problem.

Other nontrivial dynamics of the Brusselator PDE (1.1) was investigated in
[8], where the numerical evidence on the existence of 2-dimensional attracting
tori was shown. We observe, only numerically, the same phenomenon, which
exists only in PDE Brusselator model and not in the corresponding planar ODE
(1.2). First, we rigorously observe the existence of the periodic orbit when
the diffusion rates $d_{1}$ and $d_{2}$ are equal to each other. We have
proved that for parameter values $d_{1}=0.02,d_{2}=0.02,B=2,A=1$. If we
decrease the diffusion rates $d_{1}$ and $d_{2}$ (keeping them equal to each
other) we numerically observe the emergence of two dimensional attracting
torus, cf. Fig. 8. This is a numerical confirmation of the phenomenon which
was firstly observed in [8]. Note that the same paper also contains the
numerical evidence for the presence of chaos in the same system of equations.
These dynamical phenomena are further interesting challenge for computer
assisted proofs for the Brusselator system.

Figure 8. Two dimensional numerically attracting torus observed in for
$d_{1}=0.009,\;d_{2}=0.009,\;B=2.1,\;A=1.$

##  5\. Algebra

This section contains the technical results on the estimates of the
convolutions of the sine and cosine Fourier series represented with uniform
estimates on the decay of coefficients. Such operations are performed during
the calculations in the algorithm as part of the computer-assisted proofs,
where we require uniform estimates on the tails of the series. We work with
sequences $\\{u_{i}\\}_{i=1}^{\infty}$ such that first $n$ coefficients are
given by some numbers or intervals and the remainders of the sequence satisfy
$u_{i}\in\frac{[C_{u}^{-},C_{u}^{+}]}{i^{s}}$ for some $C_{u}^{-}\leq
C_{u}^{+}$ and $s\in\mathbb{R}.$

The decay of the Fourier coefficients for smooth periodic functions is related
with their regularity: if a periodic function $u:\mathbb{R}\to\mathbb{R}$ is
of class $C^{s}$, then its coefficients must decay as $\frac{1}{i^{s}}.$
Clearly, the product of two $C^{s}$ functions also has regularity $C^{s}$.
This is related with our results of this section, which state that if two
functions, represented in the sine or cosine Fourier series have some given
decay of the Fourier coefficients of the form $O\left(\frac{1}{i^{s}}\right)$
for $s>1$ then their product must have the same decay. Moreover we provide the
exact estimates for the Fourier coefficients of the product: such estimates
are needed in rigorous computations of nonlinear polynomial terms present in
the equations.

The following results are basic, so we skip the proofs. We use them several
times is the following considerations.

######  Proposition 5.1.

Let $s>1.$ The following inequality holds

$\sum_{i=1+n}^{\infty}\frac{1}{i^{s}}\leq\frac{n^{1-s}}{s-1}.$ (5.1)

######  Proposition 5.2.

Let $s>1$. If $a,b>0$ then the following inequality holds

$(a+b)^{s}\leq 2^{s-1}(a^{s}+b^{s}).$ (5.2)

The following lemma gives us the estimates on the result of multiplication of
$u$ and $v$ which are both represented in the sine Fourier series. The first
$n$ coefficients of $u$ and $v$ are given explicitly and the coefficients
indexed by numbers larger than $n$ are expressed by the polynomial decay.

######  Lemma 5.1.

Assume that

$u(x)=\sum_{i=1}^{\infty}u_{i}\sin(ix),\quad
v(x)=\sum_{i=1}^{\infty}v_{i}\sin(ix).$ (5.3)

Moreover we assume that for some $n\in\mathbb{N}$ and $s>1$ the following
estimates hold

$u_{i}\in\frac{C_{u}[-1,1]}{i^{s}}\quad\textrm{and}\quad
v_{i}\in\frac{C_{v}[-1,1]}{i^{s}}\quad\text{for $i>n,$}$ (5.4)

where $C_{v}>0,C_{u}>0$. Then

$(uv)(x)=(uv)_{0}+\sum_{k=1}^{\infty}(uv)_{k}\cos(kx),$ (5.5)

with

$(uv)_{0}=\frac{1}{2}\sum_{i=1}^{\infty}u_{i}v_{i}\quad(uv)_{k}=\frac{1}{2}\sum_{i=1}^{\infty}u_{i+k}v_{i}+\frac{1}{2}\sum_{i=1}^{\infty}u_{i}v_{i+k}-\frac{1}{2}\sum_{i=1}^{k-1}u_{i}v_{k-i}.$
(5.6)

and the following estimates hold for the coefficients of the product. For
$k=0$ we have

$(uv)_{0}\in\frac{1}{2}\sum_{i=1}^{n}u_{i}v_{i}+\frac{C_{u}C_{v}}{2}\frac{n^{1-2s}}{2s-1}[-1,1],$
(5.7)

for $1\leq k\leq 2n$ we have

$(uv)_{k}\in\frac{1}{2}\sum_{i=1}^{n}u_{i}v_{i+k}+\frac{1}{2}\sum_{i=1}^{n}u_{i+k}v_{i}-\frac{1}{2}\sum_{i=1}^{k-1}u_{i}v_{k-i}+C_{u}C_{v}\frac{n^{1-2s}}{2s-1}[-1,1],$
(5.8)

and for $k>2n$ we have

$(uv)_{k}\in\frac{D[-1,1]}{k^{s}},\quad
D=\frac{1}{2}\left(\sum_{i=1}^{n}(C_{u}|v_{i}|+C_{v}|u_{i}|)\left(1+\left(\frac{2n+1}{2n+1-i}\right)^{s}\right)+C_{u}C_{v}(2+2^{s})\frac{n^{1-s}}{s-1}\right).$
(5.9)

The above lemma provides explicit formulas for estimating the coefficients
with indexes $0$ to $2n$ of the cosine Fourier series. This is motivated by
the fact that if the sine expansions of $u$ and $w$ are finite and concentrate
on the first $n$ coefficients, then the representation of $uv$ is also finite
and only first $2n+1$ coefficients are nonzero. In the formulas in the above
lemma, whenever the terms $u_{i}$ and $v_{i}$ appear with $i>n$ in the
computation, we substitute them with the intervals $\frac{C_{u}}{i^{s}}[-1,1]$
and $\frac{C_{v}}{i^{s}}[-1,1]$, respectively. Additionally, we can use other
available estimates of $u_{i}$. For example we can use the fact that
$u_{i}\in\frac{[C_{u}^{-},C_{u}^{+}]}{i^{s}}$ (with upper and lower bound
different from each other), or that $u_{i}$ is zero for odd coefficients.

The coefficients for $k>2n$ of the cosine expansion in the above lemma are
given by the uniform polynomial decay with the same rate as the sine series
for $u$ and $v.$

###### proof of Lemma 5.2.

We have

$\displaystyle uv$
$\displaystyle=\left(\sum_{i_{1}=1}^{\infty}\sin(i_{1}x)\right)\left(\sum_{i_{2}=1}^{\infty}\sin(i_{2}x)\right)=\frac{1}{2}\sum_{i_{1},i_{2}=1}^{\infty}u_{i_{1}}v_{i_{2}}\cos((i_{1}-i_{2})x)-\frac{1}{2}\sum_{i_{1},i_{2}=1}^{\infty}u_{i_{1}}v_{i_{2}}\cos((i_{1}+i_{2})x)$
$\displaystyle=\frac{1}{2}\sum_{i_{1}=i_{2}}u_{i_{1}}v_{i_{2}}+\frac{1}{2}\sum_{k=1}^{\infty}\left(\sum_{i_{1}-i_{2}=k}u_{i_{1}}v_{i_{2}}+\sum_{i_{1}-i_{2}=-k}u_{i_{1}}v_{i_{2}}-\sum_{i_{1}-i_{2}=k}u_{i_{1}}v_{i_{2}}\right)\cos(kx).$

We express all coefficients of the resultant cosine series separately, using
the formulas

$(uv)_{0}=\frac{1}{2}\sum_{i=1}^{\infty}u_{i}v_{i},$

and

$(uv)_{k}=\frac{1}{2}\sum_{i=1}^{\infty}u_{i+k}v_{i}+\frac{1}{2}\sum_{i=1}^{\infty}u_{i}v_{i+k}-\frac{1}{2}\sum_{i=1}^{k-1}u_{i}v_{k-i},$

for $k\geq 1.$ Observe that for $i>n$ there holds

$|u_{i}v_{i}|\leq\frac{C_{v}C_{u}}{i^{2s}}.$

So, we obtain

$(uv)_{0}=\sum_{i=1}^{n}u_{i}v_{i}+\sum_{i=n+1}^{\infty}u_{i}v_{i}\in\sum_{i=1}^{n}u_{i}v_{i}+C_{u}C_{v}\sum_{i=n+1}\frac{1}{i^{2s}}[-1,1]\subset\sum_{i=1}^{n}u_{i}v_{i}+C_{u}C_{v}\frac{n^{1-2s}}{-1+2s}[-1,1].$

Now, observe that for $i>n$ the following estimates hold

$|u_{i}v_{i+k}|\leq\frac{C_{v}C_{u}}{i^{2s}},\quad|u_{i+k}v_{i}|\leq\frac{C_{v}C_{u}}{i^{2s}},$

whereas we deduce that for $1\leq k\leq 2n$ we have

$\displaystyle(uv)_{k}$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{n}u_{i+k}v_{i}+\frac{1}{2}\sum_{i=1}^{n}u_{i}v_{k+i}-\frac{1}{2}\sum_{i=1}^{k-1}u_{i}v_{k-i}+\frac{1}{2}\sum_{i=n+1}^{\infty}(u_{i+k}v_{i}+u_{i}v_{i+k})$
$\displaystyle\in\frac{1}{2}\sum_{i=1}^{n}u_{i+k}v_{i}+\frac{1}{2}\sum_{i=1}^{n}u_{i}v_{k+i}-\frac{1}{2}\sum_{i=1}^{k-1}u_{i}v_{k-i}+C_{v}C_{u}\frac{n^{1-2s}}{2s-1}[-1,1].$

Finally, for $k>2n$ we need to estimate the sums in the expression

$\displaystyle(uv)_{k}=\frac{1}{2}$
$\displaystyle\Bigg{[}\sum_{i=1}^{n}u_{i+k}v_{i}+\sum_{i=1}^{n}u_{i}v_{i+k}+\sum_{i=n+1}^{\infty}u_{i+k}v_{i}+\sum_{i=n+1}^{\infty}u_{i}v_{i+k}$
$\displaystyle+\sum_{i=1}^{n}u_{k-i}v_{i}+\sum_{i=k-n}^{k-1}u_{k-i}v_{i}+\sum_{i=n+1}^{k-n-1}u_{k-i}v_{i}\Bigg{]}.$

To this end observe, that we have

$\sum_{i=1}^{n}|u_{i+k}v_{i}|\leq\sum_{i=1}^{n}\frac{C_{u}}{(i+k)^{s}}|v_{i}|\leq\frac{C_{u}}{k^{s}}\sum_{i=1}^{n}|v_{i}|,\quad\sum_{i=1}^{n}|u_{i}v_{i+k}|\leq\frac{C_{v}}{k^{s}}\sum_{i=1}^{n}|u_{i}|,$
$\sum_{i=n+1}^{\infty}|u_{i+k}v_{i}|\leq\sum_{i=n+1}^{\infty}\frac{C_{u}C_{v}}{i^{s}(i+k)^{s}}\leq\frac{C_{u}C_{v}}{k^{s}}\sum_{i=n+1}^{\infty}\frac{1}{i^{s}}\leq\frac{C_{u}C_{v}}{k^{s}}\frac{n^{1-s}}{s-1},\
\sum_{i=n+1}^{\infty}|u_{i}v_{i+k}|\leq\frac{C_{u}C_{v}}{k^{s}}\frac{n^{1-s}}{s-1}.$

We also obtain

$\sum_{i=1}^{n}|u_{k-i}v_{i}|\leq\frac{C_{u}}{k^{s}}\sum_{i=1}^{n}|v_{i}|\left(\frac{k}{k-i}\right)^{s}=\frac{C_{u}}{k^{s}}\sum_{i=1}^{n}|v_{i}|\left(1+\frac{i}{k-i}\right)^{s}\leq\frac{C_{u}}{k^{s}}\sum_{i=1}^{n}|v_{i}|\left(1+\frac{i}{2n+1-i}\right)^{s}.$

The last inequality follows from fact that the sequence
$\\{(1+\frac{i}{k-i})\\}_{k\geq 2n+1}$ is decreasing with respect to $k.$
Similarly, we have

$\sum_{i=k-n}^{k-1}|u_{k-i}v_{i}|=\sum_{i=1}^{n}|u_{i}v_{k-i}|\leq\frac{C_{v}}{k^{s}}\sum_{i=1}^{n}|u_{i}|\left(1+\frac{i}{2n+1-i}\right)^{s}.$

The last infinite sum is estimated as follows

$\displaystyle\sum_{i=n+1}^{k-1-n}|u_{i}v_{k-i}|$ $\displaystyle\leq
C_{u}C_{v}\sum_{i=n+1}^{k-1-n}\frac{1}{i^{s}(k-i)^{s}}=\frac{C_{u}C_{v}}{k^{s}}\sum_{i=n+1}^{k-1-n}\left(\frac{1}{i}+\frac{1}{k-i}\right)^{s}$
$\displaystyle\leq\frac{2^{s-1}C_{u}C_{v}}{k^{s}}\left(\sum_{i=n+1}^{k-1-n}\frac{1}{i^{s}}+\sum_{i=n+1}^{k-1-n}\frac{1}{(k-i)^{s}}\right)\leq\frac{2^{s-1}C_{u}C_{v}}{k^{s}}\frac{2n^{-s+1}}{s-1}.$

Combining the estimates of all six sums yields directly the assertion of the
lemma. ∎

The next lemma is analogous to Lemma 5.1 and gives us the estimate on the
result of multiplication of $u$ and $v$ which are represented in cosine and
sine series, respectively. The first $n+1$ coefficients of $u$ and the first
$n$ coefficients of $v$ are given explicitly and the rest of them is expressed
by the polynomial decay.

######  Lemma 5.2.

Assume that

$u(x)=u_{0}+\sum_{i=1}^{\infty}u_{i}\cos(ix),\quad
v(x)=\sum_{i=1}^{\infty}v_{i}\sin(ix).$ (5.10)

Moreover, assume that for some $n\in\mathbb{N}$ and $s>1$ the following bounds
hold

$u_{i}\in\frac{C_{u}[-1,1]}{i^{s}}\quad\textrm{and}\quad
v_{i}\in\frac{C_{v}[-1,1]}{i^{s}}\quad\text{for $i>n,$}$ (5.11)

where $C_{v}>0,C_{u}>0.$ Then

$(uv)(x)=\sum_{k=1}^{\infty}(uv)_{k}\sin(kx),$ (5.12)

with

$(uv)_{k}=u_{0}v_{k}+\frac{1}{2}\sum_{i=1}^{\infty}u_{i}v_{i+k}-\frac{1}{2}\sum_{i=1}^{\infty}u_{i+k}v_{i}+\frac{1}{2}\sum_{i=1}^{k-1}u_{i}v_{k-i},$
(5.13)

and the following estimates on the coefficients of the product hold for $1\leq
k\leq 2n$

$(uv)_{k}\in
u_{0}v_{k}+\frac{1}{2}\left(\sum_{i=1}^{n}v_{i+k}u_{i}-\sum_{i=1}^{n}v_{i}u_{k+i}+\sum_{i=1}^{k-1}v_{i}u_{k-i}\right)+C_{u}C_{v}\frac{n^{1-2s}}{2s-1}[-1,1],$
(5.14)

while for $k>2n$ we have

$(uv)_{k}\in\frac{D[-1,1]}{k^{s}},\quad
D=|u_{0}|C_{v}+\frac{1}{2}\left(\sum_{i=1}^{n}(C_{u}|v_{i}|+C_{v}|u_{i}|)\left(1+\left(\frac{2n+1}{2n+1-i}\right)^{s}\right)+C_{u}C_{v}(2+2^{s})\frac{n^{1-s}}{s-1}\right).$
(5.15)

###### Proof.

The argument is analogous to the proof of Lemma 5.1. Namely, we have the
following representation of the product

$\displaystyle uv$
$\displaystyle=\left(u_{0}+\sum_{i_{1}=1}^{\infty}u_{i_{1}}\cos(i_{1}x)\right)\left(\sum_{i_{2}=1}^{\infty}v_{i_{2}}\sin(i_{2}x)\right)$
$\displaystyle=\sum_{k=1}^{\infty}u_{0}v_{k}\sin(kx)+\frac{1}{2}\sum_{i_{1},i_{2}=1}^{\infty}u_{i_{1}}v_{i_{2}}\sin(i_{2}-i_{1})+\frac{1}{2}\sum_{i_{1},i_{2}=1}^{\infty}u_{i_{1}}v_{i_{2}}\sin(i_{1}+i_{2})$
$\displaystyle=\sum_{k=1}^{\infty}\left(u_{0}v_{k}+\frac{1}{2}\sum_{i_{2}-i_{1}=k}u_{i_{1}}v_{i_{2}}-\frac{1}{2}\sum_{i_{2}-i_{1}=-k}u_{i_{1}}v_{i_{2}}+\frac{1}{2}\sum_{i_{2}+i_{1}=k}u_{i_{1}}v_{i_{2}}\right)\sin(kx)$

So, for natural $k\geq 1$ we represent all coefficients using the formulas

$(uv)_{k}=u_{0}v_{k}+\frac{1}{2}\sum_{i=1}^{\infty}u_{i}v_{i+k}-\frac{1}{2}\sum_{i=1}^{\infty}u_{i+k}v_{i}+\frac{1}{2}\sum_{i=1}^{k-1}u_{i}v_{k-i},$
(5.16)

For $1\leq k\leq 2n$ we write

$(uv)_{k}=u_{0}v_{k}+\frac{1}{2}\sum_{i=1}^{n}u_{i}v_{i+k}-\frac{1}{2}\sum_{i=1}^{n}u_{i+k}v_{i}+\frac{1}{2}\sum_{i=1}^{k-1}u_{i}v_{k-i}+\frac{1}{2}\sum_{i=n+1}^{\infty}(u_{i}v_{i+k}-u_{i+k}v_{i}).$

The infinite sums in the above formula are estimated in the same way as in
Lemma 5.1. For $k>2n$ we have

$\displaystyle(uv)_{k}$
$\displaystyle=u_{0}v_{k}+\frac{1}{2}\sum_{i=1}^{n}u_{i}v_{i+k}-\frac{1}{2}\sum_{i=1}^{n}u_{i+k}v_{i}+\frac{1}{2}\sum_{i=n+1}^{\infty}u_{i}v_{i+k}-\frac{1}{2}\sum_{i=n+1}^{\infty}u_{i+k}v_{i}$
$\displaystyle+\frac{1}{2}\sum_{i=1}^{n}u_{i}v_{k-i}+\frac{1}{2}\sum_{i=k-n}^{k-1}u_{i}v_{k-i}+\frac{1}{2}\sum_{i=n+1}^{k-n-1}u_{i}v_{k-i}.$

It is easy to see that the first component of the above sum can be estimated
in the following way $|u_{0}v_{k}|\leq\frac{|u_{0}|C_{v}}{k^{s}}.$ The
remaining infinite sums are estimated in the same way as in Lemma 5.1. ∎

The following simple lemmas are useful to implement operations on the infinite
interval vectors. They can be used when working both with sine and cosine
Fourier series.

######  Lemma 5.3.

Assume that the sequence $\\{u_{i}\\}_{i=1}^{\infty}$ satisfies

$u_{i}\in[u_{i}^{-},u_{i}^{+}]\ \ \text{ for }\ \ i\leq n\ \ \text{ and }\ \
u_{i}\in\frac{[C_{u}^{-},C_{u}^{+}]}{i^{s}}\ \ \text{ for }\ \ i>n.$

If $k<n$, then

$u_{i}\in[D_{u}^{-},D_{u}^{+}]\quad\text{for }i>k.$

where,

$D_{u}^{-}=\min\\{u_{k+1}^{-}(k+1)^{s},\ldots,u_{n}n^{s},C_{u}^{-}\\},\quad
D_{u}^{+}=\max\\{u_{k+1}^{+}(k+1)^{s},\ldots,u_{n}n^{s},C_{u}^{+}\\}.$

######  Lemma 5.4.

Assume that sequence $\\{u_{i}\\}_{i=1}^{\infty}$ satisfies

$u_{i}\in[u_{i}^{-},u_{i}^{+}]\ \ \text{ for }\ \ i\leq n\ \ \text{ and }\ \
u_{i}\in\frac{[C_{u}^{-},C_{u}^{+}]}{i^{s}}\ \ \text{ for }\ \ i>k.$

Then for $s_{1}<s$ there holds

$u_{i}\in[D_{u}^{-},D_{u}^{+}]\quad\text{for }i>n,$

where,

$D_{u}^{-}=\min\left\\{0,\frac{C_{u}^{-}}{(n+1)^{s-s_{1}}}\right\\},\quad
D_{u}^{+}=\max\left\\{0,\frac{C_{u}^{+}}{(n+1)^{s-s_{1}}}\right\\}.$

######  Lemma 5.5.

Assume that sequences $\\{u_{i}\\}_{i=1}^{\infty}$ and
$\\{v_{i}\\}_{i=1}^{\infty}$ satisfy

$u_{i}\in\frac{[C_{u}^{-},C_{u}^{+}]}{i^{s_{1}}},\quad
v_{i}\in\frac{[C_{v}^{-},C_{v}^{+}]}{i^{s_{2}}}\quad\text{for $i>n$},$ (5.17)

with the constants $s_{1},s_{2}.$ Then

  * •

if $s_{1}=s_{2}$ then for $i>n$ we have
$u_{i}+v_{i}\in\frac{[C_{v}^{-}+C_{v}^{-},C_{u}^{+}+C_{v}^{+}]}{i^{s}},$

  * •

if $s_{1}<s_{2}$ then for $i>n$ we have
$u_{i}+v_{i}\in\frac{[C_{u}^{-},C_{u}^{+}]+\frac{[0,1]}{(n+1)^{s_{2}-s_{1}}}[C_{v}^{-},C_{v}^{+}]}{i^{s}_{1}},$

  * •

if $s_{1}>s_{2}$ then for $i>n$ we have
$u_{i}+v_{i}\in\frac{\frac{[0,1]}{(n+1)^{s_{1}-s_{2}}}[C_{u}^{-},C_{u}^{+}]+[C_{v}^{-},C_{v}^{+}]}{i^{s_{2}}}.$

######  Lemma 5.6.

Assume that sequences $\\{u_{i}\\}_{i=1}^{\infty}$ and
$\\{v_{i}\\}_{i=1}^{\infty}$ satisfy

$u_{i}\in\frac{[C_{u}^{-},C_{u}^{+}]}{i^{s_{1}}}\quad\textrm{and}\quad
v_{i}\in\frac{[C_{v}^{-},C_{v}^{+}]}{i^{s_{2}}}\quad\text{for $i>n$},$ (5.18)

with some constants $s_{1},s_{2}.$ Then for $i>n$ we have
$u_{i}v_{i}\in\frac{[C_{u}^{-}+C_{v}^{-},C_{u}^{+}+C_{v}^{+}]}{i^{s_{1}+s_{2}}}.$

## References

  * [1] G. Arioli “Computer assisted proof of branches of stationary and periodic solutions, and Hopf bifurcations, for dissipative PDEs” In _Communications in Nonlinear Science and Numerical Simulation_ 105, 2021, pp. 106079 DOI: 10.1016/j.cnsns.2021.106079
  * [2] G. Arioli and H. Koch “A Hopf Bifurcation in the Planar Navier–Stokes Equations” In _Journal of Mathematical Fluid Mechanics_ 23, 2021, pp. 70 DOI: 10.1007/s00021-021-00592-0
  * [3] J.F.G. Auchmuty and G. Nicolis “Bifurcation analysis of nonlinear reaction-diffusion equations—I. Evolution equations and the steady state solutions” In _Bulletin of Mathematical Biology_ 37, 1975, pp. 323–365 DOI: 10.1016/S0092-8240(75)80036-X
  * [4] J.B. Berg, M. Breden, J.-P. Lessard and L. Veen “Spontaneous Periodic Orbits in the Navier–Stokes Flow” In _Journal of Nonlinear Science_ 31, 2021, pp. 41 DOI: 10.1007/s00332-021-09695-4
  * [5] K.J. Brown and F.A. Davidson “Global bifurcation in the Brusselator system” In _Nonlinear Analysis: Theory, Methods & Applications_ 24, 1995, pp. 1713–1725 DOI: https://doi.org/10.1016/0362-546X(94)00218-7
  * [6] “CAPD library” URL: http://capd.ii.uj.edu.pl
  * [7] “Code for Computer-Assisted Proofs” URL: https://github.com/Kruci-no/Brusselator-ComputerAssistedProof
  * [8] A. Dernov “Regular Dynamics and Diffusion Chaos in the Brusselator Model” In _Differential Equations_ 37, 2001, pp. 1631–1633 DOI: 10.1023/A:1017929102036
  * [9] H. Fujita “On the blowing up of solutions of the Cauchy problem for $u_{t}-\Delta u=u^{1+\alpha}$” In _J. Fac. Sci. Univ. Tokyo. Sect. 1_ 13, 1966, pp. 109–124
  * [10] G.-H. Guo, J.-H. Wu and X.-H. Ren “Hopf bifurcation in general Brusselator system with diffusion” In _Applied Mathematics and Mechanics- English Edition_ 32, 2011, pp. 1177–1186 DOI: 10.1007/s10483-011-1491-6
  * [11] A. Haraux and T. Cazenave “An Introduction to Semilinear Evolution Equations” 13, Oxford Lecture Series in Mathematics and Its Applications Clarendon Press, 1998
  * [12] P. Kalita and P. Zgliczyński “Rigorous FEM for One-Dimensional Burgers Equation” In _SIAM Journal on Applied Dynamical Systems_ 20, 2021, pp. 853–907 DOI: 10.1137/20M1338216
  * [13] T. Kapela, D. Wilczak and P. Zgliczyński “Recent advances in rigorous computation of Poincaré maps” In _Communications in Nonlinear Science and Numerical Simulation_ 110, 2022, pp. 106366 DOI: 10.1016/j.cnsns.2022.106366
  * [14] T. Kapela and P. Zgliczyński “A Lohner-type algorithm for control systems and ordinary differential inclusions” In _Discrete and Continuous Dynamical Systems - Series B_ 11, 2007, pp. 365–385 DOI: 10.3934/dcdsb.2009.11.365
  * [15] T. Kapela, M. Mrozek, D. Wilczak and P. Zgliczyński “CAPD::DynSys: A flexible C++ toolbox for rigorous numerical analysis of dynamical systems” In _Communications in Nonlinear Science and Numerical Simulation_ 101, 2021, pp. 105578 DOI: 10.1016/j.cnsns.2020.105578
  * [16] X. Li, J. Hou and Y. Shen “Slow-fast effect and generation mechanism of Brusselator based on coordinate transformation” In _Open Physics_ 14, 2016, pp. 261–268 DOI: 10.1515/phys-2016-0032
  * [17] M. Nakao, M. Plum and Y. Watanabe “Numerical Verification Methods and Computer-Assisted Proofs for Partial Differential Equations” 53, Springer Series in Computational Mathematics Springer, 2019 DOI: 10.1007/978-981-13-7669-6
  * [18] M. Nakao and Y. Watanabe “Numerical verification methods for solutions of semilinear elliptic boundary value problems” In _Nonlinear Theory and Its Applications, IEICE_ 2, 2011, pp. 2–31 DOI: 10.1587/nolta.2.2
  * [19] F. Pacella, M. Plum and D. Rütters “A computer-assisted existence proof for Emden’s equation on an unbounded L-shaped domain” In _Communications in Contemporary Mathematics_ 19, 2016, pp. 1750005 DOI: 10.1142/S0219199717500055
  * [20] B. Peña and C. Pérez-García “Stability of Turing patterns in the Brusselator model” In _Physical review. E, Statistical, nonlinear, and soft matter physics_ 64, 2001, pp. 056213 DOI: 10.1103/PhysRevE.64.056213
  * [21] E. Tadmor “The well-posedness of the Kuramoto–Sivashinsky equation” In _SIAM Journal on Mathematical Analysis_ 17, 1986, pp. 884–893 DOI: 10.1137/051706
  * [22] D. Wilczak and P. Zgliczyński “A geometric method for infinite-dimensional chaos: Symbolic dynamics for the Kuramoto–Sivashinsky PDE on the line” In _Journal of Differential Equations_ 269, 2020, pp. 8509–8548 DOI: 10.1016/j.jde.2020.06.020
  * [23] Y. You “Global dynamics of the Brusselator equations” In _Dynamics of Partial Differential Equations_ 4, 2007, pp. 167–196 DOI: 10.4310/DPDE.2007.v4.n2.a4
  * [24] P. Zgliczyński “C1-Lohner algorithm” In _Foundations of Computational Mathematics_ 2, 2002, pp. 429–465 DOI: 10.1007/s102080010025
  * [25] P. Zgliczyński “Rigorous numerics for dissipative Partial Differential Equations II. Periodic orbit for the Kuramoto–Sivashinsky PDE - a computer-assisted proof” In _Foundations of Computational Mathematics_ 4, 2004, pp. 157–185 DOI: 10.1007/s10208-002-0080-8
  * [26] P. Zgliczyński “Rigorous Numerics for Dissipative PDEs III. An effective algorithm for rigorous integration of dissipative PDEs” In _Topological Methods in Nonlinear Analysis_ 36, 2010, pp. 197–262