CODEWITHSUNDEEP
matlabgeometryellipse
John
John

Reputation: 11

Plotting an ellipse in MATLAB given in matrix form

I have an ellipse in 2 dimensions, defined by a positive definite matrix X as follows: a point x is in the ellipse if x'*X*x <= 1. How can I plot this ellipse in matlab? I've done a bit of searching while finding surprisingly little.

Figured out the answer actually: I'd post this as an answer, but it won't let me (new user):

Figured it out after a bit of tinkering. Basically, we express the points on the ellipse border (x'*X*x = 1) as a weighted combination of the eigenvectors of X, which makes some of the math to find the points easier. We can just write (au+bv)'X(au+bv)=1 and work out the relationship between a,b. Matlab code follows (sorry it's messy, just used the same notation that I was using with pen/paper):

function plot_ellipse(X, varargin)
% Plots an ellipse of the form x'*X*x <= 1

% plot vectors of the form a*u + b*v where u,v are eigenvectors of X

[V,D] = eig(X);
u = V(:,1);
v = V(:,2);
l1 = D(1,1);
l2 = D(2,2);

pts = [];

delta = .1;

for alpha = -1/sqrt(l1)-delta:delta:1/sqrt(l1)+delta
    beta = sqrt((1 - alpha^2 * l1)/l2);
    pts(:,end+1) = alpha*u + beta*v;
end
for alpha = 1/sqrt(l1)+delta:-delta:-1/sqrt(l1)-delta
    beta = -sqrt((1 - alpha^2 * l1)/l2);
    pts(:,end+1) = alpha*u + beta*v;
end

plot(pts(1,:), pts(2,:), varargin{:})

Upvotes: 1

Views: 7690

Answers (1)

fhoussiau
fhoussiau

Reputation: 61

I stumbled across this post while searching for this topic, and even though it's settled, I thought I might provide another simpler solution, if the matrix is symmetric.

Another way of doing this is to use the Cholesky decomposition of the semi-definite positive matrix E implemented in Matlab as the chol function. It computes an upper triangular matrix R such that X = R' * R. Using this, x'*X*x = (R*x)'*(R*x) = z'*z, if we define z as R*x.

The curve to plot thus becomes such that z'*z=1, and that's a circle. A simple solution is thus z = (cos(t), sin(t)), for 0<=t<=2 pi. You then multiply by the inverse of R to get the ellipse.

This is pretty straightforward to translate into the following code:

function plot_ellipse(E)
 % plots an ellipse of the form xEx = 1
 R = chol(E);
 t = linspace(0, 2*pi, 100); % or any high number to make curve smooth
 z = [cos(t); sin(t)];
 ellipse = inv(R) * z;
 plot(ellipse(1,:), ellipse(2,:))
end

Hope this might help!

Upvotes: 6

Related Questions