function [ xp, yp, maxerr ] = approx_bernstein ( f_string, a, b, n )

%*****************************************************************************80
%
%% approx_bernstein() is an interactive Bernstein approximation program.
%
%  Discussion:
%
%    This function manages an interactive computation which uses
%    n equally spaced nodes to approximate a function f(x)
%    in the interval a <= x <= b.  
%
%    The user enters a formula for f(x), the values of a and b,
%    and the number of equally spaced approximation points n.
%
%    The program returns a set of 101 pairs of values (xp,yp) that can be
%    used to plot the approximation polynomial.
%
%    The program can be invoked by a function call, in which case the
%    string specifying f(x) must be quoted:
%      [ xp, yp, maxerr ] = approx_bernstein ( 'x^2', -1, 3, 2 );
%    or, if called with no arguments, it will request them:
%      [ xp, yp, maxerr ] = approx_bernstein ( );
%      Enter function formula, like x^2: x^2
%      Enter left limit, a: -1
%      Enter right limit, b: 3
%      Enter number of interpolation points n: 2
%
%    The function is specified as a string which is either:
%    * a MATLAB expression using the argument 'x', 
%    * the name of an M-file followed by the argument '(x)'.
%
%    The string should not contain any spaces between symbols, except
%    when it is passed as a function argument in quotes.
%
%    It is not necessary to use the "dot" notation for expressions involving
%    '*', '/', or '^', but it doesn't hurt either.
%
%    Examples of function specifications:
%
%      x^2
%      x.^2
%      3/(x^4+5*x-6)
%      sin(7*x)*sqrt(x)/8
%      wiggle(x)     <-- where "wiggle.m" is a user-provided M file.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    09 August 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string f_string: the function to be approximated.
%
%    real a, b: the interval.
%
%    integer n: the number of points to use.
%
%  Output:
%
%    real xp(*), yp(*): sample values of the approximant over [a,b].
%
%    real maxerr: an estimate for the maximum error between the
%    function f(x) and the n-point Bernstein approximant over [a,b].
%

%
%  Get F, A, B, N.
%
  if ( nargin < 1 )
    f_string = input ( 'Enter function formula, like x^2: ', 's' );
  end
  f_element = vector_to_element ( f_string );
  f = str2func ( [ '@(x)', f_element ] );

  if ( nargin < 2 )
    a = input ( 'Enter left limit a: ' );
  elseif ( ischar ( a ) )
    a = str2num ( a );
  end

  if ( nargin < 3 )
    b = input ( 'Enter right limit b: ' );
  elseif ( ischar ( b ) )
    b = str2num ( b );
  end

  if ( nargin < 4 )
    n = input ( 'Enter number of Bernstein points n: ' );
  elseif ( ischar ( n ) )
    n = str2num ( n );
  end
%
%  Get the data at equally spaced points.
%
  if ( n == 1 )
    xd = ( a + b ) / 2.0;
  else
    xd = linspace ( a, b, n );
  end
  yd = f ( xd );
%
%  Estimate maximum error at 10001 points.
%
  ne = 10001;
  xe = linspace ( a, b, ne );
  xe = xe';
  ye = bernstein_approximant ( n, a, b, yd, ne, xe );
  fe = f ( xe );
  maxerr = max ( abs ( ye - fe ) );
%
%  Return a small number of points for plotting.
%
  np = 101;
  xp = linspace ( a, b, np );
  xp = xp';
  yp = bernstein_approximant ( n, a, b, yd, np, xp );
 
  return
end
function ye = bernstein_approximant ( nd, a, b, yd, ne, xe )

%*****************************************************************************80
%
%% bernstein_approximant(): Bernstein polynomial approximant to F(X) on [A,B].
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    10 August 2019
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    David Kahaner, Cleve Moler, Steven Nash,
%    Numerical Methods and Software,
%    Prentice Hall, 1989,
%    ISBN: 0-13-627258-4,
%    LC: TA345.K34.
%
%  Input:
%
%    integer nd, the number of approximation points.
%
%    real a, b, the endpoints of the approximation interval.
%
%    real yd(nd), the data values at nd equally spaced points in [A,B].
%
%    integer ne, the number of evaluation points.
%
%    real xe(ne), the evaluation points.
%
%  Output:
%
%    real ye(ne), the Bernstein approximant values.
%
  xe = xe(:);
  ye = zeros ( ne, 1 );
  be = zeros ( ne, nd );

  if ( nd == 1 )

    be(1:ne,1) = ones ( ne, 1 );
 
  else

    be(1:ne,1) = ( b - xe ) / ( b - a );
    be(1:ne,2) = ( xe - a ) / ( b - a );

    for j = 3 : nd

      k = j;
      be(1:ne,k) = ( xe - a ) ./ ( b - a ) .* be(1:ne,k-1);
      for k = j - 1 : -1 : 2
        be(1:ne,k) = ( xe - a ) ./ ( b - a ) .* be(1:ne,k-1) ...
                   + ( b - xe ) ./ ( b - a ) .* be(1:ne,k);
      end
      k = 1;
      be(1:ne,k) = ( b - xe ) ./ ( b - a ) .* be(1:ne,k);

    end

    for j = 1 : nd
      ye(1:ne) = ye(1:ne) + yd(j) .* be(1:ne,j);
    end

  end

  return
end
function f_element = vector_to_element ( f_vector )

%*****************************************************************************80
%
%% vector_to_element() rewrites a MATLAB vector function elementwise.
%
%  Discussion:
%
%    In MATLAB, by default, the operators '*', '/', and '^' are interpreted
%    as linear algebraic operations.  In many cases, a user simply desired
%    the elementwise operation instead, which MATLAB indicates using a 
%    dotted notation: '.*', './', and '.^'.
%
%    This function converts a string to a form that uses dotted operators.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    29 June 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string f_vector, a string presumably representing a MATLAB expression.
%
%  Output:
%
%    string f_element, a copy of f_vector in which all operators have been
%    replaced by their dotted (elementwise) form.
%
  n_vector = length ( f_vector );
  f_element = [''];

  c_prev = '';
  for iv = 1 : n_vector
    c = f_vector(iv);
    if ( ( c == '*' || c == '/' || c == '^' ) && ( c_prev ~= '.' ) )
      f_element = [ f_element, ' .', c, ' ' ];
    else
      f_element = [ f_element, c ];
    end
    c_prev = c;
  end

  return
end

