#! /usr/bin/env python3
#
def bvp_test3 ( ):

  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'bvp_test3():' )
  print ( '  Solve a boundary value problem for the humps function:' )
  print ( '  u'' = humps_deriv2(x)' )
  print ( '  0 <= x <= 2' )
  print ( '  u(0) = humps_fun(0)  u(2) = hump_fun(2)' )
  print ( '  Exact solution is u(x)=1.0/((x-0.3)^2+0.01 )+1.0/((x-0.9)^2 + 0.04 )-6.0' )

  a = 0.0
  b = 2.0
  n = 11
  x = np.linspace ( a, b, n )
  h = ( b - a ) / ( n - 1 )

  A = np.zeros ( [ n, n ] )
  A[0,0] = 1.0
  A[n-1,n-1] = 1.0
  for i in range ( 1, n-1 ):
    A[i,i-1] = 1.0
    A[i,i] = - 2.0
    A[i,i+1] = 1.0

  rhs = np.zeros ( n )
  rhs[0] = humps_fun ( a )
  rhs[n-1] = humps_fun ( b )
  for i in range ( 1, n-1 ):
    rhs[i] = humps_deriv2 ( x[i] ) * h**2

  u = np.linalg.solve ( A, rhs )
  v = humps_fun ( x )

  table = np.stack ( [ x, u, v ], axis = 1 )
  print ( '' )
  print ( '  X     U(X)  Uexact' )
  print ( table )

  plt.clf ( )
  plt.plot ( x, u, 'bo', linewidth = 2 )
  plt.plot ( x, v, 'r-', linewidth = 2 )
  plt.grid ( True )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- U --->' )
  plt.title ( 'u'' = humps_deriv2(x)' )
  plt.legend ( [ '11 point approximation', 'Exact solution' ] )
  filename = 'bvp_test3.png'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.show ( block = False )
  plt.close ( )
#
#  Terminate.
#
  print ( '' )
  print ( 'bvp_test3():' )
  print ( '  Normal end of execution.' )

  return

def humps_fun ( x ):

#*****************************************************************************80
#
## humps_fun() evaluates the humps function.
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    06 August 2019
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real x(): the evaluation points.
#
#  Output:
#
#    real y(): the function values.
#
  y = 1.0 / ( ( x - 0.3 )**2 + 0.01 ) \
    + 1.0 / ( ( x - 0.9 )**2 + 0.04 ) \
    - 6.0

  return y

def humps_deriv2 ( x ):

#*****************************************************************************80
#
## humps_deriv2() evaluates the second derivative of the humps function.
#
#  Discussion:
#
#    y = 1.0 / ( ( x - 0.3 )**2 + 0.01 ) \
#      + 1.0 / ( ( x - 0.9 )**2 + 0.04 ) \
#      - 6.0
#
#    ypp = - 2.0 * ( x - 0.3 ) / ( ( x - 0.3 )**2 + 0.01 )**2 \
#          - 2.0 * ( x - 0.9 ) / ( ( x - 0.9 )**2 + 0.04 )**2
#
#  Licensing:
#
#    This code is distributed under the MIT license.
#
#  Modified:
#
#    06 August 2019
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real x(): the arguments.
#
#  Output:
#
#    real ypp(): the value of the second derivative at x.
#
  u1 = - 2.0 * ( x - 0.3 )
  v1 = ( ( x - 0.3 )**2 + 0.01 )**2
  u2 = - 2.0 * ( x - 0.9 )
  v2 = ( ( x - 0.9 )**2 + 0.04 )**2

  u1p = - 2.0
  v1p = 2.0 * ( ( x - 0.3 )**2 + 0.01 ) * 2.0 * ( x - 0.3 ) 
  u2p = - 2.0
  v2p = 2.0 * ( ( x - 0.9 )**2 + 0.04 ) * 2.0 * ( x - 0.9 )

  ypp = ( u1p * v1 - u1 * v1p ) / v1**2 \
      + ( u2p * v2 - u2 * v2p ) / v2**2

  return ypp

if ( __name__ == "__main__" ):
  bvp_test3 ( )