3
$\begingroup$

I have ellipse function as below:

−3.06274 x^2 − y^2 + 1192.22 x + 152.71 y + 1.829648 x y − 196494 == 0  

This ellipse looks like: ellipse

How to get the rotate angle of this rotated ellipse?

| improve this question | | | | |
$\endgroup$
  • $\begingroup$ Have a look at this. You only have to find $\Sigma$ and $p$ by completing the square. $\endgroup$ – Henrik Schumacher Jan 22 '18 at 7:33
  • $\begingroup$ I think there's a mathematica function can help me to transform the function to complete the square.but I cannot find it. $\endgroup$ – kittygirl Jan 22 '18 at 8:09
  • 1
    $\begingroup$ You could get the EigenVectors of your ellipse and then use the dot product... $\endgroup$ – José Antonio Díaz Navas Jan 22 '18 at 9:03
1
$\begingroup$

Nasser already gave you a solution. There is a different approach though where I used nothing more than the basic knowledge, that an ellipse can be represented by

$$\frac{x^2}{a^2}+\frac{y^2}{b^2}-val=0$$

and you said in the comments, that you are seeking a quadratic form. Your example shows that the ellipse is shifted and rotated. Let us transform the x and y in the above equation and introduce a shift and a rotation. We use mx and my as translation and phi as rotation and this can be represented by a matrix multiplication

RotationMatrix[phi].{x - mx, y - my}

Now we replace our original x and y with this result

expr = x^2/a^2 + y^2/b^2 - val /. Thread[{x, y} -> RotationMatrix[phi].{x - mx, y - my}]

We get

$$\frac{((x-\text{mx}) \cos (\phi )-(y-\text{my}) \sin (\phi ))^2}{a^2}+\frac{((x-\text{mx}) \sin (\phi )+(y-\text{my}) \cos (\phi ))^2}{b^2}-\text{val}$$

You can use

expr // ExpandAll

to see all terms. What you should notice is that we have positive terms for y^2 in the form

(y^2 Sin[phi]^2)/a^2

because the sine-squared will always be positive and a^2 too. In your example those terms have a negative coefficient. Since it is an equation, we can simply multiply with -1 to turn everything around

example = (-3.06274*x^2 - y^2 + 1192.22*x + 152.71*y + 1.829648*x*y - 196494)*-1

If you look at the expanded form of expr and the form of your example, you can see that we can now compare coefficients:

eq = Thread[Flatten /@ (CoefficientList[example, {x, y}] == 
      CoefficientList[expr, {x, y}])];
eq // Column

Mathematica graphics

When I counted correctly, we have 6 valid equations for 6 parameters mx, my, a, b, phi, and val. Let us try to solve this

sol = Solve[eq, {mx, my, a, b, phi, val}];
sol // Column

Mathematica graphics

Here, you have the full set of solutions (if I'm not mistaken). You see that they mainly differ in a, b, and phi. Thinking about it, this can be easily explained. First, when there are only a^2 and b^2 in the equation, both a and -a are valid solutions. Second, you can have ellipses with larger a or larger b, they only need to be rotated differently. This mix gives the different solutions that are all valid.

Let's test the first solution by replacing the values in our original expr

expr /. sol[[1]] // Expand
(* 196494. - 1192.22 x + 3.06274 x^2 - 152.71 y - 1.82965 x y + 1. y^2 *)

As you can see, we get your original example. The last solution has a phi value of 0.362 which is in degree

phi/(2 Pi)*360 /. sol[[-1]]
(* 20.7865 *)

the same value that Nasser calculated. Finally, if you want to see the quadratic form, you can use this:

With[{e = expr},
 HoldForm[e] /. sol[[-1]]
]

$$\frac{((-350.067+y) \cos (0.362793)+(-299.196+x) \sin (0.362793))^2}{1.23774^2}+\frac{((-299.196+x) \cos (0.362793)-(-350.067+y) \sin (0.362793))^2}{0.54153^2}-8589.02$$

| improve this answer | | | | |
$\endgroup$
3
$\begingroup$

Perhaps this (with J.M.'s shortening of the code):

quad = Last@ CoefficientArrays[
    -3.06274 x^2 - y^2 + 1192.22 x + 152.71 y + 1.829648 x y - 196494,
    {x, y}, "Symmetric" -> True];
ev = Eigenvectors@quad;
Mod[ArcTan @@@ ev, Pi]

(*  {2.7788, 1.208}  *)
| improve this answer | | | | |
$\endgroup$
  • 2
    $\begingroup$ You can do Last[CoefficientArrays[-3.06274 x^2 - y^2 + 1192.22 x + 152.71 y + 1.829648 x y - 196494, {x, y}, "Symmetric" -> True]] instead to directly get a symmetric matrix. $\endgroup$ – J. M.'s technical difficulties Mar 4 '18 at 15:02
  • $\begingroup$ @J.M. Thanks. I hadn't noticed that option. Best wishes $\endgroup$ – Michael E2 Mar 4 '18 at 15:04
2
$\begingroup$

The angle is given by 1/2 ArcCot[(a - c)/(2 b)] when writing the general conic section as $a x^2+ 2 b x y + c y^2+ 2 d x+2 f y+g =0$ Hence in this case

3.06274 x^2 - y^2 + 1192.22 x + 152.71 y + 1.829648 x y - 196494=0;

The angle is -20.78 degree.

Use http://mathworld.wolfram.com/Ellipse.html

Mathematica graphics

Since in this case $c=-1$ and $a=-3.06274$. To verify, here is a manipulate, which plots the original -3.06274*x^2 - y^2 + 1192.22*x + 152.71*y + 1.829648*x*y - 196494 == 0 as ContourPlot then plots the standard ellipse equation when rotated, which is

$$ \frac{ ((x-x_0) \cos \alpha + (y-y_0) \sin\alpha)^2}{a^2} + \frac{ ((x-x_0) \sin \alpha - (y-y_0) \cos\alpha)^2}{b^2}=1 $$

Using the angle of rotation $\alpha$ found from 1/2 ArcCot[(a - c)/(2 b)] and using a and b found per

Mathematica graphics

From http://mathworld.wolfram.com/Ellipse.html

The center $x_0$ and $y_0$ are given by

Mathematica graphics

This shows exact fit. The red color is the ContourPlot using your general conic and the black ellipse is using the rotation angle.

Mathematica graphics

g1 = ContourPlot[-3.06274*x^2 - y^2 + 1192.22*x + 152.71*y + 
     1.829648*x*y - 196494 == 0, {x, 200, 500}, {y, 200, 500}, 
   ContourStyle -> Red, Frame -> False]; 

eq = -3.06274*x^2 - y^2 + 1192.22*x + 152.71*y + 1.829648*x*y - 196494; 
(*reading from the above gives *)
a = -3.06274; 
b = 1.829648/2; 
c = -1; 
d = 1192.22/2; 
f = 152.71/2; 
g = -196494; 

(*see Wolfram Mathworld for derivation*)
a0 = Sqrt[(2*(a*f^2 + c*d^2 + g*b^2 - 2*b*d*f - a*c*g))/((b^2 - 
        a*c)*(Sqrt[(a - c)^2 + 4*b^2] - (a + c)))]; 
b0 = Sqrt[(2*(a*f^2 + c*d^2 + g*b^2 - 2*b*d*f - a*c*g))/((b^2 - 
        a*c)*(-Sqrt[(a - c)^2 + 4*b^2] - (a + c)))]; 
x0 = (c d - b f)/(b^2 - a c);
y0 = (a f - b d)/(b^2 - a c);

(*compare to given conic*)
Manipulate[Module[{z = ang*(180/Pi), g2, x, y}, 

  g2 = ContourPlot[((x - x0)*Cos[z] + (y - y0)*Sin[z])^2/
       a^2 + ((x - x0)*Sin[z] - (y - y0)*Cos[z])^2/b^2 == 1, {x, 100, 
     500}, {y, 100, 600}, ContourStyle -> Black]; 

  Graphics[{First[g2], First[g1]},PlotRange -> {{100, 600}, {100, 600}}, Axes -> True]], 

 {{a, 50.18738140246613, "a"}, 50, 200, 0.01,Appearance -> "Labeled"}, 
 {{b, 114.71032896005666, "b"}, 40, 130, 0.01,Appearance -> "Labeled"}, 
 {{ang, -20.786494338349772, "angle"}, -40, 90, 0.1,Appearance -> "Labeled"}
 ]
| improve this answer | | | | |
$\endgroup$
1
$\begingroup$

You can get all parameters of that ellipse in a quite mechanical way.

The normal ellipse equation is

eli[x_, y_, a_, b_] = x^2/a^2 + y^2/b^2 - 1 == 0

to rotate the ellipse, apply this rule

rotate[phi_] := 
   Thread[{x,y} -> {{Cos[phi], -Sin[phi]}, {Sin[phi], Cos[phi]}}.{x, y}]

To move the center apply

move[c_, d_] := Thread[{x, y} -> {x - c, y - d}]

ContourPlot[
   Evaluate[eli[x, y, 50, 120] /. rotate[20 Degree] /. 
   move[300, 350]], {x, 200, 400}, {y, 200, 500}, MaxRecursion -> 8]

enter image description here

In order to find all parameters, first find y to given x in your ellipse

sol = First@
       Solve[\[Minus]3.06274 x^2 \[Minus] y^2 + 1192.22 x + 152.71 y + 
         1.829648 x y \[Minus] 196494 == 0, y, Reals] // Quiet

Generate five equations with five given x-values {250, 260, 270, 280, 290}

equ = (((First[eli[x, y, a, b]] /. rotate[phi] /. 
      move[c, d]) == -196494 + 1192.22` x - 3.06274` x^2 + 
     152.71` y + 1.829648` x y - y^2) /. sol /. x -> #) & /@ {250,
260, 270, 280, 290}

Use FindRoot to find all parameters.

FindRoot[equ, {{a, 50}, {b, 120}, {c, 300}, {d, 350}, {phi, 20}}]

(* {a -> 50.1874, b -> 114.71, c -> 299.196, d -> 350.067, phi -> 20.7865} *)
| improve this answer | | | | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.