/
Localization.m
89 lines (80 loc) · 3.03 KB
/
Localization.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
%% Low-rank approximation and localized singularities
% Nick Trefethen, April 2016
%%
% (Chebfun example approx2/Localization.m)
% [Tags: #approximation, #2D, #Chebfun2]
%% 1. Introduction
% This example is a companion to "Low-rank approximation
% and alignment with axes" [3]. Here, we look at another
% property of a function that the low-rank algorithms of
% Chebfun2 and Chebfun3 can take advantage of:
% localized singularities or near-singularities.
% Again the presentation follows the discussion of [2].
%% 2. Complex singularity near the middle of the domain
% For example, here is a bivariate Runge function with a
% broad spike, which we approximate to about
% ten digits. The rank is a bit smaller than the length,
% so we are getting some data compression.
ep = 1e-10;
f = chebfun2(@(x,y) 1./(1+(x-.2).^2+(y-.5).^2),'eps',ep);
r = rank(f), [m,n] = length(f)
%%
% If we change $1$ to $0.001$ in the denominator, the
% spike becomes much more localized. Now there is a
% dramatic difference between the rank and the length,
% showing excellent compression.
f = chebfun2(@(x,y) 1./(0.001+(x-.2).^2+(y-.5).^2),'eps',ep);
r = rank(f), [m,n] = length(f)
%%
% Let's make a plot. Each red circle shows a pivot value chosen
% by Chebfun2's approximation to Gaussian elimination with
% complete pivoting, in which
% at each step, a rank-1 cross is subtracted corresponding to
% the largest function value still remaining.
% For a full explanation, see [1].
p = f.pivotValues; a = f.pivotLocations;
x = a(:,1); y = a(:,2); n = length(p)
x = x(1:n); y = y(1:n);
for j = 1:n
plot([-1 1],y(j)*[1 1],'-k'), hold on
plot(x(j)*[1 1],[-1 1],'-k')
end
plot(x,y,'or')
set(gca,'xtick',[-1 0.2 1],'ytick',[-1 0.5 1])
axis([-1 1 -1 1]), axis square, hold off
%% 3. Real singularity outside a corner of the domain
% Here is another example. This time, the singularity is
% located outside the domain, near one corner. When the
% singularity is not very close to the corner, there is
% not much compression:
g = chebfun2(@(x,y) 1./((x+1.2).^2 + (y+1.2).^2),'eps',ep);
r = rank(g), [m,n] = length(g)
%%
% Changing 1.2 to 1.02 makes the compression striking:
g = chebfun2(@(x,y) 1./((x+1.02).^2 + (y+1.02).^2),'eps',ep);
r = rank(g), [m,n] = length(g)
%%
% Here is the analogous picture:
p = g.pivotValues; a = g.pivotLocations;
x = a(:,1); y = a(:,2); n = length(p)
x = x(1:n); y = y(1:n);
for j = 1:n
plot([-1 1],y(j)*[1 1],'-k'), hold on
plot(x(j)*[1 1],[-1 1],'-k')
end
plot(x,y,'or')
set(gca,'xtick',[-1 0.2 1],'ytick',[-1 0.5 1])
axis([-1 1 -1 1]), axis square, hold off
%% References
%
% 1. A. Townsend and L. N. Trefethen,
% An extension of Chebfun to two
% dimensions, _SIAM Journal on Scientific Computing_, 35
% (2013), C495-C518.
%
% 2. L. N. Trefethen, Cubature, approximation, and isotropy
% in the hypercube, manuscript, March 2016.
%
% 3. L. N. Trefethen, Low-rank approximation and alignment
% with axes, Chebfun example,
% http://www.chebfun.org/examples/approx2/Alignment.html.