function [ xp, yp, rmserr ] = approx_leastsquares ( f_string, a, b, m, nd )

%*****************************************************************************80
%
%% approx_leastsquares is an interactive least squares approximation program.
%
%  Discussion:
%
%    This function manages an interactive computation which constructs a
%    polynomial of degree m to approximate nd data values of a function f(x)
%    in the interval a <= x <= b.
%
%    The user enters a formula for f(x), the values of a and b,
%    the degree m, and the number of interpolation points nd.
%
%    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, rmserr ] = approx_leastsquares ( 'sin(x)', 0, 5, 5, 20 );
%    or, if called with no arguments, it will request them:
%      [ xp, yp, rmserr ] = approx_leastsquares ( );
%      Enter function formula, like x^2: sin(x)
%      Enter left limit, a: 0
%      Enter right limit, b: 5
%      Degree of polynomial: 5
%      Enter number of data points nd: 20
%
%    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 GNU LGPL license.
%
%  Modified:
%
%    11 August 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string f_string: the function to be approximated.
%
%    real a, b: the interval.
%
%    integer m: the degree of the approximating polynomial.
%
%    integer nd: the number of data points to use.
%
%  Output:
%
%    real xp(*), yp(*): sample values of the interpolant over [a,b].
%
%    real rmserr: the (discrete) root-mean-square (RMS) error between 
%    the sample values of the function f(x) and the approimating 
%    polynomial over [a,b].
%

%
%  Get f, a, b, m, nd.
%
  if ( nargin < 1 )
    f_string = input ( 'Enter function formula, like sin(x): ', '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 )
    m = input ( 'Enter degree of the approximating polynomial m: ' );
  elseif ( ischar ( m ) )
    m = str2num ( m );
  end

  if ( nargin < 5 )
    nd = input ( 'Enter number of data points nd: ' );
  elseif ( ischar ( nd ) )
    nd = str2num ( nd );
  end
%
%  Get the data values.
%
  xd = chebyspace ( a, b, nd );
  fd = f ( xd );
%
%  Evaluate the least squares approximant at the data points.
%
  yd = lagrange_approx ( m, nd, xd, fd, nd, xd );
%
%  Compute the discrete rms norm of the error.
%
  rmserr = sqrt ( sum ( ( fd(1:nd) - yd(1:nd) ).^2 ) / nd );
%
%  Evaluate the least squares approximant at the plotting points.
%
  np = 101;
  xp = linspace ( a, b, np );
  xp = xp';
  yp = lagrange_approx ( m, nd, xd, fd, np, xp );

  return
end
function x = chebyspace ( a, b, n )

%*****************************************************************************80
%
%% chebyspace creates a vector of Chebyshev spaced values.
%
%  Discussion:
%
%    A set of points with Chebyshev spacing is generated for the
%    interval [A,B].
%
%    An R8VEC is a vector of R8's.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    08 June 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries in the vector.
%
%    Input, real A, B, the interval endpoints.
%
%    Output, real X(N,1), a vector of Chebyshev spaced data.
%
  x = zeros ( n, 1 );

  if ( n == 1 )

    x(1) = ( a + b ) / 2.0;

  else

    for i = 1 : n

      theta = ( n - i ) * pi / ( n - 1 );

      c = cos ( theta );

      if ( mod ( n, 2 ) == 1 )
        if ( 2 * i - 1 == n )
          c = 0.0;
        end
      end

      x(i) = ( ( 1.0 - c ) * a  ...
             + ( 1.0 + c ) * b ) ...
             /   2.0;

    end

  end

  return
end
function yi = lagrange_approx ( m, nd, xd, yd, ni, xi )

%*****************************************************************************80
%
%% LAGRANGE_APPROX evaluates the Lagrange approximant of degree M.
%
%  Discussion:
%
%    The Lagrange approximant L(M,ND,XD,YD)(X) is a polynomial of
%    degree M which approximates the data (XD(I),YD(I)) for I = 1 to ND.
%
%    We can represent any polynomial of degree M+1 as the sum of the Lagrange 
%    basis functions at the M+1 Chebyshev points.
%
%      L(M)(X) = sum ( 1 <= I <= M+1 ) C(I) LB(M,XC)(X)
%
%    Given our data, we can seek the M+1 unknown coefficients C which minimize
%    the norm of || L(M)(XD(1:ND)) - YD(1:ND) ||.
%
%    Given the coefficients, we can then evaluate the polynomial at the
%    points XI.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 August 2012
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer M, the polynomial degree.
%
%    integer ND, the number of data points.
%    ND must be at least 1.
%
%    real XD(ND,1), the data points.
%
%    real YD(ND,1), the data values.
%
%    integer NI, the number of interpolation points.
%
%    real XI(NI,1), the interpolation points.
%
%  Output:
%
%    real YI(NI,1), the interpolated values.
%

%
%  Evaluate the Chebyshev points.
%
  xc = chebyspace ( -1.0, +1.0, m + 1 );
%
%  Evaluate the Lagrange basis functions for the Chebyshev points 
%  at the data points.
%
  ld = lagrange_basis ( m + 1, xc, nd, xd );
%
%  The value of the Lagrange approximant at each data point should
%  approximate the data value: LD * YC = YD, where YC are the unknown
%  coefficients.
%
  yc = ld \ yd;
%
%  Now we want to evaluate the Lagrange approximant at the "interpolant
%  points": LI * YC = YI;
%
  li = lagrange_basis ( m + 1, xc, ni, xi );
  yi = li * yc;

  return
end
function lb = lagrange_basis ( nd, xd, ni, xi ) 

%*****************************************************************************80
%
%% LAGRANGE_BASIS evaluates a 1D Lagrange basis.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 October 2012
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer ND, the number of data points.
%
%    real XD(ND,1), the interpolation nodes.
%
%    integer NI, the number of evaluation points.
%
%    real XI(NI,1), the evaluation points.
%
%  Output:
%
%    real LB(NI,ND), the value, at the I-th point XI, of the
%    Jth basis function.
%
  lb = zeros ( ni, nd );
  
  for i = 1 : ni
    for j = 1 : nd
      lb(i,j) = prod ( ( xi(i) - xd(1:j-1)  ) ./ ( xd(j) - xd(1:j-1)  ) ) ...
              * prod ( ( xi(i) - xd(j+1:nd) ) ./ ( xd(j) - xd(j+1:nd) ) );
    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 GNU LGPL 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

