#!/usr/bin/env python
# coding: utf-8

# In[1]:


#centered finite difference for first derivative
#second-order accurate
def fd(f, x, h):
    return (f(x+h) - f(x-h)) / (2*h)


# In[2]:


from math import sqrt
f = lambda x: sqrt(x)
fd(f, 1, 0.1)


# In[3]:


fd(f, 1, 0.01) #decreasing h 10fold decreases truncation error 100fold


# In[4]:


from numpy import sin, cos
x = 1
f = lambda x: sin(x)
df = lambda x: cos(x)
h = 0.1
fd(f, x, h), df(x)


# In[6]:


from numpy import logspace, empty
n = 40
hs = logspace(-15, 0, num=n)
actual = df(x)
errs = empty(n)
for i in range(n):
    estimate = fd(f, x, hs[i])
    errs[i] = abs(estimate - actual) / abs(actual) #fractional error


# In[9]:


from matplotlib.pyplot import plot, loglog, xlabel, ylabel
loglog(hs, errs)
xlabel('step size h')
ylabel('fractional error')


# In[18]:


#Richardson extrapolation
from numpy import zeros, nan
def richardson(f, x, h0, n):
    D = zeros((n, n))
    D[:] = nan
    for i in range(n):
        h = h0 / 2**i
        D[i, 0] = fd(f, x, h)
        for j in range(1, n):
            D[i-j, j] = (4**j/(4**j-1))*D[i-j+1, j-1] - (1/(4**j-1))*D[i-j, j-1]
    return D   


# In[20]:


from numpy import exp
f = lambda x: exp(x)
x = 1
h0 = 0.5
n = 4
D = richardson(f, x, h0, n)


# In[21]:


abs(D - f(x)) / abs(f(x)) #fractional error versus the analytic derivative of exp()


# In[30]:


#estimate derivative based on polynomial interpolation
points = [0.9, 0.95, 1, 1.05, 1.1]
y = f(points)
from numpy import polyfit, polyder, polyval
p = polyfit(points, y, 4) #degree (4) should be one less than the number of points (5)
p_der = polyder(p)
df_est = polyval(p_der, x)


# In[28]:


p, p_der


# In[31]:


df_est


# In[ ]: