-
Notifications
You must be signed in to change notification settings - Fork 188
/
lec17mldemo2.m
56 lines (45 loc) · 1.06 KB
/
lec17mldemo2.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
% MIT 18.335 - Lecture 17 MATLAB Demo 2
% The Divide-and-Conquer Algorithm
% Per-Olof Persson, November 14, 2007
format long
% Create symmetric matrix A
n=6;
randn('state',2)
A0=randn(n,n); A0=A0+A0';
A0,pause
% Reduce to tridiagonal
A=hess(A0);
A,pause
% Split
beta=A(n/2+1,n/2),pause
T1=A(1:n/2,1:n/2),pause
T2=A(n/2+1:n,n/2+1:n),pause
That1=T1;
That1(n/2,n/2)=T1(n/2,n/2)-beta;
That1,pause
That2=T2;
That2(1,1)=T2(1,1)-beta;
That2,pause
% Compute eigenvalues of That1,That2 exactly
% In the real algorithm this would be recursion
[Q1,D1]=eig(That1);
[Q2,D2]=eig(That2);
d=[diag(D1);diag(D2)],pause
z=[Q1(end,:), Q2(1,:)]',pause
neg=beta<0;
w=z*sqrt(abs(beta));
ll=-12:0.01:12;
ff=1+0*ll;
for ii=1:n
ff=ff+(-1)^neg*w(ii).^2./(d(ii)-ll);
end
plot(ll,ff)
axis([-12,12,-5,5])
line([-12,12],[0,0],'color','k')
for ii=1:n
line([d(ii),d(ii)],[-5,5],'linestyle','--','color',[.5,.5,.5])
end
pause
% Compare with true eigenvalues
l=eig(A0), pause
line(l,0*l,'linestyle','none','marker','.','markersize',16,'color',[0,.5,0])