-
Notifications
You must be signed in to change notification settings - Fork 188
/
ellipse.m
86 lines (68 loc) · 1.69 KB
/
ellipse.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
function [K,Q,D,W,l,qf] = ellipse(A)
%
% [K,Q,D,W,l,qf] = ellipse(A)
%
% Draws the ellipse
% x' * K * x = 1
% corresponding to the positive definite symmetric 2 by 2 matrix
% where K = A if A is symmetric; otherwise K = A' * A
%
% K - symmetric positive definite matrix used
% Q - orthogonal eigenvectors or principal axes of ellipse
% D - diagonal matrix with eigenvalues, so K = Q * D * Q'
% l = 1/sqrt(diag(D)) - lengths of semi-axes of ellipse
% W = Q*D^(-1/2) - linear transformation mapping unit circle to ellipse
% qf - string giving quadratic form x'*K*x
%
% See also ELLIPSOID
%
if size(A) ~= [2 2], error('argument must be 2 by 2 matrix'); end
if A == A'
K = A;
else
K = A' * A;
end
if K(1,1) == 1, q1s = ''; else q1s = [num2str(K(1,1)),' ']; end
if K(2,2) == 1, q2s = ''; else q2s = [num2str(K(2,2)),' ']; end
q12 = abs(2*K(1,2)); q12s = num2str(q12);
s12 = sign(K(1,2));
if s12 > 0, sg = [' + ',q12s,' x y'];
else if s12 < 0, sg = [' - ',q12s,' x y'];
else sg = '';
end
end
qf = [q1s,'x^2',sg,' + ',q2s,'y^2'];
disp(' ')
disp('Quadratic form')
disp(' ')
disp(qf)
disp(' ')
disp('K = ')
disp(K)
[Q,D] = eig(K);
disp('Eigenvectors - principal axes')
disp(' ')
disp(Q)
disp('Eigenvalues')
ev = diag(D);
disp(' ')
disp(ev)
disp('Lengths of principal axes')
l = 1./sqrt(ev);
disp(' ')
disp(l)
W = Q * 1/sqrt(D);
t = linspace(0,2*pi,100);
plot(W(1,1) .* cos(t) + W(1,2) .* sin(t), W(2,1) .* cos(t) + W(2,2) .* sin(t),'m')
W1 = W(1,:); W2 = W(2,:);
N = diag([1 + 1/norm(W(:,1)),1 + 1/norm(W(:,2))]);
WN = W * N;
WN1 = WN(1,:); WN2 = WN(2,:);
xa = [W1;WN1];
ya = [W2;WN2];
line(xa,ya);
xa = - xa;
ya = - ya;
line(xa,ya);
axis equal;
title(['Ellipse: ',qf,' = 1']);