Linear and Quadratic Leap-and-Land Trajectory Tracking Algorithms

Isaac Fried

Department of Mathematics, Boston University, Boston, MA, USA.

DOI: 10.4236/am.2021.127042                

Abstract

Considered in this note are numerical tracking algorithms for the accurate following of implicit curves. We start with a fixed point on the curve, and then systematically place on it additional points, one after the other. In this note we first go over the basic procedure of moving forward tangentially from an already placed point then orthogonally returning to the curve. Next, we further consider higher order forward stepping procedures for greater accuracy. We note, however, that higher order methods, desirable for greater accuracy, may harbor latent instabilities. This note suggests ways of holding such instabilities in check, to have stable and highly accurate tracing methods. The note has several supporting numerical examples, including the rounding of a dynamical “snap-through” point.

Keywords

Implicit Functions, Curve Tracing, Linear Tangential Leap, Nonlinear Leaps, Maintaining Stability, Orthogonal Landing, Cornering a Snap-Through

1. Introduction

Pursuing a smooth implicit trajectory such as an energy level line is of great interest in numerical analysis and has attracted considerable attention [1] [2] [3] [4] [5] 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 [2] 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 [6] 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 $z=f\left(x\mathrm{,}y\right)$ of which the level curve $z=0$ we desire to trace. Let $A\left(x_0\mathrm{,}{y}_0\right)$ be a point on the curve such that $f\left(A\right)=f\left(x_0,y_0\right)=0$. See Figure 1 below.

Let $f_x=\partial f/\partial x$ and $f_y=\partial f/\partial y$ be the partial derivatives of f with respect to x and y, respectively, at point $A\left(x_0\mathrm{,}{y}_0\right)$. The total differential $\text{d}z=f_x\text{d}x+f_y\text{d}y$, for differentials $\text{d}x\mathrm{,}\text{d}y$, is reduced along the level line to

$\text{d}z=0=f_x\text{d}x+f_y\text{d}y$ (1)

or in vector form

$\left(f_x\mathrm{,}{f}_y\right)\cdot \left(\text{d}x\mathrm{,}\text{d}y\right)=0$ (2)

where the first vector in the dot product is the gradient $\nabla f$ to f at point A, and where the second vector is colinear with the tangent line to $f=0$ at point A on the curve.

Thus, the equation of the tangent line to the curve at point A is

${\left(\frac{\partial f}{\partial x}\right)}_A\text{d}x+{\left(\frac{\partial f}{\partial y}\right)}_A\text{d}y=0\text{or}{\left(\frac{\partial f}{\partial x}\right)}_A\left(x-x_0\right)+{\left(\frac{\partial f}{\partial y}\right)}_A\left(y-y_0\right)=0$ (3)

and we propose to leap forward from point A and place predicted point B on the tangent line at distance $\epsilon$. Evidently, for any such $\text{d}x\mathrm{,}\text{d}y$ point $B\left(x_0+\text{d}x\mathrm{,}{y}_0+\text{d}y\right)$ is on the tangent line to $f\left(x,y\right)=f\left(x_0,y_0\right)=0$ 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 $\text{d}{x}^2+\text{d}{y}^2={\epsilon }^2$, namely, to a tangential forward leap of size $\epsilon$, results in

$\text{d}x=\epsilon \frac{f_y}{\sqrt{f_x^2+f_y^2}},\text{d}y=-\epsilon \frac{f_x}{\sqrt{f_x^2+f_y^2}}$ (4)

with signs chosen to produce a clockwise tracking.

For instance, if

$f\left(x\mathrm{,}y\right)=x^2+y^2-\mathrm{1,}$ (5)

then

$f_x=2x,f_y=2y$ (6)

and

$\text{d}x=\epsilon y_0,\text{d}y=-\epsilon x_0.$ (7)

Thus, a tangential leap originating at $A\left(x_0\mathrm{,}{y}_0\right)$ terminates here at point $B\left(x_1,y_1\right)=B\left(x_0+\epsilon y_0,y_0-\epsilon x_0\right)$, such that

$x_1^2+y_1^2=\left(x_0^2+y_0^2\right)\left(1+{\epsilon }^2\right)=1+{\epsilon }^2$ (8)

and $f\left(x_0\mathrm{,}{y}_0\right)=0$, $f\left(x_1,y_1\right)={\epsilon }^2$. See also [7].

3. Linearized Orthogonal Landing

We consider point $B\left(x_1\mathrm{,}{y}_1\right)$ as situated on the curve $f\left(x,y\right)=f\left(B\right)=f\left(x_1,y_1\right)$, having the tangent line

${\left(\frac{\partial f}{\partial x}\right)}_B\text{d}x+{\left(\frac{\partial f}{\partial y}\right)}_B\text{d}y=0\text{or}{\left(\frac{\partial f}{\partial x}\right)}_B\left(x-x_1\right)+{\left(\frac{\partial f}{\partial y}\right)}_B\left(y-y_1\right)=0.$ (9)

See Figure 1.

Point C is the intersection point of the line orthogonal to this tangent line and the curve $f=0$. To reach point C from point B we seek the corrections dx and dy such that

$f\left(x_1+\text{d}x,y_1+\text{d}y\right)=0\text{under}\text{the}\text{restriction}\text{that}{\left(f_y\right)}_B\text{d}x-{\left(f_x\right)}_B\text{d}y=0.$ (10)

This system of equations is approximately solved with the differential linearization

$f\left(x_1+\text{d}x,y_1+\text{d}y\right)=f\left(x_1,y_1\right)+{\left(f_x\right)}_B\text{d}x+{\left(f_y\right)}_B\text{d}y=0$ (11)

to yield

$\text{d}x=-\frac{f}{f_x^2+f_y^2}{f}_x,\text{d}y=-\frac{f}{f_x^2+f_y^2}{f}_y$ (12)

in which function f and its partial derivatives $f_x$ and $f_y$ are evaluated at point $B\left(x_1\mathrm{,}{y}_1\right)$.

If $\left|f\left(C\right)\right|$ 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 $x_2=x_1+\text{d}x,y_2=y_1+\text{d}y$, where $\left(x_1\mathrm{,}{y}_1\right)$ are the coordinates of predicted point B, and where $\text{d}x\mathrm{,}\text{d}y$ are from Equation (12). Now

$f\left(x_2,y_2\right)=x_2^2+y_2^2-1=\frac{1}{4\left(1+{\epsilon }^2\right)}{\epsilon }^4$ (13)

which is indeed tiny if $\epsilon \ll 1$.

To summarize: if point $A\left(x_0\mathrm{,}{y}_0\right)$ is on the curve, then point $B\left(x_1\mathrm{,}{y}_1\right)$ is the point reached by a tangential leap of size $\epsilon$ away from point A, and point $C\left(x_2\mathrm{,}{y}_2\right)$ is the point reached by a single orthogonal correction away from point B and towards the curve. For the unit circle

$f\left(A\right)=0,f\left(B\right)=O\left({\epsilon }^2\right),f\left(C\right)=O\left({\epsilon }^4\right).$ (14)

See also [7].

4. A Detailed Numerical Example

Consider the implicit function

$z=f\left(x,y\right)=y^2+xy+2x^2-4=0$ (15)

which we may turn here explicit, by a mere solution of a quadratic equation, to have

$y=-\frac{1}{2}x\pm \frac{1}{2}\sqrt{16-7x^2},-2.2857\le x\le \mathrm{2.2857.}$ (16)

At first, we look for the critical points of function f. For this we differentiate f explicitly with respect to x to have

$2y{y}^{\prime }+y+x{y}^{\prime }+4x=0$ (17)

in which $y^{\prime }=\text{d}y/\text{d}x$. Zero slope, or $y^{\prime }=0$, occurs at $y+4x=0$, which with $f\left(x\mathrm{,}y\right)=0$ yields

$x=\sqrt{\frac{2}{7}}=0.5345,y=-4\sqrt{\frac{2}{7}}=-\mathrm{2.138.}$ (18)

To determine the nature of the critical point, we implicitly differentiate Equation (17) once more with respect to x to have

$2y^{\prime }{y}^{\prime }+2y{y}^{″}+y^{\prime }+y^{\prime }+x{y}^{″}+4=0.$ (19)

At the critical point $y^{\prime }=0$, and with x and y from Equation (18), we have from Equation (19) that

$y^{″}=\frac{2}{7}\sqrt{14}>0$ (20)

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 $A\left(x_0\mathrm{,}{y}_0\right)$ is

${\left(\frac{\partial f}{\partial x}\right)}_0\text{d}x+{\left(\frac{\partial f}{\partial y}\right)}_0\text{d}y=0,\text{d}x=x-x_0,\text{d}y=y-y_0$ (21)

or here

$\left(4x_0+y_0\right)\text{d}x+\left(x_0+2y_0\right)\text{d}y=0.$ (22)

Specifically, for point $A(\; 0,2\; )$

$2\text{d}x+4\text{d}y=0,\text{or}y=-\frac{1}{2}x+2.$ (23)

A leap of $\epsilon =1.2$ brings us to point $B\left(\mathrm{1.07331,1.46334}\right)$ on the tangent line. At point B we write the function

$f\left(x,y\right)=y^2+xy+2x^2-k=0,k=y_1^2+x_1{y}_1+2x_1^2=\mathrm{6.016.}$ (24)

The tangent line to this curve at point B is

${\left(\frac{\partial f}{\partial x}\right)}_1\text{d}x+{\left(\frac{\partial f}{\partial y}\right)}_1\text{d}y=0,\text{d}x=x-x_1,\text{d}y=y-y_1$ (25)

or specifically

$5.75659\text{d}x+4\text{d}y=0.$ (26)

The line orthogonal to this second tangent line is

$4\text{d}x-5.75659\text{d}y=0.$ (27)

All as seen in Figure 2.

Next we apply the landing correction

$f\left(x_1+\text{d}x\mathrm{,}{y}_1+\text{d}y\right)=0$ (28)

by the constrained linearization

${\left(\frac{\partial f}{\partial x}\right)}_1\text{d}x+{\left(\frac{\partial f}{\partial y}\right)}_1\text{d}y=0,4\text{d}x-5.75659\text{d}y=0$ (29)

Figure 2. Level curve of an implicit function.

and obtain

$\text{d}x=-0.236176,\text{d}y=-0.164108\text{or}{x}_2=0.837137,y_2=1.29924$ (30)

for point C, which is indeed nearly on the curve. In fact $f\left(C\right)=0.177248$ instead of zero. To get point C closer to the curve, we may apply another orthogonalization and linearization. Also, a smaller leap $\epsilon$ 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 $P_1\left(x_1\mathrm{,}{y}_1\right)$, $P_2\left(x_2\mathrm{,}{y}_2\right)$, $P_3\left(x_3\mathrm{,}{y}_3\right)$. Parametrization of the interpolant to the trajectory over the three points is achieved with

$x\left(t\right)=x_1{\varphi }_1\left(t\right)+x_2{\varphi }_2\left(t\right)+x_3{\varphi }_3\left(t\right),y\left(t\right)=y_1{\varphi }_1\left(t\right)+y_2{\varphi }_2\left(t\right)+y_3{\varphi }_3\left(t\right)$ (31)

where ${\varphi }_1\mathrm{,}{\varphi }_2\mathrm{,}{\varphi }_3$ are the three Lagrange interpolation functions

${\varphi }_1\left(t\right)=\frac{1}{2}\left(t-1\right)\left(t-2\right),{\varphi }_2\left(t\right)=-t\left(t-2\right),{\varphi }_3\left(t\right)=-\frac{1}{2}t\left(t-1\right)$ (32)

with parameter $t=0,t=1$ and $t=2$, at $P_1\mathrm{,}{P}_2$ and $P_3$, respectively. At $t=3$, Equation (31) predicts the coordinates

$x_4=x\left(t=3\right)=x_1+3\left(x_3-x_2\right),y_4=y\left(t=3\right)=y_1+3\left(y_3-y_2\right).$ (33)

A higher order, cubic, parametrization of an interpolant to the trajectory over four points similarly produces

$\begin{array}{l}{x}_5=x\left(t=4\right)=-x_1+4x_2-6x_3+4x_4,\\ y_5=y\left(t=4\right)=-y_1+4y_2-6y_3+4y_4.\end{array}$ (34)

6. Bézier Curve Leaps

Use of partial derivatives allows the construction of a quadratic Bézier-like approximation to the trajectory [8] over only two points. Let $P_1\left(x_1\mathrm{,}{y}_1\right)$ and $P_2\left(x_2\mathrm{,}{y}_2\right)$ be two such points on the trajectory. The quadratic parametric approximation to both x and y between $P_1$ and $P_2$ is of the general form

$x=a_0+a_{1}t+a_2{t}^2,y=b_0+b_{1}t+b_2{t}^2,0\le t\le 1$ (35)

where the coefficients of the approximation are to be determined from the initial and end conditions at $P_1$ and $P_2$. Let $\stackrel{\dot{}}{x}$ and $\stackrel{\dot{}}{y}$ denote the derivatives of x and y with respect to parameter t. At any point t on the trajectory $x=x\left(t\right),y=y\left(t\right)$, and $f_x\stackrel{\dot{}}{x}+f_y\stackrel{\dot{}}{y}=0$. For the quadratic approximation of Equation (35), $\stackrel{\dot{}}{x}=a_1+2a_{2}t$, $\stackrel{\dot{}}{y}=b_1+2b_{2}t$. Requiring $x\mathrm{,}y\mathrm{,}\stackrel{\dot{}}{x}$, and $\stackrel{\dot{}}{y}$ at $P_1$ where $t=0$ and at $P_2$ where $t=1$ produces the equations

$x_1=a_0,y_1=b_0$

$x_2=a_0+a_1+a_2,y_2=b_0+b_1+b_2$ (36)

$f_{x1}{a}_1+f_{y1}{b}_1=0,f_{x2}\left(a_1+2a_2\right)+f_{y2}\left(b_1+2b_2\right)=0$

that are solved to yield the coefficients

$x_1=a_0,y_1=b_0$

$a_1=-\frac{2}{\Delta }f{y}_1\left(f_{x2}\delta x+f_{y2}\delta y\right),b_1=\frac{2}{\Delta }f{x}_1\left(f_{x2}\delta x+f_{y2}\delta x\right)$ (37)

$a_2=\delta x-a_1,b_2=\delta y-b_1$

where $\Delta =f_{x1}{f}_{y2}-f_{y1}{f}_{x2}$, where $\delta x=x_2-x_1$, and where $\delta y=y_2-y_1$. At $t=2$ we then have

$x_3=x\left(t=2\right)=-3x_1+4x_2-3a_1,y_3=-3y_1+4y_2-3b_1.$ (38)

7. Improved Accuracy with a Quadratic Leap

Consider the circular level line $z=f\left(x,y\right)=x^2+y^2-1=0$. A linear leap of size $\epsilon$ from point $P\left(\mathrm{0,1}\right)$ on top of the circle sends us to point $Q\left(\epsilon \mathrm{,1}\right)$ that is at distance $\frac{1}{2}{\epsilon }^2$ from the circle, provided that $\epsilon \ll 1$. A quadratic approximation through points $P_1\left(-\sin 2\theta \mathrm{,}\cos 2\theta \right)$, $P_2\left(-\sin \theta \mathrm{,}\cos \theta \right)$, $P_3\left(\mathrm{0,1}\right)$ with a relatively small $\theta$ produces out of Equation (33) the predictions

$x_4=x\left(\theta \right)=\theta +\frac{5}{6}{\theta }^3\text{and}{y}_4=y\left(\theta \right)=1-\frac{1}{2}{\theta }^2$ (39)

for the coordinates of point $P_4$, that is at the mere distance of $\left(23/24\right){\theta }^4$ from the circle.

8. Now the Landing

Let $P\left(x_0\mathrm{,}{y}_0\right)$ now be a point not on the trajectory, $f\left(x_0\mathrm{,}{y}_0\right)\ne 0$. Point P can be considered as being on the level curve $z=f\left(x,y\right)=f\left(x_0,y_0\right)$, with a gradient $\nabla f=\left(f_{x0}\mathrm{,}{f}_{y0}\right)$ that is orthogonal to the level curve $f\left(x,y\right)=f\left(x_0,y_0\right)$ and therefore also nearly orthogonal to $f\left(x\mathrm{,}y\right)=0$. Linearizing the correction $f\left(x_0+\text{d}x\mathrm{,}{y}_0+\text{d}y\right)=0$ we get $f_0+f_{x0}\text{d}x+f_{y0}\text{d}y=0$. To achieve an orthogonal aim for the landing of the corrected point we require the colinearity condition $\left(\text{d}x\mathrm{,}\text{d}y\right)=\lambda \left(f_{x0}\mathrm{,}{f}_{y0}\right)$ for scalar $\lambda$. The linearization fixes $\lambda$ as

$\lambda =-\frac{f_0}{f_{x0}^2+f_{y0}^2}\text{and}\text{d}x=\lambda f_{x0},\text{d}y=\lambda f_{y0}.$ (40)

The coordinates of point P are corrected thereby to $x_1=x_0+\text{d}x$, $y_1=y_0+\text{d}x$, and the procedure may be repeated.

Figure 3(a) below shows a fourteen-step tracking of the unit circle $f\left(x,y\right)=x^2+y^2-1=0$ by linear leaps of $\epsilon =0.49$ 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

$f\left(x,y\right)=x^2+y^2-1=0$ by linear leaps of $\epsilon =0.49$ 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

$x_{n+3}-3x_{n+2}+3x_{n+1}-x_n=0,n=1,2,3,\cdots$ (41)

we observe that it is solved by $x_n=z^n$ for number z that satisfies the cubic characteristic equation

$z^3-3z^2+3z-1={\left(z-1\right)}^3=0.$ (42)

Since this equation has three equal roots $z_1=z_2=z_3=1$,

$x_n=c_1+c_{2}n+c_3{n}^2$ (43)

where $c_1\mathrm{,}{c}_2\mathrm{,}{c}_3$ depend on the initial conditions $x_1\mathrm{,}{x}_2\mathrm{,}{x}_3$. In fact, from

$x_1=c_1+c_2+c_3,x_2=c_1+2c_2+4c_3,x_3=c_1+3c_2+9c_3$ (44)

we obtain that

$c_1=3x_1-3x_2+x_3,c_2=\frac{1}{2}\left(-5x_1+8x_2-3x_3\right),c_3=\frac{1}{2}\left(x_1-2x_2+x_3\right).$ (45)

Now, if $x_1=0,x_2=1$ and $x_3=2$, then $c_1=-1$, $c_2=1$, $c_3=0$, and $x_n=-1+n$ so that $x_{n+1}-x_n=1$ for all n. But if $x_1=0,x_2=1$ and $x_3=2+\epsilon$, then $c_1=1+\epsilon$, $c_2=1-\frac{3}{2}\epsilon$, $c_3=\frac{1}{2}\epsilon$ and $x_n=n-1+\frac{1}{2}\epsilon \left(n-1\right)\left(n-2\right)$, so that

$x_{n+1}-x_n=1+\epsilon \left(n-1\right)$ with a growing factor of $n+1$. To control and subdue the parasitic linear growth, or contraction, of the distance between $x_n$ and $x_{n+1}$, 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 $t=3$ we set $t=3+\epsilon$ in the Lagrange interpolation functions (32) to have the prediction

$x_4=x_4\left(\epsilon \right)=x_1\left(1+\frac{3}{2}\epsilon +\frac{1}{2}{\epsilon }^2\right)+x_2\left(-3-4\epsilon -{\epsilon }^2\right)+x_3\left(3+\frac{5}{2}\epsilon +\frac{1}{2}{\epsilon }^2\right)$ (46)

or

$x_4-x_3=x_1-3x_2+2x_3+\epsilon \left(\frac{3}{2}{x}_1-4x_2+\frac{5}{2}{x}_3\right)+{\epsilon }^2\left(\frac{1}{2}{x}_1-x_2+\frac{1}{2}{x}_3\right)$ (47)

and a similar expression for $y_4-y_3$. Compactly written

$x_4-x_3=a_1+a_2\epsilon +a_3{\epsilon }^2,y_4-y_3=b_1+b_2\epsilon +b_3{\epsilon }^2$ (48)

from which we obtain

$\begin{array}{l}{\left(x_4-x_3\right)}^2+{\left(y_4-y_3\right)}^2\\ =a_1^2+b_1^2+2\left(a_1{a}_2+b_1{b}_2\right)\epsilon +\left(a_2^2{b}_2^2+2a_1{a}_3+2b_1{b}_3\right){\epsilon }^2\\ +2\left(a_2{a}_3+b_2{b}_3\right){\epsilon }^3+\left(a_3^2+b_3^2\right){\epsilon }^4\end{array}$ (49)

that we equate to ${\left(x_2-x_1\right)}^2+{\left(y_2-y_1\right)}^2$ to have a quartic equation for $\epsilon$. For well placed initial points $\left|\epsilon \right|\ll 1$ we may ignore the ${\epsilon }^3$ and ${\epsilon }^4$ terms in Equation (49) to be left with a quadratic equation for $\epsilon$. Employing this stabilization procedure, we traced the unit circle as shown in Figure 3(c), using the three starting points $P_1\left(0,1\right),P_2\left(\sin \theta ,\cos \theta \right),P_2\left(\sin 2\theta ,\cos 2\theta \right)$, with $\theta =\pi /15$. The typically computed $\epsilon$ 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 $x_0$ from midpoint M on line AB. At this position the length of each rod is $l\left(x_0\right)=\sqrt{1+x_0^2}$. It is assumed that at length $l_0$ the rods are free of stress. If $l\left(x_0\right)>l_0$, 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

$\frac{2k}{m}{\left(\sqrt{1+y_0^2}-l_0\right)}^2=\frac{2k}{m}{\left(\sqrt{1+x^2}-l_0\right)}^2+{\stackrel{\dot{}}{x}}^2\mathrm{.}$ (50)

The unstressed length of the rods has a critical value

$l_{0crit}=\frac{1}{2}\frac{x_0^2}{\sqrt{1+x_0^2}-1}$ (51)

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 $k{\left(1-l_0\right)}^2$. The system comes to an unstable symmetrical standstill waiting for a disturbance to trigger a snap-up or a snap-down. For $l_0<l_{0crit}$ 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 $x=-x_0$, and then back in a periodic motion.

In fact, if $l_0=l_{0crit}\left(1-\epsilon \right),\epsilon >0$, then mass m reaches point M with a velocity

$\stackrel{\dot{}}{x}\left(x=0\right)=\frac{2k}{m}\epsilon x_0^2$ (52)

and through implicit differentiation of Equation (50) we have further that

$\ddot{x}\left(x=0\right)=0,\text{and}\stackrel{⃛}{x}\left(x=0\right)=\frac{l_0-1}{\epsilon x_0^2}$ (53)

indicating that $\stackrel{⃛}{x}\left(x=0\right)\to \infty$ as $\epsilon \to 0$. 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 $l_0\to l_{0crit}$ 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 $\epsilon =0.1$, for $x_0=\tan \left(60˚\right)$, for which $l_{0crit}=1.5$, with $l_0=1.499$, 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.

12. Conclusion

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.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	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. https://doi.org/10.1631/jzus.2006.A1115
[2]	Fried, I. (1984) Orthogonal Trajectory Accession to the Nonlinear Equilibrium Curve. Computer Methods in Applied Mechanics and Engineering, 47, 283-297. https://doi.org/10.1016/0045-7825(84)90080-X
[3]	Riks, E. (1979) An Incremental Approach to the Solution of Snapping and Buckling Problems. International Journal of Solids and Structures, 15, 329-551. https://doi.org/10.1016/0020-7683(79)90081-7
[4]	Crisfield, M.A. (1981) A Fast Incremental/Iterative Solution Procedure That Handles “Snap-Through”. Computers and Structures, 13, 55-62. https://doi.org/10.1016/0045-7949(81)90108-5
[5]	Crisfield, M.A. (1983) An Arc-Length Method Including Line Searches and Accelerations. International Journal for Numerical Methods in Engineering, 19, 1269-1289. https://doi.org/10.1002/nme.1620190902
[6]	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. https://doi.org/10.1016/0016-0032(83)90041-8
[7]	Watson, L.T. and Holzer, S.M. (1983) Quadratic Convergence of Crisfield Method. Computers and Structures, 17, 69-72. https://doi.org/10.1016/0045-7949(83)90030-5
[8]	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. https://doi.org/10.1016/j.cam.2003.10.008