#! /usr/bin/env python3
#
def back_euler_lam ( LAMBDA, xRange, yInitial, numSteps ):

#*****************************************************************************80
#
## back_euler_lam() uses Euler's backward method on a particular stiff ODE.
#
#  Input:
#    LAMBDA = value of parameter in stiff ODE.
#    xRange = [x1,x2], the solution interval.
#    yInitial = k initial values for y at x1
#    numSteps = number of equally-sized steps to take from x1 to x2
#  Output:
#    x = numSteps+1 x values.
#    y = numSteps+1 rows and k columns, with k-th row containing solution at x[k]
#
  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:
      x[k] = x[k-1] + dx
      y[k,:] = ( y[k-1,:] + dx * LAMBDA * np.sin ( x[k] ) ) / ( 1.0 + dx * LAMBDA )

  return x, y

def back_euler_lam_test ( ):

#*****************************************************************************80
#
## back_euler_lam_test() solves the stiff ODE using the backward Euler method.
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    04 October 2021
#
#  Author:
#
#    John Burkardt
#
  from stiff_ode import stiff_lambda, stiff_solution
  import matplotlib.pyplot as plt
  import numpy as np

  print ( '' )
  print ( 'back_euler_lam_test():' )
  print ( '  Solve the stiff ODE system using back_euler_lam().' )

  LAMBDA = 10000.0
  stiff_lambda ( LAMBDA )

  xRange = np.array ( [ 0.0, 2.0 * np.pi ] )
  yInit = 0.0
  numSteps = 40
  x1, y1 = back_euler_lam ( LAMBDA, xRange, yInit, numSteps )

  x2 = np.linspace ( 0.0, 2.0 * np.pi, 41 )
  y2 = stiff_solution ( x2 )

  plt.plot ( x1, y1, 'ro' )
  plt.plot ( x2, y2, 'b-', linewidth = 2 )
  plt.grid ( True )
  plt.xlabel ( '<--- X --->' )
  plt.ylabel ( '<--- Y(X) --->' )
  plt.legend ( [ 'Backward Euler', 'Exact' ] )
  plt.title ( 'Stiff ODE' )

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

  return

if ( __name__ == '__main__' ):
  back_euler_lam_test ( )