A numerical verification method to specify homoclinic orbits as application of local Lyapunov functions

We propose a verification method for specification of homoclinic orbits as application of our previous work for constructing local Lyapunov functions by verified numerics. Our goal is to specify parameters appeared in the given systems of ordinary differential equations (ODEs) which admit homoclinic orbits to equilibria. Here we restrict ourselves to cases that each equilibrium is independent of parameters. The feature of our methods consists of Lyapunov functions, integration of ODEs by verified numerics, and Brouwer’s coincidence theorem on continuous mappings. Several techniques for constructing continuous mappings from a domain of parameter vectors to a region of the phase space are shown. We present numerical examples for problems in 3 and 4-dimensional cases.


Introduction
Recently numerical verification has been widely applied to studying dynamical systems together with a lot of applications (e.g. [2, 3, 5, 7-9, 12, 16, 24-27, 29]). The second and the third authors have also proposed methods for local Lyapunov functions constructed by numerical verification as tools of analysis of dynamical systems [16]. 1 3 In study of Lyapunov functions, it is a usual manner to construct global Lyapunov functions or Lyapunov functions with wide domains, not only in theoretical studies but also in studies applying numerical verification [8,9]. Nevertheless, combining with numerical verification methods, we assert local Lyapunov functions with small domains are important tools to analyze dynamical systems according to the following reasons.
-It is possible to construct local Lyapunov functions for a wide class of dynamical systems which have hyperbolic equilibria or fixed points. -Domains of the local Lyapunov functions can be specified by numerical verification. -There are several numerical verification methods to specify trajectories in continuous dynamical systems and sequences in discrete dynamical systems, e.g. Lohner method, Taylor model and power series arithmetic [10,11,13,14,20]. It is possible to analyze global dynamics by information of trajectories or sequences together with some features around equilibria or fixed points obtained by local Lyapunov functions.
In order to show usefulness of local Lyapunov functions as tools of analysis, we propose a verification method for specification of homoclinic orbits based on combination of verified numerics and local Lyapunov functions. Early works on connecting orbits including homoclinic orbits for systems of ordinary differential equations (ODEs) by verified numerics are those of Zgliczyński and/or Wilczak [26,27], and many successive works are provided. On the other hand, Arioli, Koch, Barker, Oishi etc. gave some results on existence proof of connecting orbits and their stability [2,3,19]. Around 2011, Lessard and others succeeded in capturing connecting orbits by using analytical methods, which lead to lots of results including bifurcation of connecting orbits and further topics on partial differential equations and delay differential equations [5,7,12].
These methods are roughly sorted into two ways, topological methods and analytical methods. A representative one of topological methods is covering relation by Arioli [4] including application of mapping degree. As a typical one of analytical methods, we mention the method using parameterization and radii polynomials by Lessard and others where contraction mappings are explicitly constructed. From this point of view, our proposing method in this paper can be regarded to have the both aspects since mapping degree is adopted as a topological method and Lyapunov function is explicitly constructed in an analytical way.
Our method on proving existence of connecting orbits consists of Lyapunov functions, integration of ODEs by verified numerics, and Brouwer's coincidence theorem on continuous mappings. Some techniques in order to construct a continuous mapping from a domain of parameter vectors to a region of phase space are presented. Once such a continuous mapping is specified, one can apply both of Brouwer's fixed point theorem and Brouwer's coincidence theorem to proving the existence of a parameter vector admitting a connecting orbit. In our case we adopt Brouwer's coincidence theorem with calculation of degree of a continuous Numerical verification method for homoclinic orbits mapping by verified numerics. These ideas are based on Master's thesis of Yamano [28] supervised by the second author and advised by the third author.
We deal with autonomous systems of ODEs with a parameter vector p: where D is a closed domain in ℝ n− with a positive integer which is mentioned below. In the whole arguments, we suppose the following : -For each p ∈ D , (1.1) has a hyperbolic equilibrium x * (p) that is a saddle with (n − )-dimensional unstable manifold and -dimensional stable manifold for every p ∈ D. -The equilibrium x * (p) depends continuously on p.
In this article, we restrict ourselves to the system (1.1) with an equilibrium x * which is independent of the parameter vector p . This simplification permits us to avoid complicated mathematical arguments and to explain our method more understandable, though we believe that our method can be extended to the cases where equilibria depend on p.
Let ∶ ℝ × ℝ n × ℝ n− → ℝ n be a flow generated by (1.1), namely (t, x; p) is the solution of (1.1) with an initial point x ∈ ℝ n for a parameter vector p ∈ ℝ n− . The definition of our Lyapunov functions is as follows. We call the domain D L a Lyapunov domain.
By virtue of the numerical verification methods proposed in [16], a Lyapunov function can be defined for each fixed p with a local domain including only one equilibrium x * , and takes both positive and negative values reflecting the property that x * is a saddle point.
Fix a parameter vector p ∈ D . Since x * is a hyperbolic equilibrium of (1.1), we can construct a local Lyapunov function in a quadratic form within a domain D L (p) which is specified by numerical verification [16]. We calculate Y(p) ∈ ℝ n×n as a real symmetric matrix from the information of eigenvalues of the Jacobian matrix Df * (p) = D x f (x * , p) with respect to x at the equilibrium x * .
(1.1) dx dt = f (x;p), x ∈ ℝ n , p ∈ D, t ∈ ℝ, f ∈ C r (ℝ n × D, ℝ n ), r ≥ 1, Our method employs the local Lyapunov functions with relatively small domains in order to construct a continuous mapping F from the parameter set D to the phase space and apply Brouwer's coincidence theorem. The rest of the present paper is organized as follows. We show an outline of our method in the next section and define the mappings F 1 and F 2 which construct F. In Sect. 3, construction of Lyapunov functions and its validity based on [16] are explained. Note that the method in [16] is slightly modified since our ODEs have parameters. We state Brouwer's coincidence theorem and verification methods to compute mapping degree in Sect. 4. In Sect. 5, two lemmas on properties of our Lyapunov functions and one theorem concerning continuous mappings appeared in Sect. 8, which are proved in preparation of mathematical arguments in the successive sections. We consider a mapping called Lyapunov Tracing and prove its continuity in Sect. 6, which plays an important role in our framework. In Sect. 7, our numerical procedure with verified computation to specify homoclinic orbits is shown, which precedes mathematical arguments of construction of the continuous mappings F 1 , F 2 and F. We prove that F 1 is well-defined and continuous in Sect. 8. Unstable manifold theorem plays a central role in the proof. In Sect. 9, we describe how to construct the continuous mapping F 2 and F, and state the reason why our method verifies the parameters which admits a homoclinic orbit. Some results of numerical examples are shown in Sect. 10 including a system of 3-dimensional ODEs for illustration of our method and a 4-dimensional system of ODEs admitting the 2-dimensional unstable manifold and the 2-dimensional stable manifold of an equilibrium.

Outline of our method
We describe an outline of our approach before detailed explanation for verification of homoclinic orbits. Let p ∈ D be a parameter vector which is calculated as an approximation to a parameter vector admitting a homoclinic orbit. Take a small closed domain D p ⊂ D including p as an interior point. Note that D p is taken to be homeomorphic to a closed ball in ℝ n− in order to apply Brouwer's coincidence theorem.
Our goal is to prove that there exists a parameter vector in D p which admits a homoclinic orbit. First we construct Lyapunov functions with respect to x * . The domain D p should be sufficiently small so that we can take a common matrix Y ∈ ℝ n×n as Y(p) independently of p ∈ D p . In particular, our Lyapunov function has the form of Note that the form of the Lyapunov function does no longer depend on the parameter vector p , but the domain of the Lyapunov function is still depending on p.
As is mentioned in Sect. 3, the domain D L (p) of L(x) for each p can be defined to be compact and convex. Moreover we define the compact and convex set D L by Numerical verification method for homoclinic orbits which is called a Lyapunov domain.
The subsets L 0 , L + and L − of D L are defined as follows: In order to apply Brouwer's coincidence theorem in our frame work, we take a hyperplane Γ ⊂ ℝ n− and define the projection P ∶ ℝ n → Γ . The projection P is supposed to satisfy the following: holds, then x = x * . Construction of the projection P with the above feature is described in Sect. 7.
Let W s (x * ; p) and W u (x * ; p) be the stable and the unstable manifold for x * , respectively, which depend on a parameter vector . The set B 0 should be defined in order that there exists a certain T > 0 to attain (T, x; p) ∈ L + for any x ∈ B 0 and p ∈ D p . The mapping F 2 (x, p) is defined by using the flow y = (T, x; p) with the domain B 0 and (t 0 (y), y; p) ∈ L 0 with a positive time t 0 (y) for each y ∈ L + . Thus F 2 (x, p) ∈ L 0 holds for every x ∈ B 0 . Details are described in Sects. 7, 8 and 9. Now we define the continuous mappings F, H ∶ D p → Γ by and apply Brouwer's coincidence theorem with numerical verification of the mapping degree of F. This proves that F and H have a coincidence point, that is, F(p * ) = H(p * ) ∈ Γ holds for some p * ∈ D p . We obtain F 2 •F 1 (p * ) = x * using the feature of the projection P, which implies that there exits a flow from f 1 (p * ) ∈ W u (x * ; p * ) satisfying Namely there exists a parameter vector p = p * ∈ D p which admits a homoclinic orbit. The expression (2.1) indicates that difficultly of the computation may be caused by taking the limit t → ∞ in the definition of the continuous mapping F 2 . In order to avoid the difficulty, we adopt a tool called Lyapunov Tracing explained in Sect. 6.

Remark 2.1
There is a possibility to adopt some fixed point theorems to prove the existence of p * instead of the coincidence theorem. For example, apply a nonlinear solver of verified numerics based on a fixed point theorem to a nonlinear equation . However, there might be a difficulty in computation of F(p) for p in a neighborhood of p * which makes the verification process fail since we need long time integration to get the image F(p) for p very close to p * , which causes so much computational error. Note that we can avoid such difficulty around p * by using the coincidence theorem and Lyapunov Tracing.

Lyapunov functions
In this section, we show the procedure to construct Lyapunov functions and to specify the domains of Lyapunov functions based on [16].

Construction procedure of Lyapunov functions
Consider the autonomous system (1.1). For each p ∈ D , it is possible to construct a quadratic form as a Lyapunov function by the following steps.
1. Let Df * (p) be the Jacobian matrix of f (x; p) with respect to x at x = x * . For simplicity, assume that Df * (p) is diagonalizable by a nonsingular matrix X(p) , which is generically valid, and we have where Λ(p) =diag( 1 (p), 2 (p), … , n (p)) with eigenvalues k (p) . Note that Λ(p) and X(p) are sufficient to be calculated by floating-point arithmetic.
Note that Re( k (p) ) is not 0 for k = 1, … , n since each x * is assumed to be hyperbolic. 3. Compute a real symmetric matrix Y(p) = (Y ij (p)) n i,j=1 as follows:

3
Numerical verification method for homoclinic orbits where X −H (p) denotes the inverse matrix of the Hermitian transpose X H (p) of X(p) , which is also sufficient to be calculated by floating-point arithmetic. 4. Define the quadratic form by which is our candidate of Lyapunov function around x * . In our verification methods for homoclinic orbits, we take p ∈ D which admits an approximate homoclinic orbit and put Y ∶= Y(p) . If a parameter domain D p is a sufficiently small neighborhood of p , the quadratic form is also a candidate of Lyapunov function for each p ∈ D p . Then we define the domain D p by a closed ball centered at p with a small radius so that the above L(x) is expected to be the Lyapunov function with a common Y.

Verification of Lyapunov domains
We define the Lyapunov domain as follows. In order to ensure that D L (p) is a Lyapunov domain, a verification procedure consisting of two stages is proposed in [16]. The following steps are referred to as Stage 1 in [16], while the present argument contains the parameter dependence of objects.
-Define a matrix A(x; p) for x ∈ D L (p) , p ∈ D p by Divide the domain D L (p) into some subdomains D k , k = 1, 2, … , K and enclose each D k by an interval vector [x k ] ⊂ ℝ n . Then verify the negative definiteness of an interval matrix A([x k ]; p) . An example of such a procedure is the following: ; p)X k by interval arithmetic with verified numerics. Let [a k ] ij denote the (i, j)th entry of the matrix, which takes an interval value. 3. Verify with interval arithmetic for i = 1, 2, … , n , and apply the Gershgorin Circle theorem to the interval matrix. Note that if the above verification is succeeded for all k = 1, 2, … , K , then Define the compact and convex set D L by We also call D L a Lyapunov domain. In actual computation, the parameter vector p ∈ D p is treated by an interval vector [p] ⊂ D p with a small radius.

Brouwer's coincidence theorem and verification methods for mapping degree
We apply Brouwer's coincidence theorem to proving the existence of homoclinic orbits, and numerical verification of mapping degree of continuous mappings is the essential technique in order to apply the theorem. In this section, the definition of the mapping degree and two theorems about the mapping degree applied in our verification method are stated.
The following theorem is known as Brouwer's coincidence theorem.
Theorem 4.1 (e.g. [23]) Let B n be a unit ball in ℝ n , and S n−1 be the surface of B n . Consider a continuous mapping F : B n → B n and let F| S n−1 be the restriction of F to S n−1 . The mapping degree of F is denoted by degF.
If F| S n−1 ⊂ S n−1 and degF| S n−1 ≠ 0 holds, then for an arbitrary continuous mapping G : B n → B n , there exists a point x ∈ B n such that F(x) = G(x).
When n = 2 , the mapping degree can be defined as follows according to [23].
The mapping degree of f is then defined as degf =f (1) −f (0).
The following theorem is called the Interval Simplex Theorem which gives a verification process to check the mapping degree not to be 0 in case of n = 2 [28].
Proof Take t 0 , t 1 , t 2 as From the definition of degree of mapping, if degf = k ∈ ℤ , then Note that, e(t 0 ) = e(t 0 + k) (= e•f (1)) , and holds. From these relations, each condition of This fact implies that the order of t 0 , t 1 , t 2 should be or Now assume that k = 0 . Since e(t 0 ) = e•f (1) , we have for (4.1) and for (4.2). Therefore the case of k = 0 is in conflict with Similarly, the relations t 0 ≤ t 1 ≤ t 0 + k ≤ t 2 and t 0 ≤ t 0 + k ≤ t 1 ≤ t 2 lead to contradiction since and holds, respectively. Therefore, t 0 ≤ t 1 ≤ t 2 ≤ t 0 + k is the only possible relation com- On the other hand, if k ≥ 2 then t 0 Similar arguments hold for the cases of The assumptions of Theorem 4.3 can be checked by interval arithmetic with verified computation which is used in our numerical examples. Note that we do not have to specify the mapping f in the verification process. In the case of | deg f | > 1 or n > 2 , Aberth's method [1] can be applied as a numerical verification method to the mapping degree though it is more complicated than the Interval Simplex Theorem.
Numerical verification method for homoclinic orbits

Mathematical preliminaries
In this section, we show two lemmas on properties of our Lyapunov functions and one theorem about continuous mappings which are used in Sect. 8.
holds. ◻ Remark 5.2 Lemma 5.1 leads to the cone condition which is described in [30]. Here Q is a quadratic form and N is a given set with special properties called an h-set. Indeed, regarding Q and N as −L and D L respectively, the cone condition is satisfied from Lemma 5.1.

holds.
Proof Putting y = x 2 − x 1 + x * and we shall prove holds for any negative time t < 0 . Therefore, Lemma 5.1 ensures for any t < 0 . Since x 1 , x 2 are points on W u (x * ; p) , we have and Now assume that L(y) is 0, then some negative time t 1 < 0 exists such that holds. Let y 1 = (t 1 , x 2 ; p) − (t 1 , x 1 ; p) + x * and apply similar discussion to the above for y = y 1 , we have which is in conflict with L(y 1 ) > 0 . Therefore,

3
Numerical verification method for homoclinic orbits

holds. ◻
The following theorem is used in Sect. 8 to prove continuous dependence of the intersection point of a -dimensional ball and W u (x * ; p) with respect to p ∈ D p . Theorem 5.4 Consider bounded closed sets D ⊂ ℝ n and P ⊂ ℝ m and continuous function f ∶ D × P → ℝ n . Assume that for each p ∈ P there exists a unique Proof Consider an arbitrary sequence {p j } which satisfies p j → p (j → ∞) for each fixed p ∈ P . From the assumption, x(p j ) ∈ D uniquely exists for each p j . We should Note that the zero set of the function f is a bounded closed set. Then, the sequence {x(p j )} has convergent subsequences by virtue of Bolzano-Weierstrass theorem.

Lyapunov tracing: definition and properties
Lyapunov Tracing we shall use later is a mapping from a subset of Lyapunov domain to the surface L 0 = x ∈ D L | L(x) = 0 . It is defined using the flow (t, x; p) and the time T such that (T, x; p) ∈ L 0 holds. In this section, we refer to the arguments in [17] and show that Lyapunov Tracing is continuous with respect to x and p under the following assumptions: Assumption 6.1 There exist a closed ball B + ⊂ D L including x * as an interior point and either of the following statements holds for any p ∈ D p and any x ∈ B + ∩ L + .
1. There is a finite positive time T 0 such that (t, x; p) ∈ B + ∩ L + for all t ∈ [0, T 0 ) and Hereafter let T 0 (x, p) denote T 0 for (x, p) . Note that we treat the case when T 0 is positive in this section, but similar arguments hold for the case when T 0 is negative. See Sect. 8.

Remark 6.2
Roughly speaking, Lyapunov Tracing defined below is the mapping from points on a set onto the boundary (resp. a subset) through trajectories where the flow immediately (resp. eventually) leaves the set. Continuity of such mappings are well understood in the context of isolating blocks or index pairs (e.g. [6,22]), in which case the boundary or the subset mentioned above, so-called the exit of blocks, does contain neither inner tangential points of trajectories nor invariant sets. Continuity of mappings onto the exit in such situations is characterized by means of strong deformation retracts and the idea is based on the implicit construction of Lyapunov functions of invariant sets. The spirit of the present argument is similar to the above situation, while the exit-like boundary including equilibria is our present interest. In particular, if x is on stable or unstable manifolds of equilibria on the exit, the passage time T takes ±∞ , which is a nontrivial point for guaranteeing the continuity. Nevertheless, explicitly validated Lyapunov functions around equilibria can be applied to the continuity of mappings attaining the boundary points through solution trajectories. Here the mapping onto the boundary through trajectories is defined by means of level sets of Lyapunov functions, which is the origin of the name Lyapunov Tracing, and properties of Lyapunov functions is used to prove the continuity. When the parameter vector p ∈ D p is fixed, the continuity with respect to x is proved in [17].

Setting and the definition
We consider the autonomous systems (1.1) and assume that Assumptions 1.1 and 6.1 hold. Fix p ∈ D p and x ∈ B + ∩ L + . By definition of L + , L(x) is a positive value. By Assumption 6.1 and monotonicity of Lyapunov functions, there exist a unique time t ≥ 0 for any L ∈ (0, L(x)] such that L = L( (t, x; p)) holds. Moreover by virtue of the implicit function theorem, we regard the time t as a differentiable function of L and write (t, x; p) = (t(L), x; p).
Define the mapping R ∶ ( Note that holds for any (t, x; p) ∈ D L �{x * } . As long as L ∈ ( 0 , L(x) ] there is a positive time t > 0 with L = L( (t , x; p)) such that Numerical verification method for homoclinic orbits We write the last line of (6.2) as (t(L ), x; p). If x is outside W s (x * ; p) , then the limit of (6.2) exists since , the limit of (6.2) also exists since holds. Summarizing, for each (x, p) ∈ (B + ∩ L + ) × D p , the mapping R given in (6.1) is well-defined.
Continuity of Lyapunov Tracing is discussed in the next subsection. Our proof is based on Theorem 4.1 in [17] and consists of two parts.
Note that T l gives a continuous function T l ∶ B l × D p → ℝ + and that R l (x, p) is continuous with respect to (x, p) ∈ B l × D p , which are proved in Lemma 6.5 in the following.

Lemma 6.5 The function T l (x, p) is continuous on
Proof Take an arbitrary point (x ∞ , p ∞ ) ∈ B l × D p and any sequence of points (x i , p i ) ∈ B l × D p , which converges to (x ∞ , p ∞ ) . Let us divide the sequence into two subsequences. holds.
We prove this lemma using the continuous dependence of solutions of initial value problems for ordinary differential equations on initial values and parameters, namely for any fixed T ∈ ℝ, holds. From the boundedness of T l (x j k , p j k ) , the upper limit of T l (x j k , p j k ) and the lower limit of T l (x j k , p j k ) take finite values and . From the property of lower limit, there exists a subsequence On the other hand we have for m → ∞ by (6.3). These yield It is obvious that the continuity of T l ensures continuity of Similarly continuity of T 0 (x, p) leads to continuity of R on If we fix p ∈ D p the continuity of T 0 (x, p) with respect to x ∈ (B + ∩ L + )�W s (x * ; p) can be proved in an almost same manner as in Lemma 6.5 noting that the integrand in the right-hand side of (6.4) is still finite for l = 0 . In order to prove the continuity of T 0 (x, p) with respect to (x, p) ∈ (B + ∩ L + )� ⋃ p∈D p W s (x * ; p) × D p , we need the following lemma.

Lemma 6.6
For an arbitrary p ∈ D p and an arbitrary x ∈ (L + ∩ B + )�W s (x * ; p) there exist r > 0 and > 0 such that x � ∉ W s (x * ; p � ) holds for any x ′ and p ′ within Proof We derive contradiction assuming that, for any r > 0 and > 0 , there exist x � ∈ (L + ∩ B + ) and p � ∈ D p such that holds. Since the equilibrium point x * is a saddle, there are the stable eigenspace E = span{u 1 , … , u } and the unstable eigenspace E n− = span{u +1 , … , u n } of Df * (p) . Note that ℝ n = E ⊕ E n− holds since x * is hyperbolic, which indicates the n × n matrix U = (u 1 , u 2 , … , u n ) is invertible. We can write an arbitrary x ∈ ℝ n as [15,30]) such that holds for an arbitrary q ∈ D p . Moreover h s is Lipschitz continuous with respect to and q , respectively. Write x = ( , ) and x � = ( � , � ) . From the assumptions, ≠ h s ( , p) and Take r in (6.5) so that r < 1 4 0 holds, and we have On the other hand, from the continuity of h s , there exist r ′ > 0 and ′ > 0 for 1 4 0 such that holds. Take r and in (6.5) to be small enough so that (x, p) and (x � , p � ) satisfy (6.8). Then the left-hand side of (6.7) is estimated from below as which yields a contradiction. ◻ Proof Take an arbitrary p ∞ ∈ D p and x ∞ ∈ (B + ∩ L + )�W s (x * ; p ∞ ) . From the Lemma 6.6 there exist r > 0 and > 0 for p ∞ and Take any sequence of points (x i , p i ) ∈ B r, , i = 1, 2, … , which converges to (x ∞ , p ∞ ) . Since x i ∉ W s (x * ; p i ) , the integrand in the right-hand side of (6.4) is finite with l = 0 for each (x i , p i ) ∈ B r, . Then we can take similar argument to the proof of Lemma 6.5 and obtain the conclusion. ◻

3
Numerical verification method for homoclinic orbits Lemma 6.7 lead to the continuity of R(x, p) on (B + ∩ L + )� ⋃ p∈D p W s (x * ; p) × D p by the definition, namely we have the following proposition.

Continuity of R on
For a given r > 0 , consider a compact set D r defined by Also define From the property of Lyapunov functions, , x; p)) is negative within D r and f (x; p) ≠ 0 since x ∈ D r . Then (r) takes a positive value by the continuity of ∇ x L(x) ⋅ f (x; p) with respect to p and the compactness of D p . From the definition we have as long as (t, x; p) ∈ D r . Moreover the following lemma holds. Lemma 6.9 Fix p ∈ D p and take any l ∈ (0, L max ) and x ∈ X l . Let r = ‖x − x * ‖ and assume that (t, x) ∈ D r for 0 ≤ t ≤ T 0 (x, p) holds. Then we have Proof Note that x ∈ B + �W s (x * ; p) from the assumption (t, x; p) ∈ D r for 0 ≤ t ≤ T 0 (x, p) , which means T 0 (x, p) takes a finite positive value. Then using (6.9) we have , x; p); p)dt which implies (6.10). ◻ Now we prove the following proposition. 17]). Moreover W s (x * ; p) is represented as the graph of a Lipschitz function h s ( , p) . Then p)) is bounded and closed since B + and D p are compact. Therefore, ; p)) . Take an arbitrary (z, p) ∈ (B + ∩ W s (x * ; p)) × D p and an arbitrary sequence

Proposition 6.10 The mapping R(x, p) is continuous on
Let Since B + ∩ L 0 is bounded and closed, the sequence {y i } has a subsequence {y i j } ∞ j=1 which converges to y * ∈ B + ∩ L 0 . We prove the contradiction assuming that y * ≠ x * . Define a positive value r by r ∶= 1 2 ‖y * − x * ‖ . Take a sufficiently large finite time t l 0 such that ‖ (t l 0 , z; p) − x * ‖ < r holds, which is possible since z ∈ W s (x * ; p) . Let l 0 = L( (t l 0 , z; p); p) and note R l 0 (z, p) = (t l 0 , z; p) . By virtue of the continuity of R l 0 (z, p) from Lemma 6.5, ‖R l 0 (x i j , p i j ) − x * ‖ < r holds for sufficiently large j. Then take some large J such that ‖R l 0 (x i j , p i j ) − x * ‖ < r and ‖y * − y i j ‖ < 1 2 r hold for any j ≥ J. Let Let x l ∶= (t l , x l 0 ; p i j ) and l ∶= L(x l , p i j ) < l 0 . Applying Lemma 6.9 to x = x l , we have and This inequality and ‖y * − x * ‖ = 2r imply On the other hand, we can take an arbitrarily large time t l 0 while r is fixed. Therefore l 0 and l < l 0 can be taken arbitrarily small and it is possible to have l (r) < 1 2 r . This contradicts (6.11). 2 ◻ Numerical verification method for homoclinic orbits Theorem 6.4 follows from Propositions 6.8 and 6.10.

Numerical procedure
Before describing mathematical arguments for the existence of homoclinic orbits, we explain our numerical procedures so that readers can use our method as a tool for analysis of dynamical systems. Several figures to illustrate our procedures are shown in Sect. 10 with numerical examples. Through the procedures we use verified computation as necessary together with interval arithmetic [18]. The first procedure shows how to specify the points on the unstable manifolds and gives sufficient conditions to construct the continuous mapping F 1 mentioned in Sect. 1.
Computer programs for the verified computation of eigenvectors can be found in INTLAB, for example [21]. Put matrices and interval matrices Note that the positive definiteness of (A + ) T YA + and the negative definiteness of (A − ) T

YA − are ensured by numerical verification of [A + ] T Y[A + ] and [A − ] T Y[A
Numerical verification method for homoclinic orbits Fig. 6 for 3-dimensional case and Fig. 10 for 4-dimensional case in Sect. 10. 9. Verify that deg(P − • (T 1 (x b , p), x b ; p) ) with respect to x b ∈ B 0 (x, r) is not 0 for each fixed p ∈ D p by using Theorem 4.3 or Aberth's method and interval arithmetic with respect to [P − ] , [b] k and [D p ] . Then the ball B 0 (x, r) includes a point of W u (x * , p) for each p ∈ D p , which is proved in the next section. Note: Trial and error may be necessary to take an appropriate value of r in Step 6 in order to satisfy the conditions in Steps 7, 8 and 9.
We confirm sufficient conditions to construct the continuous mappings F 2 and F by Procedure 2 as follows, where F 2 and F are appeared in Sect. 1.

Divide the boundary D p of D p into subdomains included by interval vectors
[d] k ⊂ ℝ n , k = 1, 2, … , K � . Verify that there exists a positive time   Fig. 9 for 3-dimensional case and Fig. 11 for 4-dimensional case in Sect. 10. 6. Verify that deg(P + • (T 3 (x, p b ), x; p b )) with respect to p b ∈ D p is not 0 for each fixed x ∈ B 0 (x, r) using Theorem 4.3 or Aberth's method and interval arithmetic with respect to [P + ] , [B 0 (x, r)] and [d] k . Then D p includes a parameter vector p * ∈ D p which admits a homoclinic orbit from x * to itself.

Construction of the mapping F 1
To see that the above numerical procedures definitely provide the existence of a homoclinic orbit, we begin with constructing the mapping where F 1 (p) = (f 1 (p), p) with f 1 (p) ∈ W u (x * ; p) which are proved to be well-defined from Theorems 8.4 and 8.6. Suppose that the Lyapunov function L(x) is defined and the Lyapunov domain D L is specified. Note that for each fixed p ∈ D p , the point x * is the unique equilibrium of (1.1) in D L [16]. Furthermore, we assume the following assumption which is the counterpart of Assumption 6.1 in the unstable direction.  (x, p).
Assumption 8.1 should be checked by numerical verification. Note that Step 4 of Procedure 1 ensures Assumption 8.1.
The following steps prepare the conditions to specify the points of W u (x * ; p) for p ∈ D p .

3
Numerical verification method for homoclinic orbits Take a -dimensional ball B 0 (x, r) = {x ∈ Σ � ‖x −x‖ ≤ r} on Σ ∩ L − , which is expected to have an intersection point with W u (x * ; p) for any p ∈ D p . These are specified in Steps 5 and 6 of Procedure 1.
is a time such that (T 1 (x, p), x; p) ∈ L 0 (p) holds. 3 Note that T 1 (x b , p) is estimated for each x b ∈ B 0 (x, r) and p ∈ D p in Step 6 of Procedure 1.

Remark 8.3
Note that Q − is a retraction from Γ − to B − (y * , ) which clarifies correspondence between our method and Brouwer's coincidence theorem. Since Q − is introduced from theoretical reason, one can omit the construction of Q − in practical verification process. See examples in Sect. 10. Now we derive the conditions to prove that the ball B 0 (x, r) includes points of W u (x * ; p). Fig. 1), and assume that the following conditions hold. Combining this with the definition of Q − and Lemma 8.2, we have namely,

Numerical verification method for homoclinic orbits
This implies that x 0 is a point of W u (x * ; p) . Therefore, B 0 (x, r) includes the intersection point x 0 of Σ − and W u (x * ; p) for the fixed parameter vector p ∈ D p . ◻ To prove continuity of the mapping F 1 , we apply Theorem 5.4. In Sect. 6, we have defined the stable eigenspace E = span{u 1 , … , u } and the unstable eigenspace E n− = span{u +1 , … , u n } of Df * (p) , and pointed out that the n × n matrix U = (u 1 , u 2 , … , u n ) is invertible. By the similar discussion to Sect. 6, a Lipschitz function h u ∶ E n− × D p → E exists and the unstable manifold W u (x * ; p) is represented as Let be an (n − ) × n matrix. Then we can write = AU −1 x. We define the continuous mappings g ∶ ℝ n × ℝ n− → ℝ n by For x ∈ B 0 (x, r) , if g(x, p) = 0 , then x is in W u (x * ; p) . By virtue of Theorem 8.6, an intersection point of B 0 (x, r) and W u (x * ; p) is unique for each p ∈ D p . Therefore,

3
Numerical verification method for homoclinic orbits by taking D p , B 0 (x, r) and g for P, D and f respectively in Theorem 5.4, it is proved that intersection point of B 0 (x, r) and W u (x * ; p) is continuous with respect to p ∈ D p . This leads to the continuity of F 1 on D p .

Verification method of homoclinic orbits
Let us prove the existence of homoclinic orbits supposing that F 1 is defined in the manner of the previous section. We assume Assumption 6.1, which is verified by Step 1 in Procedure 2, and the following in addition.
Here B + is the set appeared in Assumption 6.1. Note that Step 2 of Procedure 2 ensures Assumption 9.1.
We construct the continuous mappings F 2 and F by the following steps.
where z−z * ‖z−z * ‖ is an intersection point of a line segment between z * and z ∈ Γ + and the boundary B + (z * , ) of B + (z * , ) . See the note of the step (v) in Sect. 8. Now we define the continuous mapping F as follows.

Definition 9.3
The continuous mapping F ∶ D p → B + (z * , ) ⊂ Γ + is defined by By virtue of the following theorem, we can specify the parameter vectors p ∈ D p which admit homoclinic orbits. Theorem 9.4 Let F S be the restriction of F to the boundary D p of D p (see Fig. 2) and assume that the following conditions hold: Then D p includes a parameter vector p * which admits a homoclinic orbit.
Proof From the assumptions of the theorem and Brouwer's coincidence theorem, there exists a point p * ∈ D p for an arbitrary continuous mapping H ∶ D p → B + (z * , ) such that holds. When we take H as H(p) = P + (x * ) , the point p * satisfies the following equation: Combining this with the definition of Q + and Lemma 9.2, we have From the definitions of F 1 and F 2 , this indicates Therefore p * is a parameter vector which admits a homoclinic orbit since f 1 (p * ) is a point on W u (x * ; p * ) and the left-hand side of the above expression is the integration along the trajectory from the initial point f 1 (p * ) .
Proof Put z * = P + (x * ) and consider an (n − )-dimensional ball B + (z * , ) on Γ + as B + (z * , ) = {z ∈ Γ + � ‖z − z * ‖ ≤ } with a positive constant . From the conditions of Step 5, it is possible to take a sufficiently small > 0 such that the set ⋃ p∈ D p P + •F 2 •F 1 (p) exists in the exterior of the ball B + (z * , ) for any x ∈ B(x, r) . This leads to from the definition of Q + . Moreover Step 6 of Procedure 2 verifies and all the conditions of Theorem 9.4 hold. ◻ Corollaries 8.5 and 9.5 ensure that our numerical method described in Procedures 1 and 2 leads to the proof of the existence of a homoclinic orbit.

Numerical example
In this section, we show results of two numerical examples in order to illustrate our verification method.

A 3-dimensional example
We treat modified Rössler equation: Note that -the origin is an equilibrium for any coefficients a, b, and c.
-analyzing eigenvalues of Df * , the Jacobian matrix of the right-hand side with respect to (x, y, z) at the origin, it is verified that there are 1-dimensional stable manifolds and 2-dimensional unstable manifolds when parameters take values in the following region: We consider parameters within these regions and take = 1 and n − = 2. We construct a continuous mapping F 1 by the following steps corresponding to Procedure 1 in Sect. 7.
1. An approximate solution for a = 0.38, b = 0.30, c = 4.82 with is calculated (see Fig. 3). It is observed that these coefficients admit an approximate homoclinic orbit. We fix c = 4.82 and take the coefficients a and b as parameters in p = (a, b) . The domain of parameter vectors is defined by D p = � ⟨0.38005, 2.5 × 10 −4 ⟩, ⟨0.30005, 6.5 × 10 −4 ⟩ � with p = (0.38, 0.30) . Here, ⟨c, r⟩ means the center-radius form of an interval. Note that D p is homeomorphic to a closed ball in ℝ 2 . 2. Using eigenvalues and eigenvectors of Df * (p) , we construct a Lyapunov function L around x * = 0 by the method in [16] as follows (see Fig. 4).
Simultaneously, that the domain is verified to be a Lyapunov domain for any parameters p in D p .
Numerical verification method for homoclinic orbits 3. Approximate eigenvalues ̃i and eigenvectors ṽ i of Y are obtained by floating point arithmetic : Moreover it is verified that the following interval vectors contain true eigenvectors of matrix Y by verified computation using INTLAB:   Furthermore we compute   Next we verify existence of homoclinic orbits by the following steps corresponding to Procedure 2 in Sect. 7.
1. The closed ball B + = {x ∈ D L � ‖x‖ ≤ r + } is defined with the radius r + = 0.5 and using the polar coordinates system, we divide the boundary of B + in a similar manner to Step 4 of the above procedure. It is verified that either of the following conditions of Step 1 in Procedure 2 hold.
and k ( k = 1, … , 2200 ). Moreover an interval vector [x � ] k is specified for each k ( k = 1, … , 2200 ), which satisfies (T 3 (x, p b ), x; p b ) ∈ [x � ] k for every x ∈ B 0 and p b ∈ [d] k (see Fig. 8). The verification and the specification are carried out in a similar manner to Step 6 of Procedure 1 for this example. 4. We define the subspace Γ + ⊂ ℝ 3 by Γ + = {x ∈ ℝ 3 | (A + ) T x = 0} and the projection P + ∶ ℝ 3 → Γ + by P + = I − A + (A + ) T , where I is the identity mapping on ℝ 3 and A + is defined by exact eigenvectors. In practical calculation an interval matrix  and the fact V 1 ∩ V 2 ∩ V 3 = � is confirmed by verified computation, which implies deg(P + • (T 3 (x, ⋅), x;⋅)) on D p is not 0 for each x ∈ B 0 . From this result and Corollary 9.5, interval vector D p includes a parameter vector p * ∈ D p which admits a homoclinic orbit from x * (p * ) to itself. Note that the size of the domain of parameter vectors may be taken smaller if we adopt verified computation with multiple precision arithmetic longer than double precision.
With appropriate choices of sets around approximate homoclinic orbits, trajectories through (x ∈)B 0 with p b ∈ D p (yellow curve) arrive at L 0 , the boundary of cone (blue) (color figure online)

A 4-dimensional example
We consider the following autonomous system of ODEs: Note that -the origin is an equilibrium for any coefficients a, b, and c.
-analyzing eigenvalues of Df * , the Jacobian matrix of the right-hand side with respect to (x, y, z, w) at the origin, it is observed that there are 2-dimensional stable manifolds and 2-dimensional unstable manifolds when a > 0 and b > 0 hold. Thus, we take = 2 and n − = 2.
We construct a continuous mapping F 1 by the following steps corresponding to Procedure 1 in Sect. 7.
1. An approximate solution for a = 5, b = 1.116468586, c = 0.060394705655 with is calculated and it is observed that these coefficients admit an approximate homoclinic orbit. We fix c = 0.060394705655 and take the coefficients a and b as parameters in p = (a, b) and the domain of parameter vectors with p = (5.0 , 1.116468586). 2. Using eigenvalues and eigenvectors of Df * (p) , we construct a Lyapunov function L around x * = 0 by the method in [16] as follows.
Simultaneously, that the domain 1.6448 × 10 −9 1.6289 × 10 −9 −1.7060 × 10 −11 6.4714 × 10 −11 x = r − cos 1 , y = r − sin 1 cos 2 , z = r − sin 1 sin 2 cos 3 , w = r − sin 1 sin 2 sin 3 , is confirmed in a similar manner in Step 9 of the above procedure. Thus it is verified that deg(P + • (T 3 (x, ⋅), x;⋅)) on D p is not 0 for each x ∈ B 0 . From this result and Corollary 9.5, interval vector D p includes a parameter vector p * ∈ D p which admits a homoclinic orbit from x * to itself. A numerical verification method for the existence of homoclinic orbits is proposed as an application of local Lyapunov functions with small domains which can be specified by verified numerics. We have shown numerical examples including 4-dimensional case with 2-dimensional stable manifolds and 2-dimensional unstable manifolds. This suggests our method can be applied to such problems that are not easy to treat, and we believe it indicates usefulness of local Lyapunov functions as tools for analysis of dynamical systems. Moreover Sect. 7 shows a kind of instruction manual of our method as a tool, which helps readers to apply our method to their own problems providing the computing environment has implementation of interval arithmetic and verified computation, e.g. INTLAB.
One of our future works is to treat the case where the equilibria admitting homoclinic orbits depend on the parameter vector. The mathematical arguments for the existence of homoclinic orbits become more complicated since Lyapunov functions and their level surfaces vary according to p . Nevertheless, most part of arguments are parallel to those of the present paper, and we will give comprehensive arguments on the system (1.1) without restriction to parameter-independent equilibria in a forthcoming paper.