%% Bernstein polynomials 
% Nick Trefethen, May 2012

%%
% (Chebfun example approx/BernsteinPolys.m)
% [Tags: #Bernstein, #Gibbs, #Weierstrass, #binomial]

function BernsteinPolys()

%%
% The Weierstrass Approximation Theorem asserts that a continuous function $f$
% on a bounded interval like $[0,1]$ can be approximated by polynomials (i.e.,
% approximated as closely as you like in the supremum norm). Weierstrass
% proved this in 1885 by a diffusion argument: if $f$ diffuses however little,
% it becomes an entire function, which can be approximated by truncating the
% Taylor series. (Before the diffusion, one first extends $f$ to a continuous
% with compact support on the whole real line.)

%%
% Bernstein gave a proof of the Weierstrass Approximation Theorem in 1912 that
% is a kind of discrete version of this diffusion proof: it replaces the
% continuous diffusion by a random walk on an equispaced grid in $[0,1]$.
% While this is perhaps a little more complicated conceptually, it is
% mathematically more elementary since you don't need any analysis and you
% don't need to truncate a series, for the polynomials emerge directly.

%%
% Specifically, for each positive integer $n$, the degree $n$ Bernstein
% polynomial for $f$ is
% $$ B_n(x) = \sum_{k=0}^n f(k/n) {n\choose k} x^k (1-x)^{n-k}. $$
% Note that this is basically a binomial expansion.  The formula tells us that
% to evaluate $B_n(x)$, we can imagine a biased coin that comes up heads with
% probability $x$ and tails with probability $1-x$. Then $B_n(x)$ is the
% expected result that you'll get if you start at $x=0$ and toss the coin $n$
% times, moving right on the grid if you get a heads, and evaluate $f$ when
% you finish tossing.

%%
% Let's demonstrate in Chebfun.  Here is a continuous function on $[0,1]$:
LW = 'linewidth'; lw = 1.6;
s = chebfun('s',[0 1]);
f = min(abs(s-.3),2*abs(s-.7));
f = s + max(0,1-5*f);
hold off, plot(f,LW,lw)

%%
% Since $B_n$ is a polynomial of degree $n$, we can construct it by evaluating
% it on a grid of $n+1$ points and then interpolating. For stability these
% should be Chebyshev points, not equispaced. Here is an elementary code to do
% this at least for small values of $n$. Note that it isn't really stable, so
% it turns warnings off.
function Bn = Bn(f,n)
warning off
x = chebpts(n+1,[0 1]);
Bndata = zeros(size(x));
for k = 0:n
  Bndata = Bndata + f(k/n)*nchoosek(n,k).*x.^k.*(1-x).^(n-k);
end
Bn = chebfun(Bndata,[0 1]);
warning on
end

%%
% To illustrate the behavior of Bernstein polynomials, here we see slow
% convergence as $n$ increases.
for n = [25 50 100]
  hold off, plot(f,LW,lw)
  hold on, plot(Bn(f,n),'r',LW,lw)
  title(['n = ' int2str(n)],'fontsize',14)
  snapnow
end

%%
% Note a signature feature of Bernstein polynomial approximations, their
% monotonicity in various senses.  There is never any Gibbs phenomenon.

%%
% On the other hand, since these approximations depend on the central limit
% theorem to give accuracy as $n$ gets large, they take no advantage at all
% of smoothness.  Here for example is a repetition of the last experiment
% for a far smoother function.
f = s + exp(-50*(s-.3).^2) + exp(-200*(s-.7).^2);
for n = [25 50 100]
  hold off, plot(f,LW,lw)
  hold on, plot(Bn(f,n),'r',LW,lw)
  title(['n = ' int2str(n)],'fontsize',14)
  snapnow
end

%%
% Though $f$ is now entire, the convergence is not really better than before.
% By contrast we know that $n=100$ is more than enough for Chebyshev
% interpolation to nail this function to machine precision:

length(f)

end
