 This program is the property of:

	      Hinke Osinga and James England
	      Bristol Centre for Applied Nonlinear Mathematics
	      University of Bristol
	      Queen's Building, University Walk
	      Bristol BS8 1TR, United Kingdom

     email: H.M.Osinga@bristol.ac.uk, James.England@bristol.ac.uk

 and may be used, modified and distributed freely, subject to the following
 restrictions:

1. A copy of this file (COPYRIGHT) must be distributed along with any
   copies that you redistribute; this includes copies that you have a
   modified, or copies of programs or other software products that
   include this software.

2. The header must remain in all files. Modifications of the software
   must carry prominent notices in this header, stating who changed
   the files and the date of any change. 

3. Global Manifolds 1D is distributed in the hope that it will be
   useful, but there is no warranty or other guarantee of fitness. The 
   software is provided as is without any obligation on the part of 
   the author to assist in its use, correction, modification or 
   enhancement.

4. We would appreciate an email if you use or have used this software;
   this data will remain confidential and will not be disclosed to any
   third party.

 If you use a picture produced by this software in a publication,
 please give credit with a notice such as the following: 

	Figures 1, 2 and 4 were generated with DsTool [1,2,3].

 where [1,2,3] refers to the following papers:

[1] A. Back, J. Guckenheimer, M.R. Myers, F.J. Wicklin and P.A. Worfolk,
    "DsTool: Computer assisted exploration of dynamical systems",
    Notices Amer. Math. Soc. 39(4) (1992), pp. 303-309.

[2] B. Krauskopf and H.M. Osinga,
    "Investigating torus bifurcations in the forced Van der Pol oscillator",
    in E.J. Doedel and L.S. Tuckerman (Eds.) "Numerical Methods for 
    Bifurcation Problems and Large-Scale Dynamical Systems," 
    IMA Vol. Math. Appl. 119, pp. 199-208, Springer-Verlag 2000.

[3] J. England, B. Krauskopf and H.M. Osinga, 
    "Computing one-dimensional stable manifolds of planar maps without
     the inverse", 
    Bristol Centre for Applied Nonlinear Mathematics 
    preprint #2003.02, April 2003.


======================================================================

Explanation of the parameters
------------------------------------

1) The Global Manifold 1D code gives you the option to compute either
   "both sides," the "positive side," or the "negative side" of the
   manifolds. Which side is positive or negative depends on how DsTool
   computes the corresponding eigenvector of the fixed point. 

2) The algorithm starts with a point on the respective eigendirection. 
   The initial step is specified in the next column

       initial step:			0.001
       arclength stable manifold:		10
       arclength unstable manifold:	10

   By default, the first point is found at distance 0.001 of the fixed
   point. You can think of this point as being the last point of a
   fundamental domain for one side of the manifold. The algorithm adds
   points to a list of points that approximate the manifold, until the
   specified arclength (10 is the default) is exceeded. The length of
   the computed arclength, which is then slightly larger than what you
   specified, is reported in a message.

   You can now compute the stable and unstable manifolds separately by 
   clicking on the appropriate button. Hence, in order to compute only
   one manifold, you no longer need to set the arclength of the other
   manifold to zero.

3) The parameters for the accuracy of the approximation are specified 
   in the next column:

	  alpha_min:		0.20000000000000001
	  alpha_max:		0.29999999999999999
	  (Delta alpha)_min:		9.9999999999999995e-07
	  (Delta alpha)_max:		1.0000000000000001e-05
	  Delta_min:		0.0001
	  epsilon:			0.20000000000000001
	  convergence:		1e-08

   In fact, the numbers were meant to be alpha_min = 0.2, alpha_max =
   0.3, (Delta alpha)_min = 10^-6, (Delta alpha)_max = 10^-5,
   Delta_min = 10^-4, epsilon = 0.2, and convergence = 10^-8. The
   round-offs are shown so that you see how DsTool treats these
   numbers. 

   The accuracy constraints are primarily controlled by alpha_min,
   alpha_max, (Delta alpha)_min, and (Delta alpha)_max. The parameters
   epsilon and convergence are discussed below.

   The parameter alpha denotes an approximation of the angle between
   the lines given by three successive points. If the angle is very
   small, the three points almost lie on a straight line. On the other
   hand, if the angle is large, the three points should be thought of
   as being on a sharp fold, and more points are needed to approximate
   the manifold properly. Recall that the parameter Delta is the
   distance between two successive points. The value (Delta alpha)
   denotes the product of the distance between the last two points and
   the angle between the last three points. The algorithm checks
   whether   

	     alpha_min < alpha < alpha_max,
	     (Delta alpha)_min < Delta alpha < (Delta alpha)_max.

   A point at distance Delta from the last point with alpha <
   alpha_max such that (Delta alpha) < (Delta alpha)_max, is accepted
   as being accurate enough. If also alpha < alpha_min and 
   (Delta alpha) < (Delta alpha)_min, then the next point is searched
   at distance 2 * Delta from the last point, because the computation
   was extremely accurate, i.e. more accurate than specified by the
   user. 

   If the manifold folds sharply, it may happen that the algorithm
   tries to find a point at distance Delta from the last one, with
   Delta approximately 0. It is desirable to force the algorithm to
   accept a point when Delta is smaller than Delta_min even though the
   accuracy constraints are violated, so that the computation time is 
   bounded. If this happens, an error message is given, showing you
   exactly what accuracy constraint is violated. 

   If you want to compute your manifold with higher accuracy, you
   typically decrease the initial step, the parameter Delta_min and
   the choices for (Delta alpha)_min and (Delta alpha)_max. The values
   for alpha_min and alpha_max have less effect. 

4) If the manifold converges to a fixed point, it may have a shorter 
   arclength than what you specified. The parameter convergence is
   used to detect this. The algorithm stops when two successive point
   on the manifold are only convergence (10^-8 is the default) away
   from each other, and prints the length of the computed arclength.
   Note that convergence must always be smaller than Delta_min. 

   The parameter epsilon is only used in the algorithm for computing
   the unstable manifold and for computing the stable manifold as the
   unstable manifold of the inverse map. That is, it is ignored if the
   "Way to calculate the stable manifold" is set to "Search Circle"
   and only has meaning with the "Explicit/Approx Inverse" or "Monte
   Carlo" method; see below for more details on these methods. The
   algorithm produces a sequence of points on the manifold that are
   roughly distance Delta away from each other, where Delta varies
   according to the curvature of the manifold. Since it is
   computationally expensive to find points at exactly distance Delta,
   points are accepted when the distance to the last point is in the
   interval  

	[(1 - epsilon) * Delta, (1 + epsilon) * Delta].

   The choice for Delta is automatically generated each time by the
   algorithm.

5) New to version 2 of the Global Manifold 1D code is the ability
   to select which method is used to compute the stable manifold. A
   drop-down box allows you to choose the preferred method. The options
   are "Explicit/Approximate Inverse", "Monte Carlo" or "Search Circle".
   If either of the first two options are chosen, then the stable
   manifold will be computed as the unstable manifold of the inverse
   map. The option "Explicit/Approximate Inverse" uses the inverse
   specified in the file that defines the model. If only an
   approximate inverse is specified, then this will be used as a seed
   for Newton's method. The option "Explicit/Approximate Inverse" will
   not appear in the drop-down menu if no explicit or approximate
   inverse is defined. The option "Monte Carlo" forces the algorithm
   to use random seeds for Newton's method in order to find the
   inverse, even if an inverse is specified in the model definition.

   The option "Search Circle" only appears when the map is
   two-dimensional. The method has not yet been implemented for higher 
   dimensions. If the option "Search Circle" is selected, then two
   extra parameters need to be taken into account:

   	bisection error:   9.9999999999995e-07
    	iteration max  :   1

   The Search Circle algorithm works be finding new points that map
   back onto a previously computed segment of manifold. The bisection
   error is the maximum distance that an image of the next point found
   can be away from a segment on the previously computed manifold. The
   global accuracy of the computation will not be affected by this
   bisection error if it is of the same order as the interpolation
   error. As a rule of thumb, we choose it of a similar order as 
   (Delta alpha)_max.

   The parameter iteration max is only relevant if the system does not
   have a unique inverse. For noninvertible maps the next point on a
   branch of the manifold may map to the other branch of the manifold,
   which has not been computed yet. This problem can be overcome by
   trying a higher iterate of the map; this is a similar trick as is
   used when the eigenvalue corresponding to the manifold is negative.
   The parameter iteration max specifies the maximum number of 
   iterations that are tried to find a point that maps onto the
   computed manifold. Note that iteration max refers to iterates of
   the map as defined in the model file. That is, if you are using the
   second iterate already, because the eigenvalue is negative, then a
   warning will be printed, but iteration max remains equal to 1. In
   this situation, for iteration max to have any effect, you need to
   change it to at least 3. If the eigenvalue is positive, then
   setting iteration max to 2 will already have an effect. 

6) The stable manifold in noninvertible systems may consist of several
   (sometimes infinitely many) disjoint branches. The Search Circle
   algorithm can compute such a disjoint branch if one point on it is
   known. This point must be saved in the DsTool `Selected' window and
   the "Try to extend manifold from selected point?" must be set to
   "Yes." 

   Typical first points on a disjoint branch are pre-images of the
   fixed point. These can be found by computing the backward orbit of
   the fixed point using Newton's method for finding the inverse and
   Monte Carlo points as seeds. This step must be done by hand.


======================================================================

Useful remarks
------------------

The Global Manifold 1D code should work for any map with
one-dimensional stable or unstable manifolds. The following remarks may 
help to enjoy the code.

1) The manifolds are saved in the same format as if they were
   computed with the original manifold computation of DsTool. This
   means that you can clear them by pressing "Clear manifolds" in
   either the Global Manifold 1D window, or in the Fixed Point window. 

2) As is also done in the original manifold computation of DsTool,
   Global Manifolds 1D computes each branch of the manifolds entirely
   before it is saved. Saving points may take a long time when a
   branch consists of a lot of points. 

3) DsTool does not draw lines between successive points on the
   manifold (not even for the manifold of a vector field). However,
   you can port the manifolds to Geomview (for two-dimensional systems
   just take "0" for the z-coordinate) and choose to connect the dots
   before sending them. 

4) It is possible that the Global Manifold 1D code tries to find a
   pre-image in an interval whose image gets stretched so much that
   even points that differ only in the last digit get mapped too far
   apart. When this happens, the computation stops and an error
   message is given: 

      Preimages are equal, yet no point found
      ERROR: manifold stretches too much!

   This typically happens for unbounded manifolds at a long arclength
   distance away from the fixed point. 

5) The Global Manifold 1D code is not well designed for a PERIODIC
   phase space. If the code does not behave properly, it is best to
   change the manifold type in the model definition to EUCLIDEAN and
   define auxiliary functions to divide out the periodicity. The
   manifolds will then be computed in the covering space, but can be
   projected to the proper phase space using these auxiliary functions
   as the variables.



