Predictor-Corrector Algorithm for Curves


This algorithm traces the level set of a an arbitrary smooth function f: R^n -> R^{n-1}. When n=2, it pays close attention to the existence and local resolution of singularities (points where the level curve does not have a well-defined tangent line).

How Predictor-Corrector works

The basic idea is to
  1. find one point on the curve;
  2. compute a tangent vector to the curve at that point;
  3. step out a small amount along the tangent vector;
  4. relocate the curve at a new point
  5. go to step 2
Typically, steps 1 and 4 are accomplished by using Newton's method to iterate a point near the curve onto the curve. Step 2 is accomplished by solving the linear system Df(x) T = 0 where DF(x) is the n x (n-1) Jacobian matrix of f evaluated at a point x on the curve, and where T is the tangent vector to f= 0 at the point x.

This particular implementation attempts to intelligently implement this basic strategy by

  1. Finding special points on the curve to use as initial starting points for the tracing. These "special points" include We allow the user to select which types of initial points to search for and use. Choices include the special points described above as well as points obtained by iterating onto the curve seeds corresponding to points on a uniform grid, to random points, or to points selected by the user.

  2. 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.

  3. 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.

Predictor-Corrector Panel

The panel that controls the predictor-corrector algorithm.

Controlling the Algorithm from Pisces

The Pisces panel to the predictor corrector algorithm has the following parameters which may be set from the algorithm window:
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
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: 1000

The 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.

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
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
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
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
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.

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.

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
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

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
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

Known Bugs

We are not aware of any programming bugs in the new Predictor-Corrector. However, the singularity detection and handling routines are not yet in finished form. They generally work, but sometimes fail to accurately locate or handle singularities. These routines will be improved in the future.

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


This algorithm was implemented as part of the Geometry Center's 1995 Summer Institute for Undergraduates. The Pisces team worked under the direction of Frederick J. Wicklin and with Center apprentice Nicolas Vera on this algorithm and on the related singularity algorithms. The linear algebra routines in used by this algorithm were written by Vera and Wicklin.

The summer students were

For questions about this implementation, contact Frederick J. Wicklin A good starting reference on the predictor-corrector algorithm is
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

Next: Recursive Surface Algorithm
Previous: User Interface
[Pisces] The Pisces Home Page
Comments to:
Last modified: Fri Feb 16 07:59:21 1996
Copyright © 1995 by The Geometry Center, all rights reserved.