function [ x, it ] = opt_quadratic ( f_string, x1, x2, x3, n, x_tol, y_tol )

%*****************************************************************************80
%
%% opt_quadratic uses quadratic interpolation to seek a critical point of f(x).
%
%  Discussion:
%
%    The user enters a formula for f(x), the starting values x1, x2, x3,
%    the iteration limit n, and tolerances for the size of the final x and
%    y intervals.  
%
%    The program returns an estimate for a critical point x, which might
%    be a minimum, maximum or just an inflection point.  It also returns
%    the number of iterations.
%
%    The program can be invoked by a function call, in which case the
%    string specifying f(x) must be quoted:
%
%      [ x, it ] = opt_quadratic ( 'x^2+sin(4*x)', -1, 1, 2, 25, 0.00001, 0.000001 )
%
%    or, if called with no arguments, it will request them:
%
%      [ x, it ] = opt_quadratic ( );
%      Enter function formula, like x^2: x^2+sin(4*x)
%      Enter point 1, x1: -1
%      Enter point 2, x2: 1
%      Enter point 3, x3: 2
%      Enter number of iterations, n: 25
%      Enter x tolerance: 0.00001
%      Enter y tolerance: 0.00001
%
%    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:
%
%    02 August 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string f_string: the function to be measured.
%
%    real x1, x2, x3: three distinct starting points.
%
%    integer n, the maximum number of iterations allowed.
%    n must be at least 1.
%
%    real x_tol, y_tol, tolerances for the width of the x and y intervals.
%    These must be positive values.
%
%  Output:
%
%    real x, the estimate for the location of a critical point of f(x).
%
%    integer it, the number of steps taken.

%
%  Get f, x1, x2, x3, n, x_tol, y_tol.
%
  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 )
    x1 = input ( 'Enter point 1, x1: ' );
  elseif ( ischar ( x1 ) )
    x1 = str2num ( x1 );
  end

  if ( nargin < 3 )
    x2 = input ( 'Enter point 2, x2: ' );
  elseif ( ischar ( x2 ) )
    x2 = str2num ( x2 );
  end

  if ( nargin < 4 )
    x3 = input ( 'Enter point 3, x3: ' );
  elseif ( ischar ( x3 ) )
    x3 = str2num ( x3 );
  end

  if ( nargin < 5 )
    n = input ( 'Enter iteration limit, n: ' );
  elseif ( ischar ( n ) )
    n = str2num ( n );
  end

  if ( n < 1 )
    error ( 'N must be at least 1.' );
  end

  if ( nargin < 6 )
    x_tol = input ( 'Enter x tolerance, x_tol: ' );
  elseif ( ischar ( x_tol ) )
    x_tol = str2num ( x_tol );
  end

  if ( x_tol <= 0.0 )
    error ( 'x_tol must be positive.' );
  end

  if ( nargin < 7 )
    y_tol = input ( 'Enter y tolerance, y_tol: ' );
  elseif ( ischar ( y_tol ) )
    y_tol = str2num ( y_tol );
  end

  if ( y_tol <= 0.0 )
    error ( 'y_tol must be positive.' );
  end
%
%  Carry out the procedure.
%
  x = [ x1; x2; x3 ];

  for i = 1 : n
%
%  Define the Vandermonde matrix for p, the coefficients of the 
%  interpolating quadratic polynomial.
%
    V = [ x.^2, x, ones(3,1) ];
%
%  Solve the Vandermonde system for p.
%
    p = V \ f(x);
%
%  Use p to determine the minimizer of the quadratic interpolant.
%
%  If y = p(1)*x^2+p(2)*x+p(3),
%     y' = 2*p(1)x+p(2)
%  so y' = 0 at x = -p(2)/2/p(1)
%
    x = [ x(2); x(3); - p(2) / 2.0 / p(1) ];
%
%  Accept this value if the x and y intervals are small enough.
%
    if ( ( abs (   x(2)  -   x(3) ) < x_tol ) && ...
         ( abs ( f(x(2)) - f(x(3) ) < y_tol ) ) )
      x = x(3);
      return
    end

  end

  error ( 'The method did not converge.' );

  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

