import numpy as np
from math import exp

"""
  Numerial differentiation by Richardson extrapolation method
"""

# define function to be differentiated
def func(x):
    return exp(x)

# define analytical deviation of f(x)
def dfdx(x):
    return exp(x)

#===================
# parameters
#===================
h0    = 1.0
nhmax = 15
eps   = 1.0e-6
x0 = 1.0

#===================
# main routine
#===================
print("Numerical differentiation using by Richardson extrapolation method")

fx = func(x0)
f1x = dfdx(x0)
print("")
print("eps for convergnece: {}".format(eps))
print("Analytical values:")
print("  f({})={}".format(x0, fx))
print("  df/dx({})={}".format(x0, f1x))
print("")

# create 2D array by specifing the size to be nhmax * nhmax
D = []
for i in range(nhmax):
  D.append([])
  for j in range(nhmax):
    D[i].append([])

D[0][0] = (func(x0+h0) - func(x0-h0)) / 2.0 / h0
print("D{}{}: {}".format(0, 0, D[0][0]))

for m in range(1, nhmax+1):
    h0 *= 0.5
    D[0][m] = (func(x0+h0) - func(x0-h0)) / 2.0 / h0
    print("D0{}: {}".format(m, D[0][m]))
    for k in range(m-1, -1, -1):
        m1 = m - k
        D[m1][k] = (4.0**m1 * D[m1-1][k+1] - D[m1-1][k]) / (4.0**m1 - 1)
        print("D{}{}: {}".format(m1, k, D[m1][k]))
    diff = D[m][0] - D[m-1][0]
    print("  diff = D[{}][0] - D[{}][0] = {} - {} = {}".format(m, m-1, D[m][0], D[m-1][0], diff))
    if abs(diff) < eps:
        break