Re: envelope of curves
- To: mathgroup at smc.vnet.net
- Subject: [mg118624] Re: envelope of curves
- From: Daniel Lichtblau <danl at wolfram.com>
- Date: Wed, 4 May 2011 19:47:57 -0400 (EDT)
Dan wrote:
> I would like to find the envelope of a set of curves. This evelope may
> be formed by moving the center of one curve along the path of the
> other curve. It doesn't matter which curve is held fixed and which is
> translated as can be seen in the two plots. I am pretty sure there is
> a fairly straightforward algorithm for this, but it is escaping me at
> the moment. One of the curves, drawn here in grey, is always an
> ellipse with the major and minor axes aligned with x and y. The other
> curve, drawn in red, may be like an ellipse but it may be some other
> shape. In the actual problem it is a fairly involved calculation to
> determine it, but it is always smooth and convex. If there is a closed
> form solution for the elliptical case, that would be useful. On the
> other hand, a numerical solution that can be expressed as an
> interpolated curve would probably be the most general solution for me.
>
> With[{a=5,b=1,c=0.25,d=1.5,theta=-Pi/6,n=50},
> Module[{rot,c1,c2},
> rot={{Cos[theta],-Sin[theta]},{Sin[theta],Cos[theta]}};
> c1[psi_]:={a Sin[psi],b Cos[psi]};
> c2[psi_]:=rot.{c Sin[psi],d Cos[psi]};
> Column[{
> Show[{
> ParametricPlot[c1[psi],{psi,-Pi,Pi},Axes->False,
> PlotStyle->Gray,PlotRange->All],
> ParametricPlot[c2[psi],{psi,-Pi,Pi},Axes->False,
> PlotStyle->Red],
> Table[
> ParametricPlot[c1[phi]+c2[psi],{psi,-Pi,Pi},
> Axes->False,
> PlotStyle->Directive[Red,Opacity[0.2]]
> ],{phi,-Pi,Pi,2Pi/n}
> ]
> },ImageSize->400],
> Show[{
> ParametricPlot[c1[psi],{psi,-Pi,Pi},Axes->False,
> PlotStyle->Gray,PlotRange->All],
> ParametricPlot[c2[psi],{psi,-Pi,Pi},Axes->False,
> PlotStyle->Red],
> Table[
> ParametricPlot[c1[psi]+c2[phi],{psi,-Pi,Pi},
> Axes->False,
> PlotStyle->Directive[Gray,Opacity[0.2]]
> ],{phi,-Pi,Pi,2Pi/n}
> ]
> },ImageSize->400]
> }]
> ]
> ]
>
> Thanks in advance.
> Dan
>
Not sure I understand what it is you want, but I'll try to explain what
is computed below. We take your family of parametric curves
c1[psi]+c2[phi], set up polynomials in {x,y} in terms of these, add trig
identities.
We then "extremize" via a Lagrange multiplier. This gives two new
polynomials and one new variable. We now "polynomialize" by making the
trigs into new variables (keep in mind we have identities to relate them
in pairs).
The next thing is to eliminate all variables other than {x,y}. In doing
so we obtain a polynomial. It is the implicit form of the extremal curve.
Here is the code I use for this.
grad[expr_, a_, b_] := {D[expr, a], D[expr, b]}
vals = {a -> 5, b -> 1, c -> 1/4, d -> 3/2, theta -> -Pi/6};
rot = {{Cos[theta], -Sin[theta]}, {Sin[theta], Cos[theta]}};
c1[psi_] := {a Sin[psi], b Cos[psi]};
c2[psi_] := rot.{c Sin[psi], d Cos[psi]};
coordpolys = {x, y} - (c1[phi] + c2[psi]);
reprules = {Cos[psi] -> cps, Sin[psi] -> sps, Cos[phi] -> cph,
Sin[phi] -> sph};
trigpolys = {Cos[psi]^2 + Sin[psi]^2 - 1, Cos[phi]^2 + Sin[phi]^2 - 1};
vars = {x, y};
gradpolys =
grad[coordpolys[[1]], phi, psi] -
lambda*grad[ coordpolys[[2]], phi, psi];
allpolys = Join[coords, trigpolys, gradpolys] /. vals /. reprules
elimvars = Complement[Variables[allpolys], vars];
(* Here is the polynomial system *)
Out[92]= {-((3*cps)/4) - 5*sph - (Sqrt[3]*sps)/8 +
x, -cph - (3*Sqrt[3]*cps)/4 +
sps/8 + y, -1 + cps^2 + sps^2, -1 + cph^2 + sph^2, -5*cph -
lambda*sph,
-((Sqrt[3]*cps)/8) + (3*sps)/4 - lambda*(cps/8 + (3*Sqrt[3]*sps)/4)}
(* Now compute the implicit polynomial. Can be done via GroebnerBasis,
eliminating all variables but {x,y}. *)
In[98]:= imp = First[GroebnerBasis[allpolys, vars, elimvars]]
Out[98]= -632169101025 - 15438252990*x^2 + 2327819519*x^4 +
50943632*x^6 -
3041536*x^8 + 178880770180*Sqrt[3]*x*y -
18989375300*Sqrt[3]*x^3*y +
365862560*Sqrt[3]*x^5*y + 3906560*Sqrt[3]*x^7*y +
6477954371014*y^2 -
580839892690*x^2*y^2 + 16449887888*x^4*y^2 - 158016512*x^6*y^2 +
328983054260*Sqrt[3]*x*y^3 - 17352995520*Sqrt[3]*x^3*y^3 +
196725760*Sqrt[3]*x^5*y^3 - 1341456141001*y^4 +
117801783152*x^2*y^4 -
2198334976*x^4*y^4 - 66752924000*Sqrt[3]*x*y^5 +
2511488000*Sqrt[3]*x^3*y^5 + 26694841200*y^6 - 3731788800*x^2*y^6 +
873600000*Sqrt[3]*x*y^7 - 243360000*y^8
This polynomial is irreducible (does not factor). We can plot it as below.
ct = ContourPlot[imp, {x, -8, 8}, {y, -4, 4}, Contours -> {0},
ContourShading -> None, PlotPoints -> 50]
It looks like an oblong curved rhombus, with a necktie inside (because
it is elegant, one would presume). The necktie is something of an
artifact, where the extremal equations are satisfied but all the same
other parts overlay it.
You can get an idea of how it all fits together as below.
n = 50;
plts = Table[
ParametricPlot[c1[phi] + c2[psi] /. vals, {psi, -Pi, Pi},
Axes -> False,
PlotStyle -> Directive[Red, Opacity[0.2]]], {phi, -Pi, Pi,
2 Pi/n}];
Show[{ct, plts}]
Can do similarly with the other family.
plts2 = Table[
ParametricPlot[c1[psi] + c2[phi] /. vals, {psi, -Pi, Pi},
Axes -> False,
PlotStyle -> Directive[Gray, Opacity[0.2]]], {phi, -Pi, Pi,
2 Pi/n}];
For more information about this method of computing an envelope, see
http://library.wolfram.com/infocenter/Conferences/7516/
Daniel Lichtblau
Wolfram Research