function bernstein_poly_ab_approx_test ( )

%*****************************************************************************80
%
%% BERNSTEIN_POLY_AB_APPROX_TEST examines a sequence of approximants.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 January 2016
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'BERNSTEIN_POLY_AB_APPROX_TEST\n' );
  fprintf ( 1, '  BERNSTEIN_POLY_AB_APPROX evaluates the Bernstein polynomial\n' );
  fprintf ( 1, '  approximant to a function F(X).\n' );
  fprintf ( 1, '  As N, the polynomial degree increases, the norm of the error\n' );
  fprintf ( 1, '  should decrease.\n' );

  a = 1.0;
  b = 3.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Approximate y=sin(x) over the interval [1,3]\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   N         ||F-B|| / ||F|\n' );
  fprintf ( 1, '\n' );

  for n = 0 : 20

    xdata = ( linspace ( a, b, n + 1 ) )';
    ydata = sin ( xdata );
%
%  Compare the true function and the approximant.
%
    nval = 501;
    xval = ( linspace ( a, b, nval ) )';
    ytrue = sin ( xval );
    yapprox = bernstein_poly_ab_approx ( n, a, b, ydata, nval, xval );

    f_norm = sum ( abs ( ytrue(1:nval) ) ) / ( b - a );
    e_norm = sum ( abs ( ytrue(1:nval) - yapprox(1:nval) ) ) / ( b - a );
    fprintf ( 1, '  %2d  %14.6g\n', n, e_norm / f_norm );

  end

  a = -1.0;
  b = +1.0;

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Approximate y=heaviside(x) over the interval [-1,+1]\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '   N         ||F-B|| / ||F|\n' );
  fprintf ( 1, '\n' );

  for n = 0 : 20

    xdata = ( linspace ( a, b, n + 1 ) )';
    ydata = heaviside ( xdata );
%
%  Compare the true function and the approximant.
%
    nval = 501;
    xval = ( linspace ( a, b, nval ) )';
    ytrue = heaviside ( xval );
    yapprox = bernstein_poly_ab_approx ( n, a, b, ydata, nval, xval );

    f_norm = sum ( abs ( ytrue(1:nval) ) ) / ( b - a );
    e_norm = sum ( abs ( ytrue(1:nval) - yapprox(1:nval) ) ) / ( b - a );
    fprintf ( 1, '  %2d  %14.6g\n', n, e_norm / f_norm );

  end

  return
end
