-
Notifications
You must be signed in to change notification settings - Fork 227
/
compute_laplacian_matrix.m
59 lines (50 loc) · 1.17 KB
/
compute_laplacian_matrix.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
function L = compute_laplacian_matrix(n,options)
% compute_laplacian_matrix - compute the laplacian matrix
%
% L = compute_laplacian_matrix(n,options);
%
% options.ndims set the dimensionality (1 or 2) of the operator.
% options.bound set the boundary conditions ('per' of 'sym')
%
% Copyright (c) 2007 Gabriel PeyrŽ
options.null = 0;
if isfield(options, 'bound')
bound = options.bound;
else
bound = 'per';
end
if isfield(options, 'ndims')
ndims = options.ndims;
else
ndims = 2;
end
if ndims==1
L = -2*speye(n);
for i=1:n
L(i, mod(i-2,n)+1) = 1;
L(i, mod(i,n) +1) = 1;
end
return;
end
% compute the laplacian matrix
L = -4*speye(n^2);
[Y,X] = meshgrid(1:n,1:n);
deltax = [-1 1 0 0];
deltay = [0 0 -1 1];
for i=1:4
Xi = X + deltax(i);
Yi = Y + deltay(i);
if strcmp(bound,'per')
Xi = mod(Xi-1,n)+1;
Yi = mod(Yi-1,n)+1;
else
I = find(Xi<1); Xi(I) = 2 - Xi(I);
I = find(Yi<1); Yi(I) = 2 - Yi(I);
I = find(Xi>n); Xi(I) = 2*n+1 - Xi(I);
I = find(Yi>n); Yi(I) = 2*n+1 - Yi(I);
end
k1 = X + n*(Y-1);
k2 = Xi + n*(Yi-1);
k = k1 + n^2*(k2-1);
L(k(:)) = 1;
end