A formula for the rotated ellipse can be found by taking the standard ellipse $$ \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1 \qquad\qquad\qquad\qquad {\bf{x}}(t) = (a\cos t, b\sin t) \qquad \text{ for } t\in [0, 2\pi] $$ and rotating it counter-clockwise using a rotation matrix $$ R{\bf{x}}(t) = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} a\cos t \\ b\sin t \end{pmatrix} = \begin{pmatrix} a\cos t\cos\theta - b\sin t\sin \theta\\ a\cos t\sin \theta + b \sin t \cos \theta \end{pmatrix} = {\bf{x}}_r(t) $$ I can easily find the gradient of the standard ellipse (non-rotated) by taking the derivative of the first formula. $$ g = \begin{pmatrix} \frac{2x}{a^2}, & \frac{2y}{b^2} \end{pmatrix}^\top $$
Now I tried to methods to compute the gradient at a point on the rotated ellipse.
Method 1: I take the gradient of $$ \frac{(x\cos\theta - y\sin\theta)^2}{a^2} + \frac{(x\sin\theta + y\cos\theta)^2}{b^2} = 1 $$ obtaining $$ g_{r1} = \begin{pmatrix} \frac{2(x\cos\theta - y\sin\theta)\cos\theta}{a^2} + \frac{2(x\sin\theta + y\cos\theta)\sin\theta}{b^2} \\ \frac{-2(x\cos\theta - y\sin\theta)\sin\theta}{a^2} + \frac{2(x\sin\theta + y\cos\theta)\cos\theta}{b^2} \end{pmatrix} $$
Method 2: Rotating the gradient of the non-rotated ellipse counterclockwise $$ g_{r2} = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} \frac{2x}{a^2} \\ \frac{2y}{b^2} \end{pmatrix} = \begin{pmatrix} \frac{2x\cos\theta}{a^2} - \frac{2y\sin\theta}{b^2}\\ \frac{2x\sin\theta}{a^2} + \frac{2y\cos\theta}{b^2} \end{pmatrix} $$
To me, method 1 should work while method 2 should be wrong. However, as you can see in the plot below, $g_{r1}$ seems to actually be totally wrong while $g_{r2}$ is right. How come $g_{r2}$ seems to be right but not $g_{r1}$?
Method 1 is wrong because the equation $$ \frac{(x\cos\theta - y\sin\theta)^2}{a^2} + \frac{(x\sin\theta + y\cos\theta)^2}{b^2} = 1 $$
expresses the fact that the result of rotating the point $(x,y)$ belongs to the ellipse $x^2/a^2+y^2/b^2=1$, but this is not what you need.