Continuation of homoclinic orbits in the suspension bridge equation: a computer-assisted proof

TL;DR

This paper proves the existence of symmetric homoclinic orbits in the suspension bridge equation across a range of parameters using a computer-assisted method involving Chebyshev series and contraction theorems.

Contribution

It introduces a computer-assisted proof technique for homoclinic orbits in a nonlinear differential equation, combining the uniform contraction theorem and radii polynomial approach.

Findings

Existence of symmetric homoclinic orbits for all β in [0.5, 1.9]

Effective parameterization of stable manifolds

Validated solutions using a contraction mapping approach

Abstract

In this paper, we prove existence of symmetric homoclinic orbits for the suspension bridge equation $u""+\beta u" + e^u-1=0$ for all parameter values $\beta \in [0.5,1.9]$. For each $\beta$, a parameterization of the stable manifold is computed and the symmetric homoclinic orbits are obtained by solving a projected boundary value problem using Chebyshev series. The proof is computer-assisted and combines the uniform contraction theorem and the radii polynomial approach, which provides an efficient means of determining a set, centered at a numerical approximation of a solution, on which a Newton-like operator is a contraction.

Tables

Table 1. Table 1: Coefficients S i , k ( j ) superscript subscript 𝑆 𝑖 𝑘 𝑗 S_{i,k}^{(j)} for the splitting of F ​ ( β s , x ¯ s ) 𝐹 subscript 𝛽 𝑠 subscript ¯ 𝑥 𝑠 F(\beta_{s},\bar{x}_{s}) as a polynomial in s 𝑠 s , as given in ( 4.24 ). The coefficients S i , k ( j ) superscript subscript 𝑆 𝑖 𝑘 𝑗 S_{i,k}^{(j)} for k = 0 𝑘 0 k=0 and 3 ≤ j ≤ 6 3 𝑗 6 3\leq j\leq 6 are provided in ( 4.21 ). The coefficients S i , k ( j ) superscript subscript 𝑆 𝑖 𝑘 𝑗 S_{i,k}^{(j)} for j = 1 , 2 𝑗 1 2 j=1,2 are provided in ( 4.19 ) and ( 4.20 ).

$i$	Coefficients $S_{i,k}^{\left(3\right)}$ for $k\ge 1$
$0$	$2k{\overline{x}}_{0,k}^{\left(1\right)}+{\overline{L}}_0\left[{\overline{x}}_{0,k\pm 1}^{\left(2\right)}+{\left({\overline{x}}_0^{\left(1\right)}\ast {\overline{x}}_0^{\left(2\right)}\right)}_{k\pm 1}\right]$
$1$	$2k\Delta {\overline{x}}_k^{\left(1\right)}+{\overline{L}}_0\left[\Delta {\overline{x}}_{k\pm 1}^{\left(2\right)}+{\left(\Delta {\overline{x}}^{\left(1\right)}\ast {\overline{x}}_0^{\left(2\right)}\right)}_{k\pm 1}+{\left({\overline{x}}_0^{\left(1\right)}\ast \Delta {\overline{x}}^{\left(2\right)}\right)}_{k\pm 1}\right]$
	$+\Delta \overline{L}\left[{\overline{x}}_{0,k\pm 1}^{\left(2\right)}+{\left({\overline{x}}_0^{\left(1\right)}\ast {\overline{x}}_0^{\left(2\right)}\right)}_{k\pm 1}\right]$
$2$	${\overline{L}}_0{\left(\Delta {\overline{x}}^{\left(1\right)}\ast \Delta {\overline{x}}^{\left(2\right)}\right)}_{k\pm 1}+\Delta \overline{L}\left[\Delta {\overline{x}}_{k\pm 1}^{\left(2\right)}+{\left(\Delta {\overline{x}}^{\left(1\right)}\ast {\overline{x}}_0^{\left(2\right)}\right)}_{k\pm 1}+{\left({\overline{x}}_0^{\left(1\right)}\ast \Delta {\overline{x}}^{\left(2\right)}\right)}_{k\pm 1}\right]$
$3$	$\Delta \overline{L}{\left(\Delta {\overline{x}}^{\left(1\right)}\ast \Delta {\overline{x}}^{\left(2\right)}\right)}_{k\pm 1}$
$i$	Coefficients $S_{i,k}^{\left(4\right)}$ for $k\ge 1$
$0$	$2k{\overline{x}}_{0,k}^{\left(2\right)}+{\overline{L}}_0{\overline{x}}_{0,k\pm 1}^{\left(3\right)}$
$1$	${\overline{L}}_0\Delta {\overline{x}}_{k\pm 1}^{\left(3\right)}+2k\Delta {\overline{x}}_k^{\left(2\right)}+\Delta \overline{L}{\overline{x}}_{0,k\pm 1}^{\left(3\right)}$
$2$	$\Delta \overline{L}\Delta {\overline{x}}_{k\pm 1}^{\left(3\right)}$
$3$	$0$
$i$	Coefficients $S_{i,k}^{\left(5\right)}$ for $k\ge 1$
$0$	$2k{\overline{x}}_{0,k}^{\left(3\right)}+{\overline{L}}_0{\overline{x}}_{0,k\pm 1}^{\left(4\right)}$
$1$	${\overline{L}}_0\Delta {\overline{x}}_{k\pm 1}^{\left(4\right)}+2k\Delta {\overline{x}}_k^{\left(3\right)}+\Delta \overline{L}{\overline{x}}_{0,k\pm 1}^{\left(4\right)}$
$2$	$\Delta \overline{L}\Delta {\overline{x}}_{k\pm 1}^{\left(4\right)}$
$3$	$0$
$i$	Coefficients $S_{i,k}^{\left(6\right)}$ for $k\ge 1$
$0$	$2k{\overline{x}}_{0,k}^{\left(4\right)}-{\overline{L}}_0\left[{\overline{x}}_{0,k\pm 1}^{\left(1\right)}+{\beta }_0{\overline{x}}_{0,k\pm 1}^{\left(3\right)}\right]$
$1$	$2k\Delta {\overline{x}}_k^{\left(4\right)}-\Delta \overline{L}\left[{\overline{x}}_{0,k\pm 1}^{\left(1\right)}+{\beta }_0{\overline{x}}_{0,k\pm 1}^{\left(3\right)}\right]-{\overline{L}}_0\left[\Delta {\overline{x}}_{k\pm 1}^{\left(1\right)}+{\beta }_0\Delta {\overline{x}}_{k\pm 1}^{\left(3\right)}+\Delta \beta {\overline{x}}_{0,k\pm 1}^{\left(3\right)}\right]$
$2$	$-\Delta \overline{L}\left[\Delta {\overline{x}}_{k\pm 1}^{\left(1\right)}+{\beta }_0\Delta {\overline{x}}_{k\pm 1}^{\left(3\right)}+\Delta \beta {\overline{x}}_{0,k\pm 1}^{\left(3\right)}\right]-{\overline{L}}_0\Delta \beta \Delta {\overline{x}}_{k\pm 1}^{\left(3\right)}$
$3$	$-\Delta \overline{L}\Delta \beta \Delta {\overline{x}}_{k\pm 1}^{\left(3\right)}$

Table 2. Table 2: Coefficients z ~ 2 , k ( i ) superscript subscript ~ 𝑧 2 𝑘 𝑖 \tilde{z}_{2,k}^{(i)} and z ~ 3 , k ( i ) superscript subscript ~ 𝑧 3 𝑘 𝑖 \tilde{z}_{3,k}^{(i)} for k ≠ 0 𝑘 0 k\neq 0 .

Coefficients in front of $r^2$, for $k\ge 1$
${\tilde{z}}_{2,k}^{\left(3\right)}$	$v_L\left({\left({\overline{x}}_0^{\left(1\right)}\ast u^{\left(2\right)}\right)}_{k\pm 1}+{\left(u^{\left(1\right)}\ast {\overline{x}}_0^{\left(2\right)}\right)}_{k\pm 1}+s\left[{\left(\Delta {\overline{x}}^{\left(1\right)}\ast u^{\left(2\right)}\right)}_{k\pm 1}+{\left(u^{\left(1\right)}\ast \Delta {\overline{x}}^{\left(2\right)}\right)}_{k\pm 1}\right]+u_{k\pm 1}^{\left(2\right)}\right)$
	$+u_L\left({\left({\overline{x}}_0^{\left(1\right)}\ast v^{\left(2\right)}\right)}_{k\pm 1}+{\left(v^{\left(1\right)}\ast {\overline{x}}_0^{\left(2\right)}\right)}_{k\pm 1}+s\left[{\left(\Delta {\overline{x}}^{\left(1\right)}\ast v^{\left(2\right)}\right)}_{k\pm 1}+{\left(v^{\left(1\right)}\ast \Delta {\overline{x}}^{\left(2\right)}\right)}_{k\pm 1}\right]+v_{k\pm 1}^{\left(2\right)}\right)$
	$+s\Delta \overline{L}\left({\left(u^{\left(1\right)}\ast v^{\left(2\right)}\right)}_{k\pm 1}+{\left(v^{\left(1\right)}\ast u^{\left(2\right)}\right)}_{k\pm 1}\right)+{\overline{L}}_0\left({\left(u^{\left(1\right)}\ast v^{\left(2\right)}\right)}_{k\pm 1}+{\left(v^{\left(1\right)}\ast u^{\left(2\right)}\right)}_{k\pm 1}\right)$
${\tilde{z}}_{2,k}^{\left(4\right)}$	$u_L{v}_{k\pm 1}^{\left(3\right)}+v_L{u}_{k\pm 1}^{\left(3\right)}$
${\tilde{z}}_{2,k}^{\left(5\right)}$	$u_L{v}_{k\pm 1}^{\left(4\right)}+v_L{u}_{k\pm 1}^{\left(4\right)}$
${\tilde{z}}_{2,k}^{\left(6\right)}$	$-u_L\left({\beta }_s{v}_{k\pm 1}^{\left(3\right)}+v_{k\pm 1}^{\left(1\right)}\right)-v_L\left({\beta }_s{u}_{k\pm 1}^{\left(3\right)}+u_{k\pm 1}^{\left(1\right)}\right)$
Coefficients in front of $r^3$, for $k\ge 1$
${\tilde{z}}_{3,k}^{\left(3\right)}$	$u_L{\left(u^{\left(1\right)}\ast v^{\left(2\right)}\right)}_{k\pm 1}+u_L{\left(u^{\left(2\right)}\ast v^{\left(1\right)}\right)}_{k\pm 1}+v_L{\left(u^{\left(1\right)}\ast u^{\left(2\right)}\right)}_{k\pm 1}$
${\tilde{z}}_{3,k}^{\left(4\right)}$
${\tilde{z}}_{3,k}^{\left(5\right)}$	$0$
${\tilde{z}}_{3,k}^{\left(6\right)}$

Equations

$$\frac{{\partial }^{2}U}{\partial T^2}=-\frac{{\partial }^{4}U}{\partial X^4}-e^U+1.$$

$$u^{\prime \prime \prime \prime }+c^2{u}^{\prime \prime }+e^u-1=0.$$

$$\left\{\begin{array}[]{l}u^{\prime\prime\prime\prime}+\beta u^{\prime\prime}+e^{u}-1=0\\ u(-t)=u(t)\\ \lim_{t\to\infty}u(t)=0.\end{array}\right.$$

$$\{\begin{array}{rl}& v_1^{\prime }=v_2+v_1{v}_2\\ & v_2^{\prime }=v_3\\ & v_3^{\prime }=v_4\\ & v_4^{\prime }=-\beta v_3-v_1.\end{array}$$

$$X=X_1\times X_2\times \cdots \times X_d,$$

$$\Vert x{\Vert }_X=\max \left(\Vert x^{(1)}{\Vert }_{X_1},\dots ,\Vert x^{(d)}{\Vert }_{X_d}\right),$$

$$B_r(y)=\{x\in X\mid \Vert x-y{\Vert }_X\le r\}$$

$${\bar{x}}_{\beta }\stackrel{\text{mbox}[}{=}\frac{{\beta }_1-\beta }{{\beta }_1-{\beta }_0}{\bar{x}}_{{\beta }_0}+\frac{\beta -{\beta }_0}{{\beta }_1-{\beta }_0}{\bar{x}}_{{\beta }_1}.$$

$$\underset{\beta \in [{\beta }_0,{\beta }_1]}{\sup }{\Vert T^{(j)}(\beta ,{\bar{x}}_{\beta })-{\bar{x}}_{\beta }^{(j)}\Vert }_{X_j}$$

$$\underset{\begin{array}{c}b,c\in B_r(0)\\ \beta \in [{\beta }_0,{\beta }_1]\end{array}}{\sup }{\Vert D_x{T}^{(j)}(\beta ,{\bar{x}}_{\beta }+b)c\Vert }_{X_j}$$

$$\tilde{T}:\{\begin{array}{rl}\left[{\beta }_0,{\beta }_1\right]\times B_r(0)& ⟶B_r(0)\\ (\beta ,x)& ⟼\tilde{T}(\beta ,x)\stackrel{\text{mbox}[}{=}T(\beta ,x+{\bar{x}}_{\beta })-{\bar{x}}_{\beta }\end{array}$$

$$p_j(r)\stackrel{\text{mbox}[}{=}{Y}^{(j)}+Z^{(j)}(r)-r,j=1,\dots ,d.$$

$$\mathcal{I}\stackrel{\text{mbox}[}{=}\bigcap _{j=1}^d\{r>0\mid p_j(r)<0\},$$

$${\Psi }_{\beta }(v)\stackrel{\text{mbox}[}{=}\left(\begin{array}{c}{v}_2+v_1{v}_2\\ v_3\\ v_4\\ -\beta v_3-v_1\end{array}\right),$$

$$D{\Psi }_{\beta }(0)=\left(\begin{array}{cccc}0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\\ -1& 0& -\beta & 0\end{array}\right),$$

$$\lambda (\beta )=-\frac{1}{2}\sqrt{2-\beta }+\text{i}\frac{1}{2}\sqrt{2+\beta }.$$

$$V(\beta )=\left(\begin{array}{c}1\\ \lambda (\beta )\\ \lambda (\beta {)}^2\\ \lambda (\beta {)}^3\end{array}\right).$$

$$Q_{\beta }(\theta )=\sum _{|\alpha |\ge 0}{a}_{\alpha }(\beta ){\theta }^{\alpha },\theta =\left(\begin{array}{c}{\theta }_1\\ {\theta }_2\end{array}\right)\in {\mathbb{C}}^2,a_{\alpha }(\beta )=\left(\begin{array}{c}{a}_{\alpha }^{(1)}(\beta )\\ a_{\alpha }^{(2)}(\beta )\\ a_{\alpha }^{(3)}(\beta )\\ a_{\alpha }^{(4)}(\beta )\end{array}\right)\in {\mathbb{C}}^4,$$

$$Q_{\beta }(0)=0,D{Q}_{\beta }(0)=\left(\begin{array}{cc}V(\beta )& V^{\ast }(\beta )\end{array}\right),$$

$$D{Q}_{\beta }(\theta )\left(\begin{array}{cc}\lambda (\beta )& 0\\ 0& {\lambda }^{\ast }(\beta )\end{array}\right)\theta ={\Psi }_{\beta }(Q_{\beta }(\theta )).$$

$$P_{\beta }(\theta )\stackrel{\text{mbox}[}{=}{Q}_{\beta }({\theta }_1+\mathrm{i}{\theta }_2,{\theta }_1-\mathrm{i}{\theta }_2),\text{for }\theta \in {\mathbb{R}}^2,$$

$$\sum _{|\alpha |\ge 0}({\alpha }_1\lambda (\beta )+{\alpha }_2{\lambda }^{\ast }(\beta ))a_{\alpha }(\beta ){\theta }^{\alpha }=\sum _{|\alpha |\ge 0}\left(\begin{array}{c}{a}_{\alpha }^{(2)}(\beta )+(a^{(1)}(\beta )\star a^{(2)}(\beta ){)}_{\alpha }\\ a_{\alpha }^{(3)}(\beta )\\ a_{\alpha }^{(4)}(\beta )\\ -a_{\alpha }^{(1)}(\beta )-\beta a_{\alpha }^{(3)}(\beta )\end{array}\right){\theta }^{\alpha },$$

$$(u\star v{)}_{\alpha }=\sum _{0\le \sigma \le \alpha }{u}_{\sigma }{v}_{\alpha -\sigma },\text{for all }\alpha \in {\mathbb{N}}^2,$$

$$F(\beta ,a)=\{\begin{array}{rlrl}& a_{0,0}& & \text{if }\alpha =(0,0),\\ & a_{1,0}-V(\beta )& & \text{if }\alpha =(1,0),\\ & a_{0,1}-V^{\ast }(\beta )& & \text{if }\alpha =(0,1),\\ & ({\alpha }_1\lambda (\beta )+{\alpha }_2{\lambda }^{\ast }(\beta ))a_{\alpha }-\left(\begin{array}{c}{a}_{\alpha }^{(2)}+(a^{(1)}\star a^{(2)}{)}_{\alpha }\\ a_{\alpha }^{(3)}\\ a_{\alpha }^{(4)}\\ -a_{\alpha }^{(1)}-\beta a_{\alpha }^{(3)}\end{array}\right)& & \text{for }|\alpha |\ge 2.\end{array}$$

$${\Vert u\Vert }_{1,\nu }\stackrel{\text{mbox}[}{=}\sum _{|\alpha |=0}^{\infty }|u_{\alpha }|{\nu }^{|\alpha |}<\infty .$$

$$X\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,(\ell^{1}_{\nu})^{4},\quad\text{with the norm }\left\|a\right\|_{X}\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,\max\limits_{j=1,\ldots,4}\bigl{\|}a^{(j)}\bigr{\|}_{1,\nu}.$$

$$\sum_{|\alpha|\geq 0}\bigl{|}a^{(j)}_{\alpha}\bigr{|}\nu^{|\alpha|}<\infty,\quad\text{for }j=1,\ldots,4,$$

$${\left|\theta \right|}_{\infty }\stackrel{\text{mbox}[}{=}\max \left(|{\theta }_1|,|{\theta }_2|\right)\le \nu .$$

$${\beta }_s={\beta }_0+s({\beta }_1-{\beta }_0)={\beta }_0+s\Delta \beta ,\text{for }s\in [0,1].$$

$$F({\beta }_s,a(s))=0,\text{for }s\in [0,1].$$

$$F^{[N]}({\beta }_s,\cdot )\stackrel{\text{mbox}[}{=}{\left(F_{\alpha }({\beta }_s,\cdot )\right)}_{0\le |\alpha |<N}=0,$$

Full text

Key words. Suspension bridge equation, traveling waves, contraction mapping, rigorous numerics, symmetric homoclinic orbits, stable manifolds

1 Introduction

One of the simplest models [15, 12] for a suspension bridge is the partial differential equation (PDE)

[TABLE]

Here $U(T,X)$ describes the deflection of the roadway from the rest state $U=0$ as a function of time $T$ and the spatial variable $X$ (in the direction of traffic). This paper is concerned with traveling wave solutions of (1.1), i.e., solutions $U(T,X)=u(X-cT)$ describing a disturbance with profile $u$ propagating at velocity $c$ along the surface of the bridge. In particular, we apply a computer-assisted proof method to show that there is a large range of velocities for which such a solitary wave exists.

Looking for traveling waves of (1.1) with wave speed $c$ leads to the ordinary differential equation

[TABLE]

For large positive and negative values of the independent variable $t=X-cT$ we assume the solution to converge to the equilibrium $u=0$. Due to the reversibility symmetry of the PDE in both time and space, we may restrict our attention to symmetric solutions. Hence, setting $\beta=c^{2}$, we are looking for symmetric homoclinic orbits satisfying

[TABLE]

Fourth order differential equations of the form $u^{\prime\prime\prime\prime}+\beta u^{\prime\prime}+f(u)=0$ for various nonlinearities $f$ have been studied extensively. For the bistable nonlinearity $f(u)=u^{3}-u$ the equation is a standard model in pattern formation, called the Swift-Hohenberg equation (see [18] and references therein), whereas the quadratic nonlinearity $f(u)=u^{2}-u$ appears, for example, in the study of water waves [4]. For the piecewise linear case $f(u)=\max\{u,0\}$ homoclinic solutions were obtained in [15, 8]. For the problem with the exponential nonlinearity $f(u)=e^{u}-1$ a family of periodic solutions was established in [17].

In [8] the question about existence of a symmetric homoclinic orbit of (1.3) is raised. This question was addressed by variational methods in [21], where the authors proved the result for almost all parameter values $\beta\in(0,2)$. In [20] the existence of homoclinic orbits was demonstrated for all $\beta\in(0,c_{*}^{2})\approx(0,0.5516)$, again using variational methods as well as intricate estimates on the second variation. In a different direction, using a computer-assisted proof, it was proven in [3] that (1.3) has at least 36 homoclinic solutions for the single parameter value $\beta=1.69$.

In the present paper we complement the above results by proving the following.

Theorem 1**.**

For all parameter values $\beta\in[0.5,1.9]$ there exists a symmetric homoclinic orbit of (1.3).

We remark that for $|\beta|<2$ the origin is a saddle-focus, while for $\beta>2$ it is a saddle-center. Furthermore, we note the integral identity $\int_{\mathbb{R}}|u^{\prime\prime}|^{2}-\beta|u^{\prime}|^{2}=-\int_{\mathbb{R}}(e^{u}-1)u$. Since the right hand side is non-positive, homoclinic orbits are excluded for $\beta\leq 0$. It is thus expected that the parameter range for which homoclinics exist is $\beta\in(0,2)$, or, equivalently, wave speeds $c\in(0,\sqrt{2})$. Our method for proving the result in Theorem 1 is computer-assisted. While it can certainly be extended somewhat beyond the interval $[0.5,1.9]$, it is not possible to cover the entire range $(0,2)$ in this way. Indeed, as $\beta$ decreases towards [math] the amplitude of the solution diverges ($u$ becomes very negative), whereas when $\beta$ tends to $2$ the homoclinic orbit collapses onto the trivial solution. In both limit regimes computer-assisted proofs become harder and harder. Since the result in [20] already covers the range $\beta\in(0,0.55]$, we thus focus on the parameter range $[0.5,1.9]$. We note that at $\beta=2$ a Hamiltonian-Hopf bifurcation occurs. In future work we intend to unfold this bifurcation and subsequently connect the homoclinic orbit that bifurcates to the branch covered by Theorem 1 (at that point we will know how far we have to push the current continuation technique beyond $\beta=1.9$ to connect all the way to the bifurcation point).

The rest of the paper is dedicated to the proof of Theorem 1. Our approach begins by rewriting (1.3) as a first order system for $(u_{1},u_{2},u_{3},u_{4})=(u,u^{\prime},u^{\prime\prime},u^{\prime\prime\prime})$ and then making the change of variables $(v_{1},v_{2},v_{3},v_{4})=(e^{u_{1}}-1,u_{2},u_{3},u_{4})$ to obtain

[TABLE]

There are two reasons for performing this change of variables. First, it turns the system into a polynomial vector field, which has technical advantages when performing the analysis to derive the necessary bounds. Second, while $u_{1}$ may become very negative for small values of $\beta$, the variable $v_{1}$ is always bounded from below by $-1$. Our goal is now to prove the existence of symmetric homoclinic solutions to (1.4) for all $\beta\in[0.5,1.9]$.

We split the problem into two parts. On the one hand a rigorous computational description of the local (un)stable manifold is required. On the other hand we need to solve, via a rigorous computational technique, a boundary value problem for the part of the orbit between the local invariant manifolds. We attack both parts by a continuation technique in the context of the radii polynomial approach. This parametrized Newton-Kantorovich method, adapted to a computational setting, is introduced in Section 2. In Section 3 we combine this with the parameterization method to obtain descriptions of the local (un)stable manifold of the equilibrium $0\in\mathbb{R}^{4}$. Essentially the same technique is then applied in Section 4 in a Chebyshev series setting to solve the boundary value problem. These two aspects are then combined into a rigorous computational continuation of the homoclinic solution to (1.3). We note that for smaller values of $\beta$ the boundary value problem is the more difficult part of the problem, as the orbit makes a bigger and bigger excursion away from the origin. On the other hand, for values of $\beta$ close to $2$ it is more difficult to obtain the local (un)stable manifold of the origin, as the real part of the eigenvalues tends to [math]. The algorithmic issues encountered when implementing the proof of Theorem 1 are discussed in Section 5.

Finally, let us mention that there is a growing literature on the subject of computer-assisted methods for proving existence of connecting orbits, see [22, 24, 26, 29, 30, 31, 32]. The main novel contribution of the current paper is to do rigorous continuation of a homoclinic orbit over a large range of parameter values. The method is generally applicable for connecting orbits problems in parameter dependent problems. In that sense Theorem 1, while providing a new result for traveling waves in the suspension bridge problem which complements earlier work, is an illustration.

2 The radii polynomial approach

In this section we present the functional analytic setup of our continuation method, which is formulated in terms of the radii polynomials, see Definition 5. It will be used both to find the stable manifold and to solve the boundary value problem. For more details and proofs we refer to [2, 11, 25].

Consider a sequence of Banach spaces $\left(X_{1},\|\cdot\|_{X_{1}}\right),\dots,\left(X_{d},\|\cdot\|_{X_{d}}\right)$ and the (product) Banach space

[TABLE]

with the induced norm defined by

[TABLE]

where $x=(x^{(1)},\dots,x^{(d)})\in X$ with $x^{(j)}\in X_{j}$ for $j=1,\dots,d$. Denote by

[TABLE]

the closed ball of radius $r>0$ centered at $y\in X$.

Consider an interval of parameters $[\beta_{0},\beta_{1}]\subset\mathbb{R}$ and $T:[\beta_{0},\beta_{1}]\times X\to X$ a Fréchet differentiable operator. For each $j=1,\ldots,d$, denote by $T^{(j)}:[\beta_{0},\beta_{1}]\times X\to X_{j}$ the projection of $T$ onto $X_{j}$. Let $\bar{x}_{\beta_{0}},\bar{x}_{\beta_{1}}\in X$ be approximate fixed points of $T(\beta_{0},\cdot)$ and $T(\beta_{1},\cdot)$, respectively, and define the linear interpolation

[TABLE]

Define the line of centers by $\{\bar{x}_{\beta}\mid\beta\in[\beta_{0},\beta_{1}]\}\subset X$. For each $j=1,\dots,d$, define the bounds

[TABLE]

for some $Y^{(j)}>0$ and $Z^{(j)}:\mathbb{R}^{+}\mapsto\mathbb{R}^{+}:r\to Z^{(j)}(r)$. The goal of the radii polynomial approach is to provide an efficient way to prove that an operator is a uniform contraction over a subset of $X$. This subset consists of small balls around the line of centers, provided by the linear interpolation between two numerical approximations of solutions at different parameter values.

Definition 2**.**

Let $X$ be a Banach space and $B\subset X$. Let $[\beta_{0},\beta_{1}]\subset\mathbb{R}$ be a set of parameters. A function $\tilde{T}:[\beta_{0},\beta_{1}]\times B\to B$ is a uniform contraction if there exists a constant $\kappa$ such that $0<\kappa<1$ and such that $\|\tilde{T}(\beta,x)-\tilde{T}(\beta,y)\|_{X}\leq\kappa\|x-y\|_{X}$ for all $x,y\in B$ and all $\beta\in[\beta_{0},\beta_{1}]$.

The following result is a restatement of the uniform contraction principle (e.g. see [9] for a proof).

Theorem 3** **(Uniform Contraction Principle).

Suppose there exists some $r>0$ such that

[TABLE]

is a uniform contraction, then for every $\beta\in[\beta_{0},\beta_{1}]$, there exists a unique $\tilde{x}(\beta)\in B_{r}(\bar{x}_{\beta})$ such that $T(\beta,\tilde{x}(\beta))=\tilde{x}(\beta)$. Moreover, the function $\beta\mapsto\tilde{x}(\beta)$ is of class $C^{k}$ if $(\beta,x)\mapsto T(\beta,x)$ is of class $C^{k}$.

With the bounds $Y$ and $Z$ on the residue and the derivative of $T$, see Equations (2.2) and (2.3), contractivity can be checked explicitly. This is expressed by the next theorem (we refer to [2, 11, 25] for a proof).

Theorem 4**.**

Given a set of parameters $[\beta_{0},\beta_{1}]\subset\mathbb{R}$, consider the set of centers $\{\bar{x}_{\beta}\mid\beta\in[\beta_{0},\beta_{1}]\}$ with $\bar{x}_{\beta}$ given by (2.1). Assume that $T:[\beta_{0},\beta_{1}]\times X\to X$ is an operator satisfying the bounds (2.2) and (2.3). If there exists $r>0$ such that $Y^{(j)}+Z^{(j)}(r)<r$, for each $j=1,\dots,d$, then $\tilde{T}$, defined by (2.4), is a uniform contraction (on $B_{r}(0)$).

Assuming we have determined explicit bounds $Y^{(j)}$ and $Z^{(j)}(r)$, where in practice the latter is a polynomial with positive coefficients, satisfying (2.2) and (2.3). It is convenient to introduce the radii polynomials, which provide an efficient way in verifying the hypotheses of Theorem 4.

Definition 5**.**

Let $Y=(Y^{(1)},\dots,Y^{(d)})$ and $Z=(Z^{(1)},\dots,Z^{(d)})$ be the bounds on the operator $T_{\beta}$ as given by (2.2) and (2.3) respectively. We define the radii polynomials as

[TABLE]

One can see that the radii polynomials depend on the upper bounds $Y$ and $Z$, and therefore they are not uniquely defined. But the smaller these bounds are, the higher the chances are to prove that the operator $T_{\beta}$ is a contraction over a ball around the approximation. The following result shows how the radii polynomials are used in practice to give us the value of $r$ for which we can apply Theorem 4.

Proposition 6**.**

Let

[TABLE]

and assume that $\mathcal{I}\neq\emptyset$. Then $\mathcal{I}$ is an open interval of $\mathbb{R}^{+}$, i.e., $\mathcal{I}=(r_{\rm min},r_{\rm max})$. For any $r_{0}\in(r_{\rm min},r_{\rm max})$, $\tilde{T}:B_{r_{0}}(0)\times\left[\beta_{0},\beta_{1}\right]\longrightarrow B_{r_{0}}(0)$ is a uniform contraction.

3 Parameterization of the stable manifold

In this section we compute an approximate parameterization of the (local) stable manifold at 0, and provide explicit error bounds on this parameterization. This is done by combining the ideas of the parameterization method (first introduced in [5, 6, 7], see also [13]) and of rigorous computation (following the approach of [27, 1]). Having computed the parameterization, we will be able to obtain the homoclinic connection in the next section by taking advantage of the fact that it is now enough to compute an orbit on a finite time interval, i.e., an orbit that ends up in the local stable and unstable manifolds (or rather, we compute and verify an orbit that starts from the symmetric section and ends up, after some finite time, in the local stable manifold, see (1.3)).

3.1 Looking for the stable manifold as a zero finding problem $F(\beta,a)=0$

The first step is to recast the problem of finding a parameterization as looking for a zero of a map $F$, which is the aim of this section. Setting

[TABLE]

Equation (1.4) is rewritten as $v^{\prime}=\Psi_{\beta}(v)$. The Jacobian at the origin is

[TABLE]

and one finds that for $\beta\in[0,2)$ it has two complex conjugated eigenvalues with negative real part, which we denote by $\lambda(\beta)$ and $\lambda^{*}(\beta)$:

[TABLE]

The associated eigenvectors are given by $V(\beta)$ and $V^{*}(\beta)$, where

[TABLE]

The stable manifold at [math] is thus two dimensional. Since $\Psi_{\beta}$ is analytic we may look for an analytic local parameterization of this manifold. We will look for this parameterization as a power series

[TABLE]

with standard multi-index notation: $\alpha\in\mathbb{N}^{2}$, $|\alpha|=\alpha_{1}+\alpha_{2}$, $\theta^{\alpha}=\theta_{1}^{\alpha_{1}}\theta_{2}^{\alpha_{2}}$, and satisfying

[TABLE]

together with the invariance equation

[TABLE]

Remark 7**.**

Even though the vector field $\Psi_{\beta}$ is real, the fact that we have two complex eigenvalues makes it easier to first look for a parameterization $Q_{\beta}$ of the complex manifold and then recover the real parameterization (the one which will be of interest in the next section for computing the homoclinic orbit) by considering

[TABLE]

see [16, 27] for more details. This is due to the underlying symmetry $a_{(\alpha_{2},\alpha_{1})}=a^{*}_{(\alpha_{1},\alpha_{2})}$, which is respected by the function $F$ introduced below.

Plugging the power series (3.3) into the invariance equation (3.5) we get

[TABLE]

where $\star$ stands for the Cauchy product. We recall that, given two sequences $u$ and $v$ of complex numbers (indexed over $\mathbb{N}^{2}$), their Cauchy product is the sequence defined by

[TABLE]

where $\sigma\leq\alpha$ means $\sigma_{1}\leq\alpha_{1}$ and $\sigma_{2}\leq\alpha_{2}$ (and similarly $\alpha-\sigma=(\alpha_{1}-\sigma_{1},\alpha_{2}-\sigma_{2})$).

Notice that the additional conditions (3.4) imply that the coefficients of total degree 0 and 1 are equal on both sides of Equation (3.6).

Finding an analytic parameterization of the local manifold is now equivalent to find a zero of $F(\beta,\cdot)$, defined component-wise by

[TABLE]

3.2 Getting to the fixed point formulation

Let $\nu\geq 1$ and denote by $\ell_{\nu}^{1}$ the Banach space of complex valued sequences $u=\left(u_{\alpha}\right)_{|\alpha|\geq 0}$ such that

[TABLE]

This space is a Banach algebra under the Cauchy product, which gives us control on the quadratic terms.

Lemma 8**.**

For $u,v\in\ell_{\nu}^{1}$, $\left\|u\star v\right\|_{1,\nu}\leq\left\|u\right\|_{1,\nu}\left\|v\right\|_{1,\nu}$.

Definition 9**.**

In this section we consider

[TABLE]

We are going to look for zeros $a$ of $F(\beta,\cdot)$ in the space $X$. Notice that $a\in X$ means that

[TABLE]

which ensures that the associated parameterization $Q_{\beta}$ is well defined at least for

[TABLE]

We now explain how to rigorously determine a parameterization of the manifold for all values of $\beta$ in a given interval $[\beta_{0},\beta_{1}]$. It will be more convenient to work with a rescaled parameter $s$ ranging between [math] and $1$. Therefore we define

[TABLE]

We want to get a parameterization $a(s)$ such that

[TABLE]

Since we are working on the interval $[\beta_{0},\beta_{1}]$, parameterized by $s\in[0,1]$, we have altered the notation of the parametrization of the coefficients $a$ slightly compared to Section 3.1, namely $a(s)$ instead of $a(\beta)$. We will use $a(s)$ throughout the remainder of the paper, except in Section 4.1, where the notation $a(\beta)$ is more appropriate.

We first compute approximate zeros $\bar{a}(0)$ and $\bar{a}(1)$ of $F(\beta_{0},\cdot)$ and $F(\beta_{1},\cdot)$ respectively, by solving numerically the truncated problem (for $s=0$ and $s=1$)

[TABLE]

for some $N\geq 1$, and by padding the obtained solutions with 0 to get elements of $X=(\ell^{1}_{\nu})^{4}$. We then define for $s\in[0,1]$

[TABLE]

If $\bar{a}(0)$ and $\bar{a}(1)$ are two good approximate zeros (of $F(\beta_{0},\cdot)$ and $F(\beta_{1},\cdot)$ respectively) and if $|\beta_{1}-\beta_{0}|$ is not too large, $\bar{a}(s)$ should be a good approximate zero of $F(\beta_{s},\cdot)$ for each $s\in[0,1]$. We are going to reformulate this claim into a mathematical statement and prove that in a given neighbourhood of $\bar{a}(s)$ there exist a unique zero of $F(\beta_{s},\cdot)$ for all $s\in[0,1]$. To put this in the framework described in Section 2, we consider the operator

[TABLE]

where $A$, defined below, is an approximate inverse of $D_{a}F(\beta_{0},\bar{a}(0))$. Namely, for $N$ large enough,

[TABLE]

should be a reasonably good approximation of $D_{a}F(\beta_{0},\bar{a}(0))$, where, for any $k\geq N$, $\tilde{M}_{k}$ is the $4(k+1)$ by $4(k+1)$ block diagonal matrix

[TABLE]

with $I_{4}$ the 4 by 4 identity matrix. Finally, we define $A$ as

[TABLE]

where $J$ is a numerical approximation of $\left(D_{a}F^{[N]}(\beta_{0},\bar{a}(0))\right)^{-1}$, while the $M_{k}=\tilde{M}_{k}^{-1}$ are exact inverses. The operators $A^{\dagger}$ and $A$ are then approximate inverses of each other: approximate in the finite part and exact in the infinite tail.

Remark 10**.**

To make sense of this matrix representation of $A^{\dagger}$ and $A$, as well as $\tilde{M}_{k}$ and $M_{k}$, one should think of $a_{\alpha}$ as an infinite vector where the elements are ordered according to increasing degree $|\alpha|=\alpha_{1}+\alpha_{2}$ and within fixed degree by increasing $\alpha_{2}$, while also taking into account that each $a_{\alpha}$ is a vector in $\mathbb{C}^{4}$. This means that $a$ is represented as

[TABLE]

The above representation describes the operators as infinite matrices where each element $A_{\alpha^{\prime},\alpha}$ is a linear operator on $\mathbb{C}^{4}$, i.e. a $4\times 4$ matrix that we will occasionally denote by $A_{\alpha^{\prime},\alpha}=\{A_{\alpha^{\prime},\alpha}^{(i,j)}\}_{1\leq i,j\leq 4}$.

We now follow the ideas described in Section 2, using the Banach space $X=\left(\ell^{1}_{\nu}\right)^{4}$ endowed with the norm $\|a\|_{X}=\max\limits_{j=1,\ldots,4}\|a^{(j)}\|_{1,\nu}$. In the next subsections we are going to compute the bounds $Y^{(j)}$ and $Z^{(j)}(r)$ and the associated radii polynomials, and then prove that for some positive $r$ each radii polynomial $p_{j}(r)$ is negative, which will yield (for each $s\in[0,1]$) the existence of a unique zero $a(\beta_{s})$ of $F(\beta_{s},\cdot)$ in the ball of radius $r$ around $\bar{a}(s)$. At this point we will know that $\bar{a}(s)$ defines an approximate parameterization of the stable manifold, with an error bound controlled by $r$. We will use this in Section 4 to prove the existence of a homoclinic orbit for all $\beta\in[\beta_{0},\beta_{1}]$. Moreover, derivatives of the manifold with respect to $\theta$, which will be needed in Section 4.1, can also be approximately computed with rigorous control on the error bound, see Lemma 14.

3.3 The bound $Y$

In this section we focus on the bound $Y$ defined in (2.2). Let $|A|$ denote the component-wise absolute value of $A$. In order to define the bound we are looking for, we try to bound every term $\left(T(\beta_{s},\bar{a}(s))-\bar{a}(s)\right)_{\alpha}^{(j)}$ with $|\alpha|\geq 0$ and $j=1,2,3,4$:

[TABLE]

A straightforward calculation (using that $|\lambda(\beta)|=1$ and computing the derivatives of $\lambda$ and $V$ with respect to $\beta$) yields that, for all $|\alpha|\geq 0$,

[TABLE]

where

[TABLE]

Since $\left(\bar{a}(s)\right)_{\alpha}=0$ for all $|\alpha|\geq N$ and $F(\beta,\cdot)$ is quadratic in $a$, we have that $F_{\alpha}(\beta_{s},\bar{a}(s))$ vanishes as soon as $|\alpha|\geq 2N-1$. Therefore, we define $\tilde{F}$ component-wise by

[TABLE]

and then set

[TABLE]

so that

[TABLE]

3.4 The bound $Z$

In this section we derive the bound $Z$ defined in (2.3). Let $b,c\in B_{r}(0)$. We split $D_{a}T(\beta_{s},\bar{a}(s)+b)c$ in three terms which will be easier to bound separately. For each $j=1,\ldots,4$,

[TABLE]

The bounds $Z_{i}\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,\left(Z_{i}^{(1)},Z_{i}^{(2)},Z_{i}^{(3)},Z_{i}^{(4)}\right)\in\mathbb{R}^{4}$ ($i=0,1,2$) are given in the following subsections.

3.4.1 The bound $Z_{0}$

From the definitions of $A$ and $A^{{\dagger}}$ we get

[TABLE]

The finite matrix $B=I_{2N(N+1)}-JD_{a}F^{[N]}(\beta_{0},\bar{a}(0))$ can be computed using interval arithmetic. To obtain the bound $Z_{0}$ we only need to compute the operator norm of $B$ (as acting on $(\ell^{1}_{\nu})^{4}$). This is the content of the following lemma.

Lemma 11**.**

Let $h=\left(h_{\alpha}\right)_{\alpha\in\mathbb{N}^{2}}\in\ell^{1}_{\nu}$ (with $h_{\alpha}\in\mathbb{C}$ for all $\alpha$) and $\Gamma$ a linear operator acting on $\ell^{1}_{\nu}$. Then

[TABLE]

In particular, if $\Gamma$ consists in a finite block $\Gamma^{[N]}$ of size $N(N+1)/2\times N(N+1)/2$ and a diagonal tail $\left(\gamma_{\alpha}\right)_{|\alpha|\geq N}$

[TABLE]

then

[TABLE]

Hence, we define

[TABLE]

with the notation $B_{\alpha^{\prime},\alpha}^{(i,j)}$ introduced in Remark 10, and set

[TABLE]

to obtain

[TABLE]

3.4.2 The bound $Z_{1}$

This term is the most involved one to bound tightly, so again we split it into several parts that we bound separately. For each $j=1,\ldots,4$,

[TABLE]

Let us focus first on the first term. Since

[TABLE]

we get that

[TABLE]

For the tail $|\alpha|\geq N$ we find

[TABLE]

which we estimate by

[TABLE]

Now we use Lemma 11 again and from the fact that $|(n-k)\lambda(\beta_{0})+k\lambda^{*}(\beta_{0})|\geq n|\Re(\lambda(\beta_{0}))|=n\frac{\sqrt{2-\beta_{0}}}{2}$ we infer that

[TABLE]

and we are done with the first term. For the second term

[TABLE]

Again we use Lemma 11 to obtain

[TABLE]

where we recall that $J$ is the block of $A$ corresponding to the floating point data, see (3.7), and $K^{(i,j)}$ is defined by (3.8). Finally, computing the derivative of $\lambda$ with respect to $\beta$, we get that

[TABLE]

Now we need the following lemma, which is a slightly modified version of Lemma 11.

Lemma 12**.**

Let $c=\left(c_{\alpha}\right)_{\alpha\in\mathbb{N}^{2}}\in\left(\ell^{1}_{\nu}\right)^{4}$. We denote by $e$ the vector such that for all $\alpha$, $e_{\alpha}=|\alpha|c_{\alpha}$. Then for all $1\leq j\leq 4$

[TABLE]

where

[TABLE]

Using Lemmas 11 and 12 we infer that

[TABLE]

where

[TABLE]

Finally, putting everything together, we define

[TABLE]

so that

[TABLE]

3.4.3 The bound $Z_{2}$

Since

[TABLE]

we directly use one more time Lemma 11 and set

[TABLE]

so that

[TABLE]

3.5 Use of the uniform contraction principle and error bounds

Following (2.5), we set

[TABLE]

If we find an $r>0$ such that $p^{(j)}(r)<0$ for all $j=1,\ldots,4$, then according to Proposition 6 we have validated the numerical approximation $\bar{a}(s)$ of the local stable manifold for $\beta=\beta_{s}$, for every $s\in[0,1]$.

Proposition 13**.**

For every $s\in[0,1]$, let

[TABLE]

be the approximate parameterization of the complex local stable manifold that we have computed (for $\beta=\beta_{s}$). Assume that there exists an $r>0$ such that $p^{(j)}(r)<0$ for all $j=1,\ldots,4$. Then, for each $s\in[0,1]$, there exists a parameterization $Q_{\beta_{s}}$ of the complex local stable manifold (for $\beta=\beta_{s}$) of the form

[TABLE]

which is well defined for all $\theta\in\mathbb{C}^{2}$ satisfying $|\theta|_{\infty}\leq\nu$. Let

[TABLE]

then we have the error bound $|\hat{h}_{\beta_{s}}(\theta)|_{\infty}\leq r$ for all $|\theta|_{\infty}\leq\nu$. These statements still hold true for the real (approximate and exact) local stable manifold, defined by

[TABLE]

for all $\theta\in\mathbb{R}^{2}$ satisfying $|\theta|_{2}\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,\sqrt{\theta_{1}^{2}+\theta_{2}^{2}}\leq\nu$.

Proof.

Proposition 6 yields that, for each $s\in[0,1]$, there exists a unique fixed point $a(s)$ of $T(\beta_{s},\cdot)$ in the ball of radius $r$ around $\bar{a}(s)$. The operator $A$ is injective since its non-diagonal part $J$ is invertible. The latter follows from the fact that, see (3.9),

[TABLE]

where the final inequality is implied by $p_{j}(r)<0$. Here the operator norm on $X^{[N]}\cong\mathbb{R}^{2N(N+1)}$ is induced by the one on $X=(\ell^{1}_{\nu})^{4}$. Hence the fixed point $a(s)$ of $T$ solves $F(\beta_{s},a(s))=0$. By construction $Q_{\beta_{s}}$ is a parameterization of the local stable manifold defined for $|\theta|_{\infty}\leq\nu$, and for such $\theta$,

[TABLE]

In the following section we use these approximations to rigorously prove the existence of homoclinic orbits for every parameter $\beta$ in $[0.5,1.9]$. To do so, we will also need control on the derivative of the parameterization $P_{\beta_{s}}$, which is provided by the theory of analytic functions. Define

[TABLE]

For all $s\in[0,1]$, the function $\hat{h}_{\beta_{s}}$, defined by (3.11), is analytic. Since $h_{\beta_{s}}(\theta)=\hat{h}_{\beta_{s}}(\theta_{1}+\rm{i}\theta_{2},\theta_{1}-\rm{i}\theta_{2})$, we can control the derivative of $h_{\beta_{s}}$ (on a smaller domain) by a bound on $\hat{h}_{\beta_{s}}$. This is the content of the following lemma, of which the proof can be found in [16].

Lemma 14**.**

Assume that $\hat{h}:D_{\infty,\nu}(\mathbb{C}^{2})\subset\mathbb{C}^{2}\to\mathbb{C}^{4}$ is analytic, where

[TABLE]

and $\delta>0$ is such that

[TABLE]

Consider $h:D_{2,\nu}(\mathbb{R}^{2})\subset\mathbb{R}^{2}\to\mathbb{R}^{4}$ defined by $h(\theta)=\hat{h}(\theta_{1}+\rm{i}\theta_{2},\theta_{1}-\rm{i}\theta_{2})$, where

[TABLE]

Then for any $\rho<\nu$ we have

[TABLE]

4 Parameterized families of symmetric homoclinic orbits

In this section, we apply the technique of Section 2 in a Chebyshev series setting to rigorously prove existence of parameterized families of symmetric homoclinic orbits. More precisely, we present all necessary estimates and bounds in order to demonstrate that solutions of (1.3) exist for all $\beta\in[0.5,1.9]$.

4.1 A projected boundary value problem formulation

We begin by transforming the symmetric homoclinic orbit problem (1.3) into a projected boundary value problem (BVP). In order to set up the projected BVP, we first use the symmetry of the orbit to simplify the problem and therefore solve only for “half of the orbit”. The following lemma provides a strategy to do this.

Lemma 15**.**

Let $u_{0},u_{2}$ and $t_{0}$ be arbitrary numbers, and let $u(t)$ be the solution of the initial value problem

[TABLE]

Then $u(-t+2t_{0})=u(t)$ for all $t$ for which the solution $u$ is defined.

Proof.

It is straightforward to verify that $u(-t+2t_{0})$ is also a solution of the initial value problem. By the theorem of existence and uniqueness for ODEs, it follows that $u(-t+2t_{0})=u(t)$ for all $t$ in the domain definition of $u$. ∎

Using the previous result, we fix a number $t_{0}=L>0$, and it follows that to solve (1.3), it is enough to solve

[TABLE]

The idea now is to modify the boundary value problem (4.1) in a way that the boundary value at $t=\infty$ is removed by a projected boundary value at $t=L$ where we impose at that time that $\left(u(L),u^{\prime}(L),u^{\prime\prime}(L),u^{\prime\prime\prime}(L)\right)\in W_{\rm loc}^{s}(0)$, a local stable manifold at [math]. In order to achieve this step, we use the theory of Section 3 to obtain a real-valued parameterization $P_{\beta}$ of $W_{\rm loc}^{s}(0)$ at the parameter value $\beta\in[0.5,1.9]$:

[TABLE]

which is well-defined for all $\theta\in D_{2,\tilde{\nu}}(\mathbb{R}^{2})=\left\{\theta\in\mathbb{R}^{2}:|\theta|_{2}=\sqrt{\theta_{1}^{2}+\theta_{2}^{2}}\leq\tilde{\nu}\right\}$, where the size $\tilde{\nu}=\tilde{\nu}(\beta)$ of the domain of $P_{\beta}$ changes as the parameter $\beta\in[0.5,1.9]$ varies. Using the parameterization, we impose that

[TABLE]

for some $\theta\in D_{2,\tilde{\nu}}(\mathbb{R}^{2})$, which implies that the orbit lies in the stable manifold. This introduces an indeterminacy that needs to be resolved. Namely, there is a one parameter family of pairs $(L,\theta)$ solving (4.2) while describing the same orbit. To overcome this, we impose that $\theta\in\partial D_{2,\rho}(\mathbb{R}^{2})=\left\{\theta\in\mathbb{R}^{2}:|\theta|_{2}=\sqrt{\theta_{1}^{2}+\theta_{2}^{2}}=\rho\right\}$, for some fixed $\rho<\tilde{\nu}$, and we solve for the angle $\psi$. More precisely, we consider $\theta$ such that $\sqrt{\theta_{1}^{2}+\theta_{2}^{2}}=\rho$ by setting $\theta_{1}+\text{i}\theta_{2}=\rho e^{\text{i}\psi}$ for some $\psi\in[0,2\pi)$. In this case, the evaluation of the parameterization of the local stable manifold along $\partial D_{2,\rho}(\mathbb{R}^{2})$ reduces to

[TABLE]

We slightly abuse notation by using the same notation ${P}_{\beta}$ to denote both ${P}_{\beta}(\theta)$ and ${P}_{\beta}(\psi)$. We can therefore define the projected BVP

[TABLE]

where $L>0$ and $\psi\in[0,2\pi)$ are variables. As in Section 1, we make the change of variables

[TABLE]

and set $v=(v^{(1)},v^{(2)},v^{(3)},v^{(4)})$ to obtain that $v^{\prime}=\Psi_{\beta}(v)$, where $\Psi_{\beta}:\mathbb{R}^{4}\to\mathbb{R}^{4}$ is the vector field given by the right-hand side of (1.4). We rescale time via $t\mapsto t/L$ so that (4.3) becomes

[TABLE]

A triplet $(L,\psi,v)$ satisfying (4.4) thus corresponds to a symmetric homoclinic solution of the suspension bridge equation. The rest of this section is dedicated to applying the technique of Section 2 in a Chebyshev series setting to rigorously prove existence of parameterized families of solutions of the projected BVP (4.4) for all $\beta\in[0.5,1.9]$. This begins by defining a zero finding problem $F=0$ whose solutions correspond to symmetric homoclinic solutions of the suspension bridge equation.

4.2 Setting up the zero finding problem using Chebyshev series

Now that $v^{(i)}(t)$ is defined on $[-1,1]$ and needs to solve a boundary value problem, describing $v^{(i)}(t)$ in terms of a Chebyshev series is a natural choice, see [10, 14, 24, 28]. Denote by $T_{k}:[-1,1]\to\mathbb{R}$ the $k$-th Chebyshev polynomial with $k\geq 0$, where $T_{0}(t)=1$, $T_{1}(t)=t$ and $T_{k+1}(t)=2tT_{k}(t)-T_{k-1}(t)$ for $k\geq 1$. One way to characterize the Chebyshev polynomials is through the identity $T_{k}(t)=\cos(k\arccos t)$, from which it follows that $\|T_{k}\|_{\infty}=1$, $T_{k}(1)=1$, and $T_{k}(-1)=(-1)^{k}$.

For each $i=1,2,3,4$, we expand $v^{(i)}$ using a Chebyshev series expansion, that is

[TABLE]

For each $i=1,2,3,4$, denote by $x^{(i)}\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,\{x_{k}^{(i)}\}_{k\geq 0}$ the infinite dimensional vector of Chebyshev coefficients of $v^{(i)}$. The vector field is analytic (polynomial) and therefore the solutions (if they exist) of the projected BVP (4.4) are analytic. By the Paley-Wiener theorem, this implies that the Chebyshev coefficients of each component of $v$ decay geometrically to zero. Hence, there exists a number $\nu>1$ such that $x^{(i)}\in\ell_{\nu}^{1}$ for each $i=1,2,3,4$, where

[TABLE]

We remark that throughout this section $\nu\geq 1$.

Remark 16** **(Notation).

The decay rate $\nu$ in the definition of the Banach space $\ell_{\nu}^{1}$ appears both in the current section and in Section 3. Both values need not to be the same. Therefore, to avoid confusion, we denote by $\tilde{\nu}$ the value from Section 3. Moreover, although the sequence space $\ell_{\nu}^{1}$ as considered above is slightly different from the one used in Section 3, we nevertheless use the same notation, since the spaces and norms are completely analogous to those used in Section 3.2.

The dual space can be characterized as follows.

Lemma 17**.**

The dual space $(\ell^{1}_{\nu})^{*}$ is isomorphic to

[TABLE]

For all $a\in\ell^{1}_{\nu}$ and $c\in\ell^{\infty}_{\nu^{-1}}$ we have

[TABLE]

The following lemma is analogous to Lemma 11.

Lemma 18**.**

Let $\Gamma\in B(\ell^{1}_{\nu})$, the space of bounded linear operators from $\ell^{1}_{\nu}$ to itself, acting as $(\Gamma a)_{i}=\sum_{j\geq 0}\Gamma_{i,j}a_{j}$. Define the weights $\omega=(\omega_{k})_{k\geq 0}$ by $\omega_{0}=1$ and $\omega_{k}=2\nu^{k}$ for $k\geq 1$. Then

[TABLE]

The Banach space of unknowns $x\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,(L,\psi,x^{(1)},x^{(2)},x^{(3)},x^{(4)})$ is

[TABLE]

endowed with the norm

[TABLE]

In terms of Chebyshev coefficients the differential equation $\dot{v}=L\Psi_{\beta}(v)$ becomes (see e.g. [14])

[TABLE]

for all $k\geq 1$. Here $x_{k\pm 1}^{(i)}\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,x_{k+1}^{(i)}-x_{k-1}^{(i)}$, and $*$ denotes the discrete convolution product $\ast:\ell_{\nu}^{1}\times\ell_{\nu}^{1}\to\ell_{\nu}^{1}$ defined as follows. Let $a,b\in\ell_{\nu}^{1}$, then for all $k\geq 0$ the $k$-th entry of the convolution product $a\ast b$ is given by

[TABLE]

The choice of norm and convolution product is justified by the fact $\ell_{\nu}^{1}$ is a Banach algebra, that is $\left\|a\ast b\right\|_{1,\nu}\leq\left\|a\right\|_{1,\nu}\left\|b\right\|_{1,\nu}$, for all $a,b\in\ell_{\nu}^{1}$.

The symmetry conditions $v^{(2)}(-1)=v^{(4)}(-1)=0$ reduce to

[TABLE]

and the boundary conditions $v(1)=P_{\beta}(\psi)$ become

[TABLE]

The full set of equations that we want to solve is thus $F(\beta,x)=0$, where

[TABLE]

In order to solve rigorously the problem $F(\beta,x)=0$ in the Banach space $X$, for all $\beta\in[0.5,1.9]$, we apply the radii polynomial approach of Section 2.

4.3 The finite dimensional reduction of the zero finding problem

Having identified the operator $F$ given in (4.12) whose zeros correspond to symmetric homoclinic orbits of (1.2), the next step is to compute numerical approximations, which requires considering a finite dimensional projection of the Banach space $X$ given in (4.7). Given a sequence $a=(a_{k})_{k\geq 0}\in\ell_{\nu}^{1}$, denote by $a^{[m]}=(a_{0},\dots,a_{m-1})\in\mathbb{R}^{m}$ the Galerkin projection of $a$ onto the first $m$ Chebyshev coefficients. Given an infinite dimensional vector $x=(L,\psi,x^{(1)},x^{(2)},x^{(3)},x^{(4)})\in X$, denote

[TABLE]

In this context, the finite dimensional Banach space $\mathbb{R}^{4m+2}$ is the finite dimensional projection of $X=\mathbb{R}^{2}\times(\ell_{\nu}^{1})^{4}$, and $x^{[m]}$ is the finite dimensional projection of $x$. We slightly abuse the notation by denoting $x^{[m]}\in X$ as the vector built from $x^{[m]}\in\mathbb{R}^{4m+2}$ by padding each entry $(x^{(i)})^{[m]}$ ($i=1,2,3,4$) with infinitely many zeros. The finite dimensional projection of $F$ given in (4.12) is defined as

[TABLE]

We want to compute on $F^{[m]}$, but it depends on the parameterization $P_{\beta_{s}}$, which itself depends on infinitely many Taylor coefficients. To remedy this, we consider a finite dimensional reduction of $P_{\beta_{s}}$. Recalling (3.13), for every $s$, denote by $\overline{P}_{\beta_{s}}$ the computable approximation of the stable manifold given by

[TABLE]

where $\bar{a}_{0,\alpha}$ and $\bar{a}_{1,\alpha}$ are the numerical approximations of the coefficients of the stable manifold for $\beta_{0}$ and $\beta_{1}$ respectively.

Finally, let $\overline{F}(\beta,x^{[m]})$ denote the finite dimensional projection of the operator using a Galerkin projection on the last four components and using the finite dimensional approximation $\overline{P}_{\beta}$ for the parameterization of the stable manifold. More explicitly,

[TABLE]

with $\overline{F}^{(i)}(\beta,x^{[m]})\,\stackrel{{\scriptstyle\mbox{\tiny{\raisebox{0.0pt}[0.0pt][0.0pt]{def}}}}}{{=}}\,\left\{\overline{f}_{k}^{(i)}(\beta,x^{[m]})\right\}_{k=0}^{m-1}$ for $i=1,2,3,4$, and

[TABLE]

while $\overline{f}_{k}^{(i)}(\beta,x^{[m]})=f_{k}^{(i)}(\beta,x^{[m]})$ for all $k=1,\dots,m-1$, see (4.8). Having identified $\overline{F}:\mathbb{R}\times\mathbb{R}^{4m+2}\to\mathbb{R}^{4m+2}:(\beta,x^{[m]})\mapsto\overline{F}(\beta,x^{[m]})$ defined in (4.15) as the finite dimensional reduction of $F$ given in (4.12), we can apply the finite dimensional Newton’s method to find numerical approximations. The next step is to define an infinite dimensional Newton-like operator $T:\mathbb{R}\times X\to X$ on which we apply the uniform contraction principle (via the radii polynomial approach of Section 2).

4.4 The Newton-like operator for the homoclinic orbit

Let $\beta_{0}<\beta_{1}$ be two different parameter values, and consider two numerical approximations $\bar{x}_{0}$ and $\bar{x}_{1}$ such that $F(\beta_{0},\bar{x}_{0})\approx 0$ and $F(\beta_{1},\bar{x}_{1})\approx 0$. In practice we find $\bar{x}_{i}$ by solving $\overline{F}(\beta_{i},\cdot)=0$ numerically. For every $s\in[0,1]$, set

[TABLE]

and

[TABLE]

We denote $\bar{x}_{s}=(\bar{L}_{s},\bar{\psi}_{s},\bar{x}_{s}^{(1)},\bar{x}_{s}^{(2)},\bar{x}_{s}^{(3)},\bar{x}_{s}^{(4)})\in X$ for $s\in[0,1]$, and we recall that each $\bar{x}_{s}^{(j)}$ is obtained from $(\bar{x}_{s}^{(j)})^{[m]}\in\mathbb{R}^{m}$ by padding with zeros. Similarly, we denote $\Delta\bar{x}=(\Delta\bar{L},\Delta\bar{\psi},\Delta\bar{x}^{(1)},\Delta\bar{x}^{(2)},\Delta\bar{x}^{(3)},\Delta\bar{x}^{(4)})$.

We now construct a fixed point operator $T(\beta,x)=x-AF(\beta,x)$ so that it is a uniform contraction over the interval of parameters $[\beta_{0},\beta_{1}]$, whose fixed points $x=x(\beta)$ correspond to zeros of $F(\beta,\cdot)$ at a given parameter value $\beta\in[\beta_{0},\beta_{1}]$. The operator $A$ is constructed as an approximate inverse of $DF(\beta_{0},\bar{x}_{0})$. Let $\bar{x}_{0}$ be such that $\overline{F}(\beta_{0},\bar{x}_{0})\approx 0$ and let $A^{[m]}\approx(D\overline{F}(\beta_{0},\bar{x}_{0}))^{-1}$ be a numerical approximation of the inverse of the Jacobian matrix. We decompose the $(4m+2)\times(4m+2)$ matrix $A^{[m]}$, into 36 blocks as

[TABLE]

Here $A^{[m]}_{i,j}$ is scalar for $1\leq i,j\leq 2$, $A^{[m]}_{i,j}$ is a row vector of length $m$ for $1\leq i\leq 2$, $3\leq j\leq 6$, $A^{[m]}_{i,j}$ is a column vector of length $m$ for $3\leq i\leq 6$, $1\leq j\leq 2$, and $A^{[m]}_{i,j}$ is a $m\times m$ matrix for $3\leq i,j\leq 6$.

Definition 19** **(Definition of $A$).

We extend this finite dimensional operator $A^{[m]}=\{A^{[m]}_{i,j}|1\leq i,j\leq 6\}$ to an operator $A=\{A_{i,j}|1\leq i,j\leq 6\}$ on $X$ defined block-wise as

• •

$A_{i,j}\in\mathbb{R}$* for $1\leq i,j\leq 2$, where $A_{i,j}=A^{[m]}_{i,j}$;*

• •

$A_{i,j}\in(\ell_{\nu}^{1})^{*}$* for $1\leq i\leq 2$ and $3\leq j\leq 6$, where $A_{i,j}$ is $A^{[m]}_{i,j}$ padded with zeros;*

• •

$A_{i,j}\in\ell_{\nu}^{1}$* for $3\leq i\leq 6$ and $1\leq j\leq 2$, where $A_{i,j}$ is $A^{[m]}_{i,j}$ padded with zeros;*

• •

$A_{i,j}\in B(\ell_{\nu}^{1},\ell_{\nu}^{1})$* for $3\leq i,j\leq 6$, where*

[TABLE]

with $\delta_{i,j}$ the usual Kronecker delta.

Here $(\ell_{\nu}^{1})^{*}$ is the dual of $\ell_{\nu}^{1}$. As an example, for $1\leq i\leq 2$, $3\leq j\leq 6$, we have $A_{i,j}a=A^{[m]}_{i,j}a^{[m]}$. The action of $A$ on $x=(L,\psi,x^{(1)},x^{(2)},x^{(3)},x^{(4)})\in X$ is thus

[TABLE]

where $(Ax)^{(i)}\in\mathbb{R}$ for $i=1,2$ and $(Ax)^{(i)}\in\ell^{1}_{\nu}$ for $i=3,4,5,6$.

We consider the Newton-like operator

[TABLE]

where $s\in[0,1]$ and $A$ is as in Definition 19.

Lemma 20**.**

Given the operator $A$ as in Definition 19. Then $T:\mathbb{R}\times X\to X$.

Proof.

Consider $x\in X=\mathbb{R}^{2}\times(\ell_{\nu}^{1})^{4}$ and $\beta\in\mathbb{R}$. By construction of $A$, in particular the infinite diagonal tail chosen in (4.17), it is straightforward to verify that $\left(AF(\beta,x)\right)^{(i)}\in\mathbb{R}$ for $i=1,2$ and $\left(AF(\beta,x)\right)^{(i)}\in\ell_{\nu}^{1}$ for $3\leq i\leq 6$. ∎

Showing the existence of parameterized fixed points of $T$ defined in (4.18) is done by applying the general technique of Section 2. This requires computing the bounds $Y^{(j)}$ satisfying (2.2) and the bounds $Z^{(j)}$ satisfying (2.3) for $j=1,\dots,6$. We recall that since $X=\prod_{j=1}^{6}X_{j}=\mathbb{R}^{2}\times(\ell_{\nu}^{1})^{4}$, we have that $\|\cdot\|_{X_{j}}$ denotes the absolute value for $j=1,2$ and the $\ell_{\nu}^{1}$ norm for $j=3,4,5,6$.

4.5 The $Y$ bound for the homoclinic orbit problem

We recall the definition of the bounds $Y^{(j)}$ in (2.2). In our context, $Y^{(j)}$ is a bound satisfying

[TABLE]

We begin by expanding each component of

[TABLE]

as a polynomial in $s$. Given $s\in[0,1]$ and $j=1,2,3,4$, denote $\bar{x}_{s}^{(j)}=\left(\bar{x}_{s,k}^{(j)}\right)_{k\geq 0}$.

First, $\eta^{(1)}(\beta_{s},\bar{x}_{s})=S_{0}^{(1)}+sS_{1}^{(1)}$ and $\eta^{(2)}(\beta_{s},\bar{x}_{s})=S_{0}^{(2)}+sS_{1}^{(2)}$, where

[TABLE]

Let us now expand $F^{(1)}(\beta_{s},\bar{x}_{s}),\dots,F^{(4)}(\beta_{s},\bar{x}_{s})$ as polynomials in $s$, and recall that their first component depend on the exact parameterization of the stable manifold $P_{\beta_{s}}$ which involves infinitely many Taylor coefficients. The work from Section 3 provides the existence of a function $h_{s}:D_{2,\tilde{\nu}}(\mathbb{R}^{2})\to\mathbb{R}^{4}$, see (3.14), such that

[TABLE]

As before, we slightly abuse notation by denoting $P_{\beta_{s}}(\psi)=\overline{P}_{\beta_{s}}(\psi)+h_{s}(\psi)$, where $\theta_{1}+\text{i}\theta_{2}=\rho e^{\text{i}\psi}$ for a fixed $\rho<\tilde{\nu}$. We then split the operator as

[TABLE]

where $F^{(N)}$ denotes the full infinite dimensional $F$ operator but evaluated using the (finitely computable) approximation of the parameterization $\overline{P}_{\beta_{s}}$ of order $N$, and where

[TABLE]

The size of $h_{s}^{(i)}(\psi)$ can be estimated by $r_{m}$ using Proposition 13, where $r_{m}$ is the validation radius for the manifold for $\beta_{0}\leq\beta\leq\beta_{1}$. In addition, since we know that the zeroth order term in $h_{s}^{(i)}$ vanishes, i.e. $a_{0}(s)=\bar{a}_{0}(s)=0$, we obtain a slightly sharper bound for any $\rho<\tilde{\nu}$:

[TABLE]

Hence, we can estimate $H_{s}(\bar{\psi}_{s})$ elementwise by

[TABLE]

Denoting

[TABLE]

with $F^{(j,N)}=\left\{f_{0}^{(j,N)},f_{1}^{(j)},f_{2}^{(j)},f_{3}^{(j)},\dots\right\}$ for $=1,2,3,4$, we rewrite $f_{0}^{(j,N)}(\beta_{s},\bar{x}_{s})$ as a polynomial in $s$, where we use $\Delta\overline{P}$ as defined in (4.14):

[TABLE]

for some $\xi,\zeta$ between $\bar{\psi}_{0}$ and $\bar{\psi}_{1}$ (using the mean value theorem). To obtain an explicit computable expression for $S_{1,0}^{(j+2)}$ and $S_{2,0}^{(j+2)}$, we determine

[TABLE]

by an interval arithmetic calculation, i.e., replacing $\xi$ and $\zeta$ by the interval $[\bar{\psi}_{0},\bar{\psi}_{1}]$.

For $k\geq 1$, we set ($j=1,2,3,4$)

[TABLE]

where the third order term is nonzero for $j=1$ and $j=4$ only. All terms are collected in Table 1. Then, it is possible to write the whole operator as

[TABLE]

where $S_{i}=(S_{i}^{(1)},S_{i}^{(2)},\{S_{i,k}^{(3)}\}_{k\geq 0},\{S_{i,k}^{(4)}\}_{k\geq 0},\{S_{i,k}^{(5)}\}_{k\geq 0},\{S_{i,k}^{(6)}\}_{k\geq 0})$ for $i=0,1,2,3$.

Since we evaluate $F$ using a finite dimensional approximation, $F(\beta_{s},\bar{x}_{s})$ will contain only a finite number of nonzero elements. First, we consider the entries not exceeding the dimension of the finite dimension approximation. This part gets ‘hit’ by $A^{[m]}$, and is bounded component by component:

[TABLE]

where $V=(V^{(1)},V^{(2)},\{V^{(3)}_{k}\}_{k=0}^{m-1},\{V^{(4)}_{k}\}_{k=0}^{m-1},\{V^{(5)}_{k}\}_{k=0}^{m-1},\{V^{(6)}_{k}\}_{k=0}^{m-1})^{T}$. Concerning terms that exceed the dimension of the finite dimensional projection, by using the definition of $A$, one gets for $j=3,4,5,6$

[TABLE]

The expansion of $f_{k}^{(j-2)}(\beta_{s},\bar{x}_{s})$ in powers of $s$ is given by (4.23) with the coefficients in Table 1, We note that all $S^{(j)}_{i,k}$ vanish for $k\geq 2m$. To be precise, $S^{(3)}_{i,k}$ vanishes for $k\geq 2m$, whereas when $j=4,5,6$ then $S^{(j)}_{i,k}$ vanishes already for $k\geq m+1$. Hence we define the estimates ($3\leq j\leq 6$, $m\leq k\leq 2m-1$)

[TABLE]

Having estimated all the terms appearing in the expression $\|(T(\beta_{s},\bar{x}_{s})-\bar{x}_{s})^{(j)}\|_{X_{j}}$, we set

[TABLE]

By construction, we have

[TABLE]

4.6 The $Z$ bound for the homoclinic orbit problem

We recall the definition of the bounds $Z^{(j)}$ in (2.3). In our context, $Z^{(j)}$ is a bound satisfying

[TABLE]

To simplify the manipulations of the expressions appearing in the bounds, we introduce an operator $A^{\dagger}=\{A^{\dagger}_{i,j}|1\leq i,j\leq 6\}$, where the splitting is explained in Definition 19. This operator $A^{\dagger}$ is on the one hand an ‘almost inverse’ of the operator $A$, and on the other hand it approximates $D_{x}F(\beta_{0},\bar{x}_{0})$. We define $A^{\dagger}$ piecewise, where we use the decomposition of the Jacobian $(D_{x}\overline{F}(\beta_{0},\bar{x}_{0}))=(D_{x}\overline{F}(\beta_{0},\bar{x}_{0}))_{i,j}$ into 36 blocks as in (4.16):

• •

$A_{i,j}^{\dagger}\in\mathbb{R}$ for $1\leq i,j\leq 2$, where $A_{i,j}^{\dagger}=(D_{x}\overline{F}(\beta_{0},\bar{x}_{0}))_{i,j}$;

• •

$A_{i,j}^{\dagger}\in(\ell_{\nu}^{1})^{*}$ for $1\leq i\leq 2$ and $3\leq j\leq 6$, where $A_{i,j}^{\dagger}$ is $(D_{x}\overline{F}(\beta_{0},\bar{x}_{0}))_{i,j}$ padded with zeros;

• •

$A_{i,j}^{\dagger}\in\ell_{\nu}^{1}$ for $3\leq i\leq 6$ and $1\leq j\leq 2$, where $A_{i,j}^{\dagger}$ is $(D_{x}\overline{F}(\beta_{0},\bar{x}_{0}))_{i,j}$ padded with zeros;

• •

$A_{i,j}^{\dagger}\in B(\ell_{\nu}^{1},\ell_{\nu^{\prime}}^{1})$ for $3\leq i,j\leq 6$, with $\nu^{\prime}<\nu$, where

[TABLE]

Now, we use $A^{\dagger}$ to perform the splitting

[TABLE]

As in Section 3, the bound on the first term in (4.25) can be directly computed. We set $B=I-AA^{\dagger}$, whose nonzero elements are represented by the finite matrix $I_{4m+2}-A^{[m]}D_{x}\overline{F}(\beta_{0},\bar{x}_{0})$, and we use Lemmas 17 and 18 to derive the bounds

[TABLE]

with the norms introduced in Section 4.2. This provides the desired bound on the first term of (4.25). For the second term, we set $u,v\in B_{1}(0)$ such that $b=ru$ and $c=rv$. We denote $v=(v_{L},v_{\psi},v^{(1)},v^{(2)},v^{(3)},v^{(4)})$, and similarly for $u$, $b$ and $c$. First, for $i=1,2$, we have

[TABLE]

where the final inequality follows from Lemma 17. Next we consider the $k=0$ term of the other four components. For $i=1,2,3,4$ one finds

[TABLE]

where $\zeta_{s}$ is in $[\bar{\psi}_{s}-r,\bar{\psi}_{s}+r]$, and $\xi_{s}$ is in $[\bar{\psi}_{0},\bar{\psi}_{s}]$. A direct computation shows that

[TABLE]

Combining this with (4.22) gives us a bound on the terms in (4.28):

[TABLE]

The remaining terms in (4.29) are estimated using Lemma 14. We obtain, for $i=1,2,3,4$,

[TABLE]

with, for $j=3,4,5,6$,

[TABLE]

For $k\neq 0$, we consider separately the coefficients of $r$, $r^{2}$ and $r^{3}$:

[TABLE]

The term $-A^{\dagger}v$ contributes to the $s$-independent part of $\tilde{z}_{1,k}^{(i)}$ only. Since in $\tilde{z}_{1,k}^{(i)}$ some of the terms involving $(v^{(j)})^{[m]}$ will cancel, it is useful to introduce $\widehat{v}^{(j)}$ as follows:

[TABLE]

Using this notation, for $\tilde{z}_{1,k}^{(3)}$ and $1\leq k\leq m-1$, one finds

[TABLE]

Clearly $\delta_{k,m-1}\widehat{v}_{k\pm 1}^{(2)}=\delta_{k,m-1}v_{k+1}^{(2)}=\widehat{v}_{k\pm 1}^{(2)}$, for $k\leq m-1$, and the Kronecker $\delta_{k,m-1}$ may be viewed as superfluous. For $k\geq m$ we find

[TABLE]

Once again the Kronecker $\delta_{k,m}$ may be viewed as superfluous. For $\tilde{z}_{4,k}^{(1)}$, $\tilde{z}_{5,k}^{(1)}$ and $\tilde{z}_{6,k}^{(1)}$, one finds

[TABLE]

and

[TABLE]

The $\tilde{z}_{i,k}^{(2)}$ and $\tilde{z}_{i,k}^{(3)}$ coefficients are still to be determined. For $k\neq 0$, they are given in Table 2.

Thus, we set $\tilde{z}_{1}^{(i)}=\{\tilde{z}_{1,k}^{(i)}\}_{k\geq 0}$, $\tilde{z}_{2}^{(i)}=\{\tilde{z}_{2,k}^{(i)}\}_{k\geq 0}$ and $\tilde{z}_{3}^{(i)}=\{\tilde{z}_{3,k}^{(i)}\}_{k\geq 0}$. We note that values of $\tilde{z}_{1,0}^{(i)}$ and $\tilde{z}_{2,0}^{(i)}$ are not explicitly given, but (4.30) provides bounds on these terms. We are going to abuse notation by referring to these bounds as $\tilde{z}_{1,0}^{(i)}$ and $\tilde{z}_{2,0}^{(i)}$, where we will correct for this abuse below whenever these terms get involved. We set $\tilde{z}_{3,0}^{(i)}=0$.

For $l=1,2$, one can estimate, using Equation (4.6) and the definition of $\tilde{z}^{(i)}_{j}$,

[TABLE]

and for $l=3,4,5,6$

[TABLE]

Apart from $|A_{l,i}|$ for $l,i=1,2$, which are scalars, it is not immediately obvious how to compute or estimate the terms in (4.32) and (4.33) explicitly. The norms $\|A_{l,i}\|_{\infty,\nu^{-1}}$ for $l=1,2$, $i=3,4,5,6$ and $\|A_{l,i}\|_{1,\nu}$ for $i=1,2$, $l=3,4,5,6$ can be computed directly, since they are represented by row and column vectors of length $m$. The operator norms $\left\|A_{l,i}\right\|_{B(\ell_{\nu}^{1})}$ can be computed using Lemma 18, since for $l\neq i$ they are represented by finite matrices, whereas for $l=i$ they have a decaying diagonal tail (see the analogous Lemma 11).

The norms $\|\tilde{z}^{(i)}_{2}\|_{1,\nu}$ and $\|\tilde{z}^{(i)}_{3}\|_{1,\nu}$ in the quadratic and cubic terms in $r$ can be estimated using the Banach algebra structure. Taking into account the bound on $\tilde{z}^{(i)}_{2,0}$ in (4.30), this leads to bounds

[TABLE]

with

[TABLE]

and

[TABLE]

The factor $\nu+\nu^{-1}$ in the expressions above is due to the shift in index (to the right and to the left) in $u^{(i)}_{k\pm 1}$, $v^{(i)}_{k\pm 1}$, etc.

This leaves us with estimating $|A_{l,i}\tilde{z}_{1}^{(i)}|$ and $\|A_{l,i}\tilde{z}_{1}^{(i)}\|_{1,\nu}$. Since these appear in the terms that are linear in $r$, a direct triangle inequality bound would be too rough for the method to succeed. Hence we estimate these terms more carefully below.

For the term in front of $r$ in equation (4.32), for $l=1,2$, we have

[TABLE]

Here we have corrected for our abuse of notation regarding $\tilde{z}^{(i)}_{1,0}$ by splitting it off using the triangle inequality.

Remark 21**.**

We use the bound (4.6) to estimate the convolution

[TABLE]

A similar estimate leads to

[TABLE]

Some of the terms in $\tilde{z}^{(i)}_{1,k}$ are computable directly, while others need to be estimated. To present these estimates in a structured way we introduce several computable constants. For the convolution terms involving either $v$ or $\widehat{v}$ in $\tilde{z}^{(3)}_{1,k}$ we introduce (for $k\geq 1$)

[TABLE]

Here $\mathcal{Q}_{k}(\cdot)$ and $\hat{\mathcal{Q}}_{k}(\cdot)$, defined in Remark 21, can be computed (at least finitely many of them) since $\bar{x}^{(i)}$ and $\Delta\bar{x}^{(i)}$ have only finitely many nonzero components. We now set, for $k=1,\dots,m-1$,

[TABLE]

as well as

[TABLE]

and

[TABLE]

Recall that $|A|$ denotes the component-wise absolute value. Then we have the computable estimates ($l=1,2$)

[TABLE]

For $l=3,4,5,6$, we split the estimate in three terms because of the way the $\tilde{z}_{i,0}^{(1)}$ bounds and $A$ are defined. Using (4.30), we get ($i,l=3,4,5,6$)

[TABLE]

Again we have dealt with the $\tilde{z}^{(i)}_{1,0}$ terms separately to take into account our abuse of notation.

The final two terms in (4.39) still need to be estimated. The first of these can be estimated in the same way as above, which we write (for $3\leq l,i\leq 6$) compactly as

[TABLE]

with

[TABLE]

where one should read $\hat{\hat{z}}^{(4)}_{k}=\hat{\hat{z}}^{(5)}_{k}=0$.

For the final ‘tail’ terms in (4.39), we bound these as we did for $z^{(i)}_{2}$ and $z^{(i)}_{3}$ coefficients. We obtain

[TABLE]

Therefore, recalling (4.31) and (4.34)–(4.38), for $l=1,2$, we set

[TABLE]

and for $l=3,4,5,6$, we set

[TABLE]

Finally, by construction,

[TABLE]

4.7 Use of the uniform contraction principle

Using the computable bounds $Y^{(l)}$ and $Z^{(l)}$ constructed in the previous two sections, we set

[TABLE]

If we find an $r>0$ such that $p^{(l)}(r)<0$ for all $l=1,\ldots,6$, then according to Proposition 6 we have validated the numerical approximation $\bar{x}_{s}$ of solutions to the BVP (4.4), for every $s\in[0,1]$, and hence we have proven the existence of symmetric homoclinic orbits for all $\beta\in[\beta_{0},\beta_{1}]$.

Proposition 22**.**

For every $s\in[0,1]$, let

[TABLE]

be the approximate solution of (4.4) that we have computed for $\beta=\beta_{s}$, $L=\bar{L}_{s}$ and $\psi=\bar{\psi}_{s}$. Assume that there exists an $r>0$ such that $p^{(l)}(r)<0$ for all $l=1,\ldots,6$. Then, for each $s\in[0,1]$, there exists a solution of (4.4) for $\beta=\beta_{s}$ of the form

[TABLE]

and some $L=L_{s}$ and $\psi=\psi_{s}$ satisfying $|L_{s}-\bar{L}_{s}|\leq r$ and $|\psi_{s}-\bar{\psi}_{s}|\leq r$. This solution corresponds to a (symmetric) homoclinic orbit of (1.3). Furthermore, let

[TABLE]

then we have the following uniform error bound on the (central part of) the homoclinic orbit in phase space: $|g^{(i)}_{s}(t)|\leq r$ for all $t\in[-1,1]$, $s\in[0,1]$ and $i=1,2,3,4$.

Proof.

Proposition 6 yields that, for each $s\in[0,1]$, there exists a unique fixed point $x_{s}$ of $T(\beta_{s},\cdot)$ in the ball of radius $r$ around $\bar{x}_{s}$. The operator $A$ is injective since its non-diagonal part $A^{[m]}$ is invertible. The latter follows from the fact that, see (4.26),

[TABLE]

where the final inequality is implied by $p^{(l)}(r)<0$. Here the operator norm on $X^{[m]}\cong\mathbb{R}^{4m+2}$ is induced by the one on $X$. Hence the fixed point $x_{s}$ of $T$ solves $F(\beta_{s},x_{s})=0$, and by construction $v_{s}$ is a solution of (4.4), which through the change of variables from Section 4.1 corresponds to a homoclinic solution of (1.3). The error bound follows from

[TABLE]

5 Algorithm and results

In this section we discuss some algorithmic issues. In particular, we explain how certain computational constants are chosen and how the two parts of the problem (the manifold computation and the boundary value problem) are joined together to produce the homoclinic orbit. To get the continuation started, the first thing to do is to compute the approximation of the manifold for a fixed value of $\beta$. Since the first coefficients of the parameterization depend on the steady state and the eigenvectors, which are known, one can start Newton’s method with these values in combination with zeros for all higher order coefficients. If Newton’s method does not converge, replacing the starting point with a good approximation for a slightly higher number of Taylor coefficients (which can be computed recursively) will work. Once one good approximation has been found for a particular value of the parameter $\beta$, one can use it as the starting point to find another approximation for sufficiently close values of the parameter.

Another important point for the computations is the size of the manifold that we get. Since the stable eigenvalues of the Jacobian at 0 are complex conjugates, we know that asymptotically the orbit spirals toward the origin. If the manifold we compute is large enough to contain most of the spiraling part, then we do not have to compute that part of the orbit using Chebyshev series, which is advantageous. Generally speaking, the larger the manifold is, the easier the remaining part with Chebyshev will be. Therefore we use the method developed in [1] to maximize the image of the parameterization we compute.

A natural approach to obtain a larger manifold is to try and maximize the $\tilde{\nu}$ for which we can validate the parameterization (we recall that its domain of definition is $D_{2,\tilde{\nu}}(\mathbb{R}^{2})=\left\{\theta\in\mathbb{R}^{2},\ \left|\theta\right|_{2}\leq\tilde{\nu}\right\}$). However, taking $\tilde{\nu}\gg 1$ or $\tilde{\nu}\ll 1$ leads to numerical instabilities (see for instance the quantities $K^{(i,j)}$ defined in Section 3.4.1). The key observation from [1] to avoid this phenomenon is the following. Given a parameterization

[TABLE]

and, for some $\gamma>0$, a rescaled parameterization (also with rescaled eigenvectors)

[TABLE]

the parameterization $P$ on the domain $D_{2,\gamma}(\mathbb{R}^{2})$ defines the same manifold as the rescaled parameterization $\tilde{P}$ on the domain $D_{2,1}(\mathbb{R}^{2})$. Therefore we can fix $\tilde{\nu}$ to be $1$ and instead look for the largest $\gamma$ for which we can validate the rescaled parameterization.

Another useful feature of the results of [1] is that they provide the explicit dependence of the bounds $Y$ and $Z$ with respect to the rescaling $\gamma$, enabling us to recompute bounds for any rescaling cheaply. In practice, we use the following process:

• •

Compute an approximate parameterization $P$ (that is, the coefficient $a_{\alpha}$).

• •

Compute the bounds $Y$ and $Z$ for $\beta_{0}$, without the continuation (i.e. take $\Delta a=0$ and $\Delta\beta=0$ in every estimate).

• •

Find the largest $\gamma$ for which the proof succeeds (i.e. we find an $r>0$ such that $p^{(i)}(r)<0$ for all $i=1,2,3,4$, where the four radii polynomials $p^{(i)}$ are defined in (3.10)) for the rescaled coefficients $\tilde{a}_{\alpha}=\gamma^{|\alpha|}a_{\alpha}$, while requiring the coefficients of the linear term (the one front of $r$) in each radii polynomial $p^{(i)}$ to be less than some threshold $\eta\in(0,1)$, which will be discussed below. This step yields a parameterization $\tilde{P}$ with rigorous error bounds on the domain $D_{2,1}(\mathbb{R}^{2})$.

• •

Use the parameterization $\tilde{P}$ with this $\gamma$ for the Chebyshev part and for the continuation.

Before describing in more detail the process of continuation, let us explain the role of the threshold $\eta$. Finding a positive root of a radii polynomial is impossible if its linear term is not negative, because all its other coefficients are always non-negative by construction. If the linear term is just negative enough for the proof to work at the single parameter value $\beta_{0}$, then $\Delta\beta$ has to be taken extremely small for it to remain negative for the uniform proof, since all bounds become worse monotonically in $|\Delta\beta|$. However, we want to take $\Delta\beta$ as large as possible to reduce the number of steps we have to perform to prove the existence of a symmetric homoclinic orbit for all $\beta\in[0.5,1.9]$. Hence, the addition of this threshold $\eta$ is a trade off: we get a manifold that is a bit smaller than what we could have had optimally, which makes the proof for the Chebyshev part a bit harder, but we can take larger steps in $\beta$, making the total process faster overall. In practice, we use an $\eta$ close to $0.5$ (the value we use varies slightly with $\beta$).

Once the approximation for the manifold is maximized and proven for a particular value of $\beta_{0}$, one can use it as the starting point to find the approximation for $\beta_{1}>\beta_{0}$ in order to compute an approximation for the whole interval $[\beta_{0},\beta_{1}]$. We use the same rescaling $\gamma$ for the entire interval $[\beta_{0},\beta_{1}]$. On the other hand, it is possible to use different scalings for consecutive intervals.

The value of $\Delta\beta=\beta_{1}-\beta_{0}$ that we use is not constant, and varies between $2.5\times 10^{-4}$ and $3.9\times 10^{-6}$. The smaller values are needed when $\beta_{0}\geq 1.8$. This is due to the fact that proof of the stable manifold becomes harder and harder when $\beta$ approaches $2$. Indeed, when $\beta$ goes to $2$ the real part of the stable eigenvalues (see (3.1)) goes to zero, and the problem of finding the stable manifold becomes singular (this can also be seen in the bounds derived in Section 3). Note that when the proof fails for a given interval, a smaller $\Delta\beta$ needs to be used. Thus, the algorithm needs to recompute both the manifold and the orbit for $\beta=\beta_{1}$. However, $A^{\dagger}$ and $A$ need not to be computed again for the new proofs since they both only depend on the approximation at $\beta=\beta_{0}$.

For the manifold all proofs were done using $N=30$ for the dimension of the truncated power series. For the orbit, the proof succeeds with $m=350$ for $[\beta_{0},\beta_{1}]\subset[0.5,1.8]$, and with $m=400$ otherwise. In Figure 1 one can see the profile of the solution for $\beta=0.5$, $\beta=1.2$ and $\beta=1.9$. The left part of the figure shows the decay rate of the solution using the logarithm of the absolute value of the first $50$ Chebyshev coefficients. Recall that the first component of the system is given by $v_{1}=e^{u_{1}}-1$, where $u_{1}$ is the first component of the original system, obtained after transforming the fourth order equation to a first order system. One can see that the solution for $\beta=0.5$ is really close to $-1$ for a much longer period of time than the other solutions depicted. This behaviour has an impact on the decay of the corresponding Chebyshev series. Moreover, another value affecting the decay rate of the solution is the time rescaling factor $L$ of the orbit. For $\beta=0.5$ (respectively $\beta=1.2$ and $\beta=1.9$) we have $L\approx 3.1312$ (respectively $L\approx 1.7671$ and $L\approx 2.6170$). The first three components of the solution and the local manifold can be seen in Figure 2, Figure 3 and Figure 4 for $\beta=0.5$, $\beta=1.2$ and $\beta=1.9$, respectively. The profiles of the first component $v^{(1)}$ of these three solutions can be compared in Figure 5, where half the symmetric homoclinic orbits is depicted. Furthermore, the three corresponding homoclinic solutions of the suspension bridge equation (1.2) in the original $u$-variable are presented in Figure 6.

Finally, to perform the proof successfully for the entire interval range $\beta\in[0.5,1.9]$ we had to execute the algorithm $7960$ times. Each proof took between $7$ and $10$ seconds on a laptop with an Intel Core i7 4500U processor on MATLAB R2016a. The code which was used to perform the proofs is available at [23] and uses the interval arithmetic package INTLAB [19].

Acknowledgements

Jan Bouwe van den Berg was supported by the grant NWO Vici-grant 016.123.606. Jean-Philippe Lessard was supported by Gouvernement du Canada/Natural Sciences and Engineering Research Council of Canada (NSERC), CG100747.