%% 1D subdivision and the average degree reduction
% Alex Townsend, August 2013
%%
% (Chebfun Example: roots/AverageDegreeReduction1D.m)
% [Tags: #roots, #subdivision]
function AverageDegreeReduction1D
%% 1D subdivision
% Generally, Chebfun approximates globally smooth functions $f(x)$ by global
% polynomial interpolants [1]. However, in the `roots` command, if the degree
% of a chebfun is over $49$, the algorithm recursively subdivides, effectively
% representing $f(x)$ by a piecewise smooth polynomial before finding its
% roots [1,2]. For example, if one subdivision is required, instead of solving
% $f(x) = 0$ on $[-1,1]$, the algorithm solves $f_{left}(x) = 0$ on $[-1,r]$
% and $f_{right}(x)=0$ on $[r,1]$, where $r$ is an arbitrary small constant.
% Furthermore, it continues to subdivide $f$ on smaller domains until it is of
% degree $49$ or less on each subinterval. The expectation is that if $f$ is
% of degree $n$, then $f_{left}$ and $f_{right}$ are of lower degree than $n$,
% but this is not always the case.
%% The average degree reduction in 1D
% Let $\tau$ be a parameter that measures the average degree reduction of
% $f(x)$. That is, if $f(x)$ is of degree $n$, then $f_{left}$ and $f_{right}$
% are of degree about $\tau n$. The parameter $\tau$ can be simply taken as
% $1/2$ as done in [3], but this is a simplification. We approximate $\tau$ by
% subdividing and calculating the degree reduction. The function we use to
% estimate $\tau$ is `compute_tau` (see below).
%% Oscillatory functions
% For oscillatory functions the numerical degree of a chebfun closely follows
% Nyquist's rule, i.e., resolution requires on average $\pi$ points per
% wavelength when using Chebyshev points. When the interval is halved, the
% number of oscillations is halved and hence, the number of points to resolve
% the functions on each subinterval is also halved. Therefore, we expect
% $\tau$ to be approximately $1/2$.
% For example,
M = 5e2;
f = @(x) sin(M*x);
compute_tau(f, 2^3)
%%
% Here is a diagram that shows the various degrees of the piecewise polynomial
% after the first three subdivisions.
subdivisionDiagram(f,3)
%%
% As $M\to\infty$, the function becomes more oscillatory and $\tau\to 1/2$.
%% Locally nonsmooth functions
% Functions that are locally nonsmooth at a point on the interval but are
% smooth everywhere else can have an average degree reduction significantly
% less than $1/2$. This is especially true if the subdivision almost exactly
% isolates the difficult point. For instance,
f = @(x) abs(x).^3;
compute_tau(f, 2)
subdivisionDiagram(f,1)
%%
% In this example the degree reduction after one subdivision makes
% $\tau\ll 1$. However, the subdivision within Chebfun subdivides $[-1,1]$ in
% $[-1,r]$ and $[r,1]$, where $r$ is an arbitrary constant so such small
% values for $\tau$ are rarely realised.
%%
% If the subdivision does not divide exactly on the point of discontinuity,
% then $\tau$ is just something less than $1/2$.
f = @(x) abs(x-.01).^7;
compute_tau(f, 2^3)
%% Analytic functions
% For functions that are analytic on the real line with simple poles in the
% complex plane, Elliott's method [1] can be used to understand $\tau$. For
% example, consider the function $1/(x-1.0001)$.
f = @(x) 1./(x-1.0001);
compute_tau(f, 8)
subdivisionDiagram(f,3)
%%
% This function is meromorphic in the complex $x$-plane with a simple pole at
% $x=1.0001$, and for these types of functions the Chebyshev coefficient are
% explicitly known (see [4,5]). In this case it can be shown that the
% numerical degrees of $f$ on the subinterval $[a,1]$, for $a = -1,0,1/2,1/4$,
% are the following integers:
for lvl = 1:4
c = 1.0001; % location of pole
a = 1-2.^(-lvl+1); b = 1; % interval
A = -(-2*c+b+a)/(b-a); % Elliott's formula
m(lvl) = ceil(log(-4*(a-c)*eps^(-1)/(b-a)/sqrt(A^2-1))...
./log(A+sqrt(A.^2-1)));
end
m = m'
%%
% The reason the numerical degrees found by Elliott's formula do not exactly
% match those found by Chebfun is that the Chebfun constructor has a varying
% tolerance, and the tail of the coefficients is not always truncated at the
% same level.
%% A globally difficult function
% Some functions have values of $\tau$ larger than $1/2$, at least for the
% first few subdivisions. These functions tend to be difficult at many points
% on the interval. For example,
f = chebfun(@(x) abs(prod(x-linspace(-1,1,10)).^3),'vectorise');
compute_tau(f, 2)
%% Conclusion
% The average degree reduction parameter $\tau$ can take almost any value
% between $0$ and $1$. For functions that are equally difficult throughout the
% whole interval, $\tau$ is close to $1/2$. For functions that are difficult,
% local $\tau$ is often small. Therefore, one would expect that Chebfun's
% rootfinding algorithm would be more efficient at finding the roots of
% functions that have a small $\tau$. However, this is not the case, as the
% following two examples (of roughly the same size) show:
M = 2000;
f = chebfun(@(x) sin(M*x)); % example with tau approximately 1/2.
tic, roots(f); toc
%%
c = -1.0001;
f = chebfun(@(x) 1./(x+c)); % example with tau approximately 0.375
tic, roots(f); toc
%%
% The reason is that the overhead of each subdivision (as implemented) is
% high, and many more levels of subdivision are required for $1/(x+c)$.
%%
% Here are the functions used above to compute $\tau$ and plot the subdivision
% diagrams.
function compute_tau(f, N)
% COMPUTE_TAU estimate the average degree reduction parameter
d = [-1 1]; len = length(chebfun(f,d));
newends = d(1) + (0:N)./N*(d(2)-d(1));
newlen = length(chebfun(f, newends))/N;
tau = (newlen/len).^(1/log2(N));
fprintf('tau = %1.5f',tau)
end
%%
function subdivisionDiagram(f, levels)
% SUBDIVISIONDIAGRAM draw a diagram to show subdivision and polynomial
% degrees.
g = chebfun(f);
LW = 'linewidth'; lw = 2;
for lvl = 0:levels
plot([-1 1],[lvl lvl],'k',LW,lw), hold on,
xx = linspace(-1,1,2^lvl+1);
newg = chebfun(f, xx); j = 1;
for x = xx(1:end-1)
plot([x x],[-.1 .1]+lvl,'k',LW,lw)
funs = newg.funs;
text(x+2.^(-lvl)-.09,lvl+.1,sprintf('N = %u',length(funs{j})-1));
j = j + 1;
end
plot([xx(end) xx(end)],[-.1 .1]+lvl,'k',LW,lw)
end
set(gca,'ytick',[])
set(gca,'yticklabel',[])
hold off
end
%%
end
%% References
%
% 1. Z. Battels and L. N. Trefethen, An extension of MATLAB to continuous
% functions and operators, _SIAM Journal on Scientific Computing_, 25
% (2004), pp.1743-1770.
%
% 2. J. P. Boyd, Computing zeros on a real interval through Chebyshev
% expansion and polynomial rootfinding, _SIAM Journal on Numerical
% Analysis_, 40 (2002), pp. 1666-1682.
%
% 3. J. P. Boyd, Finding the zeros of a univariate equation: Proxy
% rootfinders, Chebyshev interpolation, and the companion matrix, _SIAM
% Review_, 55 (2013), pp. 375-396.
%
% 4. D. Elliott, The evaluation and estimation of the coefficients in the
% Chebyshev series expansion of a function, _Mathematics of Computation_,
% 18 (1964), pp. 274-284.
%
% 5. M. Richardson, Exact Chebyshev expansion coefficients of a function,
% Chebfun Example, June 2012.