Computing Curves Bounding Trigonometric Planar Maps: Symbolic and Hybrid Methods

Daniel Lichtblau

A few years ago S-H Kim investigated some problems at the boundary of number theory, optimization, and geometry. One question regarded an optimal packing of certain “triangular oval” planar curves and another looked at some related transformations of (ℝ^2) to (ℝ^2). These were investigated primarily using tools from calculus but it turns out that computational algebra methods may instead be employed to particular advantage. Moreover, generalizations that are beyond the reach of such methods are still amenable to hybrid approaches using numeric and symbolic methods in tandem. We introduce some of the specific problems and generalizations, and show by detailed example how such techniques may be implemented and deployed.

1. Introduction and Overview

A common need in two and three dimensional computational geometry is to bound regions given by algebraic maps. For example, while it is typically easy to find the implicit form of a curve given parametrically by rational functions of a variable (or trigonometric functions, as discussed at length in [Hong and Schicho 1998]), it is not so well understood how to find boundary curves of maps between planar regions. Perhaps more to the point, known symbolic methods already become overwhelmed by complexity in examples of modest size. Hence there is a desire to find methods less prone to the vagaries of intermediate swell and related pitfalls of symbolic computation.

What we develop will fall into the category of computational methods. We cite theory when relevant but do not purport to develop any herein. That said, the hybrid methods to be presented draw from several areas of numeric and symbolic computation (optimization, differential equation solvers, linear algebra, lattice methods, and more), and illustrate nicely the advantages of such synergy. Moreover, as will be seen from the family of examples we study, these methods can then be applied to the study of mathematical problems.

In particular, we will investigate certain properties of trigonometric planar mappings of the form

x=cos(m a)+cos(m b)+cos(m (a-b))
y=cos(n a)+cos(n b)+cos(n (a-b))

for fixed integers m and n. We begin with a history and analysis of a particular simple case, illustrating symbolic methods that have been presented in previous venues. We proceed to develop a hybrid approach that applies when the symbolic methods are no longer up to the task.

We remark that similar maps e.g. replacing cosine with sine are also of interest and can be addressed using the same methods we will develop.

2. Background of the Problems Under Consideration

We begin with a quick review of material presented in [Lichtblau 2003b], based in part on [Kim 1998], [Kim 2003]. Given

f(n)=(Underscript[Min, (Subscript[a, j])ϵ ℝ]) Max [|(∑_(j=1)^n)(e^(i (Subscript[a, j]) ))|,| (∑_(j=1)^n)(e^(i n (Subscript[a, j]) ))|]

We want to say something useful about this function f(n).

Theorem (paraphrased).
(i) f(1)=1.
(ii) f(2)=1.
(iii) For n>3 we have f(n)=0.

Proofs are in [Kim 2000] and [Kim 2003]. One might be tempted to state, as a corollary, that the only case of interest is when n=3. We will see presently why this is not so. First we discuss the case n=3, which was of interest in Kim’s work and also lays the groundwork for the more difficult problems we will attack.

Proposition. For f(3)=(Underscript[Min, a,b∈[0,2π]])Max[|1+(e^(i a))+(e^(i b))|,|1+(e^(3i a))+(e^(3i b))|] the two moduli are equal.

Theorem. f(3) occurs at a tangential intersection of the curves

cos(a)+cos(b)+cos(a-b)=k
cos(3a)+cos(3b)+cos(3(a-b))=k

for values {a,b,k}. [These new functions are half of (the squares of the moduli - 3). Once we find k, we recover f(3) as (Sqrt[2k+3]).] Moreover at extremal values of f(3) there will be only finitely many such intersection points, up to periodicity.

Proofs again are in [Kim 2000], [Kim 2003]. The utility of this theorem is that it points the way to an algorithm because there are only finitely many values of k for which the intersections will be tangential, and we will see shortly how to find them via elimination theory.

In the theorem above we have parametrized pairs of trigonometric polynomials associated to pairs of moduli of exponential sums. As we seek a minmax it is not surprising that this would occur at parameter values that give tangential intersections of the pairs. We show some trigonometric curve pairs for k=-1.3 and k=-1.1. Computations below are done with the version 5 of Mathematica [Wolfram 2003] (Mathematica (TM) is a registered trademark of Wolfram Research, Inc.).

Needs["Graphics`"]
trigpoly[n_,a_,b_]=Cos[n a]+Cos[n b]+Cos[n (a-b)];
grad[expr_,a_,b_]:={(∂_a)expr,(∂_b)expr}
trigsubs={Cos[a]->(Subscript[c, a]),Sin[a]->(Subscript[s, a]),Cos[b]->(Subscript[c, b]),Sin[b]->(Subscript[s, b])};

pl1=ImplicitPlot[{trigpoly[1,a,b]==-1.1,trigpoly[3,a,b]==-1.1},{a,0,π},{b,π,2 π},PlotPoints->200,AspectRatio->1,PlotStyle->{{Thickness[0.012],GrayLevel[0.5],Dashing[{0.25,0.1}]},{Automatic}},DisplayFunction->Identity];
pl2=ImplicitPlot[{trigpoly[1,a,b]==-1.3,trigpoly[3,a,b]==-1.3},{a,0,π},{b,π,2 π},PlotPoints->200,AspectRatio->1,PlotStyle->{{Thickness[0.012],GrayLevel[0.5],Dashing[{0.25,0.1}]},{Automatic}},DisplayFunction->Identity];
Show[GraphicsArray[{pl1,pl2}],DisplayFunction->$DisplayFunction];

The level sets are apparently convex closed curves and for some value of -1.3<k<-1.1 we will have tangential intersections which, it will turn out, gives the optimal value we seek. As the pictures provide a hint that the level curves are smooth and convex, that of course is already useful in pointing the way to the theory. But this is just the tip of the computational iceberg in regards to the problems that will come under scrutiny. To introduce them we will next find tangential intersections of these level sets by means of Lagrange multipliers.

The gist of the computation is as follows. We begin with the two curve equations from the theorem above. As the theorem tells us the intersection of interest is tangential, we augment with the standard Lagrange multiplier relation between their gradients. We transform the equations to explicitly algebraic ones by rewriting e.g. cos(a) as (Subscript[c, a]). We must also augment by polynomials corresponding to algebraic relations such as (Subscript[c, a]^2)+(Subscript[s, a]^2)=1. From these we construct a Gröbner basis to eliminate all variables but the one of interest, using a term order efficient for this task. Specifically, the variables to be eliminated are all weighted equally, and lexicographically larger than the ones that remain. Ties are broken with graded reverse lexicographic ordering on the join of the set of elimination variables with the set of remaining variables. (This is the elimination order of Bayer and Stillman that is utilized to eliminate all but one variable. For background on construction and use of Gröbner bases, including elimination of variables from algebraic systems, see e.g. [Buchberger 1985], [Cox, Little, O’Shea 1992]. The latter reference discusses this ordering in the exercises in chapter 3, section 1.)

{e1,e2}=Map[TrigExpand[trigpoly[#,a,b]]&,{1,3}];
polys=Flatten[{e1-k,e2-k,grad[e1,a,b]-λ grad[e2,a,b],((Cos[a])^2)+((Sin[a])^2)-1,((Cos[b])^2)+((Sin[b])^2)-1}]/.trigsubs

Timing[kpoly=First[GroebnerBasis[polys,k,{λ,(Subscript[c, a]),(Subscript[s, a]),(Subscript[c, b]),(Subscript[s, b])},MonomialOrder->EliminationOrder]]]

We will factor it to see what values of k might be of interest.

Factor[kpoly]

The values -1, 1, and 3 do not in fact lead to the optimum for the original problem (at k=-1 the curves are particularly misbehaved insofar as the intersections are not even discrete). The root of interest turns out to be the real value for the cubic factor; as it is the most negative root it will give the minimal value for f(3).

solns=NSolve[kpoly==0,k]

A picture of the curve intersections for k≈-1.2 will show they intersect tangentially.

ImplicitPlot[{trigpoly[1,x,y]==k/.First[solns],trigpoly[3,x,y]==k/.First[solns]},{x,0,π},{y,π,2 π},PlotPoints->200,AspectRatio->1,PlotStyle->{{Thickness[0.012],GrayLevel[0.5],Dashing[{0.25,0.1}]},{Automatic}}];

Note that we have what appear to be both horizontal and vertical tangencies (they really are). Hence we would have had trouble had we tried to solve for either variable as an implicit function of the other. One conclusion is that Lagrange multipliers constitute a good approach to this problem. Also notice that we have a “maximal” sort of lattice packing insofar as we cannot further expand the curve components without having their interiors intersect. This link between a problem in number theory and a lattice packing was an emphasis of [Kim 2003]. We indicate pictorially the situation for the other roots in k. The case k=-1 gives intersections along line segments. For k=1 we appear to have solved a different tangential intersection problem. For k=3 we actually have degeneracies to points, which we show by plotting at k=2.9 and observing small sets of ovals.

pl3=ImplicitPlot[{trigpoly[1,a,b]==-1,trigpoly[3,a,b]==-1},{a,0,4π},{b,0,4 π},PlotPoints->200,AspectRatio->1,PlotStyle->{{Thickness[0.012],GrayLevel[0.5],Dashing[{0.25,0.1}]},{Automatic}},DisplayFunction->Identity];
pl4=ImplicitPlot[{trigpoly[1,a,b]==1,trigpoly[3,a,b]==1},{a,0,4π},{b,0,4 π},PlotPoints->200,AspectRatio->1,PlotStyle->{{Thickness[0.012],GrayLevel[0.5],Dashing[{0.25,0.1}]},{Automatic}},DisplayFunction->Identity];
pl5=ImplicitPlot[{trigpoly[1,a,b]==2.9,trigpoly[3,a,b]==2.9},{a,0,4π},{b,0,4π},PlotPoints->200,AspectRatio->1,PlotStyle->{{Thickness[0.012],GrayLevel[0.5],Dashing[{0.25,0.1}]},{Automatic}},DisplayFunction->Identity];
Show[GraphicsArray[{pl3,pl4,pl5}],DisplayFunction->$DisplayFunction];

This sets the stage for the geometry problems we will consider. The focus will be on systems of equations of the type shown above that are based on Lagrange multiplier methods. This will allow us to find curves that satisfy certain extremal value conditions.

3. The Associated Planar Map

We now investigate the function from (ℝ^2) to (ℝ^2) that is motivated by the last example.

x[a_,b_]=trigpoly[1,a,b];
y[a_,b_]=trigpoly[3,a,b];

We plot the images of around 4000 random points from the square [0,2π] × [0,2π].

pointlist[m_]:=Table[2*Pi*{Random[],Random[]},{2^m}]
points2k=pointlist[12];
data[{point__}]:={x[point],y[point]};
datalist2k=Map[data,points2k];
ListPlot[datalist2k];

We have something a bit like a Lissajous figure (insofar as those also arise from trigonometric parametrizations). At {x,y}={-1,-1} we see the image pinched to a point as boundary curves cross. This may help to explain that dimensional component of solutions to tangential intersection problem of the previous section.

The problem we now consider is to compute the boundary curve(s). Observe that our map may be regarded as rational functions of (Subscript[z, a])=(e^(i a)) and (Subscript[z, b])=(e^(i b)). Hence it is a rational map. As the boundary curves will be obtained from algebraic relations involving this map, they too are algebraic (so in a sense we are working with a bivariate flavor of [Hong and Schicho 1998]). To find these boundaries we again turn to Gröbner bases. We have two viable approaches. One is to see where the Jacobian of our map vanishes. The other possibility, similar to that used in the last section, is to find extremal values. We pursue this latter approach.

The idea is to write x and y as functions of trigonometric polynomials, again make substitutions so these become explicitly algebraic, and realize that at the extrema (that is, on the boundary curve(s)), the gradients of the two functions must be parallel. (Reason: on the boundary, for fixed x=x(a,b) we extremize y=y(a,b). This can be set up as a standard Lagrange multiplier problem with ∇y=λ∇x.) Note that one might attempt to influence the computational efficiency by converting from trigonometric to exponential polynomials, thus reducing the number of polynomials and variables in the system. This is a method advocated in [Hong and Schicho 1998]. Our experience with this particular example is that the trigonometric formulation has the advantage. But one should bear in mind that this might be dependent on both Gröbner basis implementation and problem specific details.

Here is the trigonometric formulation.

parameters={a,b};mainvars={x,y};
xpoly=x-trigpoly[1,a,b];ypoly=y-trigpoly[3,a,b];
xypolys={xpoly,ypoly};
gradientpolys=grad[xpoly,a,b]-λ grad[ypoly,a,b];
trigidentities={((Subscript[c, a])^2)+((Subscript[s, a])^2)-1,((Subscript[c, b])^2)+((Subscript[s, b])^2)-1};
elimvars={λ,(Subscript[c, a]),(Subscript[c, b]),(Subscript[s, a]),(Subscript[s, b])};
polys=TrigExpand[Join[trigidentities,xypolys,gradientpolys]]/.trigsubs

We now want to compute the polynomial in {x,y} in the elimination ideal of the system as per [Cox, Little, O’Shea 1992]. That is, we eliminate from the system above all other variables. As before, this can be attempted via a Gröbner basis computation. Note that this is in effect a generalization of implicitization a la [Kalkbrenner 1990] insofar as we have a parametrization of the region rather than just its boundary.

Timing[gb=GroebnerBasis[polys,mainvars,elimvars,Sort->True,MonomialOrder->EliminationOrder]]

Notice that we are left with one polynomial in the two variables, which means we indeed found a curve or curves. We will factor it to better understand what we have found.

envelope=Map[#[[1]]^#[[2]]&,Drop[FactorList[First[gb]],1]]

We have a cubic function of x and an implicit relation that is quadratic in y. One may observe that the latter is not exactly trivial and is unlikely to be found without recourse to computer algebra. We next plot these curves to see specifically how they bound our map.

ImplicitPlot[Evaluate[Thread[envelope==0]],{x,-1.5,3.2},{y,-2,3.2},PlotPoints->2000,AspectRatio->1,PlotStyle->{{Thickness[.01],GrayLevel[.5],Dashing[{.15,.04}]},{Automatic}}];

Before proceeding to the general case we use the picture to explain those nonminimal critical values for k derived in the previous section. The idea is to intersect with the line y=x to see what happens at level curve intersections for particular values of k. Clearly {1,1} is a local maximum. For larger values the level curves do not intersect, but they gradually shrink to (colliding) points at the global maximum of {3,3}. The remaining critical value is at {-1,-1}. From the graph we realize that when x=-1 we must have y=-1, but not vice versa. Hence the set of level curves for y=-1 contains those for x=-1 (one may wish to review the picture of level curves corresponding to k=-1; the larger dashed triangle boundaries are completely covered by their smaller solid counterparts, but not vice versa). As one set of level curves contains the other, we have a dimensional component in the intersection for this value.

4. Generalizing the Planar Maps

We have set the stage for the more general problem. While the work in [Kim 2003] shows that, for the purposes of the original minmax problem, the cases of n>3 are not relevant, in fact those give interesting planar maps from which one might wish to extract boundary curves. As an example we illustrate the n=5 case below using around 4000 points.

x[a_,b_]=trigpoly[1,a,b];y[a_,b_]=trigpoly[5,a,b];
pointlist[m_]:=Table[2*Pi*{Random[],Random[]},{2^m}]
points2k=pointlist[12];
data[{point__}]:={x[point],y[point]};
datalist2k=Map[data,points2k];
lplot=ListPlot[datalist2k];

One can discern more complexity in this plot, insofar as there appear to be “ghost” curves inside the region that perhaps correspond to curves in the domain where the Jacobian fails to have full rank. It turns out that for n=4 we are already unable to find the boundary curves in reasonable time using the method of the previous section. After conversion to exponentials a la [Hong and Schicho 1998], the program Fermat was able to do the variable elimination via a Dixon resultant computation in about 15 minutes [Lewis, private communication]. For n=5 it ran out of memory. While it is plausible that other programs might progress further, it seems unlikely that any will go far beyond the range for which results have been obtained by symbolic methods. Hence the motivation for the numeric-symbolic approach we describe below. We take for our working example the case where n=5.

xpoly=x-trigpoly[1,a,b];ypoly=y-trigpoly[5,a,b];
xypolys={xpoly,ypoly};
gradientpolys=grad[xpoly,a,b]-λ grad[ypoly,a,b];
polys=TrigExpand[Join[trigidentities,xypolys,gradientpolys]]/.trigsubs;

Now we may set values for y and solve for the rest, retaining solutions in x. The method we first show was used to advantage in a computational geometry setting in [Lichtblau 2004], in a case where standard root finding approaches were problematic.

allbutyvars={x,λ,(Subscript[c, a]),(Subscript[c, b]),(Subscript[s, a]),(Subscript[s, b])};
polyy=polys/.y->11/10;
Timing[soln=NSolve[polyy,allbutyvars,Sort->True];]

Select[Union[x/.soln,SameTest->(Abs[#1-#2]/(Abs[#1]+Abs[#2])<10^(-5)&)],Im[#]==0&]

While this does give results in modest time we perceive two difficulties. One is that it will take a very long time to get sufficiently many points to use to recover the curves. Another is that this method will not extend much further; it is quite sensitive to the degrees of the input polynomials and thus for larger n it is not likely to give us solution points in any reasonable amount of time. Hence we drop the idea of solving for several values of x for a given y, and focus instead on obtaining one such solution at a time. The most naive approach would be to use a simple root finder. These tend to be local methods and require reasonable initial guesses for solution values in order to converge to a root. Once we select a value, say for y, it is not hard to see what might be a reasonable value for x. But the best we can say for the “trigonometric” parameters is that they lie in the range between -1 and 1. Even worse, we have no information for the Lagrange multiplier. In such a situation it will turn out to be preferable to use a global method that allows for constraints on values. While one might use e.g. interval methods for constrained global root finding, in Mathematica it is more convenient to use the optimization function NMinimize. We have some latitude in how we do this. For example, we may use explicit constraints to enforce some of the polynomial equations, and take as objective function the sum of squares of the remaining polynomials (we thus hope to find solution values for which this is close to zero), perhaps plus a constant. Some trial and error experimentation indicates that it is convenient to enforce all of the equations as constraints, working with just a constant objective function. Admittedly this is not the most obvious sort of thing to set up; further information regarding use of NMinimize may be found in [Champion 2002]. Note that here it makes good sense to work with the trigonometric parametrization as we have a very good set of ranges for those variables. With a parametrization in terms of complex exponentials we would need to work separately with real and imaginary parts (in effect returning to the trigonometric parametrization), whereas with the “usual” rational parametrization of a torus our parameters would take on values in an infinite range.

Timing[{val,root}=NMinimize[{1,Flatten[{Thread[polyy==0],.9≤ x≤1.2}]},{{x,.9,1.2},{λ,-10,10},{(Subscript[c, a]),-1,1},{(Subscript[c, b]),-1,1},{(Subscript[s, a]),-1,1},{(Subscript[s, b]),-1,1}},MaxIterations->200,Method->"DifferentialEvolution"]]

At this point we will note that the choice of y=11/10 is very likely “generic”, and hence the fact that the solution values clearly indicate b=0 means that we are very likely on a boundary segment in the {a,b} domain. Thus it should not be difficult to find this part of the boundary simply by implicitizing the domain boundary curve b=0, discarding from our polynomial system the Lagrange multiplier part. We show this now.

envelopepiece1=First[GroebnerBasis[Take[polys,4]/.{(Subscript[c, b])->1,(Subscript[s, b])->0},{x,y},{ (Subscript[c, a]),(Subscript[s, a])}]]

A glance at the plot will convince us that a part of it comprises a piece of the map boundary.

envplot1=ImplicitPlot[Evaluate[envelopepiece1==0],{x,-2,3.2},{y,-2,3.2},PlotStyle->{Thickness[.01],GrayLevel[.5],Dashing[{.15,.07}]},PlotPoints->1000];

That piece was relatively easy due to a stroke of luck: we deduced that it was the image of a boundary segment. We now look at that part of the boundary in the bottom right portion of the picture. We note that for y=-11/10 we will have 1.5<x<2.

polyy=polys/.y->-11/10;
Timing[{val,root}=NMinimize[{1,Flatten[{Thread[polyy==0],1.2≤ x≤2}]},{{x,1.2,2},{λ,-10,10},{(Subscript[c, a]),-1,1},{(Subscript[c, b]),-1,1},{(Subscript[s, a]),-1,1},{(Subscript[s, b]),-1,1}},MaxIterations->200,Method->"DifferentialEvolution"]]

This time the solutions for the trigonometric parameters suggest the segment b=-a in the domain. We obtain the boundary for this piece using the same method as above.

envelopepiece2=First[GroebnerBasis[Take[polys,4]/.{(Subscript[c, b])->(Subscript[c, a]),(Subscript[s, b])->-(Subscript[s, a])},{x,y},{ (Subscript[c, a]),(Subscript[s, a])}]]

envplot2=ImplicitPlot[Evaluate[envelopepiece2==0],{x,-2,3.2},{y,-2,3.2},PlotStyle->{Thickness[.012],GrayLevel[.35],Dashing[{.12,.08}]},PlotPoints->1000];

It appears that these together give most of the boundary. We now hunt for the part in the lower middle. For this we try to find a point near the lower crossing of the boundary curve with the y axis. It turns out to be expedient to make two changes to the minimization. One is that we now fix the value of x rather than y, and the other is that we will square the polynomial constraints (this seems to have the effect of better balancing them for the problem at hand).

polyx1=polys/.x->0;
Timing[{val,root1}=NMinimize[{1,Flatten[{Thread[polyx1^2==0],-1.6≤ y≤-1.2}]},{{y,-1.6,-1.2},{λ,-10,10},{(Subscript[c, a]),-1,1},{(Subscript[c, b]),-1,1},{(Subscript[s, a]),-1,1},{(Subscript[s, b]),-1,1}},MaxIterations->200,Method->"DifferentialEvolution"]]

This time it is not obvious from the trigonometric parameter values how a domain segment or curve might give us this part of the boundary, so now we’ll have to do some real work. We proceed to find more boundary points in order to reconstruct the curve therefrom. We begin with the solution point we just found. We will want to make sure we have a good approximation. As we used a global optimizer we may suspect that the root values are accurate to but a few decimal places (though the exact values for y and the Lagrange multiplier are fairly obvious). But we now have good initial values from which to proceed further, we can refine the result to arbitrarily high precision using a local root finder.

allbutxvars={y,λ,(Subscript[c, a]),(Subscript[c, b]),(Subscript[s, a]),(Subscript[s, b])};
rootbetter1=FindRoot[polyx1==0,Evaluate[Apply[Sequence,Transpose[{allbutxvars,allbutxvars/.root1}]]],PrecisionGoal->30,AccuracyGoal->30,WorkingPrecision->40]

It turns out that by comparison with the values of root1 we actually had a good result even then, correct to about 8 decimal places. This is most likely due to use of a local polishing step in NMinimize.

To deduce the actual curve we will require a number of {x,y} pairs of points lying on it. We can obtain these in various ways. One is simply to invoke NMinimize as above but with different values for y. This can be tedious, especially as the minimization can present problems in some regions (possibly due to proximity of multiple solutions, or difficulty balancing the various constraint equations). Instead we will treat this as a homotopy continuation problem wherein the roots found above will be initial values for a differential system. As it is easier to see y as a function of x in this part of the curve, we will use the latter as independent variable. We create the system below but display it in abbreviated form due to excessive size.

varsinx=Through[allbutxvars[x]];
parampolys=polys/.Thread[allbutxvars->varsinx];
diffpolys=D[parampolys,x];
Short[odesystem1=Join[Thread[diffpolys==0],Thread[(varsinx/.x->0)==(allbutxvars/.rootbetter1)]],12]

desoln1=NDSolve[odesystem1,varsinx,{x,0,1.2}];

The message tells us we were not able to get very far. A reasonable guess as to the cause is that we may be near a tangential intersection of solution curves (which would amount to a bifurcation in the solution set to what is really a differential algebraic system). In any case, what we have already suffices to give us points that lie approximately on the boundary curve. We will use several of them as initial values to find boundary points to high precision.

solnfunctions1=First[varsinx/.desoln1];
valuelists1=Table[Join[{t->x},Thread[List[allbutxvars,solnfunctions1]]], {x,0,3/10,5/1000}]/.t->x;
highprecsolns1=Table[Flatten[{First[valuelists1[[j]]],FindRoot[Evaluate[(polys/.First[valuelists1[[j]]])==0],Evaluate[Apply[Sequence,Rest[valuelists1[[j]]]]],PrecisionGoal->100,AccuracyGoal->100,WorkingPrecision->150]}],{j,Length[valuelists1]}];

We now find more solutions by moving to another part of (what we suspect will be) the same curve. We repeat the process of finding one solution, refining it, approximating a part of the curve through it, and refining individual points thereon.

polyx2=polys/.x->-2/5;
Timing[{val,root2}=NMinimize[{1,Flatten[{Thread[polyx2^2==0],-1.2≤ y≤-.9}]},{{y,-1.2,-.9},{λ,-10,10},{(Subscript[c, a]),-1,1},{(Subscript[c, b]),-1,1},{(Subscript[s, a]),-1,1},{(Subscript[s, b]),-1,1}},MaxIterations->200,Method->"DifferentialEvolution"]]

rootbetter2=FindRoot[polyx2==0,Evaluate[Apply[Sequence,Transpose[{allbutxvars,allbutxvars/.root2}]]],PrecisionGoal->30,AccuracyGoal->30,WorkingPrecision->40];
odesystem2=Join[Thread[diffpolys==0],Thread[(varsinx/.x->-2/5)==(allbutxvars/.rootbetter2)]];
desoln2=NDSolve[odesystem2,varsinx,{x,-2/5,-1}];

Again we did not get as far as we might like, but no matter. We will take several estimated solutions and refine them as before.

solnfunctions2=First[varsinx/.desoln2];
valuelists2=Table[Join[{t->x},Thread[List[allbutxvars,solnfunctions2]]], {x,-2/5,-4/5,-5/1000}]/.t->x;
highprecsolns2=Table[Flatten[{First[valuelists2[[j]]],FindRoot[Evaluate[(polys/.First[valuelists2[[j]]])==0],Evaluate[Apply[Sequence,Rest[valuelists2[[j]]]]],PrecisionGoal->100,AccuracyGoal->100,WorkingPrecision->150]}],{j,Length[valuelists2]}];

With these values we can now attempt to reconstruct the curve. As reasonable inferences from the picture we guess that it is a polynomial of degree at most 2 in y, and 12 in x. We form the possible monomials as our “basis” set, evaluate these at the {x,y} values corresponding to our high precision approximations of points on the boundary, and find the appropriate null space. Given that we have vectors of approximate numbers, the reliability of such an undertaking depends heavily on the method. Specifically it is based on singular value decomposition rather than, say, Gaussian elimination (see e.g. [Demmel 1997]).

xyvals=Join[highprecsolns1,highprecsolns2]/.((a_/;!MatchQ[a,x|y])->b_):>Sequence[];
xypowers=Flatten[Table[x^j*y^k,{j,0,12},{k,0,2}]];
valuevectors=xypowers/.xyvals;
ns=NullSpace[valuevectors];
Length[ns]

We find that we have 16 null vectors, so we cannot formulate a unique (to scalar multiples) polynomial. But this merely indicates that we have more monomials in our basis than we require; that is, the defining polynomial is (not surprisingly) using far fewer than the full set we provided. We “squeeze” out a null vector in terms of lower degree monomial basis members by extracting the last element of the null set in row echelon form. Before we proceed we will do two things. One is to replace the null vector values that are very clearly zeros to high accuracy with exact zeros (so as not to mistakenly use them as pivots). The other is to reverse the vectors so that entries that correspond to higher degree basis members are at the beginning. Strictly speaking it would probably be better to order by total degree, but this will suffice for our purposes. After we row reduce we may again have artificial values which we replace by zero.

ns2=Chop[ns,10^(-145)];
ns3=Map[Reverse,ns2];
rred=RowReduce[ns3];
rred2=Chop[rred,10^(-145)];
N[Last[rred2]]

This looks quite promising. We rationalize, find the common denominator, remove it by multiplication, and finally form the curve by associating the resulting vector of coefficients with the monomial basis.

coeffs=Rationalize[Last[rred2]];
mult=Apply[LCM,Denominator[coeffs]];
coeffs2=coeffs*mult;
envelopepiece3=coeffs2.Reverse[xypowers]

envplot3=ImplicitPlot[Evaluate[envelopepiece3==0],{x,-2,3.2},{y,-2,3.2},PlotStyle->{Thickness[.008],GrayLevel[.2],Dashing[{.12,.08}]},PlotPoints->1000];

This certainly looks like that missing lower middle portion of the boundary. We’ll put them all in one picture to see that this is so.

ll=Show[{lplot,envplot1,envplot2,envplot3}];

We now want to verify that our curve is indeed part of the boundary of the planar map. To this end we pick specific exact values on the curve and see if they are consistent with the original polynomial system. This will happen if and only if they lie on a curve that satisfies the original system. For example, we solve for y when x=-1/2.

xycoord=Prepend[First[Solve[(envelopepiece3/.x->-1/2)==0, y]],x->-1/2]

We plug these values into the polynomial system and check whether we get a nontrivial Gröbner basis.

GroebnerBasis[polys/.xycoord,{λ,(Subscript[c, a]),(Subscript[c, b]),(Subscript[s, a]),(Subscript[s, b])},MonomialOrder->DegreeReverseLexicographic]

This indicates that the point is in fact on a boundary curve. We try another such point. For brevity we will simply ascertain that the basis is not the trivial one (containing just 1) that would arise from an inconsistent system.

xycoord2=Prepend[First[Solve[(envelopepiece3/.x->-1/4)==0, y]],x->-1/4]
GroebnerBasis[polys/.xycoord2,{λ,(Subscript[c, a]),(Subscript[c, b]),(Subscript[s, a]),(Subscript[s, b])},MonomialOrder->DegreeReverseLexicographic]=!={1}

By way of contrast, we now demonstrate that a point not on the boundary will, as claimed above, give a basis containing just 1. For this we perturb the y coordinate from the last example. (That an inconsistent polynomial system has a reduced Gröbner basis of {1} is well known from general theory, but all the same should be demonstrated from time to time in actual computation.)

GroebnerBasis[polys/.{x->-1/4,y->-261/265},{λ,(Subscript[c, a]),(Subscript[c, b]),(Subscript[s, a]),(Subscript[s, b])},MonomialOrder->DegreeReverseLexicographic]

By verifying sufficiently many points--the number being dependent on the degree of the boundary curve-- one can prove that the entire curve satisfies the boundary curve polynomial system. This might be done using interpolation methods as in [von zur Gathen and Gerhard 1999, chapter 5].

For general interest we indicate another way to recover the boundary curve, this time from one point. We will start with the exact value but it will be seen that all we actually require is a high precision approximation. For this to work it is advisable that the coordinates not be algebraic numbers (so as to avoid any algebraic dependencies not forced by the curve polynomial itself). We will take the point with x=(π/4).

xycoord=Prepend[First[Solve[(envelopepiece3/.x->-Pi/4)==0, y]],x->-Pi/4]
valuevector=xypowers/.xycoord;

This time we form a lattice with these evaluated basis monomials multiplied by a suitably large constant in the first column and an identity matrix adjoined to the right.

lat=Transpose[Prepend[IdentityMatrix[Length[valuevector]],Round[10^200*N[valuevector,210]]]];

We now reduce the lattice and take as approximate null vectors all those for which the first column is fairly small. Again we reverse the ordering so as to have later components correspond to lexicographically smaller monomials, row reduce the matrix of null vectors, and use the last vector as coefficients for our implicit polynomial.

redlat=LatticeReduce[lat];
nullvecs=Map[Drop[#,1]&,Select[redlat,Abs[#[[1]]]<20&]];
Last[RowReduce[Map[Reverse,nullvecs]]]

We see that we did in fact recover the same polynomial coefficients as with the previous approach. We remark that this method is quite similar to that of finding a minimal polynomial for a given approximate number (which, for example, is the method behind the Mathematica function NumberTheory`Recognize). This in turn may be utilized to factor a univariate polynomial as in [Lenstra 1984]. One might also use for these purposes the efficient integer relation finding method known as PSLQ [Ferguson, Bailey, and Arno 1999].

5. Discussion and Directions for Future Work

We began with a brief review of a problem that arose in a number theory setting and in some way migrated to computational geometry as an example of a sort of extremal planar packing. We showed how symbolic computation methods may be applied to determine important values in that problem. We next discussed a planar trigonometric polynomial parametrization associated in a natural way with the original problem. It gave rise to a region for which the boundary in a sense generalizes a Lissajous figure. The next step in our enquiry was to utilize symbolic algebra to determine the equation for each part of that boundary. As we did not have a direct parametrization of that boundary this task could not be done by the usual implicitization methods (which tend to be tractable to at least moderate degree).

Thus far this was a warm up to the real problem at hand. When we look at examples of such maps in higher degrees the symbolic methods become mired by computational complexity. While certainly it might be possible that a clever reformulation, or a more powerful symbolic elimination implementation, might suffice, the methods we next developed are in a sense less troubled by degree or by specifics of the problem. Hence they might be utilized when one reaches the point at which clever symbolic approaches simply fail to yield a result.

For the harder problem we employed a mix of heuristics and numeric methods. We found approximations of isolated boundary curve points using both a numeric algebraic solver and optimization software (the latter being more readily applied in harder cases). In certain cases luck was on our side and we were able to “guess” the domain segment that gave rise to the part of the boundary on which our point was found. When this was not possible we showed how to refine the approximations and generate more approximate boundary points by setting up and solving what amounts to a DAE system; we worked with it instead as an ODE system in order to maintain simplicity in the exposition and code. We then used simple numeric linear algebra to fit to a guessed curve. For this we used far more points than we had basis vectors, so as to be confident of generating the correct null space to good approximation (this seems to be all the more important given how close our sample points were clustered). We next employed row reduction to find a “minimal” null vector corresponding to the “smallest” basis monomials. This gave our solution curve polynomial.

We remark that the steps of solving an ODE (or DAE) system to obtain points on the graph, and computing a null space for a monomial basis evaluated at those points, are also utilized in [Corless, Giesbrecht, van Hoeij, Kotsireas, and Watt 2001] for the purpose of approximate factorization of bivariate polynomials. This in turn is related to work in [Corless, Giesbrecht, Kotsireas, and Watt 2001] for approximate implicitization of hypersurfaces. Our usage is actually more closely related to the former reference, in that the latter uses a calculus of variations formulation involving an integral functional, which first requires a parametrization of the curve. A big difference between our problem and that of bivariate factorization is that our hypersurface (curve, in this case) arises from an extremality condition on a surface map and hence there are added complications both to finding points on and following paths along it. One also recognizes at this point that our method might find a factor over the (approximate) reals that is a proper factor of a factor over the rationals. We will discuss this possibility in the last section.

An issue that did not arise in our particular example is as follows. The method of finding an approximate factor of the boundary curve might yield a factor over an extension field of the rationals. In this case we are faced with a few choices, and what we do will depend on what we require. If our interest is in the computational geometry aspects of the problem we might be quite happy to work with that approximate factorization over the reals. If we are working on a problem from mathematics that requires the exact factor (as was the case in [Kim 2003]) then we can attempt to find the exact form by methods that recognize algebraic numbers from numeric approximations of sufficiently high precision as discussed in [Lenstra 1984], [Ferguson, Bailey, and Arno 1999]. Alternatively we might try to find approximations of the algebraic conjugate factors, e.g. using methods from [Corless, Galligo, Kotsireas, and Watt 2002], and combine the factors to deduce the correct factor over the rationals. Finally we should note that if only a crude approximation to the curve is needed, then fast methods of edge detection from the field of image processing may suffice.

One thing to notice is that in the type of problem we discussed this issue cannot arise unless a similar factorization takes place in the domain. To be specific, recall that a necessary condition that a point be on the boundary curve is that it be the image of a point in the domain at which the Jacobian vanishes. Thus distinct factors in the boundary must come from distinct factors of the Jacobian. Since we have the Jacobian in explicit form, methods of absolute factorization might be employed to see whether there are factors over a nontrivial algebraic extension of the rationals.

It is important to realize that hybrid methods appear in several ways in this work, some lurking behind the scenes. For example, the numeric solver NSolve is based on a mix of Gröbner basis and numeric methods [Lichtblau 2000]. The function Rationalize is based on continued fraction methods (see e.g. chapter 4 of [von zur Gathen and Gerhard 1999]), which themselves can use mixed exact and approximate numeric methods when operating at high precision (though the example we showed was quite tame). Various implementations of lattice reduction make use of approximate arithmetic. Though we did not do so, a lattice reduction method may be employed to rationalize a set of approximate values collectively so as to obtain a common denominator. This is generally referred to as simultaneous diophantine approximation, and the underlying method is quite similar to that of recognizing a minimal polynomial from a root approximation; the principle difference in implementation involves whether or not a lattice is transposed. See [von zur Gathen and Gerhard 1999, chapter 17].

One obvious question to pose is just how generally useful are the methods we have discussed. To a large extent this depends on the sort of computational problem one wishes to solve. For purposes of obtaining a numeric approximation e.g. for graphing, what we have shown clearly requires too much work. It is only when one requires either an exact implicit curve or a high precision approximation that something of this sort must be done. This can arise, for example, in finding algebraic surface level curves (in effect we did that, taking as our third surface coordinate the Jacobian, and looking at the set where it vanishes). Another such computation might be finding offset curves from an algebraic planar object in cases where the complexity of a purely symbolic approach is too great. For this one might start with approximations based on fast marching methods as described in [Sethian 1999] (a simple Mathematica implementation of which may be found in [Lichtblau 2003a]).

We conclude with several questions and remarks that indicate directions for further investigation.

• Could methods from real algebraic geometry e.g. cylindrical algebraic decomposition be employed to find boundary curves efficiently?

• Can the process we described be automated? In our work it required a lot of eyeballing of the list plot to see what might be viable places to look for approximate curve points. Indeed, there were several failed attempts (not shown).

• Are there more robust ways to find individual points on the curves? While we used some global optimization software, other methods e.g. with intervals might work better. Image processing methods can be used to find approximate coordinate values for {x,y} but one also requires corresponding values from the domain space in order to obtain high precision refinements. A possible approach would be to regard {x,y} as a Taylor series in {a,b} and then “invert”. Of course on the curve of interest the inverse function theorem does not apply, because we are precisely where the Jacobian vanishes! But this is not really a concern because a “multivalued” inverse function approximation e.g. from a Puiseux series will suffice for our purposes. Alternatively we might lift the problem to three dimensions, say by adding a coordinate function z=Jacobian(x(a,b),y(a,b)), and change coordinates to invert locally the map from plane to surface.

• In the symbolic computations of boundary curves we found the extremal formulation in terms of Lagrange multipliers to be preferable to the vanishing Jacobian approach (possibly because the latter gives factors that are extraneous for our purposes). A drawback to the Lagrange multiplier method in our later numeric approach is that we really do not have a good idea of the range of values the multiplier variable might take. So it might be advantageous to revisit the vanishing Jacobian formulation in conjunction with the numeric techniques.

• Are there better ways to “track” a piece of curve to obtain a set of approximate points thereon?

• We first found candidate boundaries. We outlined a method of validating them algorithmically, by showing that individual points satisfy the boundary curve system. We chose this approach because the more direct route of computing an elimination ideal appears to be infeasible using the software available. It would be interesting to know whether there might be a feasible direct approach.

• Our method of deducing the implicit polynomial from rationalizing relied on those curves being defined by polynomials over the rationals. Is there a way to extend to curves with minimal polynomials having algebraic numbers for coefficients?

• Might this method generalize in a reasonable way to implicitization of bounding surfaces of parametrized volumetric maps?

6. Acknowledgements

I thank Ken Stolarsky (thesis advisor to Seon-Hong Kim) for prompting me to consider the many computational aspects of this family of problems. I thank Josef Schicho for reminding me of the sometimes more efficient approach to similar problems that was used in [Hong and Schicho 1998]. I thank the organizers of the ACA 2003 session on Symbolic-Numeric methods for Curves and Surfaces, John Wetzel (emeritus of the UIUC math department and organizer of the Geometric Potpourri seminar series), and the organizers of the ADG 2004 conference, all of whom provided opportunities for me to speak on various stages of this work.

7. References

B. Buchberger. Gröbner bases: An algorithmic method in polynomial ideal theory. In Multidimensional Systems Theory, chapter 6. N. K. Bose, ed. D. Reidel Publishing Company, 1985.

B. Champion. Numerical Optimization in Mathematica: An Insider’s View of NMinimize. In SCI2002, Proceedings of the 6th World Multiconference on Systemics, Cybernetics, and Informatics. Volume 16, N. Callaos, T. Ebisuzaki, B. Starr, J. M. Abe, D. Lichtblau, eds., pages 136-140. International Institute of Informatics and Systemics, 2002.
A Mathematica notebook version may be found at:
http://library.wolfram.com/infocenter/Conferences/4311/
Up to date documentation regarding NMinimize may be found at:
http://documents.wolfram.com/v5/Built-inFunctions/AdvancedDocumentation/Optimization/NMinimize/

R. M. Corless, A. Galligo, I. S. Kotsireas, S. M. Watt. A Geometric-numeric algorithm for absolute factorization of multivariate polynomials. In: Proceedings of the 2002 International Symposium on Symbolic and Algebraic Computation (ISSAC 2002), T. Mora, ed., pages 37-45. ACM Press, 2002.

R. M. Corless, M. W. Giesbrect, M. van Hoeij, I. S. Kotsireas, S. M. Watt. Towards factoring bivariate approximate polynomials. In: Proceedings of the 2001 International Symposium on Symbolic and Algebraic Computation (ISSAC 2001), B. Mourrain, ed., pages 85-92. ACM Press, 2001.

R. M. Corless, M. W. Giesbrect, I. S. Kotsireas, S. M. Watt. Numerical implicitization of parametric hypersurfaces with linear algebra. Artificial Intelligence and Symbolic Computation (AISC 2000), J.A. Campbell and E. Roanes-Lozano, eds. Lecture Notes in Artificial Intelligence 1930, pages 174-183. Springer-Verlag, 2001.

D. Cox, J. Little, D. O’Shea. Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Computer Algebra. Undergraduate Texts in Mathematics. Springer-Verlag, 1992.

J. Demmel. Applied Numerical Linear Algebra. SIAM, 1992.

H. Ferguson, D. Bailey, S. Arno. Analysis of PSLQ, an integer relation finding algorithm. Mathematics of Computation 68(225):351-369, 1999.

J. von zur Gathen, J. Gerhard. Modern Computer Algebra. Cambridge University Press, 1999.

H. Hong, J. Schicho. Algorithms for trigonometric curves (simplification, implicitization, and parameterization). Journal of Symbolic Computation 26:279-300, 1998.

M. Kalkbrenner. Implicitization of rational parametric curves and surfaces. In AAECC8, Applied Algebra, Algebraic Algorithms, and Error-Correcting Codes, S. Sakata, ed. Springer Lecture Notes in Computer Science 508, pages 249-259. Springer-Verlag, 1990.

S.-H. Kim. Packing of ovals and an extremal problem for exponential sums. Talk presented in the Geometric Potpourri Seminar series, Department of Mathematics, University of Illinois at Urbana-Champaign, September 15, 1998.
Abstract may be found at:
http://www.math.uiuc.edu/Bulletin/Abstracts/September/sep11-98geompot.html

S.-H. Kim. Sums of Polynomials, Minmax Problems and Number Theory. Ph.D thesis, Department of Mathematics, University of Illinois at Urbana-Champaign, 2000.

S.-H. Kim. Planar packings and mappings related to certain minmax problems. Mathematical Inequalities and Applications 6(2):351-374, 2003.

A. K. Lenstra. Polynomial factorization by root approximation. In EUROSAM 84, Proceedings of the International Symposium on Symbolic and Algebraic Computation, J. Fitch ed. Lecture Notes in Computer Science 174, pages 272-276, Springer-Verlag, 1984.

D. Lichtblau. Solving finite algebraic systems using numeric Gröbner bases and eigenvalues. In SCI2000, Proceedings of the World Conference on Systemics, Cybernetics, and Informatics. Volume 10, (Concepts and Applications of Systemics, Cybernetics, and Informatics). M. Torres, J. Molero, Y. Kurihara, and A. David, eds., pages 555-560. International Institute of Informatics and Systemics, 2000.

D. Lichtblau. Ordered heaps and fast marching method. 2003. Mathematica notebook available at:
http://library.wolfram.com/infocenter/Demos/4928/

D. Lichtblau. Talk presented in the session on Symbolic-Numeric methods for Curves and Surfaces at the International Conference on Applications of Computer Algebra (ACA 2003), Raleigh, North Carolina, 2003.
Abstract may be found at:
http://www.risc.uni-linz.ac.at/people/shalaby/ACA03/ 2003.

D. Lichtblau. Cylinders through five points in (ℝ^3). Submitted, 2004.

J. Sethian. Level Set Methods and Fast Marching Methods: Evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science (2nd edition). Cambridge University Press, 1999.

S. Wolfram. The Mathematica Book (5th edition). Wolfram Media, 2003.