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

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

  print ( "exercise2:" )
#
#  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 ( "             RK2" )
  print ( "  K  Y[2.0]  E[K]" )
  print ( "" )
  for k in range ( 0, 6 ):
    x, y = rk2 ( 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, y[-1,0], e[k] )

  r = np.zeros ( 5 )
  print ( "" )
  print ( "     RK2" )
  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] )
#
#  Compare Euler and RK2 errors
#
  e1 = np.zeros ( 6 )
  e2 = np.zeros ( 6 )

  print ( "" )
  print ( "     Euler  RK2" )
  print ( "  K  Error  Error" )
  print ( "" )
  for k in range ( 0, 6 ):
    x1, y1 = forward_euler ( f_ode, xRange, yInit, numSteps[k] )
    x2, y2 = rk2 ( f_ode, xRange, yInit, numSteps[k] )
    yexact = -2.0 * np.exp ( - x1 ) - 3.0 * x1 + 3.0;
    e1[k] = np.abs ( yexact[-1] - y1[-1,0] )
    e2[k] = np.abs ( yexact[-1] - y2[-1,0] )
    print ( k, e1[k], e2[k] )

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