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

  from expm_ode import expm_ode
  from forward_euler import forward_euler
  import numpy as np

  print ( "exercise1:" )
#
#  Solve the ODE.
#
  f_ode = expm_ode
  xRange = np.array ( [ 0.0, 2.0 ] )
  yInit = 1.0
  numSteps = np.array ( [ 10, 20, 40, 80, 160, 320 ] )

  e = np.zeros ( 6 )

  print ( "" )
  print ( "  K  E[K]" )
  print ( "" )
  for k in range ( 0, 6 ):
    x, y = forward_euler ( f_ode, xRange, yInit, numSteps[k] )
    yexact = -2.0 * np.exp ( - x ) - 3.0 * x + 3.0;
    e[k] = np.abs ( yexact[-1] - y[-1,0] )
    print ( k, e[k] )

  r = np.zeros ( 5 )
  print ( "" )
  print ( "  K  R[K]=E[K]/E[K+1]" )
  print ( "" )
  for k in range ( 0, 5 ):
    r[k] = e[k] / e[k+1]
    print ( k, r[k] )

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