#! /usr/bin/env python
#
def max_error_linear ( n, x, u, exact ):

#*****************************************************************************80
#
## MAX_ERROR_LINEAR estimates the max error norm of a finite element solution.
#
#  Location:
#
#    http://people.sc.fsu.edu/~jburkardt/py_src/fem1d_bvp_linear/max_error_linear.py
#
#  Discussion:
#
#    We assume the finite element method has been used, over an interval [A,B]
#    involving N nodes, with piecewise linear elements used for the basis.
#    The coefficients U(1:N) have been computed, and a formula for the
#    exact solution is known.
#
#    This function estimates the max norm of the error:
#
#      MAX_NORM = Integral ( A <= X <= B ) max ( abs ( U(X) - EXACT(X) ) ) dX
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    07 July 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, integer N, the number of nodes.
#
#    Input, real X(N), the mesh points.
#
#    Input, real U(N), the finite element coefficients.
#
#    Input, function EQ = EXACT ( X ), returns the value of the exact
#    solution at the point X.
#
#    Output, real VALUE, the estimated max norm of the error.
#
  import numpy as np

  quad_num = 8
  value = 0.0
#
#  Examine QUAD_NUM points in each element, including left node but not right.
#
  e_num = n - 1

  for e in range ( 0, e_num ):

    l = e
    xl = x[l]
    ul = u[l]

    r = e + 1
    xr = x[r]
    ur = u[r]

    for q in range ( 0, quad_num ):

      xq = ( float ( quad_num - q ) * xl   \
           + float (            q ) * xr ) \
         /   float ( quad_num     )
#
#  Use the fact that U is a linear combination of piecewise linears.
#
      uq = ( ( xr - xq      ) * ul \
           + (      xq - xl ) * ur ) \
           / ( xr      - xl )

      eq = exact ( xq )

      value = max ( value, abs ( uq - eq ) )
#
#  For completeness, check last node.
#
  xq = x[n-1]
  uq = u[n-1]
  eq = exact ( xq )

  value = max ( value, abs ( uq - eq ) )
#
#  Integral approximation requires multiplication by interval length.
#
  value = value * ( x[n-1] - x[0] )

  return value

def max_error_linear_test ( ):

#*****************************************************************************80
#
## MAX_ERROR_LINEAR_TEST tests MAX_ERROR_LINEAR.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    07 July 2015
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    None
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'MAX_ERROR_LINEAR_TEST:' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  MAX_ERROR_LINEAR computes the maximum absolute approximation error' )
  print ( '  between a function exact(x) and a piecewise linear function u()' )
  print ( '  associated with n mesh points x().'  )
  print ( '' )
  print ( '   N    Max_Error' )
  print ( '' )

  n = 3

  for test in range ( 0, 5 ):

    x_lo = 0.0
    x_hi = np.pi
    x = np.linspace ( x_lo, x_hi, n )
#
#  U is an approximation to sin(x).
#
    u = np.zeros ( n )
    for i in range ( 0, n ):
      u[i] = np.sin ( x[i] )

    e1 = max_error_linear ( n, x, u, np.sin )

    print ( '  %2d    %g' % ( n, e1 ) )

    n = 2 * ( n - 1 ) + 1
#
#  Terminate.
#
  print ( '' )
  print ( 'MAX_ERROR_LINEAR_TEST:' )
  print ( '  Normal end of execution.' )
  return

if ( __name__ == '__main__' ):
  from timestamp import timestamp
  timestamp ( )
  max_error_linear_test ( )
  timestamp ( )