Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.soft-sys.math.mathematica > #2103

Re: envelope of curves

Show all headers | View raw

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

Back to comp.soft-sys.math.mathematica | Previous | Next | Find similar | Unroll thread

Thread

Re: envelope of curves Daniel Lichtblau <danl@wolfram.com> - 2011-05-04 23:47 +0000

csiph-web