How to use an iterative method to compute y values corresponding to an x value of a rotated ellipse?

Question

I am trying to render the outline of a rotated ellipse on raster graphics.

The general form of a conic is:

H(x, y) = A x² + B xy + C y² + D x + E y + F = 0

I am currently able to render the ellipse by inputting incremental x values (starting from left most) and then solving a the resultant quadratic equation in y.

I came across this post while trying to optimise my code and found this:

Lastly, note that incremental computation saves work when evaluating H, on the line of

H(x+1, y) = H(x, y) + A.(2x+1) + B.y + D = H(x, y) + (2A).x + (A + B.y + D)

It seems like a way better idea to use this incremental computation method than calculating the roots every time I want to draw a point.

I would like to understand how I can use this method, given I have the starting (x,y) point and the ending (xmax,ymax) point (right-most point of ellipse).

Answer 1

Calling

$$ F(x,y) = a x^2+b x y + c y^2+ d x + e y + f = 0 $$

given $x = x^*$ we can calculate the corresponding $y^*$ using an iterative process like Newton's

$$ y_{k+1} = y_k - \frac{F(x^*,y_k)}{F_y(x^*,y_k)} $$

Here

$$ F_y(x,y) = b x + 2c y + e $$

NOTE

The $y_0$ initial guess should be given near the potential solution.

Attached the MATHEMATICA script which produces the previous graphics.

Clear[F, dF]
parms = {a -> 1/10, b -> 2/100, c -> 1/20, d -> 0, e -> 0, f -> -10}    
F[x_, y_] := a x^2 + b x y + c y^2 + d x + e y + f
dF[x_, y_] := b x + 2 c y + e
error = 0.001;

For[j = 1; PATH = {}, j <= 50, j++,
   x0 = RandomReal[{-10, 10}];
   y0 = RandomReal[{-12, 12}];
   path = {{x0, y0}};
   For[k = 1, k <= 10, k++, 
      y1 = (y0 - F[x0, y0]/dF[x0, y0])/. parms;
      If[Abs[y0 - y1] < error, Break[]];
      y0 = y1; 
      AppendTo[path, {x0, y0}]
   ];
   AppendTo[PATH, {path}];
]

gr1 = ContourPlot[(F[x, y] /. parms) == 0, {x, -20, 20}, {y, -20, 20}];
grpath = Table[Graphics[{Green, Line[PATH[[m]][[1]]]}], {m, 1,Length[PATH]}];
grpt = Table[Graphics[{Red, Point[PATH[[m]][[1]]]}], {m, 1, Length[PATH]}];
Show[gr1, grpath, grpt]

Answer 2

The vector equation of an ellipse with major/minor axes $a,b$ rotated by an angle $\theta$ is given by $$\begin{bmatrix}x\\y\end{bmatrix} = \begin{bmatrix}c&-s\\s&c\end{bmatrix}\begin{bmatrix}a&0\\0&b\end{bmatrix}\begin{bmatrix}\cos t\\\sin t\end{bmatrix}=RD\begin{bmatrix}\cos t\\\sin t\end{bmatrix}$$ where $c=\cos\theta$, $s=\sin\theta$, $R=\begin{bmatrix}c&-s\\s&c\end{bmatrix}$, $D=\begin{bmatrix}a&0\\0&b\end{bmatrix}$.

Then, differentiating with respect to $t$ and inverting back to the $x,y$ variables gives $$\begin{bmatrix}\dot{x}\\ \dot{y}\end{bmatrix}=RD\begin{bmatrix}-\sin t\\\cos t\end{bmatrix} = RDTD^{-1}R^{-1}\begin{bmatrix}x\\y\end{bmatrix}$$ where $T=\begin{bmatrix}0&-1\\1&0\end{bmatrix}$.

The matrix $$M=RDTD^{-1}R^{-1}=\begin{bmatrix}(\frac{a}{b}-\frac{b}{a})sc&-\frac{b}{a}s^2-\frac{a}{b}c^2\\\frac{b}{a}c^2+\frac{a}{b}s^2&(\frac{b}{a}-\frac{a}{b})sc\end{bmatrix}$$ can be calculated once before the loop.

Then $$\begin{bmatrix}\Delta x\\\Delta y\end{bmatrix}=M\begin{bmatrix}x\\y\end{bmatrix}\Delta t$$ You can choose $\Delta t$ as small as fits your resolution, for example, $0.01$. Then iterate x+=$\Delta$x, y+=$\Delta$y for $2\pi/\Delta t$ times, for example 628. The initial values of $(x,y)$ are given by $(ac,as)$.

If you want part of an ellipse, then you need to figure out the initial and final values of $t$. If you make $\Delta t$ too large or too small errors accumulate and the ellipse does not close on itself.