Pursuing a smooth implicit trajectory such as an energy level line is of great interest in numerical analysis and has attracted considerable attention      in recent years. In principle, such tracing entails marking close points on an implicit curve. Once a point is secured on the curve, a short distance forward leap is executed to place an initial point in the vicinity of the curve. Then, a returning iterative procedure is carried out to come back and land on the curve. To stay close to the curve, this stepping-out from the curve to the next initial guess is plausibly made by a tangential move at a prescribed step size. Orthogonal  iterative corrections are then carried out to bring the initially placed point ever closer to the curve.
The method, by dint of its orthogonal descent, is successful in rounding sharp bends as it follows the twists and turns of contorted trajectories. Still, the predictor leap ignores all previous information about the curving of the trajectory.
We propose here to improve the accuracy of the leap-out with a higher order parametrization of the trajectory by Bezier-like polynomial approximations. With their ability to bend, these approximations are likely to deposit an initial guess closer to the traced implicit curve for a hopefully quicker correction. We have shown that positioning an initial guess closer to the traced trajectory may well have both advantages of reducing the initial error in the following Newton-Raphson  correction, and also of improving the orthogonality of the landing approach along a gradient closer to the true one at landfall.
2. Linear Leap
Consider the implicit function of which the level curve we desire to trace. Let be a point on the curve such that . See Figure 1 below.
Let and be the partial derivatives of f with respect to x and y, respectively, at point . The total differential , for differentials , is reduced along the level line to
or in vector form
where the first vector in the dot product is the gradient to f at point A, and where the second vector is colinear with the tangent line to at point A on the curve.
Thus, the equation of the tangent line to the curve at point A is
and we propose to leap forward from point A and place predicted point B on the tangent line at distance . Evidently, for any such point is on the tangent line to at point A.
Stepping out from point A to point B is what we term a linear tangential leap.
Figure 1. Implicit level curves.
Restricting the differentials to , namely, to a tangential forward leap of size , results in
with signs chosen to produce a clockwise tracking.
For instance, if
Thus, a tangential leap originating at terminates here at point , such that
and , . See also .
3. Linearized Orthogonal Landing
We consider point as situated on the curve , having the tangent line
See Figure 1.
Point C is the intersection point of the line orthogonal to this tangent line and the curve . To reach point C from point B we seek the corrections dx and dy such that
This system of equations is approximately solved with the differential linearization
in which function f and its partial derivatives and are evaluated at point .
If is deemed not sufficiently small, then point C is taken in place of point B and the linearization is repeated.
For the unit circle the coordinates of point C are , where are the coordinates of predicted point B, and where are from Equation (12). Now
which is indeed tiny if .
To summarize: if point is on the curve, then point is the point reached by a tangential leap of size away from point A, and point is the point reached by a single orthogonal correction away from point B and towards the curve. For the unit circle
See also .
4. A Detailed Numerical Example
Consider the implicit function
which we may turn here explicit, by a mere solution of a quadratic equation, to have
At first, we look for the critical points of function f. For this we differentiate f explicitly with respect to x to have
in which . Zero slope, or , occurs at , which with yields
To determine the nature of the critical point, we implicitly differentiate Equation (17) once more with respect to x to have
At the critical point , and with x and y from Equation (18), we have from Equation (19) that
implying that the nature of the critical point is a local minimum.
Now to the leap-and-land. The equation of the tangent line to this curve at point is
Specifically, for point
A leap of brings us to point on the tangent line. At point B we write the function
The tangent line to this curve at point B is
The line orthogonal to this second tangent line is
All as seen in Figure 2.
Next we apply the landing correction
by the constrained linearization
Figure 2. Level curve of an implicit function.
for point C, which is indeed nearly on the curve. In fact instead of zero. To get point C closer to the curve, we may apply another orthogonalization and linearization. Also, a smaller leap would have resulted in point C closer to the curve.
5. Higher Order Leaps
For a quadratic leap, we start with three points already placed on the curve, say , , . Parametrization of the interpolant to the trajectory over the three points is achieved with
where are the three Lagrange interpolation functions
with parameter and , at and , respectively. At , Equation (31) predicts the coordinates
A higher order, cubic, parametrization of an interpolant to the trajectory over four points similarly produces
6. Bézier Curve Leaps
Use of partial derivatives allows the construction of a quadratic Bézier-like approximation to the trajectory  over only two points. Let and be two such points on the trajectory. The quadratic parametric approximation to both x and y between and is of the general form
where the coefficients of the approximation are to be determined from the initial and end conditions at and . Let and denote the derivatives of x and y with respect to parameter t. At any point t on the trajectory , and . For the quadratic approximation of Equation (35), , . Requiring , and at where and at where produces the equations
that are solved to yield the coefficients
where , where , and where . At we then have
7. Improved Accuracy with a Quadratic Leap
Consider the circular level line . A linear leap of size from point on top of the circle sends us to point that is at distance from the circle, provided that . A quadratic approximation through points , , with a relatively small produces out of Equation (33) the predictions
for the coordinates of point , that is at the mere distance of from the circle.
8. Now the Landing
Let now be a point not on the trajectory, . Point P can be considered as being on the level curve , with a gradient that is orthogonal to the level curve and therefore also nearly orthogonal to . Linearizing the correction we get . To achieve an orthogonal aim for the landing of the corrected point we require the colinearity condition for scalar . The linearization fixes as
The coordinates of point P are corrected thereby to , , and the procedure may be repeated.
Figure 3(a) below shows a fourteen-step tracking of the unit circle by linear leaps of and a single orthogonal Newton-Raphson landing correction.
9. A Full Numerical Experiment and a Warning
Figure 3(a) shows a fourteen-step tracking of the unit circle
by linear leaps of and a single orthogonal linear (actually Newton-Raphson) correction. The vectorial plot clearly shows the tangential leap and the orthogonal landing on the drawn circle. Figure 3 (b)
Figure 3. A leap and one orthogonal land on a unit circle.
is a vectorial plot of a nine-step quadratic leap carried out according to Equation (33), coupled with a single, linearized, orthogonal Newton-Raphson correction. Notice the relentless, undesirable growth of the step size in the quadratic method.
10. An Unstable Leap and Its Stabilization
Unfortunately, recursive formula (33) used to obtain the tracing of Figure 3(b) is unstable. Indeed, writing it in the general form
we observe that it is solved by for number z that satisfies the cubic characteristic equation
Since this equation has three equal roots ,
where depend on the initial conditions . In fact, from
we obtain that
Now, if and , then , , , and so that for all n. But if and , then , , and , so that
with a growing factor of . To control and subdue the parasitic linear growth, or contraction, of the distance between and , a proper value for the predicting t has to be chosen in Equations (31) and (32) to maintain a nearly constant step size. Instead of we set in the Lagrange interpolation functions (32) to have the prediction
and a similar expression for . Compactly written
from which we obtain
that we equate to to have a quartic equation for . For well placed initial points we may ignore the and terms in Equation (49) to be left with a quadratic equation for . Employing this stabilization procedure, we traced the unit circle as shown in Figure 3(c), using the three starting points , with . The typically computed hovered around 0.04.
11. A “Snap-Through” Oscillation of a Pre-Stressed Double Rod System
Figure 4(a) shows a system of two equal, pin-joint, supposedly massless, elastic rods of elastic constant k, fixed at pivots A and B. At their junction the rods carry a mass m that is initially held at distance from midpoint M on line AB. At this position the length of each rod is . It is assumed that at length the rods are free of stress. If , then an initial restoring pull is being exerted on m with the intention of propelling it towards point M. Assuming that the rods obey Hook's law, the energy equation of the system becomes
The unstressed length of the rods has a critical value
Figure 4. (a) A spring-mass system and (b) its computed phase portrait.
for which mass m reaches point M with zero velocity but with the rods still packing the potential energy . The system comes to an unstable symmetrical standstill waiting for a disturbance to trigger a snap-up or a snap-down. For mass m arrives at point M with a residual velocity that helps carry it over past line AB, allowing m to complete its travel to , and then back in a periodic motion.
In fact, if , then mass m reaches point M with a velocity
and through implicit differentiation of Equation (50) we have further that
indicating that as . The butterfly-shaped phase portrait of the spring-mass system in Figure 4(b) is, indeed, seen to have a waist of diminishing girth as with a turning-back that tends to become critically sharp. The computation described in Figure 4(b) was carried out with a tangential search step of size , for , for which , with , and 85 steps. Observe the deftness with which the gradient landing scheme sparsely places successive hits on flat sections of the curve, but clusters them closer around sharp bends.
We have considered in this paper the basic tangential leap followed by an orthogonal landing method for tracing implicit curves. We have then suggested a more accurate quadratic forward leap. However, such a leap may become unstable and we have shown how to restrain these possible instabilities to have a stable and accurate method. It would be of interest to consider next how to economize on the differentiations to have a more efficient method.
 Yu, Z., et al. (2006) Efficient Method for Tracing Planar Implicit Curves. Journal of Zhejiang University—Science A: Applied Physics and Engineering, 7, 1115-1123.
 Fried, I. (1984) Orthogonal Trajectory Accession to the Nonlinear Equilibrium Curve. Computer Methods in Applied Mechanics and Engineering, 47, 283-297.
 Riks, E. (1979) An Incremental Approach to the Solution of Snapping and Buckling Problems. International Journal of Solids and Structures, 15, 329-551.
 Crisfield, M.A. (1983) An Arc-Length Method Including Line Searches and Accelerations. International Journal for Numerical Methods in Engineering, 19, 1269-1289.
 Padorean, J. and Tovichackhaikul, S. (1983) Operating Characteristics of Hyperbolically and Elliptically Constrained Self-Adaptive Incremental Newton-Raphson Algorithms. Journal Franklin Institute, 316, 197-223.
 Ahn, Y.J. (2004) Approximation of Circular Arcs and Offset Curves by Bezier Curves of High Degree. Journal of Computational and Applied Mathematics, 167, 405-416.