Intersections of Surfaces

Revised  3/3/1999

This article presents an algorithm, along with a Mathematica package implementing it, that gives a parametrization for the intersection of two surfaces.

by Steven Wilkinson

If two surfaces in space intersect, their intersection is usually a curve or union of curves. How do you get your hands on a mathematical representation of those curves if you want to plot them or find out other information about them, such as their lengths? Sometimes, algebraic manipulation of the two surface equations provides an answer. For example, suppose we wish to find the intersection of the sphere centered at the origin of radius 1, given by (x^2) + (y^2) + (z^2) = 1,  and the saddle surface, given by z = (x^2) - (y^2).

Needs["Graphics`ContourPlot3D`"]

\(Block[{$DisplayFunction = Identity}, 
      sphere\  = 
        ContourPlot3D[
          x\^2 + y\^2 + z\^2, {x, \(-1\), 1}, {y, \(-1\), 1}, {z, \(-1\), 1}, 
          Contours -> {1}]; 
	
      saddle = Plot3D[x\^2 - y\^2, {x, \(-1\), 1}, {y, \(-1\), 1}]; 
	]\)

intersectingSurfaces  = Show[sphere, saddle];

To determine a representation for the curve of intersection we can simultaneously solve the equations that represent these two surfaces.

solution = 
      Solve[{(x^2) + (y^2) + (z^2) == 1, z == (x^2) - (y^2)}, {x, z}]

This produces 4 parametrizations, with parameter y,  for parts of the curve of intersection. From the solution it is evident that the values of y are restricted to be at least between ±(Sqrt[5/8]). In fact, they are even more restricted for the first two solutions. Setting the longer expression under the radical to zero and solving for y shows that the first two solutions are real-valued only for \(TraditionalForm
      \`\(-\@\(5/8\)\) ≤ y ≤ 
        \(-\@\(1\/2\ \((\(-1\) + \@5)\)\)\)\) or \(TraditionalForm
      \`\@\(1\/2\ \((\(-1\) + \@5)\)\) ≤ y ≤ 
        \@\(5/8\)\).

{r1[y_], r2[y_], r3[y_], r4[y_]} := Evaluate[{x, y, z} /. solution];

\(\(\(Block[{$DisplayFunction = Identity}, 
	
      intersection = 

        	{ParametricPlot3D[
			
            Evaluate[{r1[y], r2[y]}], {y, \(-\@\(1/2 \((\(-1\) + \@5)\)\)\), 
              \(-\@\(5/8\)\)}], 
		
          ParametricPlot3D[
			
            Evaluate[{r1[y], r2[y]}], {y, \@\(5/8\), 
              \@\(1/2 \((\(-1\) + \@5)\)\)}], 
		
          ParametricPlot3D[
            Evaluate[{r3[y], r4[y]}], {y, \(-\@\(5/8\)\), \@\(5/8\)}, 
            PlotPoints -> 150]}]; \)
	\)\)

Show[intersection];

Unfortunately, this is not a very general method because it relies heavily on our ability to algebraically  solve for two of the three coordinate variables in terms of the remaining one. This will almost never work if one or both of the surfaces are non-algebraic. For example, try it on a sphere and the trigonometric surface z = sin(x y). Even if both surfaces are algebraic as in the first example, we have already seen that certain undesirable effects can occur. In particular, that first example had four parametrizations  with different domains for one connected intersection component.
	Let me present a technique for deriving parametrizations for the curves of intersection that overcome some of these limitations. Furthermore, the algorithm gives you control over the type of parametrization you wish to obtain. The only restrictions on this technique are that it will only give the curves of intersection one connected component at a time, the surfaces must be sufficiently smooth or differentiable, and you need to know at least one point on the curve to get the algorithm started.
	The algorithm involves solving a system of ordinary differential equations with an initial condition. The method is a variation of an algorithm first presented by [Bajaj, Hoffmann,  and Lynch 1988] in which they used Taylor polynomials to approximate the curve. By setting the problem up as an initial-value problem, all the details on the accuracy and precision of the approximations can be left to the differential equations solver, here the  NDSolve command.
	I have implemented this algorithm in a Mathematica package called ParametrizeCurve`. The syntax of the one command is very similar to that of the NDSolve command since it is based on it. As an example of its use, let’s find a parametrization for the intersection of the sphere centered at the origin of radius 2 and the trigonometric surface z = sin(x y). The command requires a rough approximation to one point on the intersection, and a quick check shows that the point ((Sqrt[2]), 0, 0) lies on the curve. After finding a parametrization of the intersection we will plot it and determine its length.

Needs["ParametrizeCurve`"]

\(solution\  = \ 

      ParametrizeCurve[{x\^2 + y\^2 + z\^2 == 2, z == Sin[x\ y], x[0] == \@2, 
          y[0] == 0, z[0] == 0}, {x, y, z}, {t, 0, 11}]\)

r[t_] := Evaluate[{x[t], y[t], z[t]} /. solution];

\(Block[{$DisplayFunction\  = \ Identity}, \ 
	
      intersection\  = 
	
        ParametricPlot3D[
          Evaluate[Append[r[t], Thickness[0.015]]], {t, 0, 11}]; 
	
      sphere\  = 
        ContourPlot3D[
          x\^2 + y\^2 + z\^2, {x, \(-\@2\), \@2}, {y, \(-\@2\), \@2}, {z, 
            \(-\@2\), \@2}, Contours -> {2}, 
          ContourStyle -> EdgeForm[]]; 
	
      trigSurface = 	
        Plot3D[Sin[x\ y], {x, \(-1.5\), 1.5}, {y, \(-1.5\), 1.5}, 
          Mesh -> False]; 
	
      all = Show[intersection, sphere, trigSurface]; 
	]\)

Show[GraphicsArray[{intersection, all}]];

By default ParametrizeCurve gives parametrizations that are parametrized by arclength. Therefore, to determine the length of this intersection we only need to determine the period of the parametrization. Since the construction of this intersection was done numerically, we will determine its period by calculating the minimum value of \(TraditionalForm\`\(\(|r(t) - r(0)\)\(|\^2\)\)\) for t around 11.

FindMinimum[(r[t] - r[0]).(r[t] - r[0]), {t, 11}]

So, the approximate length of the intersection is 10.5111.

A tangential note on plotting curves on surfaces.

You may have experienced the problem of plotting curves that lie on surfaces that the above picture points out. Only parts of the curve may appear. To better illustrate this problem, redraw the above picture without the trigonometric surface and for -1 ≤ x ≤ 0.

\(\(hemisphere = 
	
      Show[intersection, sphere, 
        PlotRange -> {{\(-1.5\), 0}, Automatic, Automatic}]; \)\)

You can see that the curve appears on the outside of the sphere but not the inside. This is due to the fact that Mathematica plots surfaces as the union of polygonal pieces of planes and curves as the union of line segments. Whenever a curve is plotted with more points than a surface, as is usually the case, if the curve has nonzero curvature it will appear to lie on only one side of each polygon that makes up the surface. One way to ameliorate this situation is to slightly move the curve towards the viewer by translating the curve along the vector that points from the middle of the viewing box to the viewpoint.  This can be done by using the options BoxRatios, PlotRange and ViewPoint for the graphics. Here is how to obtain the values of those three options for hemisphere.

{boxratios, plotrange, viewpoint} = 
      Map[FullOptions[hemisphere, #]&, {BoxRatios, PlotRange, ViewPoint}];

ViewPoint is actually given in a scaled coordinate system, and to use it we will need to convert it to the coordinate system of PlotRange.

\(\(viewPtConverted\  = 
      Apply[Plus, \ plotrange, \ {1}]/2\  + \ 
\ \ \ 
        \((\(-viewpoint\)*Apply[Subtract, \ plotrange, \ {1}])\) 
          Max[boxratios]/boxratios; \)\)

Define a parametrization of the lifted curve by bringing every point toward the view point by a certain percentage, say 3%.

rLifted[t_]=
	r[t]+0.03(viewPtConverted - r[t])/√((viewPtConverted - r[t]).(viewPtConverted - r[t]));

Block[{$DisplayFunction = Identity}, 
      liftedIntersection = 
        ParametricPlot3D[
          Evaluate[Append[rLifted[t], Thickness[.015]]], {t, 0, 11}]]

\(\(hemisphere = 
	
      Show[liftedIntersection, sphere, 
        PlotRange -> {{\(-1.5\), 0}, Automatic, Automatic}]; \)\)

Here is a final picture of the sphere and the trigonometric surface with the intersection clearly delineated.

Show[liftedIntersection, sphere, trigSurface];

The Algorithm

If the two surfaces whose intersection we wish to parametrize are given by (Subscript[F, 1])(x, y, z) = 0 and (Subscript[F, 2])(x, y, z) = 0,  and we already have a parametrization, X(t) = (x(t), y(t), z(t)), of the intersection, then (Subscript[F, i])(X(t)) = 0 for i = 1, 2. Differentiate these last two equations to get

∇(Subscript[F, i]) (X (t))·(X^′) (t) = 0.

This says that X(t) is perpendicular to both the gradients ∇(Subscript[F, 1]) and ∇(Subscript[F, 2]). As long as these two gradients are linearly independent for points on the intersection, then (X^′)(t) is some multiple of the cross product \(TraditionalForm
      \`\(∇\(\(F\_1\)(X(t))\)×∇\(\(F\_2\)(X(t))\), \)\)

(X^′)(t)=λ(t)∇(Subscript[F, 1])(X(t))
                  ×∇(Subscript[F, 2])(X(t)).

For this problem we are given the two surfaces, (Subscript[F, 1])(x, y, z) = 0 and (Subscript[F, 2])(x, y, z) = 0,  and we wish to find a parametrization for the intersection, X(t). Equation (2) is a differential equation in the unknown X(t). The appearance of the scalar function λ(t) is not incidental to the discussion. For if X(t) is a parametrization for this intersection, then so is any reparametrization of it. Moreover, the derivative of the reparametrization will be a scalar multiple of the derivative of the original parametrization. That is the source of λ(t). Choose λ(t) however you want to obtain the type of parametrization you desire. For example, if you want X(t) to be parametrized by arclength, |(X^′)(t) |  = 1, then choose \(TraditionalForm
      \`λ(t) = 
        \(\((∇\(\(F\_1\)(X(t))\)×∇\(\(F\_2\)(X(t))\))\)\^\(-1
              \).\)\) One more thing is needed to solve equation (2), namely an initial condition. This corresponds to a choice of a point on the curve and a choice of a parameter value to go with it, X((Subscript[t, 0])) = ((Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])).
	As an example consider the intersection of the sphere centered at the origin of radius 2 and the trigonometric surface z = sin(x y). Here, (Subscript[F, 1])(x, y, z) = (x^2) + (y^2) + (z^2) - 1 and (Subscript[F, 2])(x, y, z) = sin(x y) - z,  so that

\(TraditionalForm
    \`∇F\_1×∇F\_2 = 
      \((\(-2\) y - 2  x\ z\ \(cos(x\ y)\), 2  x + 2  y\ z\ \(cos(x\ y)\), 
        2  x\^2\ \(cos(x\ y)\) - 2 \( y\^2\) \(cos(x\ y)\))\)\)

The point ((Sqrt[2]), 0, 0) lies on the curve of intersection, so use that as the initial point when t = 0. For the sake of simplicity take λ(t) = 1. This gives the initial condition problem:

\(TraditionalForm
    \`\((\(x\^′\)(t), \(y\^′\)(t), \(z\^′\)(t))\) = 
      \((\(-2\) y - 2  x\ z\ \(cos(x\ y)\), 2  x + 2  y\ z\ \(cos(x\ y)\), 
        2  x\^2\ \(cos(x\ y)\) - 2 \( y\^2\) \(cos(x\ y)\))\), 

    \((x(0), y(0), z(0))\) = \((\@2, 0, 0)\)\)

To solve it for 0 ≤ t ≤ 11 use NDSolve:

\(NDSolve[

      	{\(x'\)[t] == \(-2\) y[t] - 2  x[t]\ z[t]\ Cos[x[t] y[t]], 
		
        \(y'\)[t] == 2  x[t] + 2  y[t] z[t] Cos[x[t] y[t]], 
		
        \(z'\)[t] == 
          2 \( x[t]\^2\) Cos[x[t] y[t]] - 2 \( y[t]\^2\) Cos[x[t] y[t]], 
	
        x[0] == \@2, y[0] == 0, z[0] == 0}, 
	{x, y, z}, {t, 0, 11}]\)

This is essentially what ParametrizeCurve does if the two gradients ∇(Subscript[F, 1]), ∇(Subscript[F, 2]) are linearly independent along the curve.

Dealing with singular points

The algorithm has a problem if it runs into a point on the intersection where ∇(Subscript[F, 1]) and ∇(Subscript[F, 2]) are linearly dependent since equation (1) no longer implies equation (2). Such points are called singular points of the intersection. One’s first attempt to deal with such points may be to simply ignore them. Since singular points are rare (linearly dependence of two vectors in space is a degenerate case), and we are numerically solving an intial value problem, it is highly unlikely that we will ever integrate right through a singular point. However, as the next example will illustrate, difficulties can still arise. Parametrize the intersection of the sphere (x^2) + (y^2) + (z^2) = 1 and the cylinder (((2  x - 1))^2) + (((2  y))^2) = 1. With (Subscript[F, 1])(x, y, z) = (x^2) + (y^2) + (z^2) - 1 and \(TraditionalForm
      \`\(F\_2\)(x, y, z) = \((2  x - 1)\)\^2 + \((2  y)\)\^2 - 1\) a quick check shows that the point (1, 0, 0) is on the intersection and at that point ∇(Subscript[F, 1]) , ∇(Subscript[F, 2]) are both multiples of the vector (1, 0, 0).

Block[{$DisplayFunction = Identity},sphere  = 
            ContourPlot3D[
              (x^2) + (y^2) + (z^2), {x, -1, 1}, {y, -1, 1}, {z, -1, 
                1}, Contours -> {1}];
	cylinder = 
            ContourPlot3D[
              (((2  x - 1))^2) + (((2  y))^2), {x, 0, 1}, {y, -1, 1}, {
                z, -1, 1}, Contours -> {1}];
	\(solution\  = \ 
		
              ParametrizeCurve[{x\^2 + y\^2 + z\^2 == 1, 
                  \((2 x - 1)\)\^2 + \((2 y)\)\^2 == 1, x[0] == 0, y[0] == 0, 
                  z[0] == 1}, {x, y, z}, {t, 0, 8}, 
                CheckSingularPoints -> False]\);
	r[t_]:=Evaluate[{x[t], y[t], z[t]} /. solution];
	surfaces = 
              Show[sphere, cylinder, ViewPoint -> {2.6,  -1.7,  1.3}];
	intersection = 
              ParametricPlot3D[Evaluate[r[t]], {t, 0, 3}, 
                PlotRange -> {{-1, 1}, {-1, 1}, {-1, 1}}, 
                ViewPoint -> {2.6,  -1.7,  1.3}];
]

Show[GraphicsArray[{surfaces, intersection}]];

Notice that the intersection should be a “figure-eight” curve on the sphere. However, using 0 ≤ t ≤ 3 we obtain a non-smooth curve. The curve takes a sharp turn at the singular point instead of passing through to the lower part of the “figure-eight.” For a larger range of t-values, we do obtain more of the intersection, but its behavior around the singular point is not what is desired. This, in fact, is generally what happens around most singular points when you integrate system (2). For a discussion why, see [Wilkinson, 1998] where the plane case is fully investigated. The plane case can be viewed as the basis for the three-dimensional case as is shown in that same reference.
	To proceed it is useful to make an affine change of the generators (Subscript[F, 1]),  (Subscript[F, 2]) so that one of the gradients, say ∇(Subscript[F, 2]), is the zero vector at the singular point. (In the above example with the cylinder and sphere we could make the change \(TraditionalForm\`\(F\&~\_2 = F\_2 - 2  F\_1, \)\) and \(TraditionalForm\`F\&~\_1 = F\_1\)) Now differentiate equation (1) with i = 2 to get

\(TraditionalForm
    \`\((\(X\^′\)(t))\)\^T
            ·∇\^2\(\( F\_2\)(X(t))
              \)·\(\(X\^′\)(t)\) + 
        ∇\(\(F\_2\)(X(t))\)·\(\(X\^′′\)(t)\) = 
      0.\)

Here, (∇^2) (Subscript[F, 2]) is the second derivative Hessian matrix for (Subscript[F, 2]).

(∇^2) (Subscript[F, 2])=(({{(∂_(x, x))(Subscript[F, 2]), (∂_(x, y))(Subscript[F, 2]), (∂_(x, z))(Subscript[F, 2])}, {(∂_(y, x))(Subscript[F, 2]), (∂_(y, y))(Subscript[F, 2]), (∂_(y, z))(Subscript[F, 2])}, {(∂_(z, x))(Subscript[F, 2]), (∂_(z, y))(Subscript[F, 2]), (∂_(z, z))(Subscript[F, 2])}}))

Since at a singular point ∇(Subscript[F, 2]) (X(t)) = 0, equation (3) is \(TraditionalForm
      \`\((\(X\^′\)(t))\)\^T
            ·∇\^2\(\( F\_2\)(X(t))
              \)·\(\(X\^′\)(t)\) = 0.\) This along with equation (1) for i = 1 implies that (X^′)(t) is mutually perpendicular to ∇(Subscript[F, 1])(X(t)) and \(TraditionalForm
      \`\((\(X\^′\)(t))\)\^T·∇\^2 F\_2 \((X(t))\)\) at the singular point. Impose the first of two nondegeneracy assumptions on the singular points.

Nondegeneracy Assumption 1. If ((Subscript[x, 0]), (Subscript[y, 0]), (Subscript[z, 0])) is a singular point of the intersection, then the linear span of \(TraditionalForm
      \`{∇\(\(F\_1\)(x\_0, y\_0, z\_0)\), 
        ∇\(\(F\_2\)(x\_0, y\_0, z\_0)\)}\) is one-dimensional. (If the point is singular then this span is at most one-dimensional.)

With this assumption and by making ∇(Subscript[F, 2]) (X(t)) = 0 at the singular point, ∇(Subscript[F, 1]) (X(t)) ≠ 0 at the singular point. To determine X(t), solve

\(TraditionalForm\`∇F\_1·u = 0, 

    u\^T·∇\^2\( F\_2·u\) = 0\)

for u in the following way. From the first of those two equations u lies in the orthogonal space to ∇(Subscript[F, 1]); name this space ∇((Subscript[F, 1])^ ⊥) . Restrict the second of these two equations to the two-dimensional vector space ∇((Subscript[F, 1])^ ⊥) . Since this subspace is two-dimensional, there exists a transformation that rotates every vector 90° in the counterclockwise direction with respect to the normal ∇(Subscript[F, 1]). Call this transformation J. Now the equation \(TraditionalForm
      \`u\^T·∇\^2\( F\_2·u\) = 0\) implies that u is orthogonal to (∇^2) (Subscript[F, 2])·u since their dot product is zero. Using the transformation J this means that

J (∇^2) (Subscript[F, 2])·u = ρ u

for some scalar ρ. But then u is an eigenvector of \(TraditionalForm
      \`J(∇\^2 F\_2\( | 
            \_\(∇\(F\_1\^ ⊥ \)\)\))\). Since this is a transformation on a two-dimensional space, there are three cases to consider:

i. There are two real and distinct eigenvalues and thus two linearly independent one-dimensional real eigenspaces. This is the case with the intersection of the above sphere and cylinder. Can you see the two eigenspaces at the singular point? They are parallel to the two tangent lines at the singular point.

ii. There are two complex eigenvalues and thus no real eigenspaces. This would be the case if the intersection of two surfaces was an isolated point. For example, the intersection of the sphere (((x - 1))^2) + (y^2) + (z^2) = 1 and the cylinder (x^2) + (y^2) = 4 is only the point (2, 0, 0).

iii. There is only one real eigenvalue. This degenerate case is more complicated and I will not handle it here. One can check that if A is a 2×2 matrix and J is a rotation matrix, then J A has at most one real eigenvalue if and only if A is a singular matrix.

Nondegeneracy Assumption 2. If ∇(Subscript[F, 2]) (X(t)) = 0 then \(TraditionalForm
      \`\(∇\^2 F\_2\( | 
            \_\(∇\(F\_1\^ ⊥ \)\)\)\)\) is nonsingular.

Implementing the algorithm around singular points

In the implementation of the package  ParametrizeCurve` the singular points are not found in advance. To do so would be difficult even in the case where the surfaces were algebraic. Instead the package monitors the value of \(TraditionalForm
      \`\(\(|∇\(\(F\_1\)(X(t))\)×∇\(\(F\_2\)(X(t))\)\) | 
        \)\) as it integrates the initial condition problem that arises from (2). If its value becomes sufficiently close to 0, then it assumes the curve has a singular point nearby. If a singular point is near, and assuming the two nondegeracy assumptions, there are either nonreal eigenvectors or two distinct eigenvectors for equation (4). If they are not real and the package is exactly at the singular point it knows it will never move from it since it is an isolated point. If they are not real and it is close to it but not at the singular point, it can ignore the singular point since it knows it will never get to the isolated singular point. However, if the eigenvectors are real then there are four directions it can head, two opposing directions for each of the two distinct eigendirections. It chooses the direction that most closely points in the direction it has been integrating. This gives a smooth parametrization through the singular point.
	Finding these eigenvectors is accomplished by the use of the Singular Value Decomposition Theorem which states that any matrix A can be decomposed in the following way:

A = (U^T) D V

where U and V are orthogonal matrices and D is a diagonal matrix with entries made of the nonnegative eigenvalues of A (A^T) arranged in decreasing order. Mathematica has the built-in function SingularValues that computes this decomposition. Decompose the matrix whose first row is ∇(Subscript[F, 1]) (X(t)) and second row is ∇(Subscript[F, 2]) (X(t)), call that matrix A. At a singular point A is a singular matrix since its rows are linearly dependent. There are two things this decomposition accomplishes for this problems. First, we needed an affine change of generators so that ∇(Subscript[F, 2]) = 0. Since matrix U is orthogonal, (U^T) = (U^(-1)). Therefore, U A = D V. At a singular point, the last row of D is a zero row. This makes the last row of the product D V a zero row, and thus also the last row of U A. But the rows of A are ∇(Subscript[F, 1]) and ∇(Subscript[F, 2.]) Take the rows of the product U A as \(TraditionalForm\`∇F\&~\_1\) and \(TraditionalForm\`∇F\&~\_2\). (I will drop the tildes from now on. Assume ∇(Subscript[F, 2]) = 0.) The other thing to do is solve equation (4) on the two dimensional space ∇((Subscript[F, 1])^ ⊥)  where ∇(Subscript[F, 1]) is now the first row of U A.The equation U A = D V shows that the columns of orthogonal matrix V form an orthonormal basis for (ℝ^3) with the first column pointing in the same direction as ∇(Subscript[F, 1.]) Hence, the last two columns of matrix V, call those columns (Subscript[v, 2]),  (Subscript[v, 3]), form an orthonormal basis for ∇((Subscript[F, 1])^ ⊥) . The desired rotation matrix J is

J=(({{(Subscript[v, 2]), (Subscript[v, 3])}}))(({{0, -1}, {1, 0}}))(({{((Subscript[v, 2])^T)}, {((Subscript[v, 3])^T)}}))

and the projected vector \(TraditionalForm
      \`\(∇\^2 F\_2\( | 
            \_\(∇\(F\_1\^ ⊥ \)\)\)\)\) is found with the standard projection transformation

\(TraditionalForm
    \`\(proj\_\(∇\(F\_1\^ ⊥ \)\)\)(x) = 
      x\  - \ \(x
            ·∇F\_1\)\/\(∇F\_1·∇F\_1.
                \)\)

Below is the last example of the intersection of the cylinder with the sphere where this time the package checks for singular points using the above algorithm. Note the option CheckSingularPoints that is now set to True.

\(\(solution\  = \ 
		
      ParametrizeCurve[{x\^2 + y\^2 + z\^2 == 1, 
          \((2  x - 1)\)\^2 + \((2  y)\)\^2 == 1, x[0] == 0, y[0] == 0, 
          z[0] == 1}, {x, y, z}, {t, 0, 8}, 
        CheckSingularPoints -> True]; 
	
    r[t_] := Evaluate[{x[t], y[t], z[t]} /. solution]; 
	
    ParametricPlot3D[Evaluate[r[t]], {t, 0, 8}, 
      PlotRange -> {{\(-1\), 1}, {\(-1\), 1}, {\(-1\), 1}}, 
      ViewPoint -> {2.6, \ \(-1.7\), \ 1.3}]; 
\)\)

Options for the package ParametrizeCurve`

I included several options to the command  ParametrizeCurve to increase the flexibility of its use. Also, any option that NDSolve takes can be given to ParametrizeCurve and it will pass it on to NDSolve. What follows are the options to ParametrizeCurve listed in alphabetical order with their default values:

AdjustInitialPoint takes boolean values. If set to True any time the given initial point is not on the curve, the nearest point on the curve to the given point is taken as the initial point. This initial point is found by a variation of the steepest descent method. The default value is True.

CheckSingularPoints takes boolean values.If set to False the package doesn't bother to check for singular points. The package runs much faster in this mode than when it does check for singular points since NDSolve with its default settings spends a great deal of time on its adaptive scheme for choosing the stepsize whenever it nears a singular point where the package alters the system of differential equations. The default value is True.

InitialDirection expresses the general direction, as a vector, in which you wish the parametrization to initially head. There is one of two directions a parametrization can head, backwards or forwards. The package projects the given vector onto the tangent line to determine the desired direction. The default value is Automatic.

MinNorm expresses the value such that if the norm of the cross product of the gradients of the generators is less than that value, the package will automatically assume it needs to switch to the singular point algorithm. The default value is 0.

MaxNorm expresses the value such that if the norm of the cross product of the gradients of the generators is greater than that value, the package will automatically assume it need not switch to the singular point algorithm. The default value is Infinity.

ParametrizationSpeed expresses the desired speed for the parametrization. If set to UnitSpeed then the resulting parametrization is a unit speed parametrization. If set to a number, the resulting parametrization has constant speed that number. It can also be set to a function, either a pure function or an expression in terms of the parameter. If set to GradientSpeed the parametrization's speed is the same as the norm of the cross product of the gradients of the generators. The default value is UnitSpeed.

VectorNorm expresses the type of vector norm used to measure lengths. It can be Norm[p] which is the standard vector p-norm with p any whole number or ∞, or it can be a pure function. The default value is Norm[2].

In summary

The package ParametrizeCurve` implements an algorithm that allows one to numerically find parametrizations for the intersection of two surfaces. Moreover, it can parametrize through nondegenerate singular points on the intersection. Furthermore, the algorithm extends with little modifications to parametrize the curve that is the intersection of n - 1 hypersurfaces in (ℝ^n), and the package will do so. In particular, it will parametrize a curve in the plane that is given implicitly, F(x, y) = 0,  as a level curve.

References

Bajaj, C. L., Hoffmann, C. M., Hopcroft, J.E.H., Lynch, R. E. 1988.
	Tracing surface intersections. Computer Aided Geometric Design 5 (4) 285–307.
	
Wilkinson, Steven. 1998.
	Numerically parametrizing curves. College Mathematics Journal  29 (2) 104–119

About the Author

Steven Wilkinson is an associate professor in the Department of Math and CS at Northern Kentucky University. He received his Ph.D. in mathematics from Rice University in 1984.

Steven Wilkinson
Department of Mathematics and Computer Science
Northern Kentucky University
Highland Heights, Kentucky 41099-1700
wilkinson@nku.edu