Recently, while flipping through an old calculus text, I spotted an exercise that called for finding the steepest-ascent curve on a surface described by a function . Ah, thought I, that should be easy enough to do in Maple. It was. Then I thought, what if the gradient field can't be integrated analytically. Can a numeric integration generate the steepest-ascent curve? Of course. And how about if the surface itself is given digitally, not analytically by a formula. Can the steepest-ascent curve still be determined? Again, yes. The details are all given below.

Figure 1 is a graph of the surface be defined by

Steepest-ascent curves from any point on the ellipse  are integral curves of the gradient field . , itself everywhere orthogonal to the level curves for the function. With , we can use the interactive task template in Table 1 to obtain a graph of the flow lines of the gradient field. After entering the appropriate data, simply click on a point through which an integral curve is to be drawn. Constrained (one-to-one) scaling is imposed with the Context Menu.

Figure 2 shows the level curves and the orthogonal trajectories, generated by the commands hidden in the Table holding the graph.

The integral curves of the gradient field can be found by solving the system

in which case the steepest-ascent curves are given parametrically by , where  is the parameter along the curve. Alternatively, the gradient field can be integrated by solving the differential equation

in which case the steepest-ascent curves are given by .

For the surface depicted in Figure 1, the latter differential equation is , with solution . The apex of that surface is a critical point for the gradient field, that is,  at the top of the hill. Hence, it is necessary to find curves of steepest ascent, and for this, we select initial points along the level curve , that is, along the ellipse . Figure 3 shows the surface, along with a collection of steepest-ascent curves emanating from points along the base ellipse. The details are hidden behind the table containing the figure.

The function used for Figures 1 and 3 is sufficiently simple that the gradient field can be integrated analytically. Should this not be the case, numeric integration could be used. Code for this is illustrated in Table 2. The lists X and Y defining initial points on the base ellipse are defined behind Figure 3.

The first two differential equations in the list sys integrate the planar gradient field, while the third is  

The use of the parameters option in dsolve/numeric simplifies changing the initial conditions for each ascent curve. Before the resulting numerically determined curve can be graphed, the intervening line that defines a set of parameters must be executed. Surprisingly, this invocation must be in terms of floats - hence, the evalf command.

In this final phase of our investigation, we assume that the surface is given, not analytically by a formula, but digitally, by a collection of data coordinates. Because we will interpolate the data points to generate a smooth surface, we need to impose some structure on the grid supporting the heights of the surface. The grid does not have to be equispaced in the supporting xy-plane, but it has to be regular. By regular, we mean that heights are known at each point , where  is a member of an ordered list of x-coordinates in the data set, and  is a member of the corresponding list of y-coordinates. In effect, the surface is known at each point of the direct product of the lists of x- and y-coordinates.

Thinking of the function used earlier, namely

we write the coordinate lists

and generate a matrix of corresponding function values:

Next, we define , a function that interpolates the given data.

The first argument to ArrayInterpolation is a list of the lists of the independent coordinates. The second argument is the matrix of known heights. The third argument defines the "target," that is, the points at which "new" function values are to be calculated. Ordinarily, this would be an Array with the coordinates of many points. Here, it consists of the coordinates of a single point, . Thus,  is the interpolated height at the single point .

The final argument to ArrayInterpolation is the interpolating method. Since this is merely an illustrative example, we've picked "spline" for no particular reason. Figure 5 displays the surface generated by the interpolating function .

Now we must devise a "steepest-ascent algorithm" that works strictly numerically, on just data. We don't have a function whose gradient field we can integrate. We can use only the given data points and the interpolating function .

Here is our thinking. If we are at a point  where the height is , we want to move to the "nearest highest point." Of course, for a continuous surface, this is meaningless. So instead, we ask for the nearest highest point in a neighborhood of P, and define that neighborhood as a small circle around P. But how do we find the highest point on  subject to a constraint such as ? We use the Search command from Dr. Sergey Moiseev's DirectSearch package. A sequence of such optimizations generates the steepest-ascent curve as a set of points in .

For example, we start a steepest-ascent curve at the point on the base contour where , and compute the corresponding y-coordinate to be

and hence, the initial point on the steepest-ascent curve is

Additional points are then obtained as per the calculations in Table 3.

Figure 6 contains a graph of the connected points, superimposed on the surface in Figure 5.

Of course, the calculations in Table 3 would have to be repeated for a number of different initial points if we wanted to reproduce the equivalent of Figure 4. Other considerations include the "stepsize," that is, the radius of the constraint circle defining the neighborhood of each point on the surface.