- find one point on the curve;
- compute a tangent vector to the curve at that point;
- step out a small amount along the tangent vector;
- relocate the curve at a new point
- go to step 2

This particular implementation attempts to intelligently implement this basic strategy by

- Finding special points on the curve to use as initial starting
points for the tracing. These "special points" include
- turning points (where one of the derivatives to the curve vanishes)
- boundary points (where
*f=0*intersects a bounding box) - singular points (for planar curves, where
*f=0*and*grad(f)=(0,0)*)

- Detecting and locally resolving planar singularities.
By this we mean that we excise a neighborhood of the singularity
and replace that region by a local representation of
*f=0*near the singularity. This uses ideas from singularity theory, differential geometry, and algebraic geometry. - As we trace the curve, we keep track of where we are relative
to any "special points" we may have computed. If we enter a
neighborhood of a special point, then we intelligently connect
*f=0*to the point.

**Step_Length**- This field should contain a positive numerical value. It
indicates the step length for the predictor part of the algorithm.
That is, we take a step of length
`Step_Length`

in the direction tangent to the curve, before using Newton's method to correct back onto the curve. Due to this correction, the actual length of the line segments outputted by the algorithm will vary somewhat at different points along the curve. If`Moore_Penrose_2D`

is off or if the model domain is of dimension greater than 2, each segment will have length AT LEAST`Step_Length`

.

*Default value*: 0.05 **Max_Steps**- This field should contain a positive numerical value. It
indicates the maximum number of steps that the algorithm will attempt in each direction from each initial condition.

*Default value*: 1000The following parameters allow the user to choose among various methods of selecting initial points. At least one of these fields must be nonzero in order to obtain output. Several or all of them may be selected at once.

**Random_Points**- This field should contain a nonnegative integer value. This
value specifies the number of random points (all within the domain
boundaries of the model) that the algorithm will generate. Each of
these points is then used as a seed for Newton's method, which
attempts to iterate these seeds onto the curve. In a model with
domain dimension n, these iterations are confined to an n-1
dimensional subspace of the domain in order to obtain a fully
determined system. If the random point is p, this hyperplane is
chosen as the n-1 dimensional subspace orthogonal to the tangent
vector at p to the level curve f(x) = f(p). Newton's method works
somewhat differently if n = 2 and
`Moore_Penrose_2D`

is set to 1. See documentation for`Moore_Penrose_2D`

below.*Default value:*50 **Selected_Points**- This field takes a boolean value 0 (false) or 1 (true). If
the value is nonzero, the user will be prompted to enter points,
coordinate by coordinate, at the Pisces command line. The user is
allowed to enter as many points as desired. The algorithm attempts to
iterate the points onto the curve using Newton's method.
*Default value:*0 **Uniform_Points**- This field should contain a nonnegative integer value. A
value of 0 indicates that uniform grid points should not be used as
seeds. Given a positive value m for
`Uniform_Points`

in an n-dimensional model, the algorithm constructs a uniform grid with m^n points, and attempts to iterate each point onto the curve using Newton's method. Specifically, the algorithm divides the domain space into m^n uniform n-cubes, and places a grid point at the center of each cube.

*Default value:*4 **Boundary_Points**- This field should contain a nonnegative integer value. If the
value is positive, the algorithm sets up a grid along each boundary
face of the domain space, and attempts to iterate each point onto the
curve under the constraint that all points must remain along that
boundary. Given a value m for
`Boundary_Points`

, the grid used is uniform as described above, with m^(n-1) points per boundary. All initial points found by this method will be along the boundaries.

*Default value:*3 **Turning_Points**- This field should contain a nonnegative integer value. If the
value is positive, the algorithm sets up a uniform grid (where the
Turning_Points parameter specifies the fineness of the grid, as
described in #4 above), and attempts to iterate each point onto a
turning point of the curve with respect to each axis. A turning point
is defined as a point where the tangent to the curve either is not
uniquely defined (a singularity) or is perpendicular to one of the
axes. For each grid point, a separate series of iterations is
required for each axis to which a turning point might correspond.
*Default value:*0**WARNING:**When using models with a large number (n) of dimensions, be careful not to use large values (m) for`Uniform_Points`

,`Boundary_Points`

, and`Turning_Points`

. The values in these fields do not indicate the total number of seeds to use; they represent the fineness of the grid. The functions which find uniform and turning points use m^n seed points, while the one which calculates boundary points uses m^(n-1) seed points per boundary. Large values of`Turning_Points`

(say, above 10 for three or more dimensions) may cause the algorithm to be especially slow, since turning point iterations are relatively complex and since a separate set of these iterations must be executed for each coordinate axis for each point. **Singular_Points**- This field should contain a nonnegative integer value. If the
value is positive, the algorithm sets up a uniform grid (where the
`Singular_Points`

parameter specifies the fineness of the grid, as described above under uniform points), and attempts to iterate each point onto a singularity of the function. In general, we define a singularity as a point where the Jacobian of f is rank deficient. Singularity detection and handling is currently implemented only in 2 dimensions, so setting this parameter to a nonzero value with a higher dimensional model will cause an error message. In two dimensions, a singularity is simply a point where*f=0*and both partial derivatives of f are 0. Singularity searches are based upon Newton's method, but are more complicated than those described above -- both because the system is overdetermined and because the search conditions degenerate as the algorithm approaches a higher-order singularity, causing the search to become unstable.We deal with the latter problem by using an algorithm based on an extension of ideas presented in a paper by W. Govaerts, which dynamically adds search conditions as it approaches a higher-order singularity. When a singularity is found by this routine, a region surrounding it is "desingularized." That is, a local approximation to the curve in this neighborhood is computed and drawn. The singularity point itself is not used as an initial point, but near-singularity points found at the boundary of the desingularization region are used. Singularity detection is generally slow compared to other methods of finding initial points. The above warning applies, even though

*n=2*. **Relative_Sing_Size**- This field takes a real positive value. It is used to specify
the size of the desingularization region: the singular radius
is computed as the product of the
`Relative_Sing_Size`

times`Step_Length`

. For a higher order singularity, the desingularization region is a square centered at the singularity with side length of twice the singular radius. For a node, the singularity-handling routines trace a distance equal to the singular radius in each of the four tangent directions.*Default value:*1 **Moore_Penrose_2D**- This field takes a boolean value 0 (false) or 1 (true). The
value of this field has no effect unless the model is planar. If this
option is turned on and the model is planar, the predictor-corrector
algorithm uses a specialized routine for implementing Newton's
method, both for determining initial points from random, uniform, or
user-selected seeds, and for correcting back onto the curve after
taking a predictor step (the specialized functions are not used while
searching for boundary points, turning points, or singularities).
These functions are referred to as optimized because all formulas are
hard-coded for the special case of two dimensions (in which many of
the formulas become degenerate). LU decomposition and backsolve
routines are not used.
The 2D search methods also differ from the general ones (when used for a planar model) in that they use a different variant of Newton's method. These functions are called in cases when Newton's method is in fact being used to solve a system which is underdetermined by one equation. Classically, Newton's method requires a fully determined system (i.e., it involves the inversion of a square matrix). The general procedure solves this problem by introducing the additional condition that the iterations occur within a hyperplane perpendicular to a previously determined vector (the vector used is typically a tangent vector; when stepping along the curve, it is the vector indicating the direction in which the step was taken). The "optimized" functions do not add this extra condition; they implement Newton's method on an underdetermined system, using the method of the Moore-Penrose inverse described in Allgower and Georg. For the planar case, this simply involves recomputing the gradient at each iteration, and taking a Newton's step in that direction. Our experimentation shows that (after we fixed a seemingly endless string of bugs in these optimized functions) the optimized functions tend to give better results than their general counterparts, and to execute slightly faster. However, the behavior of the optimized functions is sometimes less predictable, and we have found a few cases where the general functions gave better results.

*Default value:*1 **Display_Initial_Points**- This field should contain either 0 (false) or 1 (true). If it
is 0, only the actual curve computed by the algorithm will be
displayed. When this field is set to 1, the initial points that were
used to trace the curve will also be displayed as red dots. In
addition, the algorithm will output its bounding box (separate from
Geomview's bounding box) in black when this option is set.
*Default Value:*0 **Debug**- This field should contain either 0 (false) or 1 (true). When
`Debug`

is 1, a very large quantity of unimportant information will be displayed on the standard output. This output describes whether or not Newton's method successfully converged to each initial point, and how the line segments were traced between the initial points. This option is not intended for general use, although some people find the debugging output strangely interesting.There are also some graphical effects when

`Debug`

is set to 1. After the algorithm computes all of its initial points, it selects an appropriate subset of those initial points to be used in the tracing of the curve. When both Display_Initial_Points and`Debug`

are set to 1, the initial points that were NOT selected to be used to trace the curve will also be displayed. Points that were too close to another point will be displayed as blue dots, and points that were outside of the bounding box will be displayed in green.If you want your X server to spend minutes on end outputting text debugging information, then try setting

`Debug`

to 2 (which will output everything you get with`Debug`

set to 1, and more).*Default Value:*0

You may also encounter one of the known output bugs. None of these bugs are intrinsic to the Predictor-Corrector algorithm; they are caused by the Pisces graphics code or by Geomview. Most of the glaring output bugs have been fixed, but a few remain. In particular, if you use Geomview as your output device, you will occasionally be informed that Geomview can't seek back far enough on its input pipe (and no image will be drawn).

Future improvements to predictor-corrector will include improved singularity detection and improved singularity handling (tracing within the desingularization region).

Please send bug reports to
` software@geom.umn.edu.`

The summer students were

- Erik Streed, streed@geom.umn.edu
- Ron Dror, dror@geom.umn.edu
- John Golden, jgolden@geom.umn.edu
- Brian Johnson, bwj@geom.umn.edu

Numerical Continuation Methods, An Introduction,

by Eugene L. Allgower and Kurt Georg,

published by Springer-Verlag in 1990

for the Springer Series in Computational Mathematics, Number 13

Comments to: pisces@geom.umn.edu

Last modified: Fri Feb 16 07:59:21 1996

Copyright © 1995 by The Geometry Center, all rights reserved.