#! /usr/bin/env python
"""
 CentralDiff.py is an example application of the central difference
 to numerically calculate the derivative of a function f(x) = x^3 sin(5x).

      Central Difference == Dc[f(x), h] = (f(x+0.5*h) - f(x-0.5*h)) / h

 The program compares analytical and numerical derivatives by graphing 
 results in the x range (1.0: 3.0).   The program then graphically explores 
 the values of the central difference as a function of h. 

 Paul Eugenio
 Florida State University
 Mar 22, 2019
"""

from __future__ import division, print_function

import numpy as np
import matplotlib.pyplot as plt

# use latex rendering in lables (on MacOS comment out the line below)
plt.rc('text', usetex=True)



def numerical_dervative(func, x, h=1e-8):
    """ 
     Numerical derivative using Central Difference

     Central Difference == Dc[f(x), h] = (f(x+0.5*h) - f(x-0.5*h)) / h

     Usage: dFuncDx = numerical_dervative(func, x)  
     ## use default value of h for optimal result.
    """
    return (func(x + 0.5*h) - func(x - 0.5*h)) / h





#
# main program
#

# program specific functions

def myFunction(x):
    """ f(x) = x^3 sin(5x) """
    return x**3 * np.sin(5*x)

dFdxNumerical = lambda x : numerical_dervative(myFunction, x)


def dFdx(x):
    """ 
     analytic derivative 
        of f(x)  -> df/dx =  3x^3 sin(5x) + 5x^3 cos(5x)
    """
    return 3*x**2 * np.sin(5*x) + 5*x**3 * np.cos(5*x)


# generatting plot data
xLow, xHigh = 1.0, 3.0
nPoints = 100

x = np.linspace(xLow,xHigh, nPoints)
functionPoints =   myFunction(x)
#cDiffPoints =  centralDiff(myFunction, x)       # central difference (numerical derivative)
cDiffPoints =  dFdxNumerical(x)       # could also have used =numerical_dervative(myFunction, x) rather than the lambda function
dFdxPoints =  dFdx(x)    

# plotting F(x) & dFdx as a continuous line & plotting the
# central difference of F(x) as triangular points
fig1 = plt.figure(1)
plt.plot(x, functionPoints, "-")
plt.plot(x, dFdxPoints, "-")
plt.plot(x, cDiffPoints, "v")

# dressing up plot
plt.grid(True)
plt.title("Central Difference")
plt.xlabel("x")
plt.legend([r"$f(x) = x^3 sin(5x)$", r"$\frac{df}{dx}$",r"$D^c _h [f(x)]$"])
fig1.savefig("centralDiff-a.png")
fig1.show()


#
# Now fixing x=2 and calculating the central difference of F(x=2) as
# a function of h where Dc[F(x)] = (f(x+0.5h) - f(x-0.5h)) / h
#

# initial values
x = 2
h = 1

hPoints, cDiffPoints = [], []            # container for plotting points 

# calculate points h vs Dc[f(x),h]
for i in range(20):
    h *= 0.15
    cDiffPoints += [numerical_dervative(myFunction, x, h)]   

# generate second figure and dress up plot
fig2 = plt.figure(2)
plt.semilogx(hPoints, cDiffPoints, "o-")
plt.xlabel("h")
plt.ylabel(r"$D^c _h [f(x=2)]$")
fig2.savefig("centralDiff-b.png")
fig2.show()


input("Enter [return] to exit")