function q = quad_gauss ( f_string, a, b, n )

%*****************************************************************************80
%
%% quad_gauss() is an interactive Gauss quadrature program.
%
%  Discussion:
%
%    This function manages an interactive computation which uses
%    an n-point Gauss quadrature rule to estimate the integral of a function
%    f(x) over the interval [a,b].
%
%    The user enters a formula for f(x), the values of a and b,
%    and the number of Gauss points n.
%
%    The program can be invoked by a function call, in which case the
%    string specifying f(x) must be quoted:
%      q = quad_gauss ( 'x^2', -1, 3, 2 )
%    or, if called with no arguments, it will request them:
%      q = quad_gauss ( );
%      Enter function formula, like x^2: x^2
%      Enter left limit, a: -1
%      Enter right limit, b: 3
%      Enter number of Gauss 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:
%
%    29 June 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    string f_string: the function whose integral is desired.
%
%    real a, b: the interval.
%
%    integer n: the number of Gauss points to use.
%
%  Output:
%
%    real q: the integral estimate.
%

%
%  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 Gauss points n: ' );
  elseif ( ischar ( n ) )
    n = str2num ( n );
  end
%
%  Get the N-point Gauss rule for [-1,+1].
%
  [ x, w ] = legendre_ek_compute ( n );
%
%  Convert rule to [a,b].
%
  x = ( ( 1.0 - x ) * a + ( x + 1.0 ) * b ) / 2.0;
  w = w * ( b - a ) / 2.0;
%
%  Evaluate the integrand.
%
  fx = f ( x );
%
%  Apply the rule.
%
  fx = fx(:);
  q = w' * fx;

  fprintf ( 1, '  quad_gauss: n = %d, integral estimate q = %24.16g\n', n, q );

  return
end
function [ x, w ] = legendre_ek_compute ( n )

%*****************************************************************************80
%
%% legendre_ek_compute(): Legendre quadrature rule by the Elhay-Kautsky method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    19 April 2011
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    This version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%  Input:
%
%    integer N, the order.
%
%  Output:
%
%    real X(N), the abscissas.
%
%    real W(N), the weights.
%

%
%  Define the zero-th moment.
%
  zemu = 2.0;
%
%  Define the Jacobi matrix.
%
  bj = zeros ( n, 1 );

  for i = 1 : n
    bj(i) = ( i * i ) / ( 4 * i * i - 1 );
  end
  bj(1:n) = sqrt ( bj(1:n) );

  x = zeros ( n, 1 );

  w = zeros ( n, 1 );
  w(1) = sqrt ( zemu );
%
%  Diagonalize the Jacobi matrix.
%
  [ x, w ] = imtqlx ( n, x, bj, w );

  w(1:n) = w(1:n).^2;

  return
end
function [ d, z ] = imtqlx ( n, d, e, z )

%*****************************************************************************80
%
%% imtqlx() diagonalizes a symmetric tridiagonal matrix.
%
%  Discussion:
%
%    This routine is a slightly modified version of the EISPACK routine to
%    perform the implicit QL algorithm on a symmetric tridiagonal matrix.
%
%    The authors thank the authors of EISPACK for permission to use this
%    routine.
%
%    It has been modified to produce the product Q' * Z, where Z is an input
%    vector and Q is the orthogonal matrix diagonalizing the input matrix.
%    The changes consist (essentially) of applying the orthogonal 
%    transformations directly to Z as they are generated.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 February 2013
%
%  Author:
%
%    Original FORTRAN77 version by Sylvan Elhay, Jaroslav Kautsky.
%    This version by John Burkardt.
%
%  Reference:
%
%    Sylvan Elhay, Jaroslav Kautsky,
%    Algorithm 655: IQPACK, FORTRAN Subroutines for the Weights of
%    Interpolatory Quadrature,
%    ACM Transactions on Mathematical Software,
%    Volume 13, Number 4, December 1987, pages 399-415.
%
%    Roger Martin, James Wilkinson,
%    The Implicit QL Algorithm,
%    Numerische Mathematik,
%    Volume 12, Number 5, December 1968, pages 377-383.
%
%  Input:
%
%    integer N, the order of the matrix.
%
%    real D(N), the diagonal entries of the matrix.
%
%    real E(N), the subdiagonal entries of the
%    matrix, in entries E(1) through E(N-1). 
%
%    real Z(N), a vector to be operated on.
%
%  Output:
%
%    real D(N), the diagonal entries of the diagonalized matrix.
%
%    real Z(N), the value of Q' * Z, where Q is the matrix that 
%    diagonalizes the input symmetric tridiagonal matrix.
%
  itn = 30;

  prec = eps;

  if ( n == 1 )
    return
  end

  e(n) = 0.0;

  for l = 1 : n

    j = 0;

    while ( true )

      for m = l : n

        if ( m == n )
          break
        end

        if ( abs ( e(m) ) <= prec * ( abs ( d(m) ) + abs ( d(m+1) ) ) )
          break
        end

      end

      p = d(l);

      if ( m == l )
        break
      end

      if ( j == itn )
        fprintf ( 1, '\n' );
        fprintf ( 1, 'IMTQLX - Fatal error!\n' );
        fprintf ( 1, '  Iteration limit exceeded.\n' );
        error ( 'IMTQLX - Fatal error!' );
      end

      j = j + 1;
      g = ( d(l+1) - p ) / ( 2.0 * e(l) );
      r = sqrt ( g * g + 1.0 );

      if ( g < 0.0 )
        t = g - r;
      else
        t = g + r;
      end

      g = d(m) - p + e(l) / ( g + t );
 
      s = 1.0;
      c = 1.0;
      p = 0.0;
      mml = m - l;

      for ii = 1 : mml

        i = m - ii;
        f = s * e(i);
        b = c * e(i);

        if ( abs ( g ) <= abs ( f ) )
          c = g / f;
          r =  sqrt ( c * c + 1.0 );
          e(i+1) = f * r;
          s = 1.0 / r;
          c = c * s;
        else
          s = f / g;
          r =  sqrt ( s * s + 1.0 );
          e(i+1) = g * r;
          c = 1.0 / r;
          s = s * c;
        end

        g = d(i+1) - p;
        r = ( d(i) - g ) * s + 2.0 * c * b;
        p = s * r;
        d(i+1) = g + p;
        g = c * r - b;
        f = z(i+1);
        z(i+1) = s * z(i) + c * f;
        z(i) = c * z(i) - s * f;

      end

      d(l) = d(l) - p;
      e(l) = g;
      e(m) = 0.0;

    end

  end

  for ii = 2 : n

     i = ii - 1;
     k = i;
     p = d(i);

     for j = ii : n
       if ( d(j) < p )
         k = j;
         p = d(j);
       end
     end

     if ( k ~= i )
       d(k) = d(i);
       d(i) = p;
       p = z(i);
       z(i) = z(k);
       z(k) = p;
     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

