Osinga, Hinke M

Rokni, Reza

English

In this paper we study the numerical accuracy of computing one-dimensional manifolds of (non-hyperbolic) equilibria in a planar vector field for which the manifolds are known explicitly. We consider the (strong) stable manifolds of a saddle, a sink, and a centre-stable equilibrium

Explore Bristol Research

                          Osinga, H. M., & Rokni, R. (2003). Numerical study of manifoldcomputations.Link to publication record in Explore Bristol ResearchPDF-documentUniversity of Bristol - Explore Bristol ResearchGeneral rightsThis document is made available in accordance with publisher policies. Please cite only the publishedversion using the reference above. Full terms of use are available:http://www.bristol.ac.uk/pure/about/ebr-terms.htmlTake down policyExplore Bristol Research is a digital archive and the intention is that deposited content should not beremoved. However, if you believe that this version of the work breaches copyright law please contactopen-access@bristol.ac.uk and include the following information in your message:• Your contact details• Bibliographic details for the item, including a URL• An outline of the nature of the complaintOn receipt of your message the Open Access Team will immediately investigate your claim, make aninitial judgement of the validity of the claim and, where appropriate, withdraw the item in questionfrom public view.Numerical study of manifold computationsHinke OsingaDepartment of Engineering MathematicsUniversity of BristolBristol BS8 1TR, UKH.M.Osinga@bristol.ac.ukReza RokniSchool of Mathematical SciencesUniversity of ExeterExeter EX4 4QE, UKReza@maths.ex.ac.ukAbstractIn this paper we study the numerical accuracy of computing one-dimensional man-ifolds of (non-hyperbolic) equilibria in a planar vector field for which the manifoldsare known explicitly. We consider the (strong) stable manifolds of a saddle, a sink,and a centre-stable equilibrium.IntroductionThe numerical approximation of stable or unstable manifolds is often used as a tool toinvestigate the global behaviour of a dynamical system. Most methods are designed formanifolds of hyperbolic equilibria and utilize the fact that integration can be done suchthat the manifold is (locally) an attractor. This means that a bounded approximation ofthe manifold can be computed as close to the true manifold as one likes. In this paper,we consider the cases where a one-dimensional (strong) stable manifold of a vector field iscomputed. This necessarily implies that it is the manifold of an equilibrium. We considerboth hyperbolic and non-hyperbolic equilibria.The standard method to compute the one-dimensional stable manifold of an equilib-rium in a vector field is to choose a point on the local manifold close to the equilibriumand integrate it backward in time; for example, see [3]. Most methods use the linearmanifold as an approximation for the local manifold, because the (strong) stable mani-fold is tangent at the equilibrium to the corresponding (strong) stable eigenspace [2]. Ourinvestigations use this method as well.If the equilibrium is hyperbolic and has only a one-dimensional stable manifold, thenthis manifold is an attractor in backward time, at least locally near the equilibrium. Thismeans that initially the distance of points on the computed trajectory to the true manifolddecays. This property is utilized to show that the computational error decays to zero asδ → 0, even though the integration time that is needed to reach a specific finite arclengthevidently goes to infinity in the process; see [4].If one considers a vector field with an equilibrium that is not hyperbolic, or if oneis interested in finding the strong stable manifold of an equilibrium with more than onestable eigenvalue, then standard error bounds no longer apply. In this paper we investigatehow the maximum error varies with δ using a two-dimensional vector field for which themanifolds are known explicitly. We consider the hyperbolic saddle, the hyperbolic sinkand the case where the equilibrium has one zero and one stable eigenvalue.1Explicit manifoldsWe consider an example for which the (strong) stable manifold is known explicitly, namely{x˙ = αx,y˙ = β y + γ x2,(1)with the explicit solution{x(t) = x0 eαt,y(t) = γ2α−βx20 e2αt +(y0 −γ2α−βx20)eβt.(2)Here we assume that α < 0 and β > α, so that the origin has at least a one-dimensionalstable manifold. Note that this automatically ensures that 2α − β 6= 0. If β > 0, thenthe origin is a saddle and we can immediately see from (2) that trajectories with initialcondition (x0, y0) ∈ R2 go to the origin if and only ify0 =γ2α− βx20. (3)Hence, Eq. (3) defines the stable manifold W s(0) of the origin. Similarly, if β = 0, theorigin is non-hyperbolic with a strong stable manifold defined by Eq. (3). In fact, here anypoint on the y-axis is an equilibrium and they all have one stable and one zero eigenvalue.A system of this type is of interest, for example, in adaptive control theory [5].Let us now assume that α < β < 0, so that the origin is stable with a one-dimensionalstrong stable manifold. Since both α and β are negative, all solutions go to the origin.However, since the attraction in the x-direction is stronger than in the y-direction, mostsolutions converge to the origin tangent to the y-axis. The only exception are the solutionson the strong stable manifold W ss(0), which is tangent to the x-axis at the origin. Hence,we are looking for solutions with dydx(t)→ 0 as t→∞. For x0 6= 0, we havedydx(t) =y˙(t)x˙(t)=β y(t) + γ x(t)2αx(t)=βγ2α−βx20 e2αt + β(y0 −γ2α−βx20)eβt + γ x20 e2αtαx0 eαt=2γ2α− βx0 eαt +βα x0(y0 −γ2α− βx20)e(β−α)t.Since α < β < 0 and thus e(β−α)t → ∞ as t → ∞, the strong stable manifold is againdefined by Eq. (3).Throughout this paper, we choose γ = β − 2α such that the (strong) stable manifoldis always given by the graph of the function y = −x2.2Numerical approximationsOur aim is to compute the (strong) stable manifold of the origin for our example (1) withappropriate choices of α < 0, β > α and γ = β−2α. The linear (strong) stable manifold isthe x-axis. Hence, we approximate the branch {(x, γ2α−βx2) | x > 0} up to arclength L bythe horizontal segment from the origin to (δ, 0) followed by the backward-time trajectorythrough (δ, 0), where δ > 0 is small. We use fourth-order Runge-Kutta integration withfixed step size h = 10−5 to approximate such a trajectory. This obviously introduces anintegration error, but it appears to be negligible in our computations. Note that the totalintegration time needed to reach arclength L increases as δ decreases.The numerical approximation consists of a finite set Mδ(L) = {(xi, yi)}Ni=0 of integra-tion points, where (x0, y0) is the origin, (x1, y1) = (δ, 0) and N is the total number ofintegration steps such that (xN , yN) is the first point along the trajectory for which theapproximate arclength distance to the origin is at least L. To keep the total number ofpoints in Mδ(L) small, we consider only a subset of the integration mesh points such thatwe maintain a good resolution of the curve. We ensure that the angle formed by any threeconsecutive points does not exceed a prescribed maximum αmax = 0.3. Furthermore, theproduct of this angle with the distance between the last two of these three points neverexceeds (∆α)max = 10−9. In practice, this means that points will be distributed accordingto curvature, while the bound (∆α)max ensures a small interpolation error. For ease ofnotation, this reduced set of mesh points is still denoted by Mδ(L), and it contains NLpoints, where typically, NL ¿ N .The piecewise linear curve through the mesh points in Mδ(L) forms the approximate(strong) stable manifold. The distance of this approximation to the true manifold isconsidered pointwise. For convenience, we only consider the difference in the y-coordinate.Hence, we define the error of our computation asdist(Mδ(L),Ws(0)) = max0≤i<NL∣∣∣∣yi − γ2α− βx2i∣∣∣∣ .Given how we compute the approximate manifold, this error depends on the initial point(δ, 0) and on the arclength L of the computation.We consider three different cases, namely the case of a hyperbolic saddle, a hyperbolicsink and a non-hyperbolic centre-stable equilibrium. For each case we compute approxi-mations of the (strong) stable manifold up to arclength L = 3 with δ = 0.001. Figure 1shows the behaviour of the pointwise error as a function of the arclength distance to theorigin for all three cases. The stable manifold of the hyperbolic saddle is an attractor inbackward time. Hence, the first point (δ, 0) of the backward trajectory is furthest awayfrom the true manifold (namely 10−6 = |−δ2|), because the trajectory converges to themanifold. In contrast, the strong stable manifold of a hyperbolic sink is a repellor inbackward time. The largest error now occurs at the last integration point. The erroralong the stable manifold of the centre-stable equilibrium does not change at all.Let us now study how the maximum error depends on δ. For each case we calculatemanifolds Mδ(L) with L = 3 fixed throughout, while δ varies. For L = 3, the reducedmesh of Mδ(L) consists of NL ≈ 130 000 points, even though the required integration timemay grow as large as one million.30 0.5 1 1.5 2 2.5 310−1010−910−810−710−610−510−4errorarclengthsaddlecentre-stablesinkFigure 1: The pointwise error versus arclength distance to the origin of approximate(strong) stable manifolds Mδ(L) for the hyperbolic saddle, the hyperbolic sink and thecentre-stable equilibrium; for all cases δ = 0.001 and L = 3.δ 10−1 10−2 10−3 10−4 10−5error saddle 10−2 10−4 10−6 10−8 10−10error centre 10−2 10−4 10−6 10−8 10−10error sink 3.95 10−2 1.25 10−3 3.94 10−5 1.25 10−6 3.94 10−8Table 1: Approximations of the (strong) stable manifold of the origin in (1) up to arclengthL = 3. For the hyperbolic saddle (α, β, γ) = (−1, 1, 3), for the non-hyperbolic centre-stableequilibrium (α, β, γ) = (−1, 0, 2), and for the hyperbolic sink (α, β, γ) = (−2,−1, 3).Hyperbolic saddleThe origin in example (1) is a hyperbolic saddle for α = −1, β = 1 and γ = 3. Note thatthe parameters are such that W s(0) is defined as the graph of Eq. (3).The second row in Table 1 shows the results of our computations when varying δ from10−1 to 10−5. The error decreases quadratically with δ; see also Fig. 2. This behaviour isas expected, because the largest error is always at the initial point.Centre-stable equilibriumThe origin in example (1) has one zero and one stable eigenvalue for α = −1, β = 0and γ = 2. Note that this example is degenerate: the y-axis is invariant and consistsentirely of non-hyperbolic equilibria. Hence, it is not suprising that the computational410−6 10−5 10−4 10−3 10−2 10−1 10010−1010−810−610−410−2100errorδsinksaddle and centre-stableFigure 2: The maximal error dist(Mδ(L),Ws(0)) as a function of δ with L = 3 forthe hyperbolic saddle, the hyperbolic sink and the centre-stable equilibrium (same as thesaddle).error remains the same along the entire manifold. (This also shows that the integrationerror is small enough so that it is negligible.)The maximal error behaves exactly the same as for the case of the hyperbolic saddle;see Table 1. Hence, again, the error decreases quadratically with δ.Hyperbolic sinkExample (1) has a hyperbolic sink at the origin for α = −2, β = −1 and γ = 3.Since the strong stable manifold is a repellor in backward time, we now find that thelargest error is at the last computed point. The maximal error depends on δ as shownin the last row of Table 1. We see that, as before, the error decreases with δ. Moreprecisely, when δ decreases by a factor 10, the error decreases by a factor 3.16 10−2. Thequantitative differences with the hyperbolic saddle are best seen in a log log-plot of theerror versus δ as visualised in Fig. 2.For this example it is also of interest to investigate how the error depends on thearclength L, because the largest error is always at the last integration point. We computedapproximations of the strong stable manifold using the same variations for δ, but nowalso varying L from 100 to 104. The maximal error, i.e. the distance in the y-direction ofthe last computed point to the true manifold, is shown in Table 2.As before, for fixed arclength the error decreases by a factor 3.16 10−2 when δ isdecreased by a factor 10. That is, the contraction rate across the rows in Table 2 isconstant. The contraction rate across the columns in Table 2 is also fairly constant.Namely, when the arclength is reduced by a factor 10, the error decreases by a factor0.56. This reduction rate is slightly smaller (no less than 0.5) when L is relatively small.5δ = 10−1 δ = 10−2 δ = 10−3 δ = 10−4 δ = 10−5L = 100 2.77 10−2 8.74 10−4 2.76 10−5 8.74 10−7 2.76 10−8L = 101 5.52 10−2 1.74 10−3 5.51 10−5 1.74 10−6 5.51 10−8L = 102 9.98 10−2 3.15 10−3 9.97 10−5 3.15 10−6 9.97 10−8L = 103 17.78 10−2 5.62 10−3 17.78 10−5 5.62 10−6 17.78 10−8L = 104 31.62 10−2 10.00 10−3 31.62 10−5 10.00 10−6 31.64 10−8Table 2: Example (1) with (α, β, γ) = (−2,−1, 3): Errors of approximations of the strongstable manifold of the origin for varying arclengths L and initial points (δ, 0).DiscussionThe example that we consider in this paper is constructed so that the manifolds areknown explicitly. This necessarily means that the vector field is relatively simple and thecomputational error of finding the (strong) stable manifolds behaves extremely regular.In general, one should compute two or more approximations while varying δ andconsider the contraction rates of these successive approximations. Such relative errorsand contraction rates give a good idea of the behaviour of the true error and even giveinformation about the size of this error. It is our experience that the contraction rate ofthe error for more general C2 vector fields is still of order 10−2; see [1].Standard error bounds on integration routines predict that the error increases expo-nentially with integration time. Indeed, when taking into account integration time in thedistance calculations, this may well be true, but we are interested in the geometric dis-tance between the manifolds. Hence, in practice, the standard error bounds are nowherenear the actual error. For extremely small δ, i.e. very large integration time, one expectsthat the integration error will start to have an effect and the approximation error willprobably not converge to 0. However, our computations show that approximations witherrors of the order 10−8 are certainly reasonable to expect.References[1] H.M. Osinga, G.R. Rokni-Lamooki, S. Townley, “Numerical approximations of strong(un)stable manifolds,” Applied Nonlinear Mathematics Preprint 2002.19, University of Bris-tol, 2002.[2] J. Palis and W. de Melo, Geometric Theory of Dynamical Systems, Springer-Verlag, 1982.[3] T.S. Parker and L.O. Chua, Practical Numerical Algorithms for Chaotic Systems, Springer-Verlag, 1989.[4] A.M. Stuart and A.R. Humphries, Dynamical Systems and Numerical Analysis, CambridgeUniversity Press, 1996.[5] S. Townley, “An example of a globally stabilizing adaptive controller with a genericallydestabilizing parameter estimate,” IEEE Trans. Automat. Control 44(11): 2238–2241, 1999.6

Numerical study of manifold computations

Abstract

Similar works

Full text

Available Versions

Explore Bristol Research