# A framework for rigorous computational methods using Haar wavelets for
differential equations

Guilherme Nakassima , Corresponding authorInstituto de Ciências Matemáticas e
de Computação, Universidade de São Paulo. 400 Trabalhador São-Carlense Avenue,
São Carlos, 13566-590, São Paulo, Brazil. E-mail<EMAIL_ADDRESS>Marcio
Gameiro Department of Mathematics, Rutgers, The State University of New
Jersey, 57 US Highway 1, Piscataway, 08854, NJ, USA. E-mail:
<EMAIL_ADDRESS>
###### Abstract

This work presents a framework for a-posteriori error-estimating algorithms
for differential equations which combines the radii polynomial approach with
Haar wavelets. By using Haar wavelets, we obtain recursive structures for the
matrix representations of the differential operators and quadratic
nonlinearities, which can be exploited for the radii polynomial method in
order to get error estimates in the $L^{2}$ sense. This allows the method to
be applicable when the system or solution is not continuous, which is a
limitation of other radii-polynomial-based methods. Numerical examples show
how the method is implemented in practice.

Keywords: Rigorous computation, Computer-assisted proofs, Haar wavelets,
Nonlinear dynamical systems.

2020 MSC: 34A34, 34L30, 65G20, 65H10, 65T60.

##  1 Introduction

Rigorous computation is an area under active development since the 1980s [18].
With an steady increase of computing power, numerical methods became viable
tools for analyzing differential equations and gaining insight on structures
such as invariant objects. However, standard numerical methods provide only
approximations; the results are non-rigorous and cannot be used in formal
proofs. They can only be used to gain insights on the true structures of the
system. Moreover, some structures, such as bifurcations, may still be hidden
even when using very accurate numerical methods.

Rigorous computational methods try to fill these gaps, providing
mathematically valid estimates and bounds for truncation and rounding errors,
and rigorously proving the existence of such hidden structures. Over the
years, a number of such methods were developed, such as rigorous integration
[29, 8, 17], Conley index methods [4, 14], self-consistent bounds [30, 28] and
discretization methods [19, 9]. A more thorough review can be seen in [20] and
references therein.

Of particular interest to us are the _radii polynomials_ approach [5, 7, 27,
13, 11, 6, 2, 22, 26]. These methods recast the problem of investigating the
existence of structures as finding solutions to functional equations, usually
expanding the solutions in terms of a basis. Then, usual numerical methods are
employed to find an approximate solution to these equations. Finally, using
fixed point theorems, we can guarantee the existence of a true, rigorous
solution of the functional equations within certain bounds of the numerical
solution. The hypotheses of the fixed point theorems are in turn proven to be
satisfied with the aid of the so-called radii polynomials.

In this work we present a new radii polynomial method employing Haar wavelets.
While many other bases were already employed by this approach, such as Taylor
series [22], Fourier series [27] and Chebyshev polynomials [13], to the best
of our knowledge, no attempt has been made to combine wavelet methods and
rigorous computations. We also believe that this method can be can be a
framework to build other methods upon, such as rigorous continuation methods
and methods for partial differential equations.

Wavelets are functions that form an orthonormal basis for the $L^{2}$ function
space. While wavelet theory was only relatively recently formalized, their
special properties — such as time and frequency localization — made them
widely applicable in many fields, such as signal processing and compression
algorithms. This poses an interesting case, because most of the aforementioned
radii polynomial methods were proposed to work with smooth functions, while
our wavelet-based radii polynomial method works in more general settings.

The radii polynomial method presented in this paper is based on the ideas of a
numerical method for solving differential equations, the _Haar wavelet method_
, which was first proposed in [3]. It is assumed that the highest-order
derivative is expressed in terms of wavelets, and the solutions are given by
the integral of the series, essentially rewriting the differential equation in
its integral form. This allows the method to work with only Haar wavelets, the
simplest wavelet available, and leads to a matrix representation of the
integral operator. The simplicity of the Haar wavelet allows this matrix to be
easily and recursively calculated. The original Haar wavelet method was
further analyzed, developed and applied in several publications [12, 21, 1,
24, 16, 15].

In our method, we find the functional equations for the radii polynomial
method using the same expansions of the Haar wavelet method. This essentially
transforms the differential equation into an integral one, and allow us to use
a radii polynomial theorem similar to [25]. Also, by using Haar wavelets, the
integral operator and nonlinearities can be represented using infinite but
recursive matrices, allowing us to make the estimates needed in the radii
polynomial approach.

The work is organized in the following way. In Section 2 we introduce the Haar
wavelet and its integral, and review some of their properties. In Section 3 we
introduce the radii polynomial method and prove the theorems that guarantee
the existence of a true solution, provided that certain estimates are
satisfied. In Section 4 we study quadratic nonlinearities in order to prove
estimates needed in the method. In Section 5 we present some examples
illustrating the applications of the proposed method.

##  2 The Haar wavelets and their integral

Here we introduce the Haar wavelet and its integral, which are one of the
pillars of the proposed methods. Consider the space $L^{2}([0,1])$ of the
square-integrable functions $f:[0,1]\to\mathbb{R}$ with respect to the usual
Lebesgue measure. The Haar wavelets are a family of functions
$\\{\phi,\psi_{j,k}\\}\subset L^{2}([0,1])$ defined, for $j=0,1,\ldots$ and
$k=0,1,\ldots,2^{j}-1$, by

$\phi(t):=\begin{cases}1&,\quad 0\leq t<1\\\
0&,\quad\text{otherwise}\end{cases}\qquad\text{,}\qquad\psi_{j,k}(t):=\begin{cases}{2^{j/2}}&,\quad\frac{k}{2^{j}}\leq
t<\frac{k+0.5}{2^{j}}\\\\[3.0pt] -{2^{j/2}}&,\quad\frac{k+0.5}{2^{j}}\leq
t<\frac{k+1}{2^{j}}\\\\[3.0pt] 0&,\quad\text{otherwise}\end{cases}$ (1)

The Haar wavelets form an orthonormal basis for $L^{2}([0,1])$; the proof is
in many standard texts in wavelet theory, see e.g. [10]. Hence any function
$y(t)\in L^{2}([0,1])$ can be expanded into a unique Haar wavelet series

$y(t)=c_{1}\phi(t)+\sum_{j=0}^{\infty}\sum_{k=0}^{2^{j}-1}c_{j,k}\psi_{j,k}(t),$

where $c_{1}=\int_{0}^{1}y(t)\phi(t)\,dt$,
$c_{j,k}=\int_{0}^{1}y(t)\psi_{j,k}(t)\,dt$, and the sum converges in
$L^{2}([0,1])$. If we make $i=2^{j}+k+1$, then the sequence
$(c_{i})_{i=1}^{\infty}\in\ell^{2}(\mathbb{R})$.

######  Notation.

We can change between the “one-index” and “two-indices” notations, depending
on which is more convenient in each case. One can be converted to the other by
making, for all $i>2$,

$i=2^{j}+k+1\iff\begin{aligned} j&=\lfloor\text{log}_{2}i\rfloor\\\
k&=i-2^{j}-1\end{aligned}$

where $\lfloor\cdot\rfloor$ is the floor function. The index $i=1$ is reserved
for the scaling function $\phi$.

Conversely, any $(c_{i})_{i=1}^{\infty}\in\ell^{2}(\mathbb{R})$ defines a
unique $f(t)\in L^{2}([0,1])$ by making

$f(t)=\sum_{i=1}^{\infty}c_{i}\psi_{i}(t)=c_{1}\phi(t)+\sum_{j=0}^{\infty}\sum_{k=0}^{2^{j}-1}c_{j,k}\psi_{j,k}(t).$

This means that we can define an invertible operator
$\mathcal{H}:L^{2}([0,1])\to\ell^{2}(\mathbb{R})$ defined element-wise as

$(\mathcal{H}f)_{i}:=\int_{0}^{1}f(t)\psi_{i}(t)\,dt.$

The inverse Haar transform $\mathcal{H}^{-1}$ is given by

$(\mathcal{H}^{-1}\text{{c}})(t):=\text{{h}}^{T}(t)\text{{c}}=\sum_{i=1}^{\infty}c_{i}\psi_{i}(t)\quad,\quad\text{{h}}(t):=(\psi_{1}(t),\psi_{2}(t),...)$

It is worth noting that the Haar transform is an isometry, due to the fact
that the Haar wavelet system is an orthonormal basis of $L^{2}([0,1])$.

Another interesting and useful property of the Haar wavelets is what we call
the “nesting property”:

######  Proposition 2.1 (Nesting property).

Let $\psi_{j,k}$ and $\psi_{m,n}$ be two Haar wavelets such that $j<m$. If
$\text{supp }\psi_{j,k}\cap\text{supp }\psi_{m,n}\not=\emptyset$, then either
$\text{supp }\psi_{m,n}\subseteq[\frac{k}{2^{j}},\frac{k+0.5}{2^{j}})$ or
$\text{supp }\psi_{m,n}\subseteq[\frac{k+0.5}{2^{j}},\frac{k+1}{2^{j}})$.

The proof is simple and will be omitted; it consists in comparing the supports
of the wavelets $\phi_{j,k}$, which are dyadic intervals of length $2^{-j}$,
for different $j$. Intuitively, it means that, if the supports of two wavelets
at different resolutions overlap, then the support of the “finer” wavelet
(that is, the higher-resolution one) is entirely nested within either the
positive or the negative part of the “coarser” wavelet.

For this work, we are also interested in the integral of the Haar wavelets,
and how the integral relates to the wavelet themselves. The integrals of the
Haar wavelet family in the interval $[0,1]$ are the triangular functions given
by

$w_{1}(t)=t\qquad\text{,}\qquad
w_{j,k}(t)=\begin{cases}{2^{j/2}}\left(t-\frac{k}{2^{j}}\right)&,\quad\frac{k}{2^{j}}\leq
t\leq\frac{k+0.5}{2^{j}}\\\\[3.0pt]
{2^{j/2}}\left(\frac{k+1}{2^{j}}-t\right)&,\quad\frac{k+0.5}{2^{j}}\leq
t\leq\frac{k+1}{2^{j}}\\\\[3.0pt] 0&,\quad\text{otherwise}.\end{cases}$ (2)

We also extend the one-index notation to the Haar wavelet integrals $w_{i}$.

The Haar wavelet integrals are continuous functions in $[0,1]$; thus they are
square-integrable in that interval, and can be expanded in Haar wavelet series
themselves:

$w_{i}(t)=\sum_{l=1}^{\infty}P_{i,l}\psi_{l}(t)=P\text{{h}}^{T}(t)\ ,\
P_{i,l}:=\int_{0}^{1}\psi_{l}(t)w_{i}(t)\,dt.$ (3)

While $P$ is expressed by an infinite matrix, there is a recursive formula to
compute it. We must first define the Haar matrix of our wavelet system:

######  Definition 2.1.

For a given resolution $J$, the Haar matrix $H_{M}$ of order $M=2^{J+1}$ is
given element-wise by

$(H_{M})_{p,q}:=\psi_{p}(t_{q})$ (4)

where $\displaystyle t_{q}=\frac{q-0.5}{M}$, for $q=1,\dots,M$.

The discrete Haar transform matrix $HT_{M}$ of order $M=2^{J+1}$ is defined by

$HT_{M}:=\frac{1}{\sqrt{M}}H_{M}$ (5)

A particularly important fact is that $HT_{M}$ is unitary for all $J$, and
hence $H_{M}$ is invertible and $(H_{M})^{-1}=\frac{1}{M}H_{M}^{T}$.

######  Theorem 2.1.

The infinite matrix $P$ can be recursively calculated as

$P_{1}=\frac{1}{2}\qquad,\qquad
P_{2m}=\begin{bmatrix}P_{m}&-\displaystyle\frac{1}{4\sqrt{m^{3}}}H_{m}\\\\[7.22743pt]
\displaystyle\frac{1}{4\sqrt{m^{3}}}H_{m}^{T}&0_{m\times m}\end{bmatrix}$ (6)

for $m=2^{j}$ and $j=0,1,2,\dots$.

The proof for this formula is in [3], with some modifications to account for
the fact that we are using normalized wavelets.

Let us prove that $P^{T}\text{{c}}\in\ell^{2}(\mathbb{R})$ for all
$\text{{c}}\in\ell^{2}(\mathbb{R})$. We first need to define some projections.
Given a resolution level $J$, let $M=2^{J+1}$ as before, and define the
projection $\Pi_{M}$ as

$\displaystyle\Pi_{M}:\ell^{2}(\mathbb{R})$ $\displaystyle\to\mathbb{R}^{M}$
(7) c $\displaystyle\mapsto(c_{1},...,c_{M}).$

We identify the vector $(c_{1},...,c_{M})\in\mathbb{R}^{M}$ with its infinite-
dimensional counterpart $(c_{1},...,c_{M},0,0,...)\in\ell^{2}(\mathbb{R})$. We
also define the projection $\Pi_{\infty}\in B(\ell^{2}(\mathbb{R}))$ as

$\Pi_{\infty}\text{{c}}:=(I-\Pi_{M})\text{{c}}.$

######  Notation.

Given $\text{{c}}\in\ell^{2}(\mathbb{R})$, we sometimes divide it into blocks
of length $2^{n}$, $n=0,1,2,...$ as

$\text{{c}}=\left(\text{{c}}^{*}_{0},\text{{c}}^{*}_{1},\text{{c}}^{*}_{2},\text{{c}}^{*}_{4},...\right)^{T}\
,\text{ where }\text{{c}}^{*}_{0}=c_{1}\text{ ,
}\text{{c}}^{*}_{2^{n}}=\left(c_{2^{n-1}+1},c_{2^{n-1}+2},...,c_{2^{n}}\right)^{T}.$
(8)

Also, given a matrix $A$, we denote $A_{m_{1}:m_{2},n_{1}:n_{2}}$ the
submatrix of $A$ given by

$A_{m_{1}:m_{2},n_{1}:n_{2}}=\begin{bmatrix}A_{m_{1},n_{1}}&\cdots&A_{m_{1},n_{2}}\\\
\vdots&&\vdots\\\ A_{m_{2},n_{1}}&\cdots&A_{m_{2},n_{2}}\end{bmatrix}$

If we wish to take all rows or all columns of $A$, we denote
$A_{*,n_{1}:n_{2}}$ and $A_{m_{1}:m_{2},*}$, respectively. Lastly, to reduce
notation clutter, we denote $\text{{c}}_{M}=\Pi_{M}\text{{c}}$ and
$\text{{c}}_{\infty}=\Pi_{\infty}\text{{c}}$.

######  Proposition 2.2.

For $\text{{{c}}}\in\ell^{2}(\mathbb{R})$, $M=2^{J+1}$, $J=0,1,2,...$ and the
projections as defined before,

$\displaystyle\left\|\Pi_{M}P^{T}\text{{{c}}}\right\|_{\ell^{2}}\leq\frac{1}{\sqrt{3}}\left(4-\frac{1}{2^{2J+2}}\right)^{\frac{1}{2}}\|\text{{{c}}}\|_{\ell^{2}}\quad\quad,\quad\quad\left\|\Pi_{\infty}P^{T}\text{{{c}}}\right\|_{\ell^{2}}\leq\frac{1}{\sqrt{3}}\frac{\|\text{{{c}}}\|_{\ell^{2}}^{2}}{2^{J+1}}$

In particular, letting $J\to\infty$,
$P^{T}\text{{{c}}}\in\ell^{2}(\mathbb{R})$ and
$P^{T}:\ell^{2}(\mathbb{R})\to\ell^{2}(\mathbb{R})$ is a bounded linear
operator.

###### Proof.

From the structure of $P$ in Theorem 2.1 and the block representation of c, we
have

$\left(P^{T}\text{{c}}\right)_{2^{j}}^{*}=-\frac{1}{2^{\frac{3j}{2}+2}}H_{2^{j}}^{T}(\Pi_{2^{j}}\text{{c}})+\sum_{q=j+1}^{\infty}\frac{1}{2^{\frac{3q}{2}+2}}H_{2^{q}}\text{{c}}^{*}_{2^{q}}$
(9)

Recalling that $H_{M}=\sqrt{M}HT_{M}$ and $HT_{M}$ is a unitary matrix, we can
bound the above term by

$\displaystyle\left\|\left(P^{T}\text{{c}}\right)_{2^{j}}^{*}\right\|_{\ell^{2}}$
$\displaystyle\leq\frac{1}{2^{j+2}}\left\|HT_{2^{j}}^{T}(\Pi_{2^{j}}\text{{c}})\right\|_{\ell^{2}}+\sum_{q=j+1}^{\infty}\frac{1}{2^{q+2}}\left\|HT_{2^{q}}\text{{c}}^{*}_{2^{q}}\right\|_{\ell^{2}}$
$\displaystyle\leq\frac{\|\text{{c}}\|_{\ell^{2}}}{2^{j+2}}+\sum_{q=j+1}^{\infty}\frac{1}{2^{q+2}}\left\|\text{{c}}^{*}_{2^{q}}\right\|_{\ell^{2}}\leq\frac{\|\text{{c}}\|_{\ell^{2}}}{2^{j+1}}$
(10)

The only term left is $\left(P^{T}\bar{\text{{c}}}\right)_{1}$. We can bound
it with

$\displaystyle\left|\left(P^{T}\bar{\text{{c}}}\right)_{1}\right|$
$\displaystyle=\left|\frac{1}{2}c_{1}-\sum_{q=0}^{\infty}\frac{1}{2^{\frac{3q}{2}+2}}H_{2^{q}}\text{{c}}^{*}_{2^{q}}\right|\leq\frac{1}{2}\left|c_{1}\right|+\sum_{q=0}^{\infty}\frac{1}{2^{q+2}}\left\|HT_{2^{q}}\text{{c}}^{*}_{2^{q}}\right\|_{\ell^{2}}\leq\left\|\text{{c}}\right\|_{\ell^{2}}$
(11)

Thus, the norm of $\Pi_{M}P^{T}\text{{{c}}}$ can be estimated by

$\displaystyle\left\|\Pi_{M}P^{T}\text{{{c}}}\right\|_{\ell^{2}}^{2}$
$\displaystyle\leq\left|\left(P^{T}\text{{c}}\right)_{1}\right|^{2}+\sum_{j=0}^{J}\left\|\left(P^{T}\text{{c}}\right)_{2^{j}}^{*}\right\|_{\ell^{2}}^{2}=\left(4-\frac{1}{2^{2J+2}}\right)\frac{\|\text{{c}}\|_{\ell^{2}}^{2}}{3}$

Analogously, the norm of $\Pi_{\infty}P^{T}\text{{{c}}}$ is bounded by

$\displaystyle\left\|\Pi_{\infty}P^{T}\text{{{c}}}\right\|_{\ell^{2}}^{2}$
$\displaystyle\leq\sum_{j=J+1}^{\infty}\left\|\left(P^{T}\bar{\text{{c}}}\right)_{2^{j}}^{*}\right\|_{\ell^{2}}^{2}\leq\sum_{j=J+1}^{\infty}\frac{\|\text{{c}}\|_{\ell^{2}}^{2}}{2^{2j+2}}=\frac{1}{3}\frac{\|\text{{c}}\|_{\ell^{2}}^{2}}{2^{2J+2}}$

∎

##  3 The radii polynomial approach

In this section, we introduce the radii polynomial [5, 7, 27, 13, 11, 6]
approach for rigorous computation. Consider an initial value problem

$\begin{cases}\dot{u}=f(u,t)\\\ u(0)=u_{0}\end{cases}$ (12)

and suppose we find a numerical, approximate solution $\bar{u}(t)$. Our aim is
to prove the existence of a true solution $\tilde{u}(t)$ in some neighborhood
of $\bar{u}(t)$. This is done using the radii polynomial method.

For our work, suppose that $\dot{u}(t)\in L^{2}([0,1])$. Then, we can write
$\dot{u}$ and $u$ using the Haar wavelet and its integral as

$\dot{u}(t)=\sum_{i=1}^{\infty}c_{i}\psi_{i}(t)\quad,\quad
u(t)=\sum_{i=1}^{\infty}c_{i}w_{i}(t)+u_{0}$ (13)

Substituting back into the differential equation (12) and taking the Haar
transform,

$F(\text{{c}}):=\text{{c}}-\mathcal{H}\left(f\left(t,\sum_{i=1}^{\infty}c_{i}w_{i}(t)+u_{0}\right)\right)=0$
(14)

Thus, we have a map $F:\ell^{2}(\mathbb{R})\to\ell^{2}(\mathbb{R})$ such that
finding a solution of the differential equation implies finding a zero of $F$.
Conversely, due to the uniqueness of the wavelet series that represents
$\dot{u}$ — and consequently of the series that represents $u(t)$ — finding a
zero of $F$ is equivalent to finding a solution to (12).

Now we recast the problem of finding the zeros of $F$ to finding a fixed point
of a map $T:\ell^{2}(\mathbb{R})\to\ell^{2}(\mathbb{R})$ near the numerical
solution $\bar{u}$. This is done by showing that $T$ is a contraction near
$\bar{u}$. First, define the operator
$A:\ell^{2}(\mathbb{R})\to\ell^{2}(\mathbb{R})$ as

$A\text{{x}}=A_{M}\Pi_{M}\text{{x}}+\Pi_{\infty}\text{{x}},$

where $A_{M}:\mathbb{R}^{M}\to\mathbb{R}^{M}$ is a finite-dimensional,
computational approximation for the inverse
$D(\Pi_{M}F(\bar{\text{{{x}}}}))^{-1}$, with $\bar{\text{{{x}}}}$ the solution
numerically obtained. Then we can define the fixed-point map
$T:\ell^{2}(\mathbb{R})\to\ell^{2}(\mathbb{R})$ by

$T(\text{{x}}):=\text{{x}}-AF(\text{{x}})=(\Pi_{M}-A_{M}\Pi_{M}F)(\text{{x}})+\Pi_{\infty}(\text{{x}}-F(\text{{x}}))$
(15)

Its derivative, which is used for the radii polynomial method, is given by

$DT(\text{{x}})=\Pi_{M}-A_{M}\Pi_{M}DF(\text{{x}})+\Pi_{\infty}(I-DF(\text{{x}}))$
(16)

######  Remark 3.1.

Since many of the matrices in this work are block matrices, one can compute
their inverse as

$\begin{bmatrix}A&B\\\
C&D\end{bmatrix}^{-1}=\begin{bmatrix}A^{-1}+A^{-1}B(D-CA^{-1}B)^{-1}CA^{-1}&-A^{-1}B(D-CA^{-1}B)^{-1}\\\
-(D-CA^{-1}B)^{-1}CA^{-1}&(D-CA^{-1}B)^{-1}\end{bmatrix}$

provided that the matrices $A$ and $D-CA^{-1}B$ are invertible. Calculating
the inverse with this formula can be faster than directly inverting the full
matrix $D(\Pi_{M}F(\bar{\text{{{x}}}}))$.

######  Notation.

In order to help visualize the operators and reduce clutter in notation, we
employ the following “block matrix” notation for an operator
$C:\ell^{2}(\mathbb{R})\to\ell^{2}(\mathbb{R})$ whenever it is convenient:

$C\text{{x}}=\begin{bmatrix}C_{M}&C_{M,\infty}\\\
C_{\infty,M}&C_{\infty}\end{bmatrix}\begin{pmatrix}\Pi_{M}\text{{x}}\\\
\Pi_{\infty}\text{{x}}\end{pmatrix}.$ (17)

We refer to this as the finite-infinite decomposition.

Intuitively, the first term of (15) is a Newton-like map for the finite terms
we computed numerically. In the second term, we hope the “tail” of $F$ will
contract to zero by itself — which is what happens with quadratic
nonlinearities. This is all motivated by the fact that, as $M$ increases, the
new elements of the matrix $P_{M}$ become smaller.

We now formally prove that a fixed point of $T$ corresponds to a zero of $F$:

######  Proposition 3.1.

Suppose the map $T$ as defined above in (15) is a contraction in some closed
neighborhood of $\ell^{2}(\mathbb{R})$. Then $T$ has a unique fixed point
$\tilde{\text{{{c}}}}$ in that neighborhood. Moreover, $\tilde{\text{{{c}}}}$
is a fixed point of $T$ if and only if it is a zero of $F$ as defined earlier
in (14).

###### Proof.

Since $T$ is a contraction in a closed neighborhood of $\ell^{2}(\mathbb{R})$,
the Banach Fixed Point Theorem guarantees that it has a fixed point
$\tilde{\text{{c}}}$ in the same neighborhood. Also, if $\tilde{\text{{c}}}$
is a zero of $F$, then a straightforward calculation shows that it is a fixed
point of $T$.

It remains to prove that the fixed point $\tilde{\text{{c}}}$ is a zero of
$F$. By (15)

$T(\tilde{\text{{c}}})-\tilde{\text{{c}}}=0=A_{M}\Pi_{M}F(\tilde{\text{{c}}})+\Pi_{\infty}F(\tilde{\text{{c}}}).$

Since $A_{M}\Pi_{M}F(\tilde{\text{{c}}})\in\Pi_{M}(\ell^{2}(\mathbb{R}))$, we
have

$\displaystyle\Pi_{\infty}F(\tilde{\text{{c}}})=0\quad,\quad
A_{M}\Pi_{M}F(\tilde{\text{{c}}})=0$

and since $A_{M}$ is invertible, then $\Pi_{M}F(\tilde{\text{{c}}})=0$ as
well. Thus $F(\tilde{\text{{c}}})=0$. ∎

To prove that $T$ is actually a contraction near of our numerical solution
$\bar{\text{{c}}}$, we use the radii polynomials. First, we define the closed
neighborhood

$\overline{B_{\omega}(\bar{\text{{c}}},r)}=\left\\{\text{{y}}\in\ell^{2}(\mathbb{R}):\|\Pi_{M}(\text{{y}}-\bar{\text{{c}}})\|\leq\omega
r\text{ and
}\|\Pi_{\infty}(\text{{y}}-\bar{\text{{c}}})\|\leq(1-\omega)r\right\\}$ (18)

in which $T$ will be a contraction. $\omega\in(0,1)$ is a “trade-off
parameter”: we can loosen the radius in the infinite part, at the cost of
tightening the radius in the finite part, and vice-versa. For the next
calculations, we assume $\omega$ is fixed, though in practice it is chosen
later.

Next, we need bounds $Y_{M}$ and $Y_{\infty}$, and polynomials $Z_{M}(r)$ and
$Z_{\infty}(r)$ such that

$\displaystyle\|\Pi_{M}(T(\bar{\text{{c}}})-\bar{\text{{c}}})\|_{\ell^{2}}\leq
Y_{M}$ (19)
$\displaystyle\|\Pi_{\infty}(T(\bar{\text{{c}}})-\bar{\text{{c}}})\|_{\ell^{2}}\leq
Y_{\infty}$ (20)
$\displaystyle\sup_{x_{1},x_{2}\in\overline{B(r)}}\|\Pi_{M}(DT(\bar{\text{{c}}}+x_{1})x_{2})\|\leq
Z_{M}(r)r$ (21)
$\displaystyle\sup_{x_{1},x_{2}\in\overline{B(r)}}\|\Pi_{\infty}(DT(\bar{\text{{c}}}+x_{1})x_{2})\|\leq
Z_{\infty}(r)r$ (22)

Then, we can define the radii polynomials as

$\displaystyle p_{M}(r)$ $\displaystyle:=Z_{M}(r)r-\omega r+Y_{M}\quad,\quad
p_{\infty}(r)$ $\displaystyle:=Z_{\infty}(r)r-(1-\omega)r+Y_{\infty}$ (23)

######  Theorem 3.1.

Consider the radii polynomials as defined in (23). If there exists an
$r_{0}>0$ such that $p_{M}(r_{0})<0$ and $p_{\infty}(r_{0})<0$, then there
exists a unique $\tilde{\text{{{c}}}}\in B_{\omega}(\bar{\text{{{c}}}},r_{0})$
such that $T(\tilde{\text{{{c}}}})=\tilde{\text{{{c}}}}$.

###### Proof.

Due to the Banach Fixed Point Theorem, we only need to prove that
$T\left(\overline{B_{\omega}(\bar{\text{{c}}},r_{0})}\right)\subseteq\overline{B_{\omega}(\bar{\text{{c}}},r_{0})}$
and that $T$ is a contraction when restricted to
$\overline{B_{\omega}(\bar{\text{{c}}},r_{0})}$.

We first prove that $T$ restricted to
$\overline{B_{\omega}(\bar{\text{{c}}},r_{0})}$ is a contraction. If
$y_{1},y_{2}\in\overline{B_{\omega}(\bar{\text{{c}}},r_{0})}$, then by the
Mean Value Theorem

$\displaystyle\|T(y_{1})-T(y_{2})\|_{\ell^{2}}$
$\displaystyle\leq\sup_{x\in\overline{B_{\omega}(\bar{\text{{c}}},r_{0})}}\|DT(x)\|\|y_{1}-y_{2}\|_{\ell^{2}}=\sup_{x_{1}\in\overline{B(r_{0})}}\|DT(\bar{\text{{c}}}+x_{1})\|\|y_{1}-y_{2}\|_{\ell^{2}}$

Hence, we must show that $\|DT(\bar{\text{{c}}}+x_{1})\|<1$ for
$x_{1}\in\overline{B(r_{0})}$. Observe that

$\displaystyle\sup_{x_{1}\in\overline{B(r_{0})}}\|DT(\bar{\text{{c}}}+x_{1})\|$
$\displaystyle=\frac{1}{r_{0}}\sup_{x_{1},x_{2}\in\overline{B(r_{0})}}\|DT(\bar{\text{{c}}}+x_{1})x_{2}\|_{\ell^{2}}$
$\displaystyle\leq\frac{1}{r_{0}}\bigg{(}\sup_{x_{1},x_{2}\in\overline{B(r_{0})}}\|\Pi_{M}DT(\bar{\text{{c}}}+x_{1})x_{2}\|_{\ell^{2}}+\sup_{x_{1},x_{2}\in\overline{B(r_{0})}}\|\Pi_{\infty}DT(\bar{\text{{c}}}+x_{1})x_{2}\|_{\ell^{2}}\bigg{)}$
$\displaystyle\leq Z_{M}(r_{0})+Z_{\infty}(r_{0})$

So we have

$\|T(y_{1})-T(y_{2})\|_{\ell^{2}}\leq(Z_{M}(r_{0})+Z_{\infty}(r_{0}))\|y_{1}-y_{2}\|_{\ell^{2}}$

But since $p_{M}(r_{0})<0$ and $p_{\infty}(r_{0})<0$,

$\displaystyle Z_{M}(r_{0})r_{0}+Z_{\infty}(r_{0})r_{0}-r_{0}$
$\displaystyle\leq(Z_{M}(r_{0})r_{0}-\omega
r_{0}+Y_{M})+(Z_{\infty}(r_{0})r_{0}-(1-\omega)r_{0}+Y_{\infty})$
$\displaystyle\hskip 20.00003pt=p_{M}(r_{0})+p_{\infty}(r_{0})<0$

Thus $Z_{M}(r_{0})+Z_{\infty}(r_{0})<1$, and $T$ restricted to
$\overline{B_{\omega}(\bar{\text{{c}}},r_{0})}$ is a contraction.

Now we must prove that
$T\left(\overline{B_{\omega}(\bar{\text{{c}}},r_{0})}\right)\subseteq\overline{B_{\omega}(\bar{\text{{c}}},r_{0})}$.
If $y\in\overline{B_{\omega}(\bar{\text{{c}}},r_{0})}$, then

$\displaystyle\|\Pi_{M}(T(y)-\bar{\text{{c}}})\|_{\ell^{2}}$
$\displaystyle\leq\|\Pi_{M}(T(y)-T(\bar{\text{{c}}}))\|_{\ell^{2}}+\|\Pi_{M}(T(\bar{\text{{c}}})-\bar{\text{{c}}})\|_{\ell^{2}}$
$\displaystyle\leq Z_{M}(r_{0})\|y-\bar{\text{{c}}}\|_{\ell^{2}}+Y_{M}\leq
Z_{M}(r_{0})r_{0}+Y_{M}<\omega r_{0}$

and similarly for $\Pi_{\infty}(T(y)-\bar{\text{{c}}})$

$\displaystyle\|\Pi_{\infty}(T(y)-\bar{\text{{c}}})\|_{\ell^{2}}$
$\displaystyle\leq\|\Pi_{\infty}(T(y)-T(\bar{\text{{c}}}))\|_{\ell^{2}}+\|\Pi_{\infty}(T(\bar{\text{{c}}})-\bar{\text{{c}}})\|_{\ell^{2}}$
$\displaystyle\leq
Z_{\infty}(r_{0})\|y-\bar{\text{{c}}}\|_{\ell^{2}}+Y_{\infty}\leq
Z_{\infty}(r_{0})r_{0}+Y_{\infty}<(1-\omega)r_{0}$

and hence $T(y)\in\overline{B_{\omega}(\bar{\text{{c}}},r_{0})}$. ∎

Thus, if the radii polynomial method is successful in finding an $r_{0}$, then
the solution $\bar{\text{{c}}}$ found by the numerical method is “close” to
the wavelet coefficients of true solution $\tilde{\text{{c}}}$ in the
$\ell^{2}(\mathbb{R})$ sense, that is,
$\|\bar{\text{{c}}}-\tilde{\text{{c}}}\|_{\ell^{2}}\leq r_{0}$. Or
equivalently, the numerical approximation
$\bar{u}(t)=\bar{\text{{c}}}^{T}\text{{h}}(t)$ is “close” to the true solution
$\tilde{u}(t)$ in the $L^{2}([0,1])$ sense, that is,
$\|\bar{u}-\tilde{u}\|_{L^{2}}\leq r_{0}$.

##  4 Nonlinear terms

In this section we study quadratic nonlinearities in more depth. This may seem
restrictive, but there are many interesting systems involving those, such as
the Lorenz system. Furthermore, we believe that estimates for higher
nonlinearities can be computed with similar techniques.

Consider two functions $u$ and $v$ such that $\dot{u},\dot{v}\in
L^{2}([0,1])$. Their expansions into Haar wavelet integrals as in (13) are

$\displaystyle
u(t)=u_{0}+\sum_{i=1}^{\infty}c_{i}w_{i}(t)=u_{0}+\text{{w}}^{T}(t)\text{{c}}\
,\
v(t)=v_{0}+\sum_{i=1}^{\infty}d_{i}w_{i}(t)=v_{0}+\text{{w}}^{T}(t)\text{{d}}$

Thus, considering the product $u(t)v(t)$, we have

$\displaystyle
u(t)v(t)=u_{0}v_{0}+u_{0}\text{{w}}^{T}(t)\text{{d}}+v_{0}\text{{w}}^{T}(t)\text{{c}}+(W(t))(\text{{c}},\text{{d}})$

where, for $t\in[0,1]$,
$(W(t))(\text{{c}},\text{{d}}):=\text{{c}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{d}}$.
A crucial observation is that, for any given $t\in[0,1]$, $W(t)$ _is a
symmetric bilinear form_.

The next theorem shows that $W(t)(\text{{c}},\text{{d}})\in L^{2}([0,1])$ for
any pair $\text{{c}},\text{{d}}\in\ell^{2}(\mathbb{R})$, which allows us to
calculate its Haar transform and use the radii polynomial methods developed in
Section 3. However, its proof is lengthy and will be left to A for clarity.

######  Theorem 4.1.

The bilinear form
$\left(W(\text{{{c}}},\text{{{d}}})\right)(t):=\text{{{c}}}^{T}\text{{{w}}}(t)\text{{{w}(t)}}\text{{{d}}}$
is bounded in $L^{2}([0,1])$ for all
$\text{{{c}}},\text{{{d}}}\in\ell^{2}(\mathbb{R})$, that is, there exists
$C>0$ such that

$\left\|W(\text{{{c}}},\text{{{d}}})\right\|_{L^{2}}\leq
C\left\|\text{{{c}}}\right\|_{\ell^{2}}\left\|\text{{{d}}}\right\|_{\ell^{2}}\text{
for all }\text{{{c}}},\text{{{d}}}\in\ell^{2}(\mathbb{R}).$

Additionally, making $\text{{{a}}}:=P^{T}\text{{{c}}}$ and
$\text{{{b}}}:=P^{T}\text{{{d}}}$, we can write

$\text{{{c}}}^{T}\text{{{w}}}(t)\text{{{w}(t)}}\text{{{d}}}=\text{{{a}}}^{T}\Omega(t)\text{{{b}}}+\text{{{a}}}^{T}\Omega^{T}(t)\text{{{b}}}+\text{{{a}}}^{T}\Theta(t)\text{{{b}}}$
(24)

where the operator $\Omega(t)$ can be recursively defined as

$\Omega_{1}(t)=0\quad,\quad\Omega_{2m}(t)=\begin{bmatrix}\Omega_{m}(t)&\Upsilon_{m}(t)\\\\[7.22743pt]
0_{m\times m}&0_{m\times m}\end{bmatrix}$ (25)

with $\Upsilon_{m}(t)$ being a $m\times m$ matrix defined element-wise as

$(\Upsilon_{m})_{i,l}(t)=(H_{m})_{i,l}\psi_{m+l}(t),$

and the operator $\Theta(t)$ can be represented as an infinite diagonal matrix
given by

$\Theta(t)=\left[\begin{array}[]{cccc}\psi_{1}^{2}(t)&\phantom{0}&\phantom{0}&\phantom{0}\\\
\phantom{0}&\psi_{2}^{2}(t)&\phantom{0}&\phantom{0}\\\
\phantom{0}&\phantom{0}&\psi_{3}^{2}(t)&\phantom{0}\\\
\phantom{0}&\phantom{0}&\phantom{0}&\ddots\\\ \end{array}\right]$ (26)

with zeros omitted for clarity. With the terms defined as above,

$\mathcal{H}(\text{{{c}}}^{T}\text{{{w}}}(t)\text{{{w}(t)}}\text{{{d}}})=\mathcal{H}(\text{{{a}}}^{T}\Omega(t)\text{{{b}}})+\mathcal{H}(\text{{{a}}}^{T}\Omega(t)\text{{{b}}})+\mathcal{H}(\text{{{a}}}^{T}\Omega(t)\text{{{b}}})$

and

$\displaystyle\mathcal{H}(\text{{{a}}}^{T}\Omega(t)\text{{{b}}})$
$\displaystyle=(\tilde{\Omega}^{T}\text{{{a}}})\odot\text{{{b}}}$ (27)
$\displaystyle\mathcal{H}(\text{{{a}}}^{T}\Omega^{T}(t)\text{{{b}}})$
$\displaystyle=(\tilde{\Omega}^{T}\text{{{b}}})\odot\text{{{a}}}$
$\displaystyle\mathcal{H}(\text{{{a}}}^{T}\Theta(t)\text{{{b}}})$
$\displaystyle=\Gamma^{T}(\text{{{a}}}\odot\text{{{b}}})$

where $\tilde{\Omega}^{T}$ and $\tilde{\Gamma}^{T}$ are recursively defined as

$\begin{aligned}
\tilde{\Omega}_{1}=0\quad,\quad\tilde{\Omega}_{2m}=\begin{bmatrix}\tilde{\Omega}_{m}&H_{m}\\\
0_{m}&0_{m}\end{bmatrix}\\\
\Gamma_{1}=1\quad,\quad\Gamma_{2m}=\begin{bmatrix}\Gamma_{m}&0_{m}\\\
H_{m}^{T}&0_{m}\end{bmatrix}\end{aligned}\quad,\text{ for }m=2^{j}\text{ and
}j=0,1,2,3,...$

We now present some estimates required for the radii polynomial method. The
full proof for those estimates are lengthy and left B. The main strategy
consists in employing both the recursive block structures of the matrices from
Theorems 2.1 and 4.1 and the finite-infinite decomposition from (17). These
estimates provide tighter bounds which increase the likelihood of finding an
$r_{0}$ which satisfies Theorem 3.1. We believe that similar estimates may be
applied for higher-degree polynomial nonlinearities.

######  Proposition 4.1.

Given $\bar{\text{{{c}}}},\bar{\text{{{d}}}}\in\mathbb{R}^{M}$, $M=2^{J+1}$
for some $J\geq 0$, and $\bar{\text{{{a}}}}=P^{T}\bar{\text{{{c}}}}$,
$\bar{\text{{{b}}}}=P^{T}\bar{\text{{{d}}}}$, the following estimates are
valid:

  1. i)

$\Pi_{M}P^{T}\bar{\text{{{c}}}}=P_{M}^{T}\bar{\text{{{c}}}}$

  2. ii)

$\Pi_{M}\mathcal{H}(\bar{\text{{{a}}}}^{T}\Omega(t)\bar{\text{{{b}}}})=(\tilde{\Omega}_{M}^{T}P_{M}^{T}\bar{\text{{{c}}}})\odot(P_{M}^{T}\bar{\text{{{d}}}})$

  3. iii)

$\Pi_{M}\mathcal{H}(\bar{\text{{{a}}}}^{T}\Theta(t)\bar{\text{{{b}}}})=\Gamma^{T}_{M}(\bar{\text{{{a}}}}_{M}\odot\bar{\text{{{b}}}}_{M})+\Gamma^{T}_{\infty,M}(\bar{\text{{{a}}}}_{\infty}\odot\bar{\text{{{b}}}}_{\infty})$;
moreover,

$\left\|\Gamma^{T}_{\infty}(\bar{\text{{{a}}}}_{\infty}\odot\bar{\text{{{b}}}}_{\infty})\right\|_{\ell^{2}}\leq\frac{\sqrt{2}\,\|\bar{\text{{{c}}}}\|_{\ell^{2}}\|\bar{\text{{{d}}}}\|_{\ell^{2}}}{\left(4-\sqrt{2}\right)2^{\frac{3J}{2}+3}}.$

######  Proposition 4.2.

Given $\bar{\text{{{c}}}},\bar{\text{{{d}}}}\in\mathbb{R}^{M}$, $M=2^{J+1}$
for some $J\geq 0$, and $\bar{\text{{{a}}}}=P^{T}\bar{\text{{{c}}}}$,
$\bar{\text{{{b}}}}=P^{T}\bar{\text{{{d}}}}$, the following estimates are
valid:

  1. i)

$\displaystyle\left\|\Pi_{\infty}P^{T}\bar{\text{{{c}}}}\right\|_{\ell^{2}}\leq\frac{1}{\sqrt{3}}\frac{\left\|\bar{\text{{{c}}}}\right\|_{\ell^{2}}}{2^{J+2}}$

  2. ii)

$\displaystyle\left\|\Pi_{\infty}\mathcal{H}(\bar{\text{{{a}}}}^{T}\Omega(t)\bar{\text{{{b}}}})\right\|_{\ell^{2}}\leq\frac{1}{\sqrt{3}}\frac{\left\|\bar{\text{{{c}}}}\right\|_{\ell^{2}}\left\|\bar{\text{{{d}}}}\right\|_{\ell^{2}}}{2^{2J+4}}$

  3. iii)

$\displaystyle\left\|\Pi_{\infty}\mathcal{H}(\bar{\text{{{a}}}}^{T}\Theta(t)\bar{\text{{{b}}}})\right\|_{\ell^{2}}\leq\frac{1}{21\sqrt{7}}\frac{\left\|\bar{\text{{{c}}}}\right\|_{\ell^{2}}\left\|\bar{\text{{{d}}}}\right\|_{\ell^{2}}}{2^{3J+6}}.$

######  Proposition 4.3.

Given $\text{{{x}}},\text{{{y}}}\in\ell^{2}(\mathbb{R})$ and
$\bar{\text{{{c}}}}\in\mathbb{R}^{M}$,

  1. i)

$\displaystyle\left\|(P_{\infty,M}^{T}\,\text{{{y}}}_{\infty})\right\|_{\ell^{2}}\leq\frac{\left\|\text{{{y}}}\right\|_{\ell^{2}}}{2^{J+2}}$

  2. ii)

$\displaystyle\left\|\left(\tilde{\Omega}^{T}P^{T}\right)_{M,\infty}\text{{{y}}}_{\infty}\right\|_{\ell^{2}}\leq\left(1+\sqrt{2}\right)\frac{\left\|\text{{{y}}}\right\|_{\ell^{2}}}{2^{\frac{J+3}{2}}}$

  3. iii)

$\displaystyle\left\|\Gamma_{\infty,M}^{T}(\bar{\text{{{a}}}}_{\infty}\odot\Pi_{\infty}P^{T}\text{{{y}}})\right\|_{\ell^{2}}\leq\frac{\sqrt{2}}{4-\sqrt{2}}\frac{\left\|\bar{\text{{{c}}}}\right\|_{\ell^{2}}\left\|\text{{{y}}}\right\|_{\ell^{2}}}{2^{\frac{3J}{2}+3}}$

  4. iv)

$\displaystyle\left\|A_{M}\Pi_{M}\mathcal{H}(\text{{{x}}}^{T}P\,\Omega(t)P^{T}\text{{{y}}})\right\|_{\ell^{2}}\leq
K_{1}\|\text{{{x}}}\|_{\ell^{2}}\|\text{{{y}}}\|_{\ell^{2}}$

  5. v)

$\displaystyle\left\|A_{M}\Pi_{M}\mathcal{H}(\text{{{x}}}^{T}P\Theta(t)P^{T}\text{{{y}}})\right\|_{\ell^{2}}\leq
K_{2}\|\text{{{x}}}\|_{\ell^{2}}\|\text{{{y}}}\|_{\ell^{2}}$

where

$\displaystyle K_{1}$
$\displaystyle=\left\|A_{M}\,\text{{diag}}\left(\|(P_{M}^{T})_{i,*}\|_{\ell^{2}}\right)\,\tilde{\Omega}_{M}^{T}P_{M}^{T}\right\|+\frac{\|A_{M}\tilde{\Omega}_{M}^{T}P_{M}^{T}\|}{2^{J+2}}+\frac{(1+\sqrt{2})\|A_{M}P_{M}^{T}\|}{2^{\frac{J+3}{2}}}+\frac{(1+\sqrt{2})\|A_{M}\|}{2^{\frac{3J+7}{2}}}$
$\displaystyle K_{2}$
$\displaystyle=\left\|A_{M}\Gamma_{M}^{T}\text{{diag}}(\|(P_{M}^{T})_{i,*}\|_{\ell}^{2})P_{M}^{T}\right\|+\frac{\left\|A_{M}\Gamma_{M}^{T}P_{M}^{T}\right\|}{2^{J+1}}+\frac{\left\|A_{M}\Gamma_{M}^{T}\right\|}{2^{2J+4}}+\frac{\sqrt{2}\|A_{M}\|}{(4-\sqrt{2})2^{\frac{3J}{2}+4}}$

######  Proposition 4.4.

Given $\text{{{x}}},\text{{{y}}}\in\ell^{2}(\mathbb{R})$ and
$\bar{\text{{{c}}}}\in\mathbb{R}^{M}$, the following estimates are valid:

  1. i)

$\displaystyle\left\|\Pi_{\infty}\mathcal{H}(\bar{\text{{{c}}}}^{T}\text{{{w}}}(t)\text{{{w}}}^{T}(t)\text{{{y}}})\right\|_{\ell^{2}}\leq
D_{1}\frac{\left\|\bar{\text{{{c}}}}\right\|_{\ell^{2}}\left\|\text{{{y}}}\right\|_{\ell^{2}}}{2^{\frac{3J}{2}}}\quad,\quad
D_{1}=\frac{1}{8\sqrt{7}}\left(3+\frac{\sqrt{2}}{4}+\frac{4}{4-\sqrt{2}}\right);$

  2. ii)

$\displaystyle\left\|\Pi_{\infty}\mathcal{H}(\text{{{x}}}^{T}\text{{{w}}}(t)\text{{{w}}}^{T}(t)\text{{{y}}})\right\|_{\ell^{2}}\leq
D_{2}\frac{\left\|\text{{{x}}}\right\|_{\ell^{2}}\left\|\text{{{y}}}\right\|_{\ell^{2}}}{2^{\frac{3J}{2}+2}}\quad,\quad
D_{2}=\frac{8+6\sqrt{2}}{\sqrt{7}\left(4-\sqrt{2}\right)}$.

##  5 Examples

In this section we illustrate the implementation of our method by means of
three examples: the logistic equation, the logistic equation with a
discontinuous forcing term, and the Lorenz system. The method was implemented
in MATLAB R2021b using the INTLAB package for interval arithmetic [23]. The
files for these examples are available at
https://github.com/gknakassima/RigComp-HaarWavelet.

###  5.1 Logistic equation

As a first simple example we look at the logistic equation

$\dot{u}=\lambda u(1-u)\quad,\quad u(0)=u_{0}$ (28)

since it has a polynomial nonlinearity and its analytical solution is given by
$u(t)=\frac{u_{0}e^{\lambda t}}{1-u_{0}+u_{0}e^{\lambda t}}$. Using the
expansions in (13), we have the functional equation equivalent to (28):

$\text{{c}}^{T}\text{{h}}(t)-\lambda\left(\text{{c}}^{T}\text{{w}}(t)+u_{0}-u_{0}^{2}-2u_{0}\text{{c}}^{T}\text{{w}}(t)-(W(t))(\text{{c}},\text{{c}})\right)=0.$
(29)

####  5.1.1 Obtaining a numerical approximation

In order to obtain a finite-dimensional approximation $\bar{\text{{c}}}$ of
the solution, we first consider a truncated version of our matrix equation.
Given a resolution level $J>0$ and making $M=2^{J+1}$, we apply the projection
$\Pi_{M}$ to all sequences of (29), obtaining

$\text{{c}}_{M}^{T}\text{{h}}_{M}(t)-\lambda\left(\text{{c}}_{M}^{T}\text{{w}}_{M}(t)+u_{0}-u_{0}^{2}-2u_{0}\text{{c}}_{M}^{T}\text{{w}}_{M}(t)-(W_{M}(t))(\text{{c}}_{M},\text{{c}}_{M})\right)=0$

where, for $\text{{a}}_{M},\text{{b}}_{M}\in\mathbb{R}^{M}$,
$(W_{M}(t))(\text{{a}}_{M},\text{{b}}_{M}):=\text{{a}}_{M}^{T}\text{{w}}_{M}(t)\text{{w}}_{M}^{T}(t)\text{{b}}_{M}$.
Since the equation holds for all $t\in[0,1]$, we sample it at the times
$t_{l}=\frac{l-0.5}{2^{J}}$ for $l=1,...,M$. Recalling that

$H_{M}=[\text{{h}}_{M}(t_{1}),\text{{h}}_{M}(t_{2}),...,\text{{h}}_{M}(t_{l}),...,\text{{h}}_{M}(t_{2^{J+1}})],$

we can organize the time samples in matrix form as

$\text{{c}}_{M}^{T}H_{M}+\lambda(2u_{0}-1)\text{{c}}_{M}^{T}P_{M}H_{M}+\lambda(u_{0}^{2}-u_{0})\text{{e}}^{T}+\lambda(\text{{W}}_{M}(\text{{c}}_{M},\text{{c}}_{M}))^{T}=0$

with e and $\text{{W}}_{M}(\text{{c}}_{M},\text{{c}}_{M})$ being $M\times 1$
vectors given by

$\text{{e}}:=(1,1,...,1)^{T}\quad,\quad\big{(}\text{{W}}_{M}(\text{{c}}_{M},\text{{c}}_{M})\big{)}_{i}:=(W_{M}(t_{i}))(\text{{c}}_{M},\text{{c}}_{M}).$

By transposing this system and multiplying by
$(H_{M}^{T})^{-1}=\frac{1}{M}H_{M}$ we finally arrive at the equation to solve
numerically:

$\text{{c}}_{M}+\lambda(2u_{0}-1)P_{M}^{T}\text{{c}}_{M}+\frac{\lambda(u_{0}^{2}-u_{0})}{M}H_{M}\text{{e}}+\frac{\lambda}{M}H_{M}\text{{W}}_{M}(\text{{c}}_{M},\text{{c}}_{M})=0.$
(30)

As this is a nonlinear equation, we use Newton’s method. Define
$F_{M}:\mathbb{R}^{M}\to\mathbb{R}^{M}$ as

$F_{M}(\text{{c}}_{M}):=\text{{c}}_{M}+\lambda(2u_{0}-1)P_{M}^{T}\text{{c}}_{M}+\frac{\lambda(u_{0}^{2}-u_{0})}{M}H_{M}\text{{e}}+\frac{\lambda}{M}H_{M}\text{{W}}_{M}(\text{{c}}_{M},\text{{c}}_{M})$
(31)

Thus, we apply Newton’s method by iteratively calculating

$\text{{c}}^{p+1}_{M}=\text{{c}}^{p}_{M}-A_{M}(\text{{c}}^{p}_{M})F_{M}(\text{{c}}^{p}_{M})$
(32)

where $A_{M}(\text{{c}}^{p}_{M})$ is a numerical approximation for
$\Big{(}DF_{M}(\text{{c}}^{p}_{M})\Big{)}^{-1}$ and $\text{{c}}^{p}_{M}$ is
the result from the $p$-th iteration.

####  5.1.2 Estimates for the radii polynomials

Here, we provide the bounds for Theorem 3.1. The maps $T$ and $DT$ for the
fixed point theorem are given by (15) and (16), respectively. For this
example, the map $F:\ell^{2}(\mathbb{R})\to\ell^{2}(\mathbb{R})$ is given as

$F(\text{{x}}):=\lambda(u_{0}^{2}-u_{0})e_{1}+\text{{x}}+\lambda(2u_{0}-1)P^{T}\text{{x}}+\lambda\mathcal{H}\left(\text{{x}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{x}}\right).$
(33)

and its derivative $DF(\text{{x}})$ applied to
$\text{{y}}\in\ell^{2}(\mathbb{R})$ is given by

$(DF(\text{{x}}))\text{{y}}=\text{{y}}+\lambda(2u_{0}-1)P^{T}\text{{y}}+2\lambda\mathcal{H}(\text{{x}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{y}}).$
(34)

where the last equality comes from the symmetry of the bilinear form.

Before proceeding, it is worth outlining the general strategy for the
estimates. We separate the operator matrices according to the finite-infinite
decomposition in (17). Then, all the finite-dimensional parts are collected
together and left for the computer to calculate, while we use the analytic
estimates from Section 4 for the infinite parts.

  * •

$Y_{M}$: Using the decomposition of
$\lambda\mathcal{H}\left(\bar{\text{{c}}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\bar{\text{{c}}}\right)$
and Proposition 4.1,

$\displaystyle\Pi_{M}(T(\bar{\text{{c}}})-\bar{\text{{c}}})$
$\displaystyle=-A_{M}\Pi_{M}\mathcal{H}F(\bar{\text{{c}}})$
$\displaystyle=-A_{M}\Big{[}\lambda(u_{0}^{2}-u_{0})\text{{e}}+\bar{\text{{c}}}+\lambda(2u_{0}-1)P_{M}^{T}\bar{\text{{c}}}+2\lambda\big{(}(\tilde{\Omega}_{M}^{T}P_{M}^{T}\bar{\text{{c}}})\odot(P_{M}^{T}\bar{\text{{c}}})\big{)}$
$\displaystyle\hskip
40.00006pt+\lambda\Gamma^{T}_{M}(\bar{\text{{{a}}}}_{M}\odot\bar{\text{{{a}}}}_{M})\Big{]}-\lambda
A_{M}\Pi_{M}\Gamma^{T}_{\infty,M}(\bar{\text{{{a}}}}_{\infty}\odot\bar{\text{{{a}}}}_{\infty}),$

with $\bar{\text{{{a}}}}=P^{T}\bar{\text{{{c}}}}$. Note that the term in the
brackets can be computationally evaluated. Using the bound from Proposition
4.1 (iii) for the last term, we can define $Y_{M}$ as

$\displaystyle Y_{M}$
$\displaystyle:=\Big{\|}A_{M}\Big{[}\lambda(u_{0}^{2}-u_{0})\text{{e}}+\bar{\text{{c}}}+\lambda(2u_{0}-1)P_{M}^{T}\bar{\text{{c}}}+2\lambda\big{(}(\tilde{\Omega}_{M}^{T}P_{M}^{T}\bar{\text{{c}}})\odot(P_{M}^{T}\bar{\text{{c}}})\big{)}+\lambda\Gamma^{T}_{M}(\bar{\text{{{a}}}}_{M}\odot\bar{\text{{{a}}}}_{M})\Big{]}\Big{\|}_{\ell^{2}}$
(35) $\displaystyle\hskip
20.00003pt+\frac{|\lambda|\|A_{M}\|\sqrt{2}}{\left(4-\sqrt{2}\right)2^{\frac{3J}{2}+3}}\|\bar{\text{{c}}}\|_{\ell^{2}}^{2}$

  * •

$Y_{\infty}$: Observe that

$\Pi_{\infty}(T(\bar{\text{{c}}})-\bar{\text{{c}}})=\Pi_{\infty}\left(\lambda(2u_{0}-1)P^{T}\bar{\text{{c}}}+\lambda\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)\bar{\text{{a}}}+\bar{\text{{a}}}^{T}\Omega^{T}(t)\bar{\text{{a}}}+\bar{\text{{a}}}^{T}\Theta(t)\bar{\text{{a}}}\right).$

With the estimates from Proposition 4.2, we can make $Y_{\infty}$ as

$Y_{\infty}:=\frac{\left|\lambda(2u_{0}-1)\right|}{\sqrt{3}}\frac{\left\|\bar{\text{{c}}}\right\|_{\ell^{2}}}{2^{J+2}}+\frac{|\lambda|}{\sqrt{3}}\frac{\left\|\bar{\text{{c}}}\right\|_{\ell^{2}}^{2}}{2^{2J+3}}+\frac{|\lambda|}{21\sqrt{7}}\frac{\left\|\bar{\text{{c}}}\right\|_{\ell^{2}}^{2}}{2^{3J+6}}.$
(36)

  * •

$Z_{M}$: Using the finite-infinite decomposition and the fact that
$\bar{\text{{c}}}\in\mathbb{R}^{M}$,

$\displaystyle\Pi_{M}(DT(\bar{\text{{c}}}+\text{{x}}))\text{{y}}$
$\displaystyle=\left(I_{M}-A_{M}B_{1}\right)\text{{y}}_{M}-A_{M}B_{2}P_{\infty,M}^{T}\text{{y}}_{\infty}-2\lambda
A_{M}\text{diag}(P_{M}^{T}\bar{\text{{c}}})\big{(}\tilde{\Omega}^{T}P^{T}\big{)}_{M,\infty}\text{{y}}_{\infty}$
$\displaystyle\hskip 30.00005pt-2\lambda
A_{M}\Gamma_{\infty,M}^{T}\Pi_{\infty}\left(P^{T}\bar{\text{{c}}}\odot
P^{T}\text{{y}}\right)-2\lambda
A_{M}\Pi_{M}\mathcal{H}(\text{{x}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{y}})$

where

$\displaystyle B_{1}$
$\displaystyle:=I_{M}+\lambda(2u_{0}-1)P_{M}^{T}+2\lambda\text{diag}\left(\bar{\text{{c}}}^{T}P_{M}\tilde{\Omega}_{M}\right)P_{M}^{T}$
$\displaystyle\hskip
50.00008pt+2\lambda\text{diag}(P_{M}^{T}\bar{\text{{c}}})\tilde{\Omega}_{M}^{T}P_{M}^{T}+2\lambda\Gamma_{M}^{T}\text{diag}\left(P_{M}^{T}\bar{\text{{c}}}\right)P_{M}^{T}$
$\displaystyle B_{2}$
$\displaystyle:=\lambda(2u_{0}-1)I_{M}+2\lambda\text{diag}\left(\bar{\text{{c}}}^{T}P_{M}\tilde{\Omega}_{M}\right)+2\lambda\text{diag}(P_{M}^{T}\bar{\text{{c}}})+2\lambda\Gamma_{M}^{T}\text{diag}\left(P_{M}^{T}\bar{\text{{c}}}\right)$

and $I_{M}$ is the $M\times M$ identity matrix. While the expression of
$B_{1}$ and $B_{2}$ seem complicated, all terms are finite-dimensional and
hence their norms can be calculated computationally. Thus, using the estimates
from Proposition 4.3 with the expression for $DT$ to bound the terms which
cannot be easily estimated computationally,

$\left\|\Pi_{M}(DT(\bar{\text{{c}}}+\text{{x}}))\text{{y}}\right\|_{\ell^{2}}\leq\left(C_{1}+C_{2}\left\|\text{{x}}\right\|_{\ell^{2}}\right)\left\|\text{{y}}\right\|_{\ell^{2}}$

where

$\displaystyle C_{1}$
$\displaystyle:=\left\|I_{M}-A_{M}B_{1}\right\|+\frac{\left\|A_{M}B_{2}\right\|}{2^{J+2}}+\frac{(1+\sqrt{2})\left\|\lambda
A_{M}\text{diag}(P_{M}^{T}\bar{\text{{c}}})\right\|}{2^{\frac{J+1}{2}}}+\frac{\sqrt{2}\,\|\lambda
A_{M}\|\,\|\bar{\text{{c}}}\|_{\ell^{2}}}{(4-\sqrt{2})\,2^{\frac{3J}{2}+2}}$
$\displaystyle C_{2}$ $\displaystyle:=2|\lambda|\left(2K_{1}+K_{2}\right)$

and $K_{1}$ and $K_{2}$ are as in Proposition 4.3. Hence, we can make
$Z_{M}(r)$ as

$Z_{M}(r):=C_{1}+C_{2}r$ (37)

  * •

$Z_{\infty}$: We have that

$\Pi_{\infty}(DT(\bar{\text{{c}}}+\text{{x}}))\text{{y}}=\lambda(2u_{0}-1)\Pi_{\infty}P^{T}\text{{y}}+2\lambda\Pi_{\infty}\mathcal{H}(\bar{\text{{c}}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{y}})+2\lambda\Pi_{\infty}\mathcal{H}(\text{{x}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{y}})$

All terms are infinite-dimensional and need to be analitically estimated.
Using the estimates from Propositions 4.4 and 2.2, we can make $Z_{\infty}(r)$
as

$Z_{\infty}(r):=\left(\frac{\left|\lambda(2u_{0}-1)\right|}{2^{J+1}\sqrt{3}}+\frac{|\lambda|D_{1}\left\|\bar{\text{{c}}}\right\|_{\ell^{2}}}{2^{\frac{3J}{2}}}\right)+\frac{|\lambda|D_{2}}{2^{\frac{3J}{2}+2}}\,r$
(38)

with $D_{1}$ and $D_{2}$ as in Proposition 4.4.

####  5.1.3 Results

Figure 1 shows the numerical solutions for $J=6$ and $J=10$ compared to the
true solution, using $\lambda=6$ and $u_{0}=0.2$ for both cases. Visually, the
numerical solutions agrees with the true one.

(a) $J=6$

(b) $J=10$

Figure 1: Numerical and true solutions for the logistic equation.

Figure 2 shows the radius $r_{0}$ obtained as $J$ increases for different
$\omega$. It is clear that, as $J$ increases, the radius $r_{0}$ decreases;
this is due to more terms being calculated more accurately, instead of only
being bounded by analytical estimates. Also, smaller values of $\omega$ yield
tighter radii; however, if $\omega$ is too small the method will not work, as
there will not be a true solution within
$\overline{B_{\omega}(\bar{\text{{c}}},r)}$. Figure 3(a) shows the radii
obtained with $\omega$ optimized up to two significant digits.

$J$ |  $r_{0}$ ($\omega=0.6$) |  $r_{0}$ ($\omega=0.75$) |  $r_{0}$ ($\omega=0.85$)
---|---|---|---
6 | $2.1677704\times 10^{-2}$ | $3.4976922\times 10^{-2}$ | $5.9222878\times 10^{-2}$
7 | $1.0690405\times 10^{-2}$ | $1.7163569\times 10^{-2}$ | $2.8785441\times 10^{-2}$
8 | $5.3120948\times 10^{-3}$ | $8.5120420\times 10^{-3}$ | $1.4224769\times 10^{-2}$
9 | $2.6483940\times 10^{-3}$ | $4.2402808\times 10^{-3}$ | $7.0756199\times 10^{-3}$
10 | $1.3223867\times 10^{-3}$ | $2.1164777\times 10^{-3}$ | $3.5294192\times 10^{-3}$

Figure 2: Radius $r_{0}$ obtained for different values of $\omega$

Lastly, Figure 3(b) shows the radii and computation time. As $J$ increases,
the computation time is expected to increase; however, Figure 3(b) shows that
after a certain point it increases more rapidly than $r_{0}$ decreases. This
is expected as the size of the matrices quadruples for every increase of $J$,
and so one must carefully balance the needed precision with computing time.

(a)

(b)

$J$ | $\omega$ | $r_{0}$ | Time
---|---|---|---
6 | 0.51 | $1.7651807\times 10^{-2}$ | $0.283$
7 | 0.31 | $6.1826325\times 10^{-3}$ | $0.377$
8 | 0.21 | $2.6863885\times 10^{-3}$ | $1.499$
9 | 0.14 | $1.4138969\times 10^{-3}$ | $9.447$
10 | 0.089 | $1.2789917\times 10^{-3}$ | $70.958$
11 | 0.060 | $7.4420035\times 10^{-4}$ | $474.220$

Figure 3: (a) Radius $r_{0}$ obtained with more optimized $\omega$ for each
resolution level $J$. (b) Comparison between $r_{0}$ and time elapsed.

###  5.2 Logistic equation with a discontinuous forcing term

The next example is again the logistic equation, but with a discontinuous
forcing term

$\begin{aligned} &\dot{u}=\lambda u(1-u)+g(t)\\\
&u(0)=u_{0}\end{aligned}\quad,\quad g(t)=\begin{cases}1&,\text{if
}t\leq\frac{1}{2}\\\ 0&,\text{if }t>\frac{1}{2}\\\ \end{cases}$ (39)

While this is a Ricatti equation which can be explicitly solved, we can see
from the equation itself that the solution should not be smooth, since $g$ is
discontinuous. Nonetheless, we can find a verification radius in the $L^{2}$
sense.

The functional equation for this case is similar to (28); making
$\text{{g}}=\mathcal{H}(g)$,

$\text{{c}}^{T}\text{{h}}(t)-\lambda\left(\text{{c}}^{T}\text{{w}}(t)+u_{0}-u_{0}^{2}-2u_{0}\text{{c}}^{T}\text{{w}}(t)-(W(t))(\text{{c}},\text{{c}})\right)-\text{{g}}^{T}\text{{h}}(t)=0.$
(40)

####  5.2.1 Numerical approximation

Applying the same method as in the previous example, we get the equation to
obtain the numerical, finite-dimensional approximation $\bar{\text{{c}}}$,
which is similar to before:

$\text{{c}}_{M}+\lambda(2u_{0}-1)P_{M}^{T}\text{{c}}_{M}+\frac{\lambda(u_{0}^{2}-u_{0})}{M}H_{M}\text{{e}}+\frac{\lambda}{M}H_{M}\text{{W}}_{M}(\text{{c}}_{M},\text{{c}}_{M})-\frac{1}{M}H_{M}\text{{g}}_{M}=0$
(41)

where $\text{{g}}_{M}:=(g(t_{1})\,,\,...\,,\,g(t_{M}))^{T}$. Again, since this
is a nonlinear equation, we will use Newton’s method.

####  5.2.2 Estimates for the radii polynomials

In order to apply the radii polynomial method, we use the following functional
equation $F:\ell^{2}(\mathbb{R})\to\ell^{2}(\mathbb{R})$:

$F(\text{{x}}):=\lambda(u_{0}^{2}-u_{0})e_{1}+\text{{x}}+\lambda(2u_{0}-1)P^{T}\text{{x}}+\lambda\mathcal{H}\left(\text{{x}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{x}}\right)-\mathcal{H}(g).$

Observe that $g=\frac{1}{2}(\psi_{1}+\psi_{2})$, and thus
$\mathcal{H}(g)=\left(\frac{1}{2},\frac{1}{2},0,0,...\right)^{T}=:\text{{g}}$;
in particular, $\text{{g}}\in\mathbb{R}^{M}$. Hence, the map used is given by

$F(\text{{x}})=\lambda(u_{0}^{2}-u_{0})e_{1}+\text{{x}}+\lambda(2u_{0}-1)P^{T}\text{{x}}+\lambda\mathcal{H}\left(\text{{x}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{x}}\right)-\text{{g}}.$
(42)

The maps $T$ and $DT(\text{{x}})$ for the radii polynomials are the same as
(15) and (16), respectively. Actually, since g does not depend on x, the
derivative $DF$ is the same as in the non-forced logistic equation from (34).

Using the same methods as before, we have the following bounds for the radii
polynomial method:

$\displaystyle Y_{M}$
$\displaystyle:=C_{0}+\frac{|\lambda|\|A_{M}\|\sqrt{2}}{\left(4-\sqrt{2}\right)2^{\frac{3J}{2}+3}}\|\bar{\text{{c}}}\|_{\ell^{2}}^{2}$
(43) $\displaystyle Y_{\infty}$
$\displaystyle:=\frac{\left|\lambda(2u_{0}-1)\right|}{\sqrt{3}}\frac{\left\|\bar{\text{{c}}}\right\|_{\ell^{2}}}{2^{J+2}}+\frac{|\lambda|}{\sqrt{3}}\frac{\left\|\bar{\text{{c}}}\right\|_{\ell^{2}}^{2}}{2^{2J+3}}+\frac{|\lambda|}{21\sqrt{7}}\frac{\left\|\bar{\text{{c}}}\right\|_{\ell^{2}}^{2}}{2^{3J+6}}$
(44) $\displaystyle Z_{M}(r)$ $\displaystyle:=C_{1}+C_{2}r$ (45)
$\displaystyle Z_{\infty}(r)$
$\displaystyle:=\left(\frac{\left|\lambda(2u_{0}-1)\right|}{2^{J+1}\sqrt{3}}+\frac{|\lambda|D_{1}\left\|\bar{\text{{c}}}\right\|_{\ell^{2}}}{2^{\frac{3J}{2}}}\right)+\frac{|\lambda|D_{2}}{2^{\frac{3J}{2}+2}}\,r$
(46)

where $C_{1}$, $C_{2}$, $D_{1}$ and $D_{2}$ are as in the estimates for the
non-forced logistic equation, and

$\displaystyle C_{0}$
$\displaystyle=\Big{\|}A_{M}\Big{[}\lambda(u_{0}^{2}-u_{0})\text{{e}}+\bar{\text{{c}}}+\lambda(2u_{0}-1)P_{M}^{T}\bar{\text{{c}}}$
$\displaystyle\hskip
50.00008pt+2\lambda\Big{(}(\tilde{\Omega}_{M}^{T}P_{M}^{T}\bar{\text{{c}}})\odot(P_{M}^{T}\bar{\text{{c}}})\Big{)}+\lambda\Gamma^{T}_{M}(\bar{\text{{{a}}}}_{M}\odot\bar{\text{{{a}}}}_{M})-\text{{g}}\Big{]}\Big{\|}_{\ell^{2}}.$

####  5.2.3 Results

Figure 4 shows the results using the Haar wavelet method compared to numerical
integration, using $\lambda=6$ and $u_{0}=0.2$. For the numerical integration,
we used the same amount of points as the Haar wavelet method, that is,
$2^{J+1}$ points. It can be seen that the numerical integration tends to
smooth the graph at $t=0.5$, while our method preserves the original shape.

(a) $J=6$

(b) $J=10$

Figure 4: Numerical and true solutions for the logistic equation with forcing
term

Figure 5 shows the verification radius as $J$ increases. It is worth noting
that the solution is not smooth; nonetheless, our method returned verification
radii similar to the non-forced logistic equation.

$J$ | $\omega$ | $r_{0}$
---|---|---
6 | 0.53 | $2.6161420\times 10^{-2}$
7 | 0.31 | $9.2508029\times 10^{-3}$
8 | 0.20 | $5.1495598\times 10^{-3}$
9 | 0.13 | $3.1382945\times 10^{-3}$
10 | 0.086 | $1.7710909\times 10^{-3}$
11 | 0.057 | $1.1107730\times 10^{-3}$

Figure 5: Verification radius for the forced logistic equation

###  5.3 Lorenz system

The Lorenz system is given by

$\displaystyle\dot{x}$ $\displaystyle=\sigma(y-x)$ (47) $\displaystyle\dot{y}$
$\displaystyle=x(\rho-z)-y$ $\displaystyle\dot{z}$ $\displaystyle=xy-\beta z$

where $\sigma$, $\rho$ and $\beta$ are positive parameters, usually taken as
$\sigma=10$, $\beta=\frac{8}{3}$ and $\rho=28$. This is a well-studied system,
and with these parameters the system exhibits chaotic behavior with a strange
attractor.

Since this is a system of equations, we must first some of our definitions in
order to apply our method. First, we define the spaces
$X_{s}:=\ell^{2}(\mathbb{R})\times\ell^{2}(\mathbb{R})\times\ell^{2}(\mathbb{R})$
and $X_{f}:=L^{2}([0,1])\times L^{2}([0,1])\times L^{2}([0,1])$, and endow
them with the norms

$\displaystyle\|(\text{{c}}_{1},\text{{c}}_{2},\text{{c}}_{3})\|_{X_{s}}$
$\displaystyle:=\max\left\\{\|\text{{c}}_{1}\|_{\ell^{2}},\|\text{{c}}_{2}\|_{\ell^{2}},\|\text{{c}}_{3}\|_{\ell^{2}}\right\\}$
(48) $\displaystyle\|(f_{1},f_{2},f_{3})\|_{X_{f}}$
$\displaystyle:=\max\left\\{\|f_{1}\|_{L^{2}},\|f_{2}\|_{L^{2}},\|f_{3}\|_{L^{2}}\right\\}$

With these norms, $X_{s}$ and $X_{f}$ are still Banach spaces; though they are
no longer Hilbert spaces, the methods of Section 3 are still applicable. Also,
an operator $A:X_{s}\to X_{s}$ can be expressed using block matrix notation as

$A=\begin{bmatrix}A_{1,1}&A_{1,2}&A_{1,3}\\\ A_{2,1}&A_{2,2}&A_{2,3}\\\
A_{3,1}&A_{3,2}&A_{3,3}\end{bmatrix}$

where $A_{i,j}:\ell^{2}(\mathbb{R})\to\ell^{2}(\mathbb{R})$ for $i,j=1,2,3$.
If those are bounded, $A$ is bounded and

$\|A\|_{B(X_{s})}=\max_{1\leq i\leq
3}\sum_{j=1}^{3}\|A_{i,j}\|_{B(\ell^{2})}.$ (49)

Similar notation will be used when $A_{i,j}\in B(\mathbb{R}^{M})$. Lastly, we
make a small abuse of notation and extend the notation for the operators in
$\ell^{2}(\mathbb{R})$ such as the projections $\Pi_{M}$ and $\Pi_{\infty}$ to
$X_{s}$ by applying them element-wise:

$\displaystyle\Pi_{M}(\text{{c}}_{1},\text{{c}}_{2},\text{{c}}_{3})$
$\displaystyle:=(\Pi_{M}\text{{c}}_{1},\Pi_{M}\text{{c}}_{2},\Pi_{M}\text{{c}}_{3})$
$\displaystyle\Pi_{\infty}(\text{{c}}_{1},\text{{c}}_{2},\text{{c}}_{3})$
$\displaystyle:=(\Pi_{\infty}\text{{c}}_{1},\Pi_{\infty}\text{{c}}_{2},\Pi_{\infty}\text{{c}}_{3}).$

Applying the same methods used to obtain (29) to each equation in (47), we
obtain the functional equation

$\displaystyle\text{{c}}_{x}^{T}\text{{h}}(t)-\sigma\text{{c}}_{y}^{T}\text{{w}}(t)+\sigma\text{{c}}_{x}^{T}\text{{w}}(t)-\sigma(y_{0}-x_{0})=0$
(50)
$\displaystyle\text{{c}}_{y}^{T}\text{{h}}(t)+x_{0}\text{{c}}_{z}^{T}\text{{w}}(t)-(\rho-
z_{0})\text{{c}}_{x}^{T}\text{{w}}(t)+\text{{c}}_{y}^{T}\text{{w}}(t)-(x_{0}(\rho-
z_{0})-y_{0})+W(t)(\text{{c}}_{x},\text{{c}}_{z})=0$
$\displaystyle\text{{c}}_{z}^{T}\text{{h}}(t)-y_{0}\text{{c}}_{x}^{T}\text{{w}}(t)-x_{0}\text{{c}}_{y}^{T}\text{{w}}(t)+\beta\text{{c}}_{z}^{T}\text{{w}}(t)-(x_{0}y_{0}-\beta
z_{0})-W(t)(\text{{c}}_{x},\text{{c}}_{y})=0$

where $\text{{c}}_{x}=\mathcal{H}(\dot{x})$,
$\text{{c}}_{y}=\mathcal{H}(\dot{y})$ and
$\text{{c}}_{z}=\mathcal{H}(\dot{z})$.

####  5.3.1 Numerical approximation

Using the same techniques used to obtain (30) and (41) to each equation in
(50), we obtain the system to be solved numerically with Newton’s method:

$\displaystyle\bar{\text{{c}}}_{x}-\frac{\sigma}{M}P_{M}^{T}\bar{\text{{c}}}_{y}+\frac{\sigma}{M}P_{M}^{T}\bar{\text{{c}}}_{x}-\frac{\sigma}{M}(y_{0}-x_{0})H_{M}\text{{e}}=0$
(51)
$\displaystyle\bar{\text{{c}}}_{y}+\frac{1}{M}H_{M}\text{{W}}_{M}(\bar{\text{{c}}}_{x},\bar{\text{{c}}}_{z})+x_{0}P_{M}^{T}\bar{\text{{c}}}_{z}-(\rho-
z_{0})P_{M}^{T}\bar{\text{{c}}}_{x}+P_{M}^{T}\bar{\text{{c}}}_{y}-\frac{1}{M}(x_{0}(\rho-
z_{0})-y_{0})H_{M}\text{{e}}=0$
$\displaystyle\bar{\text{{c}}}_{z}-\frac{1}{M}H_{M}\text{{W}}_{M}(\bar{\text{{c}}}_{x},\bar{\text{{c}}}_{y})-y_{0}P_{M}^{T}\bar{\text{{c}}}_{x}-x_{0}P_{M}^{T}\bar{\text{{c}}}_{y}+\beta
P_{M}^{T}\bar{\text{{c}}}_{z}-\frac{1}{M}(x_{0}y_{0}-\beta
z_{0})H_{M}\text{{e}}=0$

######  Remark 5.1.

One interesting remark in [12] is that one can use results from lower
resolutions as initial guesses for the Newton’s method for higher resolution
levels, instead of using a high resolution level right from the start. For the
Lorenz system, this can reduce convergence problems and overall calculation
time.

####  5.3.2 Estimates for the radii polynomials

For the Lorenz system, the maps $T$ and $DT$ are as in (15) and (16)
respectively, with
$\text{{c}}=(\text{{c}}_{x},\text{{c}}_{y},\text{{c}}_{z})\in X_{s}$. The map
$F:X_{s}\to X_{s}$ is given by

$F(\text{{c}}):=\left(\begin{array}[]{l}\text{{c}}_{x}-\sigma
P^{T}\text{{c}}_{y}+\sigma
P^{T}\text{{c}}_{x}-\sigma(y_{0}-x_{0})\text{{e}}_{1}\\\\[5.0pt]
\text{{c}}_{y}+\mathcal{H}\left(\text{{c}}_{x}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{c}}_{z}\right)+x_{0}P^{T}\text{{c}}_{z}-(\rho-
z_{0})P^{T}\text{{c}}_{x}+P^{T}\text{{c}}_{y}-(x_{0}(\rho-
z_{0})-y_{0})\text{{e}}_{1}\\\\[5.0pt]
\text{{c}}_{z}-\mathcal{H}\left(\text{{c}}_{x}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{c}}_{y}\right)-y_{0}P^{T}\text{{c}}_{x}-x_{0}P^{T}\text{{c}}_{y}+\beta
P^{T}\text{{c}}_{z}-(x_{0}y_{0}-\beta z_{0})\text{{e}}_{1}\end{array}\right)$

and for $\text{{v}}=(\text{{v}}_{x},\text{{v}}_{y},\text{{v}}_{z})\in X_{s}$,

$DF(\text{{c}})\text{{v}}:=\left(\begin{array}[]{l}\text{{v}}_{x}-\sigma
P^{T}\text{{v}}_{y}+\sigma P^{T}\text{{v}}_{x}\\\\[5.0pt]
\text{{v}}_{y}+\mathcal{H}\left(\text{{c}}_{x}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{v}}_{z}\right)+\mathcal{H}\left(\text{{c}}_{z}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{v}}_{x}\right)\\\\[2.0pt]
\hskip 70.0001pt+x_{0}P^{T}\text{{v}}_{z}-(\rho-
z_{0})P^{T}\text{{v}}_{x}+P^{T}\text{{v}}_{y}\\\\[5.0pt]
\text{{c}}_{z}-\mathcal{H}\left(\text{{c}}_{x}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{v}}_{y}\right)-\mathcal{H}\left(\text{{c}}_{y}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{v}}_{x}\right)\\\\[2.0pt]
\hskip 70.0001pt-y_{0}P^{T}\text{{v}}_{x}-x_{0}P^{T}\text{{v}}_{y}+\beta
P^{T}\text{{v}}_{z}\end{array}\right).$

For the Lorenz system, we have the following estimates for the radii
polynomial method. Using the norms defined in (48) and (49) and applying the
same techniques and techniques as before to each equation, we find the
following bounds:

$\displaystyle Y_{M}$
$\displaystyle=\|A_{M}F_{M}(\bar{\text{{c}}})\|_{X_{s}}+\frac{\|A_{M}\|_{B(X_{s})}\|\bar{\text{{c}}}_{x}\|_{\ell^{2}}}{(4-\sqrt{2})2^{\frac{3J+5}{2}}}\max\big{\\{}\|\bar{\text{{c}}}_{y}\|_{\ell^{2}},\|\bar{\text{{c}}}_{z}\|_{\ell^{2}}\big{\\}}$
(52) $\displaystyle Y_{\infty}$
$\displaystyle=\max\\{Y_{\infty}^{1},Y_{\infty}^{2},Y_{\infty}^{3}\\}$
$\displaystyle Z_{M}(r)$ $\displaystyle=\alpha_{1}+\alpha_{2}r$ $\displaystyle
Z_{\infty}(r)$ $\displaystyle=\gamma_{1}+\gamma_{2}r$

where the constants above are given by

$\displaystyle Y_{\infty}^{1}$
$\displaystyle:=\frac{|\sigma|}{\sqrt{3}}\frac{\|\bar{\text{{c}}}_{x}-\bar{\text{{c}}}_{y}\|_{\ell^{2}}}{2^{J+2}}$
$\displaystyle Y_{\infty}^{2}$
$\displaystyle:=\frac{\|\bar{\text{{c}}}_{x}\|_{\ell^{2}}\|\bar{\text{{c}}}_{z}\|_{\ell^{2}}}{\sqrt{3}\
2^{2J+5}}+\frac{\|\bar{\text{{c}}}_{x}\|_{\ell^{2}}\|\bar{\text{{c}}}_{z}\|_{\ell^{2}}}{(21\sqrt{7})\,2^{3J+6}}+\frac{\|x_{0}\bar{\text{{c}}}_{z}-(\rho-
z_{0})\bar{\text{{c}}}_{x}+\bar{\text{{c}}}_{y}\|_{\ell^{2}}}{\sqrt{3}\
2^{J+2}}$ $\displaystyle Y_{\infty}^{3}$
$\displaystyle:=\frac{\|\bar{\text{{c}}}_{x}\|_{\ell^{2}}\|\bar{\text{{c}}}_{y}\|_{\ell^{2}}}{2^{2J+5}\
\sqrt{3}}+\frac{\|\bar{\text{{c}}}_{x}\|_{\ell^{2}}\|\bar{\text{{c}}}_{y}\|_{\ell^{2}}}{(21\sqrt{7})\,2^{3J+6}}+\frac{\|\beta\bar{\text{{c}}}_{z}-y_{0}\bar{\text{{c}}}_{x}-x_{0}\bar{\text{{c}}}_{y}\|_{\ell^{2}}}{\sqrt{3}\
2^{J+2}}$ $\displaystyle\alpha_{1}$
$\displaystyle:=\left\|I_{M}-A_{M}B_{1}\right\|_{B(X_{s})}+\frac{\left\|A_{M}B_{2}\right\|_{B(X_{s})}}{2^{J+2}}$
$\displaystyle\hskip
40.00006pt+\frac{1+\sqrt{2}}{2^{\frac{J+3}{2}}}\left\|A_{M}B_{2}\right\|_{B(X_{s})}+\frac{\sqrt{2}}{4-\sqrt{2}}\left\|A_{M}^{\dagger}\begin{pmatrix}0\\\
\|\bar{\text{{c}}}_{x}\|_{\ell^{2}}+\|\bar{\text{{c}}}_{z}\|_{\ell^{2}}\\\
\|\bar{\text{{c}}}_{x}\|_{\ell^{2}}+\|\bar{\text{{c}}}_{y}\|_{\ell^{2}}\end{pmatrix}\right\|_{X_{s}}$
$\displaystyle\alpha_{2}$ $\displaystyle:=\|4C+2D\|_{B(X_{s})}$
$\displaystyle\gamma_{1}$
$\displaystyle:=\max\bigg{\\{}\frac{|\sigma|}{\sqrt{3}\ 2^{J+1}}\ ,\
\frac{D_{1}\big{(}\|\bar{\text{{c}}}_{z}\|_{\ell^{2}}+\|\bar{\text{{c}}}_{x}\|_{\ell^{2}}\big{)}}{2^{\frac{3J}{2}}}+\frac{1+|\rho-
z_{0}|+|x_{0}|}{\sqrt{3}\ 2^{J+2}}\ ,$ $\displaystyle\hskip
80.00012pt\frac{D_{1}\big{(}\|\bar{\text{{c}}}_{z}\|_{\ell^{2}}+\|\bar{\text{{c}}}_{x}\|_{\ell^{2}}\big{)}}{2^{\frac{3J}{2}}}+\frac{\beta+|x_{0}|+|y_{0}|}{\sqrt{3}\
2^{J+2}}\bigg{\\}}$ $\displaystyle\gamma_{2}$
$\displaystyle:=\frac{D_{2}}{2^{\frac{3J}{2}+1}}.$

and the auxiliary quantities to calculate the constants are given by

$A_{M}^{\dagger}:=\begin{bmatrix}\|A_{M_{x,x}}\|_{B(\ell^{2})}&\|A_{M_{x,y}}\|_{B(\ell^{2})}&\|A_{M_{x,z}}\|_{B(\ell^{2})}\\\
\|A_{M_{y,x}}\|_{B(\ell^{2})}&\|A_{M_{y,y}}\|_{B(\ell^{2})}&\|A_{M_{y,z}}\|_{B(\ell^{2})}\\\
\|A_{M_{z,x}}\|_{B(\ell^{2})}&\|A_{M_{z,y}}\|_{B(\ell^{2})}&\|A_{M_{z,z}}\|_{B(\ell^{2})}\end{bmatrix}$
$B_{1}:=\begin{bmatrix}I_{M}+\sigma P_{M}^{T}&-\sigma P_{M}^{T}&0\\\
(B_{M}(\bar{\text{{c}}}_{z})-(\rho-
z_{0})I_{M})P_{M}^{T}&I_{M}+P_{M}^{T}&(B_{M}(\bar{\text{{c}}}_{x})+x_{0}I_{M})P_{M}^{T}\\\
(-B_{M}(\bar{\text{{c}}}_{y})-y_{0}I_{M})P_{M}^{T}&-(B_{M}(\bar{\text{{c}}}_{x})+x_{0}I_{M})P_{M}^{T}&I_{M}+P_{M}^{T}\end{bmatrix}$
$B_{M}(\bar{\text{{c}}}_{i}):=\text{diag}(\tilde{\Omega}_{M}^{T}P_{M}^{T}\bar{\text{{c}}}_{i})+\Gamma_{M}^{T}\text{diag}(P_{M}^{T}\bar{\text{{c}}}_{i})$
$B_{2}:=\begin{bmatrix}\sigma I_{M}&-\sigma I_{M}&0\\\
B_{M}(\bar{\text{{c}}}_{z})&I_{M}&B_{M}(\bar{\text{{c}}}_{x})\\\
-B_{M}(\bar{\text{{c}}}_{y})&-B_{M}(\bar{\text{{c}}}_{x})&I_{M}\end{bmatrix}$
$C=\begin{bmatrix}0&C_{x,y}&C_{x,z}\\\ 0&C_{y,y}&C_{y,z}\\\
0&C_{z,y}&C_{z,z}\end{bmatrix}\quad,\quad
D=\begin{bmatrix}0&D_{x,y}&D_{x,z}\\\ 0&D_{y,y}&D_{y,z}\\\
0&D_{z,y}&D_{z,z}\end{bmatrix}$

and for $i,j\in\\{x,y,z\\}$

$\displaystyle C_{i,j}$
$\displaystyle:=\|A_{M_{i,j}}\,\text{diag}\left(\|(P_{M}^{T})_{i,*}\|_{\ell^{2}}\right)\tilde{\Omega}_{M}^{T}P_{M}^{T}\|+\frac{\left\|A_{M_{i,j}}\tilde{\Omega}_{M}^{T}P_{M}^{T}\right\|}{2^{J+2}}$
$\displaystyle\hskip
80.00012pt+\frac{(1+\sqrt{2})\left\|A_{M_{i,j}}P_{M}^{T}\right\|}{2^{\frac{J+3}{2}}}+\frac{(1+\sqrt{2})\left\|A_{M_{i,j}}\right\|}{2^{\frac{3J+7}{2}}}$
$\displaystyle D_{i,j}$
$\displaystyle:=\left\|A_{M_{i,j}}\Gamma_{M_{i,j}}^{T}\text{diag}(\|(P_{M}^{T})_{i,*}\|_{\ell}^{2})P_{M}^{T}\right\|+\frac{\left\|A_{M_{i,j}}\Gamma_{M_{i,j}}^{T}P_{M}^{T}\right\|}{2^{J+1}}$
$\displaystyle\hskip
80.00012pt+\frac{\left\|A_{M_{i,j}}\Gamma_{M_{i,j}}^{T}\right\|}{2^{2J+4}}+\frac{\sqrt{2}}{\left(4-\sqrt{2}\right)}\,\frac{\|A_{M_{i,j}}\|}{2^{\frac{3J}{2}+4}}$

and $D_{1}$ and $D_{2}$ are as in Proposition 4.4.

####  5.3.3 Results

Figure 6 shows the approximation obtained with the Haar wavelet method using
$J=10$ and numerical integration (fewer points from the Haar wavelet method
are displayed for clarity); for the latter, we used $M^{2}=2^{2J+2}$ points in
order to obtain a precise result. Visually, there seems to be good agreement
between both results.

Figure 6: Results for the Lorenz system

For the same resolution level, the numerical result was rigorously verified by
the radii polynomial method using the estimates described above, obtaining an
$r_{0}=3.9868504\times 10^{-2}$ for $\omega=0.45$ in which the true solution
lies in $X_{s}$.

##  6 Conclusions and future work

We developed a radii polynomial method using the Haar wavelet approach for
differential equations, and illustrated the method by applying it to three
differential equations. One advantage of our method over previous methods
based on the radii polynomials approach is that, due to the use of the Haar
wavelets, our method does not require the solutions to be smooth.

In the future we plan to develop the estimates for higher-order derivatives.
While a higher-order differential equation can be transformed into a system of
first-order equations, this increases the size of the matrices. Thus it might
be interesting to use operators that directly represent higher-order
derivatives. While some of those have already been used for usual numerical
methods, we need to compute the estimates needed for our radii polynomial
method.

Furthermore, we only presented the estimates needed for the radii polynomials
for quadratic nonlinearities, since the main goal of this paper is to present
the general method and illustrate how to compute the estimates and apply the
method. In the future we plan to extend these estimates to include higher-
order polynomial nonlinearities, as this would greatly expand the
applicability of our method.

Lastly, we believe this method can be a basis to build other methods using
similar techniques, such as continuation methods.

### Acknowledgments

The work of Guilherme K. Nakassima was supported by the Coordenação de
Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code
001. The work of Marcio Gameiro was partially supported by the National
Science Foundation under awards DMS-1839294 and HDR TRIPODS award CCF-1934924,
DARPA contract HR0011-16-2-0033, National Institutes of Health award R01
GM126555, and by CNPq grant 309073/2019-7.

## Appendices

###  A Proof of Theorem 4.1

Here we present the full proof of Theorem 4.1. While lengthy, this proof not
only validates the theorem that allows us to deal with quadratic
nonlinearities, but also hints at how to find bounds for the radii
polynomials.

Let us outline the general strategy. We first make heuristical calculations to
find an expression for the Haar transforms; then, by proving that they are
indeed in $\ell^{2}(\mathbb{R})$, the uniqueness of the Haar wavelet series
justifies the calculations.

Let $\text{{c}},\text{{d}}\in\ell^{2}(\mathbb{R})$, and assume that
$\text{{c}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{d}}\in L^{2}([0,1])$. We
have

$\text{{c}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{d}}=\text{{a}}^{T}\text{{h}}(t)\text{{h}}^{T}(t)\text{{b}}$

where we denote $\text{{a}}=P^{T}\text{{c}}$ and $\text{{b}}=P^{T}\text{{d}}$
for simplicity. Let us study the matrix $\text{{h}}(t)\text{{h}}^{T}(t)$ in
more depth. First, by definition, we have $\psi_{j,k}(t)\psi_{j,q}(t)\not=0$
if and only if $k=q$. For $\psi_{j,k}(t)\psi_{p,q}(t)$ when $j\not=p$, without
loss of generality, consider $j<p$ for a fixed $p$. Because of the nesting
property, $\psi_{j,k}$ is constant (possibly zero) in
$\operatorname{supp}\psi_{p,q}$, so we can assume
$\psi_{j,k}(t)=\psi_{j,k}\left(\frac{q+0.5}{2^{p}}\right)$ for
$t\in\operatorname{supp}\psi_{p,q}$. Thus, we can write
$\psi_{j,k}(t)\psi_{p,q}(t)$ as

$\psi_{j,k}(t)\psi_{p,q}(t)=\psi_{j,k}\left(\frac{q+0.5}{2^{p}}\right)\psi_{p,q}(t)=(H_{m/2})_{i,l}\psi_{l}(t).$

Hence, after adjusting indices and taking into account the symmetry of the
matrix, we have proved the following:

######  Lemma .1.

The product $\text{{{h}}}(t)\text{{{h}}}^{T}(t)$ can be recursively calculated
for $m=2^{j}$, $j=0,1,2,\dots$ as

$\text{{{h}}}_{1}(t)\text{{{h}}}_{1}^{T}(t)=\psi_{1}^{2}(t)=\phi^{2}(t)\quad,\quad\text{{{h}}}_{2m}(t)\text{{{h}}}_{2m}^{T}(t)=\begin{bmatrix}\text{{{h}}}_{m}(t)\text{{{h}}}_{m}^{T}(t)&\Upsilon_{m}(t)\\\\[7.22743pt]
\Upsilon_{m}^{T}(t)&\Delta_{m}(t)\end{bmatrix}$

where, for $i,l=1,...,2^{m}$, $\Gamma_{m}$ and $\Delta_{m}$ are $m\times m$
matrices defined element-wise as

$(\Upsilon_{m})_{i,l}(t)=(H_{m})_{i,l}\psi_{m+l}(t)\quad,\quad(\Delta_{m})_{i,l}(t)=\begin{cases}\psi_{i}^{2}(t)&,\text{
if }i=l\\\ 0&,\text{ otherwise}.\end{cases}$

This justifies the decomposition of $\text{{{h}}}(t)\text{{{h}}}^{T}(t)$ as in
Theorem 4.1, that is,
$\text{{h}}(t)\text{{h}}^{T}(t)=\Omega(t)+\Omega^{T}(t)+\Theta(t)$. We treat
each term separately. First, for $\text{{a}}^{T}\Omega(t)\text{{b}}$, dividing
b in $2^{j}$ vector blocks, the product $\Omega(t)\text{{b}}$ is given
element-wise by

$\displaystyle(\Omega(t)\text{{b}})_{i}$
$\displaystyle=\sum_{p=j}^{\infty}(\Upsilon_{2^{p}})_{i,*}\text{{b}}^{*}_{2^{p}}=\sum_{p=j}^{\infty}\sum_{q=1}^{2^{p}}(H_{2^{p}})_{i,q}\,\psi_{2^{p}+q}(t)\,b_{2^{p}+q}$
$\displaystyle=\sum_{p=j}^{\infty}(H_{2^{p}})_{i,*}\,(\text{{b}}_{2^{p}}^{*}\odot\text{{h}}_{2^{p}}^{*}(t))=:(\tilde{\Omega}(\text{{b}}\odot\text{{h}}(t)))_{i}$

where $\odot$ denotes the Hadamard (i.e. element-wise) product, and the matrix
$\tilde{\Omega}$ can be recursively constructed as

$\displaystyle\tilde{\Omega}_{1}=0\quad,\quad\tilde{\Omega}_{2m}=\begin{bmatrix}\tilde{\Omega}_{m}&H_{m}\\\
0_{m}&0_{m}\end{bmatrix}\quad,\text{ for }m=2^{j}\text{ and }j=0,1,2,...$

Thus, the product $\text{{a}}^{T}\Omega(t)\text{{b}}$ is given by

$\displaystyle\text{{a}}^{T}\Omega(t)\text{{b}}$
$\displaystyle=(\text{{c}}^{T}P)\tilde{\Omega}\left(\text{{b}}\odot\text{{h}}(t)\right)=\sum_{i=1}^{\infty}\sum_{l=1}^{\infty}c_{i}(P\tilde{\Omega})_{i,l}b_{l}\psi_{l}(t).$

In particular, if we assume $\text{{a}}^{T}\Omega(t)\text{{b}}\in
L^{2}([0,1])$, then

$\left(\mathcal{H}\left(\text{{a}}^{T}\Omega(t)\text{{b}}\right)\right)_{l}=b_{l}\sum_{i=1}^{\infty}c_{i}(P\tilde{\Omega})_{i,l}=(\text{{c}}^{T}P\tilde{\Omega})_{l}b_{l}.$
(53)

or written in another way,
$\mathcal{H}\left(\text{{a}}^{T}\Omega(t)\text{{b}}\right)=\big{(}\tilde{\Omega}^{T}P^{T}\text{{c}}\big{)}\odot\text{{b}}$.

Also, since the output of $\text{{a}}^{T}\Omega^{T}(t)\text{{b}}$ is a real
number, then
$\text{{a}}^{T}\Omega^{T}(t)\text{{b}}=\left(\text{{a}}^{T}\Omega^{T}(t)\text{{b}}\right)^{T}=\text{{b}}^{T}\Omega(t)\text{{a}}$,
and thus
$\mathcal{H}\left(\text{{a}}^{T}\Omega^{T}(t)\text{{b}}\right)=\big{(}\tilde{\Omega}^{T}P^{T}\text{{d}}\big{)}\odot\text{{a}}$
as well.

For the term $\mathcal{H}(\text{{a}}^{T}\Theta(t)\text{{b}})$, we have

$\text{{a}}^{T}\Theta(t)\text{{b}}=\sum_{i=1}^{\infty}a_{i}b_{i}\psi_{i}^{2}(t)\quad,\quad\psi_{i}^{2}(t)=\begin{cases}2^{j}&,\displaystyle\text{
if }\frac{k}{2^{j}}\leq t\leq\frac{k+1}{2^{j}};\\\\[4.30554pt] 0&,\text{
otherwise.}\end{cases}$

Then $\psi_{i}^{2}\in L^{2}([0,1])$, and we can write it as a Haar wavelet
series:

$\psi_{i}^{2}(t)=\sum_{l=1}^{\infty}\gamma_{i,l}\psi_{l}(t)=\Gamma\text{{h}}(t)\quad,\quad\gamma_{i,l}:=\int_{0}^{1}\psi_{i}^{2}(t)\psi_{l}(t)\,dt.$

Similar to $P$, the matrix $\Gamma$ also has a recursive structure:

######  Lemma .2.

The matrix $\Gamma$ can be recursively calculated by

$\Gamma_{1}=1\quad,\quad\Gamma_{2m}=\begin{bmatrix}\Gamma_{m}&0_{m}\\\
H_{m}^{T}&0_{m}\end{bmatrix}\quad,\text{ for }m=2^{j}\text{ and
}j=0,1,2,3,...$ (54)

###### Proof.

First, calculating $\Gamma_{1}=\gamma_{11}$ is straightforward. By definition,
for $i,l=1,2,...,m$, we have

$(\Gamma_{2m})_{i,l}=\int_{0}^{1}\psi_{i}^{2}(t)\psi_{l}(t)\,dt=(\Gamma_{m})_{i,l}$

Suppose now that $m+1\leq l\leq 2m$, and denote using the two-index notation
$\psi_{l}=\psi_{p,q}$ and $\psi_{i}=\psi_{j,k}$. For $2\leq i\leq 2m$, due to
the nesting property, $\psi_{j,k}$ is constant (possibly zero) in
$\operatorname{supp}\psi_{p,q}$. However, even if $\psi_{j,k}$ is non-zero in
$\operatorname{supp}\psi_{p,q}$,

$(\Gamma_{2m})_{i,l}=\int_{0}^{1}\psi_{i}^{2}(t)\psi_{l}(t)\,dt=2^{j}\int_{\frac{k}{2^{j}}}^{\frac{k+1}{2^{j}}}\psi_{p,q}(t)\,dt=2^{j}\int_{\frac{q}{2^{p}}}^{\frac{q+1}{2^{p}}}\psi_{p,q}(t)\,dt=0.$

Similar reasoning applies for $i=1$ (in which case $\psi_{1}(t)=\phi(t)\equiv
1$ in $[0,1]$).

Suppose now that $m+1\leq i\leq 2m$ and $l\leq m$. For $l=1$ a straightforward
calculation show that $(\Gamma_{2m})_{i,l}=1=(H_{m})_{1,i}$. For $l>1$, since
$p<j$, we can apply the same reasoning as in the proof of Lemma .1, yielding

$\gamma_{i,l}=2^{j}\int_{\frac{k}{2^{j}}}^{\frac{k+1}{2^{j}}}\psi_{p,q}(t)\,dt=2^{j}\int_{\frac{k}{2^{j}}}^{\frac{k+1}{2^{j}}}\psi_{p,q}\left(\frac{k+0.5}{2^{j}}\right)\,dt=\psi_{l}(t_{i})=(H_{m})_{l,i}.$

∎

Hence,

$\text{{a}}^{T}\Theta(t)\text{{b}}=\sum_{i=1}^{\infty}a_{i}b_{i}\psi_{i}^{2}(t)=\sum_{i=1}^{\infty}a_{i}b_{i}\sum_{l=1}^{\infty}\gamma_{i,l}\psi_{l}(t)=(\text{{a}}\odot\text{{b}})^{T}\Gamma\,\text{{h}}(t)$

Thus, if $\Gamma^{T}(\text{{a}}\odot\text{{b}})\in\ell^{2}(\mathbb{R})$, then
$\mathcal{H}(\text{{a}}^{T}\Theta(t)\text{{b}})=\Gamma^{T}(\text{{a}}\odot\text{{b}})$.

Now we must prove that our tentative Haar transforms are indeed elements of
$\ell^{2}(\mathbb{R})$; Theorem 4.1 follows then from the uniqueness of the
Haar series. Before that, we prove a few lemmas:

######  Lemma .3.

The matrix $P\tilde{\Omega}$ is recursively given by

$P_{2m}\tilde{\Omega}_{2m}=\begin{bmatrix}P_{m}\tilde{\Omega}_{m}&P_{m}H_{m}\\\\[10.0pt]
\displaystyle\frac{1}{4\sqrt{m^{3}}}H_{m}^{T}\tilde{\Omega}_{m}&\displaystyle\frac{1}{4\sqrt{m}}I_{m}\end{bmatrix}.$
(55)

This is proven by multiplying the recursive formulas for $P$ and
$\tilde{\Omega}$.

######  Lemma .4.

The matrix $H_{m}^{T}\tilde{\Omega}_{m}$ is given element-wise as

$(H_{m}^{T}\tilde{\Omega}_{m})_{i,l}=\begin{cases}0&\text{, if }l=1;\\\
\psi_{l}^{2}(t_{i})&\text{, otherwise.}\end{cases}$ (56)

###### Proof.

First, observe that, since the first column of the matrix $\tilde{\Omega}_{m}$
is zero, then the first column of $H_{m}^{T}\tilde{\Omega}_{m}$ is also zero,
proving the case $l=1$.

For $l\geq 2$, fix an element $(H_{m}^{T}\tilde{\Omega}_{m})_{i,l}$ and make
$\psi_{i}=\psi_{j,k}$ and $\psi_{l}=\psi_{p,q}$ using the two-index notation.
Due to the structure of $H_{m}^{T}$ and $\tilde{\Omega}_{m}$,

$\displaystyle(H_{m}^{T}\tilde{\Omega}_{m})_{i,l}=\text{{h}}_{2^{p}}^{T}(t_{i})\text{{h}}_{2^{p}}(t_{q+1})=1+\sum_{r=0}^{p-1}\sum_{s=0}^{2^{r}-1}\psi_{r,s}(t_{i})\psi_{r,s}(t_{q+1}).$

where $t_{i}=2^{-j}(k-0.5)$ and $\displaystyle t_{q+1}=2^{-p}(q+0.5)$. Since
$t_{i}$ and $t_{q+1}$ are fixed, for each $r$ there is at most a single
wavelet $\psi_{r,s_{r}}$ whose support contains both $t_{i}$ and $t_{q+1}$,
because the intervals where wavelets at the same resolution level are non-zero
do not overlap. Thus,

$(H_{m}^{T}\tilde{\Omega}_{m})_{i,l}=1+\sum_{r=0}^{p-1}\psi_{r,s_{r}}(t_{i})\psi_{r,s_{r}}(t_{q+1})$

If neither $\psi_{r,s_{r}}(t_{i})$ nor $\psi_{r,s_{r}}(t_{q+1})$ are zero,
only two cases may occur: either
$\psi_{r,s_{r}}(t_{i})=\psi_{r,s_{r}}(t_{q+1})$ or
$\psi_{r,s_{r}}(t_{i})=-\psi_{r,s_{r}}(t_{q+1})$. The possible situations are
depicted in Figure 7, supposing without loss of generality that $t_{i}\leq
t_{q+1}$.

$t_{i}$$t_{q+1}$ (a) $\psi_{r,s_{r}}(t_{i})=\psi_{r,s_{r}}(t_{q+1})$

$t_{i}$$t_{q+1}$ (b) $\psi_{r,s_{r}}(t_{i})=-\psi_{r,s_{r}}(t_{q+1})$

$t_{i}$$t_{q+1}$ (c) $\psi_{r,s_{r}}(t_{i})=\psi_{r,s_{r}}(t_{q+1})$

Figure 7: All possible situations for the product
$\psi_{r,s_{r}}(t_{i})\psi_{r,s_{r}}(t_{q+1})$

Let us study what happens when we change the resolution level $r$:

  1. a)

Suppose that $\psi_{r,s_{r}}(t_{i})=\psi_{r,s_{r}}(t_{q+1})$ (situations (a)
or (c) in Figure 7) for every $r=0,...,p-1$. Then $t_{i}=t_{q+1}$, since for
$r=p-1$ we only sample $\psi_{r,s_{r}}$ at the times
$t_{1}=\frac{s_{r}+0.25}{2^{r}}$ and $t_{2}=\frac{s_{r}+0.75}{2^{r}}$, and we
have $\psi_{r,s_{r}}(t_{1})=-\psi_{r,s_{r}}(t_{2})$. Thus

$(H_{m}^{T}\tilde{\Omega}_{m})_{i,l}=1+\sum_{r=0}^{p-1}\psi_{r,s_{r}}^{2}(t_{i})=1+\sum_{r=0}^{p-1}2^{r}=2^{p}=\psi_{l}^{2}(t_{i})$

  2. b)

Suppose that $\psi_{r,s_{r}}(t_{i})=-\psi_{r,s_{r}}(t_{q+1})$ (situation (b)
in Figure 7) happens for some $r$ for some $r\leq p-1$. Then
$t_{i}\in[\frac{s}{2^{r}},\frac{s+0.5}{2^{r}}]$ and
$t_{q+1}\in[\frac{s+0.5}{2^{r}},\frac{s+1}{2^{r}}]$. Due to the nesting
property of the Haar wavelets, no finer wavelet has both $t_{i}$ and $t_{q+1}$
in its support, and for all coarser resolutions $\rho<r$ we have
$\psi_{\rho,s_{\rho}}(t_{i})=\psi_{\rho,s_{\rho}}(t_{q+1})$. Hence

$(H_{m}^{T}\tilde{\Omega}_{m})_{i,l}=1+\sum_{\rho=0}^{r-1}\psi_{\rho,s_{\rho}}^{2}(t_{i})-\psi_{r,s_{r}}^{2}(t_{i})=1+\sum_{\rho=0}^{r-1}2^{\rho}-2^{r}=0$

Lastly, since $p>r$ and $t_{q}\in\text{supp }\psi_{p,q}$, then
$t_{i}\not\in\text{supp}\psi_{p,q}$. Thus
$\psi_{l}^{2}(t_{i})=\psi_{p,q}^{2}(t_{i})=0=(H_{m}^{T}\tilde{\Omega}_{m})_{i,l}$.

∎

Now we finally prove that
$(\tilde{\Omega}^{T}P^{T}\text{{{c}}})\odot\text{{{b}}}$ and
$\Gamma^{T}(\text{{{a}}}\odot\text{{{b}}})$ are indeed in
$\ell^{2}(\mathbb{R})$.

######  Proposition .1.

The sequences $(\tilde{\Omega}^{T}P^{T}\text{{{c}}})\odot\text{{{b}}}$ and
$\Gamma^{T}(\text{{{a}}}\odot\text{{{b}}})$ are in $\ell^{2}(\mathbb{R})$ and
satisfy, for some $C_{1},C_{2}>0$,

$\displaystyle\big{\|}(\tilde{\Omega}^{T}P^{T}\text{{{c}}})\odot\text{{{b}}}\big{\|}_{\ell^{2}}$
$\displaystyle\leq
C_{1}\|\text{{{c}}}\|_{\ell^{2}}\|\text{{{d}}}\|_{\ell^{2}}$
$\displaystyle\left\|\Gamma^{T}(\text{{{a}}}\odot\text{{{b}}})\right\|_{\ell^{2}(\mathbb{R})}$
$\displaystyle\leq
C_{2}\|\text{{{c}}}\|_{\ell^{2}}\|\text{{{d}}}\|_{\ell^{2}}$

###### Proof.

For $(\tilde{\Omega}^{T}P^{T}\text{{{c}}})\odot\text{{{b}}}$, if we divide c
as in (8) and using the block structure from Lemma .3, each block of
$(\tilde{\Omega}^{T}P^{T}\text{{{c}}})$ is given by

$\displaystyle(\tilde{\Omega}^{T}P^{T}\text{{{c}}})_{1}$
$\displaystyle=\sum_{r=0}^{\infty}\frac{1}{2^{\frac{3r}{2}+2}}(\text{{{c}}}^{*}_{2^{r}})^{T}H^{T}_{2^{r}}\tilde{\Omega}_{2^{r}}$
$\displaystyle(\tilde{\Omega}^{T}P^{T}\text{{{c}}})^{*}_{2^{j}}$
$\displaystyle=(c_{1},...,c_{2^{j}})P_{2^{j}}H_{2^{j}}+\frac{1}{2^{\frac{j}{2}+2}}(\text{{{c}}}^{*}_{2^{j}})^{T}+\sum_{r=j+1}^{\infty}\frac{1}{2^{\frac{3r}{2}+2}}(\text{{{c}}}^{*}_{2^{r}})^{T}H^{T}_{2^{r}}\tilde{\Omega}_{2^{r}}$

We bound each term in the right-hand side:

  * •

$\|(c_{1},...,c_{2^{j}})P_{2^{j}}H_{2^{j}}\|_{\ell^{2}}\leq\big{\|}H_{2^{j}}^{T}\big{\|}\big{\|}P_{2^{j}}^{T}(c_{1},...,c_{2^{j}})\big{\|}_{\ell^{2}}\leq\frac{\|\text{{{c}}}\|_{\ell^{2}}}{2^{\frac{j}{2}+1}}$

  * •

$\left\|\frac{1}{2^{\frac{j}{2}+2}}(\text{{{c}}}^{*}_{2^{j}})^{T}\right\|_{\ell^{2}}\leq\frac{\left\|\text{{{c}}}\right\|_{\ell^{2}}}{2^{\frac{j}{2}+2}}$

  * •

From Lemma .4, $\max H^{T}_{2^{j}}\Omega^{T}_{2^{j}}=2^{j}$; thus

$\left\|\sum_{r=j+1}^{\infty}\frac{1}{2^{\frac{3r}{2}+2}}(\text{{{c}}}^{*}_{2^{r}})^{T}H^{T}_{2^{r}}\Omega^{T}_{2^{r}}\right\|_{\ell^{2}}\leq\sum_{r=j+1}^{\infty}\frac{1}{2^{\frac{3r}{2}+2}}\left\|\text{{{c}}}\right\|_{\ell^{2}}\big{\|}H^{T}_{2^{r}}\Omega^{T}_{2^{r}}\big{\|}_{\ell^{2}}\leq\frac{1+\sqrt{2}}{2^{\frac{j}{2}+2}}\left\|\text{{{c}}}\right\|_{\ell^{2}}$

Hence,

$\big{\|}(\tilde{\Omega}^{T}P^{T}\text{{{c}}})^{*}_{2^{j}}\big{\|}_{\ell^{2}}\leq\left(1+\frac{\sqrt{2}}{8}\right)\frac{\left\|\text{{{c}}}\right\|_{\ell^{2}}}{2^{\frac{j}{2}}}\quad,\quad|(\tilde{\Omega}^{T}P^{T}\text{{{c}}})_{1}|\leq\frac{1+\sqrt{2}}{4}\left\|\text{{{c}}}\right\|_{\ell^{2}}$
(57)

Also, for $i=2^{j},2^{j}+1,...,2^{2j}-1$,
$\big{|}(P^{T}\text{{{d}}})_{i}\big{|}\leq\big{\|}(P^{T}\text{{{d}}})^{*}_{2^{j}}\big{\|}_{\ell^{2}}\leq\frac{\left\|\text{{{d}}}\right\|_{\ell^{2}}}{2^{j+1}}$,
and thus

$\displaystyle\left\|\big{(}(\tilde{\Omega}^{T}P^{T}\text{{{c}}})\odot\text{{{b}}}\big{)}^{*}_{2^{j}}\right\|_{\ell^{2}}$
$\displaystyle\leq\big{\|}(\tilde{\Omega}^{T}P^{T}\text{{{c}}})^{*}_{2^{j}}\big{\|}_{\ell^{2}}\max_{2^{j}\leq
i\leq
2^{2j}-1}\big{|}(P^{T}\text{{{d}}})_{i}\big{|}\leq\left(1+\frac{\sqrt{2}}{8}\right)\frac{\left\|\text{{{c}}}\right\|_{\ell^{2}}\left\|\text{{{d}}}\right\|_{\ell^{2}}}{2^{\frac{3j}{2}+1}}$

Therefore, we can bound
$(\tilde{\Omega}^{T}P^{T}\text{{{c}}})\odot\text{{{b}}}$ by

$\displaystyle\left\|(\tilde{\Omega}^{T}P^{T}\text{{{c}}})\odot\text{{{b}}}\right\|_{\ell^{2}}^{2}$
$\displaystyle=\big{|}\big{(}(\tilde{\Omega}^{T}P^{T}\text{{{c}}})\odot\text{{{b}}}\big{)}_{1}\big{|}^{2}+\sum_{j=0}^{\infty}\left\|\big{(}(\tilde{\Omega}^{T}P^{T}\text{{{c}}})\odot\text{{{b}}}\big{)}^{*}_{2^{j}}\right\|_{\ell^{2}}^{2}$
$\displaystyle\leq\frac{1}{112}(169+79\sqrt{2})\left\|\text{{{c}}}\right\|_{\ell^{2}}^{2}\left\|\text{{{d}}}\right\|_{\ell^{2}}^{2}=C_{1}^{2}\left\|\text{{{c}}}\right\|_{\ell^{2}}^{2}\left\|\text{{{d}}}\right\|_{\ell^{2}}^{2}.$

For $\Gamma^{T}(\text{{{a}}}\odot\text{{{b}}})$, from (54) we have

$\displaystyle\left|\left(\Gamma^{T}(\text{{{a}}}\odot\text{{{b}}})\right)_{1}\right|$
$\displaystyle=\left|a_{1}b_{1}+\sum_{r=0}^{\infty}H_{2^{r}}\left(\text{{{a}}}^{*}_{2^{r}}\odot\text{{{b}}}^{*}_{2^{r}}\right)\right|\leq|a_{1}b_{1}|+\sum_{r=0}^{\infty}\|H_{2^{r}}\|\left\|\text{{{a}}}^{*}_{2^{r}}\right\|_{\ell^{2}}\max\left|\text{{{b}}}^{*}_{2^{r}}\right|$
$\displaystyle\leq\left\|\text{{{c}}}\right\|_{\ell^{2}}\left\|\text{{{d}}}\right\|_{\ell^{2}}+\sum_{r=0}^{\infty}\frac{\left\|\text{{{c}}}\right\|_{\ell^{2}}\left\|\text{{{d}}}\right\|_{\ell^{2}}}{2^{\frac{3r}{2}+2}}\leq\frac{18+\sqrt{2}}{14}\left\|\text{{{c}}}\right\|_{\ell^{2}}\left\|\text{{{d}}}\right\|_{\ell^{2}}$
$\displaystyle\left\|\left(\Gamma^{T}(\text{{{a}}}\odot\text{{{b}}})\right)^{*}_{2^{j}}\right\|_{\ell^{2}}$
$\displaystyle=\left\|\sum_{r=j+1}^{\infty}H_{2^{r}}\left(\text{{{a}}}^{*}_{2^{r}}\odot\text{{{b}}}^{*}_{2^{r}}\right)\right\|_{\ell^{2}}\leq\sum_{r=j+1}^{\infty}\|H_{2^{r}}\|\left\|\text{{{a}}}^{*}_{2^{r}}\right\|_{\ell^{2}}\max\left|\text{{{b}}}^{*}_{2^{r}}\right|$
$\displaystyle\leq\sum_{r=j+1}^{\infty}\frac{\left\|\text{{{c}}}\right\|_{\ell^{2}}\left\|\text{{{d}}}\right\|_{\ell^{2}}}{2^{\frac{3r}{2}+2}}=\frac{1+2\sqrt{2}}{4}\frac{\left\|\text{{{c}}}\right\|_{\ell^{2}}\left\|\text{{{d}}}\right\|_{\ell^{2}}}{2^{\frac{3r}{2}}}$

Finally, we bound the sequence $\Gamma^{T}(\text{{{a}}}\odot\text{{{b}}})$
with

$\displaystyle\left\|\Gamma^{T}(\text{{{a}}}\odot\text{{{b}}})\right\|_{\ell^{2}}^{2}$
$\displaystyle=\left|\left(\Gamma^{T}(\text{{{a}}}\odot\text{{{b}}})\right)_{1}\right|^{2}+\sum_{j=0}^{\infty}\left\|\left(\Gamma^{T}(\text{{{a}}}\odot\text{{{b}}})\right)^{*}_{2^{j}}\right\|_{\ell^{2}}^{2}$
$\displaystyle\leq\frac{1}{49}(113+23\sqrt{2})\left\|\text{{{c}}}\right\|_{\ell^{2}}\left\|\text{{{d}}}\right\|_{\ell^{2}}=C_{2}^{2}\left\|\text{{{c}}}\right\|_{\ell^{2}}\left\|\text{{{d}}}\right\|_{\ell^{2}}$

∎

###  B Proofs of quadratic estimates from Section 4

Here we prove the quadratic estimates from Propositions 4.1—4.4. As stated in
the paper, the main strategy is to employ both the recursive block structures
of the matrices from Theorems 2.1 and 4.1 and the finite-infinite
decomposition from (17). For clarity, Figure 8 shows how they overlap for the
operator $P$; the other matrices follow a similar pattern. We also draw
insights from A to bound the sums that appear in the proof. We believe that
similar procedures may be applied for higher-degree polynomial nonlinearities.

Figure 8: Overlay of the block structure and finite-infinite decomposition for
$P$

####  B.1 Proof of Proposition 4.1

  1. i)

For $i\leq M$, since $\bar{\text{{{c}}}}\in\mathbb{R}^{M}$, the $i$-th element
of $\Pi_{M}P^{T}\bar{\text{{{c}}}}$ is given by

$\left(P^{T}\bar{\text{{{c}}}}\right)_{i}=\sum_{l=1}^{M}\left(P^{T}\right)_{i,l}c_{l}=\left(P_{M}^{T}\bar{\text{{{c}}}}\right)_{i}.$

  2. ii)

For $i\leq M$,

$\left(\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)\bar{\text{{b}}})\right)_{i}=(\bar{\text{{{c}}}}^{T}P\,\tilde{\Omega})_{i}\bar{\text{{b}}}_{i}=\sum_{l=1}^{M}c_{l}(P\,\tilde{\Omega})_{l,i}\bar{\text{{b}}}_{i}=(\bar{\text{{{c}}}}^{T}P_{M}\tilde{\Omega}_{M})_{i}\bar{\text{{b}}}_{i}.$

with the last equality due to Lemma .3. Thus, from item _(i)_

$\displaystyle\Pi_{M}\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)\bar{\text{{b}}})$
$\displaystyle=(\bar{\text{{{c}}}}^{T}P_{M}\tilde{\Omega}_{M})^{T}\odot\Pi_{M}\bar{\text{{b}}}=(\tilde{\Omega}_{M}^{T}P_{M}^{T}\bar{\text{{{c}}}})\odot(P_{M}^{T}\bar{\text{{{d}}}}).$

  3. iii)

Using the finite-infinite decomposition, we can separate
$\Pi_{M}\mathcal{H}(\bar{\text{{a}}}^{T}\Theta(t)\bar{\text{{b}}})$ as

$\Pi_{M}\mathcal{H}(\bar{\text{{a}}}^{T}\Theta(t)\bar{\text{{b}}})=\Gamma^{T}_{M}(\bar{\text{{a}}}_{M}\odot\bar{\text{{b}}}_{M})+\Gamma^{T}_{\infty,M}(\bar{\text{{a}}}_{\infty}\odot\bar{\text{{b}}}_{\infty})$

We now need to bound the infinite sum in the second term. Adapting the
expression of $\bar{\text{{b}}}=P^{T}\bar{\text{{{d}}}}$ in the proof of
Proposition 2.2 for $2^{j}+l>M$ and $\bar{\text{{{d}}}}\in\mathbb{R}^{M}$,

$\bar{b}_{2^{j}+l}=-\frac{1}{2^{\frac{3j}{2}+2}}\sum_{q=1}^{M}\left(H_{2^{j}}^{T}\right)_{l,q}\bar{c}_{q},$

and thus, using the structure of $\Gamma$, (58) and the fact that
$H_{2^{j}}H_{2^{j}}^{T}=2^{j}I_{2^{j}}$,

$\displaystyle\left|\left(\Gamma^{T}_{\infty}(\Pi_{\infty}\bar{\text{{a}}}\odot\Pi_{\infty}\bar{\text{{b}}})\right)_{i}\right|$
$\displaystyle=\left|\sum_{j=J+1}^{\infty}\sum_{l=1}^{2^{j}}\left(H_{2^{j}}\right)_{i,l}\bar{a}_{2^{j}+l}\bar{b}_{2^{j}+l}\right|$
$\displaystyle\leq\sum_{j=J+1}^{\infty}\left|\sum_{l=1}^{2^{j}}-\left(H_{2^{j}}\right)_{i,l}\left(\frac{1}{2^{\frac{3j}{2}+2}}\sum_{q=1}^{M}\left(H_{2^{j}}^{T}\right)_{l,q}c_{q}\right)\bar{b}_{2^{j}+l}\right|$
$\displaystyle\leq\sum_{j=J+1}^{\infty}\frac{1}{2^{\frac{3j}{2}+2}}\left|\left(\sum_{q=1}^{M}\sum_{l=1}^{2^{j}}\left(H_{2^{j}}\right)_{i,l}\left(H_{2^{j}}^{T}\right)_{l,q}c_{q}\right)\right|\frac{\|\bar{\text{{{d}}}}\|_{\ell^{2}}}{2^{j+1}}=\frac{\sqrt{2}\,\|\bar{\text{{{d}}}}\|_{\ell^{2}}}{\left(4-\sqrt{2}\right)2^{\frac{3J}{2}+3}}|c_{i}|$

and thus

$\|\Gamma^{T}_{\infty}(\Pi_{\infty}\bar{\text{{a}}}\odot\Pi_{\infty}\bar{\text{{b}}})\|_{\ell^{2}}\leq\frac{\sqrt{2}\,\|\bar{\text{{{c}}}}\|_{\ell^{2}}\|\bar{\text{{{d}}}}\|_{\ell^{2}}}{\left(4-\sqrt{2}\right)2^{\frac{3J}{2}+3}}.$

####  B.2 Proof of Proposition 4.2

  1. i)

Using the finite-infinite decomposition, the recursive block structure of $P$
and the fact that $\bar{\text{{c}}}\in\mathbb{R}^{M}$, we have for $j>J$

$\left\|\left(P^{T}\bar{\text{{{c}}}}\right)_{2^{j}}^{*}\right\|_{\ell^{2}}=\left\|-\frac{1}{2^{\frac{3j}{2}+2}}H_{2^{j}}^{T}\bar{\text{{{c}}}}\right\|_{\ell^{2}}\leq\frac{1}{2^{\frac{3j}{2}+2}}\left\|H_{2^{j}}^{T}\right\|\left\|\bar{\text{{{c}}}}\right\|_{\ell^{2}}=\frac{\left\|\bar{\text{{{c}}}}\right\|_{\ell^{2}}}{2^{j+2}}$
(58)

Estimating as we did in Theorem .1,

$\left\|\Pi_{\infty}\left(P^{T}\bar{\text{{{c}}}}\right)\right\|_{\ell^{2}}=\sqrt{\sum_{j=J+1}^{\infty}\left\|\left(P^{T}\bar{\text{{{c}}}}\right)_{2^{j}}^{*}\right\|_{\ell^{2}}^{2}}\leq\frac{1}{\sqrt{3}}\frac{\left\|\bar{\text{{{c}}}}\right\|_{\ell^{2}}}{2^{J+2}}.$

  2. ii)

For $j>J$, dividing
$\Pi_{\infty}\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)\bar{\text{{b}}})$ in
blocks,

$\displaystyle\left\|\left(\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)\bar{\text{{b}}})\right)_{2^{j}}^{*}\right\|_{\ell^{2}}$
$\displaystyle=\left\|-\frac{1}{2^{\frac{3j}{2}+2}}\left(H_{2^{j}}^{T}P_{2^{j}}^{T}\bar{\text{{{c}}}}\right)\odot\left(H_{2^{j}}^{T}\bar{\text{{{d}}}}\right)\right\|_{\ell^{2}}\leq\frac{1}{2^{\frac{3j}{2}+2}}\left(\max_{1\leq
k\leq
2^{j}}\left|\left(H_{2^{j}}^{T}\bar{\text{{{c}}}}\right)_{k}\right|\right)\left\|H_{2^{j}}^{T}P_{2^{j}}^{T}\bar{\text{{{d}}}}\right\|_{\ell^{2}}$
$\displaystyle\leq\frac{1}{2^{\frac{3j}{2}+2}}\left\|H_{2^{j}}^{T}\bar{\text{{{c}}}}\right\|_{\ell^{2}}\left\|H_{2^{j}}^{T}P_{2^{j}}^{T}\bar{\text{{{d}}}}\right\|_{\ell^{2}}\leq\frac{\left\|\bar{\text{{{c}}}}\right\|_{\ell^{2}}\left\|\bar{\text{{{d}}}}\right\|_{\ell^{2}}}{2^{j+4}}$

and thus

$\displaystyle\left\|\Pi_{\infty}\left(\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)\bar{\text{{b}}})\right)\right\|_{\ell^{2}}$
$\displaystyle=\sqrt{\sum_{j=J+1}^{\infty}\left\|\left(\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)\bar{\text{{b}}})\right)_{2^{j}}^{*}\right\|_{\ell^{2}}^{2}}\leq\frac{1}{\sqrt{3}}\frac{\left\|\bar{\text{{{c}}}}\right\|_{\ell^{2}}\left\|\bar{\text{{{d}}}}\right\|_{\ell^{2}}}{2^{2J+4}}.$

  3. iii)

To estimate
$\Pi_{\infty}\mathcal{H}(\bar{\text{{a}}}^{T}\Theta(t)\bar{\text{{a}}})$, with
the block matrix structure from (54) for $i>M$,

$\displaystyle\left|\left(\mathcal{H}(\bar{\text{{a}}}^{T}\Theta(t)\bar{\text{{b}}})\right)_{i}\right|$
$\displaystyle=\left|\Gamma^{T}_{i,*}(\bar{\text{{a}}}\odot\bar{\text{{b}}})\right|=\left|\sum_{p=j+1}^{\infty}\frac{1}{2^{3p+4}}\sum_{q=1}^{2^{p}}(H_{2^{p}})_{i,q}\left((H^{T}_{2^{p}})_{q,*}\bar{\text{{{c}}}}\right)\left((H^{T}_{2^{p}})_{q,*}\bar{\text{{{d}}}}\right)\right|$
$\displaystyle\leq\sum_{p=j+1}^{\infty}\frac{1}{2^{3p+4}}\max_{1\leq q\leq
2^{p}}\left|\left((H^{T}_{2^{p}})_{q,*}\bar{\text{{{d}}}}\right)\right|\sum_{q=1}^{2^{p}}\left|(H_{2^{p}})_{i,q}\left((H^{T}_{2^{p}})_{q,*}\bar{\text{{{c}}}}\right)\right|$
$\displaystyle\leq\sum_{p=j+1}^{\infty}\left|(I_{2^{p}})_{i,*}\bar{\text{{{c}}}}\right|\frac{\|\bar{\text{{{d}}}}\|_{\ell^{2}}}{2^{3p+6}}=\frac{|c_{i}|}{2^{3j+6}}\frac{\|\bar{\text{{{d}}}}\|_{\ell^{2}}}{7}$

and hence

$\left\|\Pi_{\infty}\mathcal{H}(\bar{\text{{a}}}^{T}\Theta(t)\bar{\text{{b}}})\right\|_{\ell^{2}}\leq\frac{1}{21\sqrt{7}}\frac{\|\bar{\text{{{c}}}}\|_{\ell^{2}}\|\bar{\text{{{d}}}}\|_{\ell^{2}}}{2^{3J+6}}.$

####  B.3 Proof of Proposition 4.3

The following Lemma helps estimating terms of the type
$A_{M}(B_{M}\bar{\text{{{x}}}}\odot C_{M}\bar{\text{{{y}}}})$ for arbitrary
$\bar{\text{{{x}}}},\bar{\text{{{y}}}}\in\mathbb{R}^{M}$, which would be
tricky otherwise.

######  Lemma .5.

Given $\bar{\text{{{x}}}},\bar{\text{{{y}}}}\in\mathbb{R}^{M}$ and $M\times M$
matrices $A_{M},B_{M},C_{M}$, then

$\left\|A_{M}(B_{M}\bar{\text{{{x}}}}\odot
C_{M}\bar{\text{{{y}}}})\right\|_{\ell^{2}}\leq\left\|A_{M}\,\text{{diag}}(\|(C_{M})_{l,*}\|_{\ell^{2}})B_{M}\right\|\|\bar{\text{{{x}}}}\|_{\ell^{2}}\|\bar{\text{{{y}}}}\|_{\ell^{2}},$
(59)

where

$\text{{diag}}(\bar{\text{{{x}}}}):=\left[\begin{array}[]{cccc}x_{1}&\lx@intercol\hfil\phantom{00}\hfil\lx@intercol&\lx<EMAIL_ADDRESS>\lx@intercol\hfil\phantom{00}\hfil\lx@intercol&x_{2}&\lx@intercol\hfil\phantom{00}\hfil\lx@intercol&\lx<EMAIL_ADDRESS>\lx@intercol\hfil\hbox{\multirowsetup$0$}\hfil\lx@intercol&\ddots&\lx@intercol\hfil\phantom{00}\hfil\lx@intercol\\\
\lx@intercol\hfil\phantom{00}\hfil\lx@intercol&\lx@intercol\hfil\phantom{00}\hfil\lx@intercol&\lx@intercol\hfil\phantom{00}\hfil\lx@intercol&x_{n}\end{array}\right]\quad,\quad\left(\|(C_{M})_{l,*}\|_{\ell^{2}}\right)=\begin{pmatrix}\|(C_{M})_{1,*}\|_{\ell^{2}}\\\
\|(C_{M})_{2,*}\|_{\ell^{2}}\\\ \vdots\\\ \|(C_{M})_{M,*}\|_{\ell^{2}}\\\
\end{pmatrix}.$

###### Proof.

For each element of $A_{M}(B_{M}\bar{\text{{{x}}}}\odot
C_{M}\bar{\text{{{y}}}})$

$\displaystyle\left|\left(A_{M}(B_{M}\bar{\text{{{x}}}}\odot
C_{M}\bar{\text{{{y}}}})\right)_{i}\right|$
$\displaystyle=\left|\sum_{l=1}^{M}(A_{M})_{i,l}(B_{M}\bar{\text{{{x}}}})_{l}(C_{M}\bar{\text{{{y}}}})_{l}\right|$
$\displaystyle\leq\left|\sum_{l=1}^{M}(A_{M})_{i,l}(B_{M}\bar{\text{{{x}}}})_{l}\|(C_{M})_{l,*}\|_{\ell^{2}}\|\bar{\text{{{y}}}}\|_{\ell^{2}}\right|=\left|(A_{M})_{i,*}\text{diag}\left(\|(C_{M})_{l,*}\|_{\ell^{2}}\right)B_{M}\bar{\text{{{x}}}}\right|\|\bar{\text{{{y}}}}\|_{\ell^{2}}.$

Taking the $\ell^{2}$ norm, the result follows. ∎

  1. i)

Applying the finite-infinite decomposition to $P$ and observing its block
structure as in Figure 8, we have that

$\displaystyle\left\|(P_{\infty,M}^{T}\,\text{{y}}_{\infty})\right\|_{\ell^{2}}$
$\displaystyle=\left\|\sum_{q=J+1}^{\infty}\frac{1}{2^{\frac{3q}{2}+2}}\left(H_{2^{q}}\right)_{1:M,*}\text{{y}}_{2^{q}}^{*}\right\|_{\ell^{2}}\leq\sum_{q=J+1}^{\infty}\frac{1}{2^{\frac{3q}{2}+2}}\left\|\left(H_{2^{q}}\right)_{1:M,*}\right\|\left\|\text{{y}}_{2^{q}}^{*}\right\|_{\ell^{2}}$
$\displaystyle\leq\sum_{q=J+1}^{\infty}\frac{1}{2^{\frac{3q}{2}+2}}\left\|H_{2^{q}}\right\|\left\|\text{{y}}_{2^{q}}^{*}\right\|_{\ell^{2}}\leq\frac{\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{J+2}}.$
(60)

  2. ii)

Fix $j>J$ and make $\text{{z}}_{2^{j}}:=H_{2^{j}}\text{{y}}_{2^{j}}^{*}$. From
Lemma .3,

$\displaystyle\left(\tilde{\Omega}^{T}P^{T}\right)_{M,\infty}\text{{y}}_{\infty}$
$\displaystyle=\sum_{j=J+1}^{\infty}\frac{1}{2^{\frac{3j}{2}+2}}\left(\tilde{\Omega}_{2^{j}}^{T}H_{2^{j}}\right)_{1:M,*}\text{{y}}_{2^{j}}^{*}$
$\displaystyle=\sum_{j=J+1}^{\infty}\frac{1}{2^{\frac{3j}{2}+2}}\left(\tilde{\Omega}_{2^{j}}^{T}\right)_{1:M,*}H_{2^{j}}\text{{y}}_{2^{j}}^{*}=\sum_{j=J+1}^{\infty}\frac{1}{2^{\frac{3j}{2}+2}}\left(\tilde{\Omega}_{2^{j}}^{T}\right)_{1:M,*}\text{{z}}_{2^{j}}$

Observe that, for $p=0,1,...,J$,

$\displaystyle\left(\tilde{\Omega}_{2^{j}}^{T}\text{{z}}_{2^{j}}\right)_{1}$
$\displaystyle=0\quad,\quad\left\|\left(\tilde{\Omega}_{2^{j}}^{T}\text{{z}}_{2^{j}}\right)_{2^{p}}^{*}\right\|_{\ell^{2}}$
$\displaystyle=\left\|H_{2^{p}}^{T}(z_{1},...,z_{2^{p}})^{T}\right\|_{\ell^{2}}\leq
2^{\frac{p}{2}}\left\|\text{{z}}_{2^{j}}\right\|_{\ell^{2}}$

and therefore, for each $j>J$,

$\displaystyle\left\|\tilde{\Omega}_{2^{j}}^{T}\text{{z}}_{2^{j}}\right\|_{\ell^{2}}$
$\displaystyle=\sqrt{\sum_{p=0}^{j}\left\|\left(\tilde{\Omega}_{2^{j}}^{T}\text{{z}}_{2^{j}}\right)_{2^{p}}^{*}\right\|_{\ell^{2}}^{2}}\leq
2^{\frac{j+1}{2}}\left\|\text{{z}}_{2^{j}}\right\|_{\ell^{2}}\leq
2^{j+\frac{1}{2}}\left\|\text{{y}}_{2^{j}}^{*}\right\|_{\ell^{2}}.$

Hence, we can bound the norm of
$\left(\tilde{\Omega}^{T}P^{T}\right)_{M,\infty}\text{{y}}_{\infty}$ with

$\displaystyle\left\|\left(\tilde{\Omega}^{T}P^{T}\right)_{M,\infty}\text{{y}}_{\infty}\right\|_{\ell^{2}}$
$\displaystyle\leq\sum_{j=J+1}^{\infty}\frac{1}{2^{\frac{3j}{2}+2}}\left\|\left(\tilde{\Omega}_{2^{j}}^{T}\right)_{1:M,*}\text{{z}}_{2^{j}}\right\|_{\ell^{2}}$
$\displaystyle\leq\sum_{j=J+1}^{\infty}\frac{1}{2^{\frac{j+3}{2}}}\left\|\text{{y}}_{2^{j}}^{*}\right\|_{\ell^{2}}\leq\left(1+\sqrt{2}\right)\frac{\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{\frac{J+3}{2}}}$

  3. iii)

From the block structure of (54), and using (2) and (58),

$\displaystyle\left\|\Gamma_{\infty,M}^{T}(\bar{\text{{a}}}_{\infty}\odot\Pi_{\infty}P^{T}\text{{y}})\right\|_{\ell^{2}}$
$\displaystyle=\left\|\sum_{q=J+1}^{\infty}(H_{2^{q}})_{1:M,*}\left(\bar{\text{{a}}}^{*}_{2^{q}}\odot(P^{T}\text{{y}})^{*}_{2^{q}}\right)\right\|_{\ell^{2}}$
$\displaystyle\leq\sum_{q=J+1}^{\infty}\left\|(H_{2^{q}})\right\|\left(\max_{1\leq
r\leq
2^{q}}\left|\left(\bar{\text{{a}}}^{*}_{2^{q}}\right)_{r}\right|\right)\left\|(P^{T}\text{{y}})^{*}_{2^{q}}\right\|_{\ell^{2}}$
$\displaystyle\leq\sum_{q=J+1}^{\infty}2^{\frac{q}{2}}\frac{\left\|\bar{\text{{c}}}\right\|_{\ell^{2}}}{2^{q+2}}\frac{\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{q+1}}=\frac{\sqrt{2}}{4-\sqrt{2}}\frac{\left\|\bar{\text{{c}}}\right\|_{\ell^{2}}\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{\frac{3J}{2}+3}}$

  4. iv)

Applying the finite-infinite decomposition to
$A_{M}\Pi_{M}\mathcal{H}(\text{{x}}^{T}P\,\Omega(t)P^{T}\text{{y}})$, we have

$\displaystyle\Pi_{M}\mathcal{H}(\text{{x}}^{T}P\,\Omega(t)P^{T}\text{{y}})$
$\displaystyle=\big{(}\text{{x}}_{M}^{T}P_{M}\tilde{\Omega}_{M}\big{)}^{T}\odot\big{(}P_{M}^{T}\text{{y}}_{M}\big{)}+\left(\text{{x}}_{M}^{T}P_{M}\tilde{\Omega}_{M}\right)^{T}\odot\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right)$
$\displaystyle\hskip
20.00003pt+\big{(}\text{{x}}_{\infty}^{T}(P\,\tilde{\Omega})_{M,\infty}\big{)}^{T}\odot\left(P_{M}^{T}\text{{y}}_{M}\right)+\big{(}\text{{x}}_{\infty}^{T}(P\,\tilde{\Omega})_{M,\infty}\big{)}^{T}\odot\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right).$

Using (59), we can bound the first term with Lemma .5:

$\displaystyle\left\|A_{M}\left(\left(\text{{x}}_{M}^{T}P_{M}\tilde{\Omega}_{M}\right)^{T}\odot\left(P_{M}^{T}\text{{y}}_{M}\right)\right)\right\|_{\ell^{2}}\leq\left\|A_{M}\,\text{diag}\left(\|(P_{M}^{T})_{i,*}\|_{\ell^{2}}\right)\,\tilde{\Omega}_{M}^{T}P_{M}^{T}\right\|\
\|\text{{x}}\|_{\ell^{2}}\|\text{{y}}\|_{\ell^{2}}.$

For the second term, using the bound from (i)),

$\displaystyle\left|\left(A_{M}\left(\left(\text{{x}}_{M}^{T}P_{M}\tilde{\Omega}_{M}\right)^{T}\odot\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right)\right)\right)_{i}\right|$
$\displaystyle=\left|\sum_{l=1}^{M}\left(A_{M}\right)_{i,l}\left(\text{{x}}_{M}^{T}P_{M}\tilde{\Omega}_{M}\right)_{l}\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right)_{l}\right|$
$\displaystyle\leq\left|\sum_{l=1}^{M}\left(A_{M}\right)_{i,l}\left(\text{{x}}_{M}^{T}P_{M}\tilde{\Omega}_{M}\right)_{l}\right|\frac{\|\text{{y}}\|_{\ell^{2}}}{2^{J+2}}=\left|\left(A_{M}\tilde{\Omega}_{M}^{T}P_{M}^{T}\text{{x}}_{M}\right)_{i}\right|\frac{\|\text{{y}}\|_{\ell^{2}}}{2^{J+2}}$

and thus

$\left\|A_{M}\left(\left(\text{{x}}_{M}^{T}P_{M}\tilde{\Omega}_{M}\right)^{T}\odot\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right)\right)\right\|_{\ell^{2}}\leq\frac{\left\|A_{M}\tilde{\Omega}_{M}^{T}P_{M}^{T}\right\|}{2^{J+2}}\|\text{{x}}\|_{\ell^{2}}\|\text{{y}}\|_{\ell^{2}}.$

For the third term, with a similar method to the previous term and the
estimate for $\text{{x}}_{\infty}^{T}(P\,\tilde{\Omega})_{M,\infty}$,

$\displaystyle\left|\left(A_{M}\left(\left(\text{{x}}_{\infty}^{T}(P\,\tilde{\Omega})_{M,\infty}\right)^{T}\odot\left(P_{M}^{T}\text{{y}}_{M}\right)\right)\right)_{i}\right|$
$\displaystyle=\left|\sum_{q=1}^{M}(A_{M})_{i,q}\left(P_{M}^{T}\text{{y}}_{M}\right)_{q}\left(\text{{x}}_{\infty}^{T}(P\,\tilde{\Omega})_{M,\infty}\right)_{q}\right|$
$\displaystyle\leq\left(\max_{1\leq q\leq
M}\left|\left(\text{{x}}_{\infty}^{T}(P\,\tilde{\Omega})_{M,\infty}\right)_{q}\right|\right)\left|\sum_{q=1}^{M}(A_{M})_{i,q}\left(P_{M}^{T}\text{{y}}_{M}\right)_{q}\right|$
$\displaystyle\leq\left\|\text{{x}}_{\infty}^{T}(P\,\tilde{\Omega})_{M,\infty}\right\|_{\ell^{2}}\left|\left(A_{M}P_{M}^{T}\text{{y}}_{M}\right)_{i}\right|\leq\left(1+\sqrt{2}\right)\frac{\left\|\text{{x}}\right\|_{\ell^{2}}}{2^{\frac{J+3}{2}}}\left|\left(A_{M}P_{M}^{T}\text{{y}}_{M}\right)_{i}\right|$

and thus

$\displaystyle\left\|A_{M}\left(\left(\text{{x}}_{\infty}^{T}(P\,\tilde{\Omega})_{M,\infty}\right)^{T}\odot\left(P_{M}^{T}\text{{y}}_{M}\right)\right)\right\|_{\ell^{2}}$
$\displaystyle\leq\left(1+\sqrt{2}\right)\frac{\left\|\text{{x}}\right\|_{\ell^{2}}}{2^{\frac{J+3}{2}}}\left\|A_{M}P_{M}^{T}\text{{y}}_{M}\right\|_{\ell^{2}}$
$\displaystyle\leq\left(1+\sqrt{2}\right)\left\|A_{M}P_{M}^{T}\right\|\frac{\left\|\text{{x}}\right\|_{\ell^{2}}\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{\frac{J+3}{2}}}$

For the fourth term, with the estimates for
$\text{{x}}_{\infty}^{T}\left((P\,\tilde{\Omega})_{M,\infty}\right)$ and
$P_{\infty,M}^{T}\text{{y}}_{\infty}$, we obtain

$\displaystyle\left\|A_{M}\left(\left(\text{{x}}_{\infty}^{T}(P\,\tilde{\Omega})_{M,\infty}\right)^{T}\odot\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right)\right)\right\|_{\ell^{2}}$
$\displaystyle=\left\|A_{M}\,\text{diag}\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right)\left(\text{{x}}_{\infty}^{T}(P\,\tilde{\Omega})_{M,\infty}\right)^{T}\right\|_{\ell^{2}}$
$\displaystyle\leq\left\|A_{M}\right\|\left\|P_{\infty,M}^{T}\text{{y}}_{\infty}\right\|_{\ell^{2}}\left\|\text{{x}}_{\infty}^{T}(P\,\tilde{\Omega})_{M,\infty}\right\|_{\ell^{2}}\leq\left\|A_{M}\right\|\left(1+\sqrt{2}\right)\frac{\left\|\text{{x}}\right\|_{\ell^{2}}\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{\frac{3J+7}{2}}}$

Finally, adding the four bounds, we obtain the desired bound.

  5. v)

Applying the finite-infinite decomposition to both $P^{T}$ and $\Gamma^{T}$,
we have

$\displaystyle
A_{M}\Pi_{M}\mathcal{H}(\text{{x}}^{T}P\Theta(t)P^{T}\text{{y}})$
$\displaystyle=A_{M}\Gamma_{M}^{T}\Big{[}\left(P_{M}^{T}\text{{x}}_{M}\right)\odot\left(P_{M}^{T}\text{{y}}_{M}\right)+\left(P_{M}^{T}\text{{x}}_{M}\right)\odot\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right)$
$\displaystyle\hskip
30.00005pt+\left(P_{\infty,M}^{T}\text{{x}}_{\infty}\right)\odot\left(P_{M}^{T}\text{{y}}_{M}\right)+\left(P_{\infty,M}^{T}\text{{x}}_{\infty}\right)\odot\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right)\Big{]}$
$\displaystyle\hskip
30.00005pt+A_{M}\Gamma_{\infty,M}^{T}\Pi_{\infty}\left(P^{T}\text{{x}}\odot
P^{T}\text{{y}}\right).$

We estimate it term by term again. For the first one, using Lemma .5,

$\displaystyle\left\|A_{M}\Gamma_{M}^{T}\left(\left(P_{M}^{T}\text{{x}}_{M}\right)\odot\left(P_{M}^{T}\text{{y}}_{M}\right)\right)\right\|_{\ell^{2}}$
$\displaystyle\leq\left\|A_{M}\Gamma_{M}^{T}\text{diag}(\|(P_{M}^{T})_{i,*}\|_{\ell^{2}})P_{M}^{T}\right\|\|\text{{x}}\|_{\ell^{2}}\|\text{{y}}\|_{\ell^{2}}.$

For the second term, we can use (i)) to bound it element-wise by

$\displaystyle\left|\Big{(}A_{M}\Gamma_{M}^{T}\left(\left(P_{M}^{T}\text{{x}}_{M}\right)\odot\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right)\right)\Big{)}_{i}\right|$
$\displaystyle=\left|\sum_{l=1}^{M}(A_{M}\Gamma_{M}^{T})_{i,l}\left(P_{M}^{T}\text{{x}}_{M}\right)_{l}\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right)_{l}\right|$
$\displaystyle\leq\left|\sum_{l=1}^{M}(A_{M}\Gamma_{M}^{T})_{i,l}\left(P_{M}^{T}\text{{x}}_{M}\right)_{l}\right|\frac{1}{2^{J+2}}\left\|\text{{y}}\right\|_{\ell^{2}}=\frac{\left|(A_{M}\Gamma_{M}^{T}P_{M}^{T}\text{{x}}_{M})_{i}\right|}{2^{J+2}}\left\|\text{{y}}\right\|_{\ell^{2}}$

and thus

$\displaystyle\left\|A_{M}\Gamma_{M}^{T}\left(\left(P_{M}^{T}\text{{x}}_{M}\right)\odot\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right)\right)\right\|_{\ell^{2}}$
$\displaystyle=\left(\sum_{i=1}^{M}\left|\Big{(}A_{M}\Gamma_{M}^{T}\left(\left(P_{M}^{T}\text{{x}}_{M}\right)\odot\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right)\right)\Big{)}_{i}\right|^{2}\right)^{\frac{1}{2}}$
$\displaystyle\leq\left(\sum_{i=1}^{M}\left|(A_{M}\Gamma_{M}^{T}P_{M}^{T}\text{{x}}_{M})_{i}\right|^{2}\right)^{\frac{1}{2}}\frac{\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{J+2}}=\left\|A_{M}\Gamma_{M}^{T}P_{M}^{T}\text{{x}}_{M}\right\|_{\ell^{2}}\frac{\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{J+2}}$
$\displaystyle\leq\frac{\left\|A_{M}\Gamma_{M}^{T}P_{M}^{T}\right\|}{2^{J+2}}\left\|\text{{x}}\right\|_{\ell^{2}}\left\|\text{{y}}\right\|_{\ell^{2}}$

The same procedure applied to the third term yields the same bound as above.
For the fourth term, using previous estimates,

$\displaystyle\left\|A_{M}\Gamma_{M}^{T}\left(\left(P_{\infty,M}^{T}\text{{x}}_{\infty}\right)\odot\left(P_{\infty,M}^{T}\text{{y}}_{\infty}\right)\right)\right\|_{\ell^{2}}$
$\displaystyle\leq\left\|A_{M}\Gamma_{M}^{T}\right\|\left\|P_{\infty,M}^{T}\text{{x}}_{\infty}\right\|_{\ell^{2}}\left\|P_{\infty,M}^{T}\text{{y}}_{\infty}\right\|_{\ell^{2}}$
$\displaystyle\leq\left\|A_{M}\Gamma_{M}^{T}\right\|\frac{\left\|\text{{x}}\right\|_{\ell^{2}}\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{2J+4}}.$

For the fifth term, using the block structure of $\Gamma^{T}$ in (54) and
previous estimates,

$\displaystyle\left\|A_{M}\Gamma_{\infty,M}^{T}\Pi_{\infty}\left(P^{T}\text{{x}}\odot
P^{T}\text{{y}}\right)\right\|_{\ell^{2}}$
$\displaystyle=\left\|A_{M}\sum_{p=J+1}^{\infty}(H_{2^{p}})_{1:M,*}\left(\text{diag}\left(P^{T}\text{{x}}\right)_{2^{p}}^{*}\right)\left(P^{T}\text{{y}}\right)_{2^{p}}^{*}\right\|_{\ell^{2}}$
$\displaystyle\leq\|A_{M}\|\sum_{p=J+1}^{\infty}\|H_{2^{p}}\|\left\|\left(P^{T}\text{{x}}\right)_{2^{p}}^{*}\right\|_{\ell^{2}}\left\|\left(P^{T}\text{{y}}\right)_{2^{p}}^{*}\right\|_{\ell^{2}}$
$\displaystyle\leq\|A_{M}\|\sum_{p=J+1}^{\infty}2^{\frac{p}{2}}\frac{\left\|\text{{x}}\right\|_{\ell^{2}}}{2^{p+2}}\frac{\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{p+2}}=\frac{\sqrt{2}\,\|A_{M}\|}{4-\sqrt{2}}\,\frac{\left\|\text{{x}}\right\|_{\ell^{2}}\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{\frac{3J}{2}+4}}$

Hence, adding the estimates, the result follows.

####  B.4 Proof of Proposition 4.4

  1. i)

To estimate
$\mathcal{H}(\bar{\text{{c}}}^{T}\text{{w}}(t)\text{{w}}^{T}\text{{y}})$, we
again decompose in three parts and estimate each one separately:

$\mathcal{H}(\bar{\text{{c}}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{y}})=\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)P^{T}\text{{y}})+\mathcal{H}(\bar{\text{{a}}}^{T}\Omega^{T}(t)P^{T}\text{{y}})+\mathcal{H}(\bar{\text{{a}}}^{T}\Theta(t)P^{T}\text{{y}}).$

     * •

$\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)P^{T}\text{{y}})$: Since
$\bar{\text{{c}}}\in\mathbb{R}^{M}$,
$\left(\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)P^{T}\text{{y}})\right)_{i}=(\bar{\text{{c}}}^{T}P_{2^{j}}H_{2^{j}})_{i}(P^{T}\text{{y}})_{i}$.
Therefore, applying the block decomposition to
$\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)P^{T}\text{{y}})$ and some
previously calculated estimates from (2) and (58),

$\displaystyle\left\|\left(\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)P^{T}\text{{y}})\right)_{2^{j}}^{*}\right\|_{\ell^{2}}$
$\displaystyle=\left\|\left(H_{2^{j}}^{T}P_{2^{j}}^{T}\bar{\text{{c}}}\right)_{2^{j}}^{*}\odot\left(P^{T}\text{{y}}\right)_{2^{j}}^{*}\right\|_{\ell^{2}}\leq\left\|\text{diag}\left((P^{T}\text{{y}})_{2^{j}}^{*}\right)\right\|\,\left\|\left(H_{2^{j}}^{T}P_{2^{j}}^{T}\bar{\text{{c}}}\right)_{2^{j}}^{*}\right\|_{\ell^{2}}$
$\displaystyle\hskip 20.00003pt=\max_{2^{j}+1\leq i\leq
2^{j+1}}\left|\left(P^{T}\text{{y}}\right)_{i}\right|\,\left\|\left(H_{2^{j}}^{T}P_{2^{j}}^{T}\bar{\text{{c}}}\right)_{2^{j}}^{*}\right\|_{\ell^{2}}\leq\frac{\|\bar{\text{{c}}}\|_{\ell^{2}}\|\text{{y}}\|_{\ell^{2}}}{2^{\frac{3j}{2}+3}}$

and thus

$\left\|\Pi_{\infty}\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)P^{T}\text{{y}})\right\|=\sqrt{\sum_{j=J+1}^{\infty}\left\|\left(\mathcal{H}(\bar{\text{{a}}}^{T}\Omega(t)P^{T}\text{{y}})\right)_{2^{j}}^{*}\right\|^{2}}\leq\frac{1}{\sqrt{7}}\frac{\|\bar{\text{{c}}}\|_{\ell^{2}}\|\text{{y}}\|_{\ell^{2}}}{2^{\frac{3J}{2}+3}}.$

     * •

$\mathcal{H}(\bar{\text{{a}}}^{T}\Omega^{T}(t)P^{T}\text{{y}})$: Since
$\bar{\text{{a}}}^{T}\Omega^{T}(t)P^{T}\text{{y}}=(\bar{\text{{a}}}^{T}\Omega^{T}(t)P^{T}\text{{y}})^{T}=(P^{T}\text{{y}})^{T}\Omega(t)\bar{\text{{a}}}$,
then using the estimate from (57),

$\displaystyle\left\|\left(\mathcal{H}(\bar{\text{{a}}}^{T}\Omega^{T}(t)P^{T}\text{{y}})\right)_{2^{j}}^{*}\right\|_{\ell^{2}}$
$\displaystyle=\left\|\text{diag}\left((P^{T}\bar{\text{{c}}})_{2^{j}}^{*}\right)\big{(}\tilde{\Omega}^{T}P^{T}\text{{y}}\big{)}_{2^{j}}^{*}\right\|_{\ell^{2}}$
$\displaystyle\leq\max_{2^{j}+1\leq i\leq
2^{j+1}}\left|(P^{T}\bar{\text{{c}}})_{i}\right|\left(1+\frac{\sqrt{2}}{8}\right)\frac{\left\|\text{{{y}}}\right\|_{\ell^{2}}}{2^{\frac{j}{2}}}\leq\left(1+\frac{\sqrt{2}}{8}\right)\frac{\|\bar{\text{{c}}}\|_{\ell^{2}}\|\text{{y}}\|_{\ell^{2}}}{2^{\frac{3j}{2}+2}},$

and the bound for
$\Pi_{\infty}\mathcal{H}(\bar{\text{{a}}}^{T}\Omega^{T}(t)P^{T}\text{{y}})$
becomes

$\left\|\Pi_{\infty}\mathcal{H}(\bar{\text{{a}}}^{T}\Omega^{T}(t)P^{T}\text{{y}})\right\|_{\ell^{2}}\leq\frac{1}{\sqrt{7}}\left(1+\frac{\sqrt{2}}{8}\right)\frac{\|\bar{\text{{c}}}\|_{\ell^{2}}\|\text{{y}}\|_{\ell^{2}}}{2^{\frac{3J}{2}+2}}.$

     * •

$\mathcal{H}(\bar{\text{{a}}}^{T}\Theta(t)P^{T}\text{{y}})$: Using the block
structure of $\Gamma^{T}$ from (54) and the estimates (58) and (2), we have

$\displaystyle\left\|\left(\Gamma^{T}\left(\bar{\text{{a}}}\odot
P^{T}\text{{y}}\right)\right)_{2^{j}}^{*}\right\|_{\ell^{2}}$
$\displaystyle=\left\|\sum_{p=j+1}^{\infty}(H_{2^{p}})_{2^{j}+1:2^{j+1},*}\left(\bar{\text{{a}}}_{2^{p}}^{*}\odot(P^{T}\text{{y}})_{2^{p}}^{*}\right)\right\|_{\ell^{2}}$
$\displaystyle\leq\sum_{p=j+1}^{\infty}\left\|(H_{2^{p}})_{2^{j}+1:2^{j+1},*}\right\|\left\|\text{diag}(\bar{\text{{a}}}_{2^{p}}^{*})\right\|\left\|(P^{T}\text{{y}})_{2^{p}}^{*}\right\|_{\ell^{2}}$
$\displaystyle\leq\sum_{p=j+1}^{\infty}\left\|H_{2^{p}}\right\|\max_{2^{p}+1\leq
i\leq
2^{p+1}}\left|(P^{T}\bar{\text{{c}}})_{i}\right|\left\|(P^{T}\text{{y}})_{2^{p}}^{*}\right\|_{\ell^{2}}$
$\displaystyle\leq\sum_{p=j+1}^{\infty}2^{\frac{p}{2}}\frac{\left\|\bar{\text{{c}}}\right\|_{\ell^{2}}}{2^{p+2}}\frac{\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{p+1}}=\frac{\sqrt{2}}{4-\sqrt{2}}\frac{\left\|\bar{\text{{c}}}\right\|_{\ell^{2}}\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{\frac{3j}{2}+3}},$

and thus

$\displaystyle\left\|\Pi_{\infty}\mathcal{H}(\bar{\text{{a}}}^{T}\Theta(t)P^{T}\text{{y}})\right\|_{\ell^{2}}\leq\frac{\sqrt{2}}{\sqrt{7}(4-\sqrt{2})}\frac{\left\|\bar{\text{{c}}}\right\|_{\ell^{2}}\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{\frac{3J}{2}+3}}.$

The result follows from summing the estimates.

  2. ii)

We proceed in the same way as we did for
$\Pi_{\infty}\mathcal{H}(\bar{\text{{c}}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{y}})$,
that is, separate it in three parts

$\mathcal{H}(\text{{x}}^{T}\text{{w}}(t)\text{{w}}^{T}(t)\text{{y}})=\mathcal{H}(\text{{x}}^{T}\Omega(t)P^{T}\text{{y}})+\mathcal{H}(\text{{x}}^{T}\Omega^{T}(t)P^{T}\text{{y}})+\mathcal{H}(\text{{x}}^{T}\Theta(t)P^{T}\text{{y}})$

and estimate each one.

     * •

$\mathcal{H}(\text{{x}}^{T}\Omega(t)P^{T}\text{{y}})$: Applying the block
decomposition to the expression of
$\mathcal{H}(\text{{x}}^{T}\Omega(t)P^{T}\text{{y}})$ and the estimates from
(2) and (58),

$\displaystyle\left\|\left(\mathcal{H}(\text{{x}}^{T}\Omega(t)P^{T}\text{{y}})\right)_{2^{j}}^{*}\right\|_{\ell^{2}}$
$\displaystyle=\left\|(\tilde{\Omega}^{T}P^{T}\text{{x}})_{2^{j}}^{*}\odot(P^{T}\text{{y}})_{2^{j}}^{*}\right\|_{\ell^{2}}\leq\left\|\text{diag}\left(P^{T}\text{{y}}\right)_{2^{j}}^{*}\right\|\left\|(\tilde{\Omega}^{T}P^{T}\text{{x}})_{2^{j}}^{*}\right\|_{\ell^{2}}$
$\displaystyle\leq\frac{\|\text{{y}}\|_{\ell^{2}}}{2^{j+1}}\left(3+2\sqrt{2}\right)\frac{\|\text{{x}}\|_{\ell^{2}}}{2^{\frac{j}{2}+2}}=\frac{\left(3+2\sqrt{2}\right)}{2^{\frac{3j}{2}+3}}\|\text{{x}}\|_{\ell^{2}}\|\text{{y}}\|_{\ell^{2}}$

and then

$\left\|\Pi_{\infty}\mathcal{H}(\text{{x}}^{T}\Omega(t)P^{T}\text{{y}})\right\|_{\ell^{2}}\leq\frac{\left(3+2\sqrt{2}\right)}{\sqrt{7}}\frac{\|\text{{x}}\|_{\ell^{2}}\|\text{{y}}\|_{\ell^{2}}}{2^{\frac{3J}{2}+3}}$

     * •

$\mathcal{H}(\text{{x}}^{T}\Omega^{T}(t)P^{T}\text{{y}})$: Since
$\text{{x}}^{T}\Omega^{T}(t)P^{T}\text{{y}}=(\text{{x}}^{T}\Omega^{T}(t)P^{T}\text{{y}})^{T}=\text{{y}}^{T}\Omega^{T}(t)P^{T}\text{{x}}$,
repeating the same process for the previous item, we also have

$\left\|\Pi_{\infty}\mathcal{H}(\text{{x}}^{T}\Omega^{T}(t)P^{T}\text{{y}})\right\|_{\ell^{2}}\leq\frac{\left(3+2\sqrt{2}\right)}{\sqrt{7}}\frac{\|\text{{x}}\|_{\ell^{2}}\|\text{{y}}\|_{\ell^{2}}}{2^{\frac{3J}{2}+3}}$

     * •

$\mathcal{H}(\text{{x}}^{T}\Theta(t)P^{T}\text{{y}})$: Using the block
decomposition and the block structure of $\Gamma$,

$\displaystyle\left\|\left(\Gamma^{T}\left(P^{T}\text{{x}}\odot
P^{T}\text{{y}}\right)\right)_{2^{j}}^{*}\right\|_{\ell^{2}}$
$\displaystyle=\left\|\sum_{p=j+1}^{\infty}(H_{2^{p}})_{2^{j}+1:2^{j+1},*}\left((P^{T}\text{{x}})_{2^{p}}^{*}\odot(P^{T}\text{{y}})_{2^{p}}^{*}\right)\right\|_{\ell^{2}}$
$\displaystyle\leq\sum_{p=j+1}^{\infty}\left\|(H_{2^{p}})\right\|\,\left\|\text{diag}((P^{T}\text{{x}})_{2^{p}}^{*})\right\|\left\|(P^{T}\text{{y}})_{2^{p}}^{*}\right\|_{\ell^{2}}$
$\displaystyle\leq\sum_{p=j+1}^{\infty}2^{\frac{p}{2}}\frac{\left\|\text{{x}}\right\|_{\ell^{2}}}{2^{p+1}}\frac{\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{p+1}}=\frac{\sqrt{2}}{4-\sqrt{2}}\frac{\left\|\text{{x}}\right\|_{\ell^{2}}\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{\frac{3j}{2}+2}}$

and thus

$\left\|\Pi_{\infty}\mathcal{H}(\text{{x}}^{T}\Theta(t)P^{T}\text{{y}})\right\|_{\ell^{2}}\leq\frac{\sqrt{2}}{\sqrt{7}\left(4-\sqrt{2}\right)}\frac{\left\|\text{{x}}\right\|_{\ell^{2}}\left\|\text{{y}}\right\|_{\ell^{2}}}{2^{\frac{3J}{2}+2}}\\\
$

Putting together the estimates, we have the desired result.

## References

  * [1] I. Aziz and Siraj ul Islam.  New algorithms for the numerical solution of nonlinear Fredholm and Volterra integral equations using Haar wavelets.  Journal of Computational and Applied Mathematics, 239:333–345, 2013\.
  * [2] M. Breden and C. Kuehn.  Rigorous validation of stochastic transition paths.  Journal des Mathematiques Pures et Appliquees, 131:88–129, 2019\.
  * [3] C. F. Chen and C. H. Hsiao.  Haar wavelet method for solving lumped and distributed-parameter systems.  IEE Proceedings - Control Theory and Applications, 144(1):87–94, 1997.
  * [4] S. Day, Y. Hiraoka, K. Mischaikow, and T. Ogawa.  Rigorous numerics for global dynamics: a study of the Swift-Hohenberg equation.  SIAM Journal on Applied Dynamical Systems, 4(1):1–31, 2005.
  * [5] S. Day, J. Lessard, and K. Mischaikow.  Validated continuation for equilibria of PDEs.  SIAM Journal on Numerical Analysis, 45(4):1398–1424, 2007.
  * [6] J.-L. Figueras, M. Gameiro, J.-P. Lessard, and R. de la Llave.  A framework for the numerical computation and a posteriori verification of invariant objects of evolution equations.  SIAM Journal on Applied Dynamical Systems, 16(2):1070–1088, 2017\.
  * [7] M. Gameiro, J.-P. Lessard, and K. Mischaikow.  Validated continuation over large parameter ranges for equilibria of PDEs.  Mathematics and Computers in Simulation, 79(4):1368–1382, 2008\.
  * [8] The CAPD Group.  Coomputer assisted proofs in dynamics group.  http://capd.ii.uj.edu.pl.  Accessed: April 13, 2022.
  * [9] K. Hashimoto, T. Kinoshita, and M. T Nakao.  Numerical verification of solutions for nonlinear parabolic problems.  Numerical Functional Analysis and Optimization, 41(12):1495–1514, 2020.
  * [10] Eugenio Hernández and Guido Weiss.  A First Course on Wavelets.  Studies in Advanced Mathematics. CRC Press, 2000 N.W. Corporate Blvd., Boca Raton, Florida 33431, 1 edition, 1996.
  * [11] A. Hungria, J.-P. Lessard, and J. D. Mireles James.  Rigorous numerics for analytic solutions of differential equations: The radii polynomial approach.  Mathematics of Computation, 85(299):1427–1459, 2016.
  * [12] Ü. Lepik.  Solving fractional integral equations by the Haar wavelet method.  Applied Mathematics and Computation, 214(2):468–478, 2009.
  * [13] J.-P. Lessard and C. Reinhardt.  Rigorous numerics for nonlinear differential equations using Chebyshev series.  SIAM Journal on Numerical Analysis, 52(1):1–22, 2014.
  * [14] S. Maier-paape, U. Miller, K. Mischaikow, and T. Wanner.  Rigorous numerics for the Cahn-Hilliard equation on the unit square.  Revista Matemática Complutense, 21(2):351–426, 2008.
  * [15] J. Majak, B. Shvartsman, K. Karjust, M. Mikola, A. Haavajõe, and M. Pohlak.  On the accuracy of the Haar wavelet discretization method.  Composites Part B, 80:321–327, 2015.
  * [16] V. Mehandiratta, M. Mehra, and G. Leugering.  An approach based on Haar wavelet for the approximation of fractional calculus with application to initial and boundary value problems.  Special issue paper of Mathematical Methods in Applied Sciences, pages 1–20, 2020.
  * [17] T. Miyaji, P. Pilarczyk, M. Gameiro, H. Kokubu, and K. Mischaikow.  A study of rigorous ODE integrators for multi-scale set-oriented computations.  Applied Numerical Mathematics, 107:34–47, 2016.
  * [18] M. T. Nakao.  A computational verification method of existence of solutions for nonlinear elliptic equations.  North-Holland Mathematics Studies, 160:101–120, 1989.
  * [19] M. T. Nakao, T. Kimura, and T. Kinoshita.  Constructive a priori error estimates for a full discrete approximation of the heat equation.  SIAM Journal on Numerical Analysis, 51(3):1525–1541, 2013.
  * [20] M. T. Nakao, M Plum, and Y. Watanabe.  Numerical Verification Methods and Computer-Assisted Proofs for Partial Differential Equations, volume 53 of Springer Series in Computational Mathematics.  Springer Nature Singapore Pte Ltd., 152 Beach Road, #21-01/04 Gateway East, Singapore 189721, 1 edition, 2019.
  * [21] Ö. Oruç.  A non-uniform Haar wavelet method for numerically solving two-dimensional convection-dominated equations and two-dimensional near singular elliptic equations.  Computers & Mathematics with Applications, 77(7):1799–1820, 2019\.
  * [22] C. Reinhardt and J. D. Mireles James.  Fourier–Taylor parameterization of unstable manifolds for parabolic partial differential equations: Formalism, implementation and rigorous validation.  Indagationes Mathematicae, 30(1):39–80, 2019.
  * [23] S.M. Rump.  INTLAB - INTerval LABoratory.  In Tibor Csendes, editor, Developments in Reliable Computing, pages 77–104. Kluwer Academic Publishers, Dordrecht, 1999.  http://www.ti3.tuhh.de/rump/.
  * [24] Siraj ul Islam, I. Aziz, and M. Ahmad.  Numerical solution of two-dimensional elliptic PDEs with nonlocal boundary conditions.  Computers & Mathematics with Applications, 69(1):180–205, 2015\.
  * [25] J. B. Van Den Berg, C. M. Groothedde, and J. F. Williams.  Rigorous computation of a radially symmetric localized solution in a Ginzburg-Landau problem.  SIAM Journal on Applied Dynamical Systems, 14(1):423–447, 2015\.
  * [26] J. B. van den Berg and R. Sheombarsing.  Validated computations for connecting orbits in polynomial vector fields.  Indagationes Mathematicae, 31(2):310–373, 2020.
  * [27] J. B. van der Berg, J.-P. Lessard, and K. Mischaikow.  Global smooth solution curves using rigorous branch following.  Mathematics of Computation, 79(271):1565–1584, 2010.
  * [28] P. Zgliczyński.  Attracting fixed points for the Kuramoto-Sivashinsky equation: A computer assisted proof.  SIAM Journal on Applied Dynamical Systems, 1(2):215–235, 2002.
  * [29] P. Zgliczyński.  $C^{1}$ Lohner algorithm.  Foundations of Computational Mathematics, 2(4):429–465, 2002.
  * [30] P. Zgliczyński and K. Mischaikow.  Rigorous numerics for partial differential equations: The Kuramoto-Sivashinsky equation.  Foundations of Computational Mathematics, 1(3):255–288, 2001.