-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathdevilliers.m
174 lines (157 loc) · 6.14 KB
/
devilliers.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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
function [D,V,g,c,N,A,D2,V2,g2]=devilliers(m,N,K,L)
% [D,V,g,c,N,A,D2,V2,g2]=DEVILLIERS(m,N,K,L)
%
% Builds the tridiagonal matrix that we need to diagonalize in order to
% find the Jacobi/Bessel expansion coefficients of the radial part of the
% Slepian functions concentrating to a disk-shaped region.
%
% INPUT:
%
% m The fixed order of the concentration eigenfunctions [default: 3]
% N The Shannon number of the full concentration problem [default: 3]
% K The number of functions required [default: N]
% L The maximum degree of the Jacobi/Bessel expansion [default: 2N]
%
% OUTPUT:
%
% D The expansion coefficients, which is an (L+1)x(L+1) matrix with
% the row index the degree and column the eigenfunction rank
% g The integral equation eigenvalues
% V The concentration eigenvalues
% c The concentration factor c=2*sqrt(N)
% N The Shannon number of the concentration problem
% A The tridiagonal matrix whose eigenvectors are D
% D2 An alternative orthogonal set of coefficients from a symmetric A
% V2 An alternative set of concentration eigenvalues belonging to D2
% g2 An alternative set of integral equation eigenvalues
%
% SEE ALSO: SWDISK
%
% Last modified by dongwang-at-princeton.edu, 08/06/2008
% Last modified by fjsimons-at-alum.mit.edu, 06/30/2018
%
% NOTES:
%
% @Article{DeVilliers+2003,
% author = {G. D. de Villiers and F. B. T. Marchaud and E. R. Pike},
% title = {Generalized {G}aussian quadrature applied to an
% inverse problem in antenna theory: {II.} {T}he
% two-dimensional case with circular symmetry},
% journal = IP,
% year = 2003,
% volume = 19,
% pages = {755--778, doi: 10.1088/0266-5611/19/3/317}
% }
%
% @Article{Shkolnisky2007,
% author = {Yoel Shkolnisky},
% title = {Prolate spheroidal wave functions on a disc ---
% {I}ntegration and approximation of two-dimensional
% bandlimited functions},
% journal = ACHA,
% year = 2007,
% volume = 22,
% pages = {235--256, doi: 10.1016/j.acha.2006.07.002}
% }
%
% @Article{Slepian64,
% author = {David Slepian},
% title = {Prolate Spheroidal Wave Functions, {F}ourier
% analysis and Uncertainty --- {IV}: {E}xtensions to
% many dimensions; generalized prolate spheroidal
% functions},
% journal = {Bell Syst.~Tech.~J.},
% year = 1964,
% volume = 43,
% number = 6,
% pages = {3009--3057, doi: 10.1002/j.1538-7305.1964.tb01037.x}
% }
%
% This is probably equivalent to the three-term recursion of Slepian
% (1964). See also Shkolnisky (2007) and Xiao (2001) for the m=1/2
% case, i.e. for the "regular" one-dimensional prolate spheroidals.
% Supply defaults just in case, though usually this is a called routine
defval('m',3)
defval('N',3)
defval('K',ceil(N))
% Some heuristic defaults - make identical in SWDISK
defval('L',min(max([2*K ceil(2*N) 2*m 10]),84))
% Stick to the notation of Slepian (1964)
c=2*sqrt(N);
disp(sprintf('c = %4.2g ; m = %i ; N = %4.2g ; K = %i ; L = %i ',...
c,m,N,K,L))
% Where does the range start? From zero.
k=0:L; kk=0:L-1;
% Make the sub, main and superdiagonal
warning off MATLAB:divideByZero
% Eq. (52) in doi: 10.1088/0266-5611/19/3/317
ldiag=-c^2*(m+kk+1).^2./(2*kk+m+1)./(2*kk+m+2);
% Eq. (50) in doi: 10.1088/0266-5611/19/3/317
ondiag=(2*k+m+1/2).*(2*k+m+3/2)+...
c^2/2*(1+m^2./(2*k+m)./(2*k+m+2));
% Eq. (51) in doi: 10.1088/0266-5611/19/3/317
udiag=-c^2*( kk+1).^2./(2*kk+m+2)./(2*kk+m+3);
warning on MATLAB:divideByZero
% Form the tridiagonal non-symmetric matrix
A=tridiag(ldiag,ondiag,udiag);
% Now for the symmetric variety. Common form for both diagonals
% doi: 10.1016/j.acha.2006.07.002
uldiag=-c^2*(kk+1).*(m+kk+1)./sqrt(2*kk+m+1)./(2*kk+m+2)./sqrt(2*kk+m+3);
% And notice that A and AA and AAA have the same eigenvalues! But now, how do the
% eigenvectors scale?
AA=tridiag(uldiag,ondiag,uldiag);
% Supply the correct entry when the order was zero
if m==0
A(1,1)=1/2*3/2+c^2/2;
AA(1,1)=A(1,1);
end
% Now find the eigenvectors of this thing
if K>=L & K<=L+1
% Find them all - somewhat slower
[D,v]=eig(A);
% But when you solve the symmetric system you get orthogonal eigenvectors!
[D2,v2]=eig(AA);
difer(sort(v(:))-sort(v2(:)),[],[],NaN)
elseif K<L
% Return the smallest eigenvalues of the Sturm-Liouville system; these
% correspond to the largest eigenvalues of the concentration problem
OPTS.disp=0;
[D,v]=eigs(A,K,'SM',OPTS);
[D2,v2]=eigs(AA,K,'SM',OPTS);
difer(sort(v(:))-sort(v2(:)),[],[],NaN)
elseif K>L+1
error(sprintf('You must increase truncation degree to at least %i',K))
end
% Note that the columns of D are the eigenvectors; that they are not
% necessarily an orthogonal set; that each column has been normalized to
% energy one; so when using them in the expansion, we will still need to
% normalize the functions according to our particular wishes of the day.
% Note that the sorting may be non-existent or wrong
[v,i]=sort(diag(v)); D=D(:,i);
% The rows of D are the degrees l, the columns the rank of the function
% Now find the concentration eigenvalues as per
% Eq. (46) in 10.1002/j.1538-7305.1964.tb01037.x
g=c^(m+1/2)*D(1,1:K)./2^(m+1)./factorial(m+1)./sum(D(:,1:K),1);
V=g.^2*c;
g2=c^(m+1/2)*D2(1,1:K)./2^(m+1)./factorial(m+1)./sum(D2(:,1:K),1);
V2=g2.^2*c;
% Check that the eigenvalues are properly contained
if ~isnan(nanmean([V(V>1) V(V<0)]));
V1=V-1; V1=V1(V1>0); V0=V+1; V0=V0(V0<1);
if nanmean([V1 V0])>1000*eps
error('Eigenvalues exceeding 0 to 1 - definitely increase L');
end
end
% Build in a check on the convergence - we want the last term in the row
% sum to not matter at all for the first K requested functions (noting
% that technically, we may have requested K=L+1).
avco=abs(D(L+1,1:K))./sum(abs(D(1:L,1:K)),1);
% The tough criterion is when the sum converges, period
tough=mean(abs(avco(1:K)));
%disp(sprintf('Tougher convergence estimate %6.3e',tough))
% The lenient criterion is when the sum converges for the big eigenvalues
lenient=mean(abs(avco(g(1:K)>eps)));
%disp(sprintf('Lenient convergence estimate %6.3e',lenient))
% Formal error reporting on these two criteria
difer(tough,[],[])
difer(lenient,[],[])