%% Bernstein polynomials 
% Nick Trefethen, 15th 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.

%%
% 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$.
% $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