#! /usr/bin/env python3
#
def midpoint ( f_ode, xRange, yInitial, numSteps ):

#*****************************************************************************80
#
## midpoint() applies the midpoint method to an ODE.
#
#  Input:
#    f_ode(x,y) evaluates the right hand side;
#    xRange is a two dimensional vector of beginning and final values for x
#    yInitial is a column vector for the initial value of y
#    numSteps is the number of evenly-spaced steps to divide up the interval xRange
#  Output:
#    x returns the values of the independent variable.
#    y returns the values of the dependent variable.
#
  from scipy.optimize import fsolve
  import numpy as np

  x = np.zeros ( numSteps + 1 )
  y = np.zeros ( ( numSteps + 1, np.size ( yInitial ) ) )

  dx = ( xRange[1] - xRange[0] ) / numSteps

  for k in range ( 0, numSteps + 1 ):

    if ( k == 0 ):
      x[0] = xRange[0]
      y[0,:] = yInitial
    else:

      xm = x[k-1] + 0.5 * dx
      ym = y[k-1,:] + 0.5 * dx * f_ode ( x[k-1], y[k-1,:] )
      ym = fsolve ( bef, ym, args = ( f_ode, x[k-1], y[k-1,:], xm ) )

      x[k] = x[k-1] + dx
      y[k,:] = 2.0 * ym - y[k-1,:] 

  return x, y

def bef ( yp, f, to, yo, tp ):

#*****************************************************************************80
#
## bef() evaluates the backward Euler residual.
#
#  Discussion:
#
#    We are seeking a value YP defined by the implicit equation:
#
#      YP = YO + ( TP - TO ) * F ( TP, YP )
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    20 October 2020
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real yp: the estimated solution value at the new time.
#
#    function f: evaluates the right hand side of the ODE.  
#
#    real to, yo: the old time and solution value.
#
#    real tp: the new time.
#
#  Output:
#
#    real value: the residual.
#
  value = yp - yo - ( tp - to ) * f ( tp, yp )

  return value

def midpoint_test ( ):

#*****************************************************************************80
#
## midpoint_test() tests midpoint().
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    20 October 2020
#
#  Author:
#
#    John Burkardt
#
  import numpy as np
  import platform

  print ( '' )
  print ( 'midpoint_test():' )
  print ( '  Python version: %s' % ( platform.python_version ( ) ) )
  print ( '  Test midpoint.' )

  tspan = np.array ( [ 0.0, 2.0 ] )
  y0 = 5.1765
  n = 100
  humps_midpoint ( tspan, y0, n )

  tspan = np.array ( [ 0.0, 5.0 ] )
  y0 = np.array ( [ 5000, 100 ] )
  n = 200
  predator_prey_midpoint ( tspan, y0, n )
#
#  Terminate.
#
  print ( '' )
  print ( 'midpoint_test:' )
  print ( '  Normal end of execution.' )
  return

def humps_midpoint ( tspan, y0, n ):

#*****************************************************************************80
#
## humps_midpoint(): humps ODE using midpoint.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    20 October 2020
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real tspan[2]: the time span
#
#    real y0[2]: the initial condition.
#
#    integer n: the number of steps to take.
#
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'humps_midpoint():' )
  print ( '  Solve the humps ODE system using midpoint().' )

  t, y = midpoint ( humps_deriv, tspan, y0, n )

  plt.clf ( )

  plt.plot ( t, y, 'r-', linewidth = 2 )

  a = tspan[0]
  b = tspan[1]
  if ( a <= 0.0 and 0.0 <= b ):
    plt.plot ( [a,b], [0,0], 'k-', linewidth = 2 )

  ymin = min ( y )
  ymax = max ( y )
  if ( ymin <= 0.0 and 0.0 <= ymax ):
    plt.plot ( [0, 0], [ymin,ymax], 'k-', linewidth = 2 )

  plt.grid ( True )
  plt.xlabel ( '<--- T --->' )
  plt.ylabel ( '<--- Y(T) --->' )
  plt.title ( 'humps ODE solved by midpoint()' )

  filename = 'humps_midpoint.jpg'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.show ( block = False )
  plt.close ( )

  return

def humps_deriv ( x, y ):

#*****************************************************************************80
#
## humps_deriv() evaluates the 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
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 April 2020
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real x[:], y[:]: the arguments.
#
#  Output:
#
#    real yp[:]: the value of the derivative at x.
#
  yp = - 2.0 * ( x - 0.3 ) / ( ( x - 0.3 )**2 + 0.01 )**2 \
       - 2.0 * ( x - 0.9 ) / ( ( x - 0.9 )**2 + 0.04 )**2

  return yp

def predator_prey_midpoint ( tspan, y0, n ):

#*****************************************************************************80
#
## predator_prey_midpoint(): predator ODE using midpoint_fsolve.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    20 October 2020
#
#  Author:
#
#    John Burkardt
#
#  Input:
#
#    real tspan[2]: the time span
#
#    real y0[2]: the initial condition.
#
#    integer n: the number of steps to take.
#
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'predator_prey_midpoint():' )
  print ( '  Solve the predator prey ODE system using midpoint().' )

  t, y = midpoint ( predator_prey_deriv, tspan, y0, n )

  plt.clf ( )

  plt.plot ( y[:,0], y[:,1], 'r-', linewidth = 2 )

  plt.grid ( True )
  plt.xlabel ( '<--- Prey --->' )
  plt.ylabel ( '<--- Predators --->' )
  plt.title ( 'predator prey ODE solved by midpoint' )

  filename = 'predator_prey_midpoint.jpg'
  plt.savefig ( filename )
  print ( '  Graphics saved as "%s"' % ( filename ) )
  plt.show ( block = False )
  plt.close ( )

  return

def predator_prey_deriv ( t, rf ):

#*****************************************************************************80
#
## predator_prey_deriv() evaluates the right hand side of the system.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    22 April 2020
#
#  Author:
#
#    John Burkardt
#
#  Reference:
#
#    George Lindfield, John Penny,
#    Numerical Methods Using MATLAB,
#    Second Edition,
#    Prentice Hall, 1999,
#    ISBN: 0-13-012641-1,
#    LC: QA297.P45.
#
#  Input:
#
#    real T, the current time.
#
#    real RF[2], the current solution variables, rabbits and foxes.
#
#  Output:
#
#    real DRFDT[2], the right hand side of the 2 ODE's.
#
  import numpy as np

  r = rf[0]
  f = rf[1]

  drdt =    2.0 * r - 0.001 * r * f
  dfdt = - 10.0 * f + 0.002 * r * f

  drfdt = np.array ( [ drdt, dfdt ] )

  return drfdt

def timestamp ( ):

#*****************************************************************************80
#
## timestamp() prints the date as a timestamp.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license. 
#
#  Modified:
#
#    21 August 2019
#
#  Author:
#
#    John Burkardt
#
  import time

  t = time.time ( )
  print ( time.ctime ( t ) )

  return

if ( __name__ == '__main__' ):
  timestamp ( )
  midpoint_test ( )
  timestamp ( )