<jupyter_start><jupyter_text>Lecture 8 Numerical Differentiation  The following was adpated from: http://gilgamesh.cheme.cmu.edu/doc/software/jacapo/9-numerics/9.1-numpy/9.2-integration.htmlnumerical-differentiation numpy has a function called numpy.diff() that is similar to the one found in matlab. It calculates the differences between the elements in your list, and returns a list that is one element shorter, which makes it unsuitable for plotting the derivative of a function.there is also a function in scipy called scipy.misc.derivate, which calculates a derivative to arbitrary order, however, it takes as input a function rather than a set of data, which makes its use limited.  Using loops to find finite difference derivatives Loops in python are pretty slow (relatively speaking) but they are usually trivial to understand. In this script we show some simple ways to construct derivative vectors using loops. It is implied in these formulas that the data points are equally spaced.<jupyter_code>import numpy as np
import matplotlib.pylab as plt
from pylab import *
import time
%matplotlib inline

'''
These are quick simple ways to calculate numerical derivatives. They
work well for very smooth data. they are surprisingly fast even up to
10000 points in the vector.
'''

x = np.linspace(0,3,3000)   # steps of 0.001
y = x**2+2*x-5
dy_analytical = 2*x+2
'''
lets use a forward difference method:
that works up until the last point, where there is not
a forward difference to use. there, we use a backward difference.
'''

tf1 = time.time()
dyf = [0.0]*len(x)
for i in range(len(y)-1):
    dyf[i] = (y[i+1] - y[i])/(x[i+1]-x[i])
#set last element by backwards difference
dyf[-1] = (y[-1] - y[-2])/(x[-1] - x[-2])  # backward difference for the last point in the array

print (' Forward difference took %1.16f seconds' % (time.time() - tf1))

'''and now a backwards difference'''
tb1 = time.time()
dyb = [0.0]*len(x)
#set first element by forward difference
dyb[0] = (y[0] - y[1])/(x[0] - x[1])
for i in range(1,len(y)):
    dyb[i] = (y[i] - y[i-1])/(x[i]-x[i-1])

print (' Backward difference took %1.16f seconds' % (time.time() - tb1))

'''and now, a centered formula'''
tc1 = time.time()
dyc = [0.0]*len(x)
dyc[0] = (y[0] - y[1])/(x[0] - x[1])
for i in range(1,len(y)-1):
    dyc[i] = (y[i+1] - y[i-1])/(x[i+1]-x[i-1])
dyc[-1] = (y[-1] - y[-2])/(x[-1] - x[-2])

print (' Centered difference took %1.16f seconds' % (time.time() - tc1))

#plot(x,y)
plt.plot(x,dy_analytical,label='analytical derivative')
plt.plot(x,dyf,'--',label='forward')
plt.plot(x,dyb,'--',label='backward')
plt.plot(x,dyc,'--',label='centered')
plt.title('Quadratic function')

plt.legend(loc='upper left')
#plt.savefig('9.2.1-simple-diffs.png')
plt.show()<jupyter_output>Forward difference took 0.0025918483734131 seconds
 Backward difference took 0.0029859542846680 seconds
 Centered difference took 0.0026650428771973 seconds<jupyter_text>Let's compute the errors of each algorithm<jupyter_code>diff_dyf = (dyf-dy_analytical)/dy_analytical
diff_dyb = (dyb-dy_analytical)/dy_analytical
diff_dyc = (dyc-dy_analytical)/dy_analytical

plt.plot(x[:-1],diff_dyf[:-1],'--',label='forward')
plt.plot(x[1:],diff_dyb[1:],'--',label='backward')
plt.plot(x[1:-1],diff_dyc[1:-1],'--',label='centered')
plt.title('Fractional errors')
plt.xlabel('x')
plt.ylabel('(numerical-analytical)/analytical')

legend(loc='upper right')
plt.show()
plt.plot(x[1:-1],diff_dyc[1:-1],'g--')
'''
These are quick simple ways to calculate numerical derivatives. They
work well for very smooth data. they are surprisingly fast even up to
10000 points in the vector.
'''

x = np.linspace(0,np.pi,300)   # step of 0.01
y = np.sin(x)
dy_analytical = np.cos(x)
'''
lets use a forward difference method:
that works up until the last point, where there is not
a forward difference to use. there, we use a backward difference.
'''

tf1 = time.time()
dyf = [0.0]*len(x)
for i in range(len(y)-1):
    dyf[i] = (y[i+1] - y[i])/(x[i+1]-x[i])
#set last element by backwards difference
dyf[-1] = (y[-1] - y[-2])/(x[-1] - x[-2])

print (' Forward difference took %1.16f seconds' % (time.time() - tf1))

'''and now a backwards difference'''
tb1 = time.time()
dyb = [0.0]*len(x)
#set first element by forward difference
dyb[0] = (y[0] - y[1])/(x[0] - x[1])
for i in range(1,len(y)):
    dyb[i] = (y[i] - y[i-1])/(x[i]-x[i-1])

print (' Backward difference took %1.16f seconds' % (time.time() - tb1))

'''and now, a centered formula'''
tc1 = time.time()
dyc = [0.0]*len(x)
dyc[0] = (y[0] - y[1])/(x[0] - x[1])
for i in range(1,len(y)-1):
    dyc[i] = (y[i+1] - y[i-1])/(x[i+1]-x[i-1])
dyc[-1] = (y[-1] - y[-2])/(x[-1] - x[-2])

print (' Centered difference took %1.16f seconds' % (time.time() - tc1))


# the centered formula is the most accurate formula here
#plot(x,y)
plt.plot(x,dy_analytical,label='analytical derivative')
plt.plot(x,dyf,'--',label='forward')
plt.plot(x,dyb,'--',label='backward')
plt.plot(x,dyc,'--',label='centered')

plt.legend(loc='lower left')
#savefig('9.2.1-simple-diffs.png')
plt.show()<jupyter_output>Forward difference took 0.0003221035003662 seconds
 Backward difference took 0.0003881454467773 seconds
 Centered difference took 0.0006759166717529 seconds<jupyter_text>Let's compute the errors of each algorithm<jupyter_code>diff_dyf = dyf-dy_analytical
diff_dyb = dyb-dy_analytical
diff_dyc = dyc-dy_analytical

plt.plot(x[:-1],diff_dyf[:-1],'--',label='forward')
plt.plot(x[1:],diff_dyb[1:],'--',label='backward')
plt.plot(x[1:-1],diff_dyc[1:-1],'--',label='centered')
plt.title('Differences')
plt.xlabel('x')
plt.ylabel('(numerical-analytical)')

legend(loc='upper right')
plt.show()
plt.plot(x[1:-1],diff_dyc[1:-1],'g--')
'''
These are quick simple ways to calculate numerical derivatives. They
work well for very smooth data. they are surprisingly fast even up to
10000 points in the vector.
'''

x = np.linspace(0.1,3,300)
y = np.exp(x**2)
dy_analytical = 2*x*np.exp(x**2)
'''
lets use a forward difference method:
that works up until the last point, where there is not
a forward difference to use. there, we use a backward difference.
'''

tf1 = time.time()
dyf = [0.0]*len(x)
for i in range(len(y)-1):
    dyf[i] = (y[i+1] - y[i])/(x[i+1]-x[i])
#set last element by backwards difference
dyf[-1] = (y[-1] - y[-2])/(x[-1] - x[-2])

print (' Forward difference took %1.16f seconds' % (time.time() - tf1))

'''and now a backwards difference'''
tb1 = time.time()
dyb = [0.0]*len(x)
#set first element by forward difference
dyb[0] = (y[0] - y[1])/(x[0] - x[1])
for i in range(1,len(y)):
    dyb[i] = (y[i] - y[i-1])/(x[i]-x[i-1])

print (' Backward difference took %1.16f seconds' % (time.time() - tb1))

'''and now, a centered formula'''
tc1 = time.time()
dyc = [0.0]*len(x)
dyc[0] = (y[0] - y[1])/(x[0] - x[1])
for i in range(1,len(y)-1):
    dyc[i] = (y[i+1] - y[i-1])/(x[i+1]-x[i-1])
dyc[-1] = (y[-1] - y[-2])/(x[-1] - x[-2])

print (' Centered difference took %1.16f seconds' % (time.time() - tc1))


#plot(x,y)
plt.plot(x,dy_analytical,label='analytical derivative')
plt.plot(x,dyf,'--',label='forward')
plt.plot(x,dyb,'--',label='backward')
plt.plot(x,dyc,'--',label='centered')
plt.title('Exponential function')

legend(loc='upper left')
#savefig('9.2.1-simple-diffs.png')
#show()<jupyter_output>Forward difference took 0.0004148483276367 seconds
 Backward difference took 0.0003969669342041 seconds
 Centered difference took 0.0006060600280762 seconds<jupyter_text>Plot the fractional error of the derivatives<jupyter_code>r_dyf = dyf/dy_analytical-1
r_dyb = dyb/dy_analytical-1
r_dyc = dyc/dy_analytical-1

plt.plot(x[:-1],r_dyf[:-1],'--',label='forward')
plt.plot(x[1:],r_dyb[1:],'--',label='backward')
plt.plot(x[1:-1],r_dyc[1:-1],'--',label='centered')
plt.title('Fractional errors')
plt.xlabel('x')
plt.ylabel('(numerical-analytical)/analytical')

legend(loc='upper center')
plt.show()

plt.plot(x[1:-1],r_dyc[1:-1],'g--')<jupyter_output><empty_output><jupyter_text>Derivatives in the presense of noise<jupyter_code>x = np.linspace(0,2*np.pi,1000)
y = np.sin(x) + 0.01*np.random.random(size=x.shape)   # units are "Volts", amplitude is 1V, noise is 0.01 V
dy_analytical = np.cos(x)

#2-point formula, forward diff
dyf = [0.0]*len(x)
for i in range(len(y)-1):
    dyf[i] = (y[i+1] - y[i])/(x[i+1]-x[i])
#set last element by backwards difference
dyf[-1] = (y[-1] - y[-2])/(x[-1] - x[-2])

'''
calculate dy by 4-point center differencing using array slices

\frac{y[i-2] - 8y[i-1] + 8[i+1] - y[i+2]}{12h}

y[0] and y[1] must be defined by lower order methods
and y[-1] and y[-2] must be defined by lower order methods
'''

dy = np.zeros(y.shape,float) #we know it will be this size
h = x[1]-x[0] #this assumes the points are evenely spaced!
dy[2:-2] = (y[0:-4] - 8*y[1:-3] + 8*y[3:-1] - y[4:])/(12.*h)

dy[0] = (y[1]-y[0])/(x[1]-x[0])
dy[1] = (y[2]-y[1])/(x[2]-x[1])
dy[-2] = (y[-2] - y[-3])/(x[-2] - x[-3])
dy[-1] = (y[-1] - y[-2])/(x[-1] - x[-2])

fg3 = plt.figure(figsize=(10,8))
plt.plot(x,y,label='function')
plt.plot(x,dy_analytical,label='analytical derivative')
plt.plot(x,dyf,'r-',label='2pt-forward diff')
plt.plot(x,dy,'k--',lw=2,label='4pt-centered diff')

plt.legend(loc='upper center')
# compute a running average of the derivative
Nave = 50    # number = 2*Nave+1

dyave = [np.mean(dy[i-Nave:i+Nave]) for i in range(Nave,len(dy)-Nave)]
plt.plot(x,dy_analytical,label='analytical derivative')
plt.plot(x[Nave:-Nave],dyave,'k--',lw=2,label='running average of 4pt diff')

plt.legend(loc='upper center')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1985726280436608, 'fraction_numerical': 0.052162048698572626, 'mean_word_length': 4.761185006045949, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '112441', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4097, 'n_tokens_neox': 3806, 'n_words': 937, 'path': 'Week10/L08/Lecture08a.ipynb', 'repo_name': 'ds-connectors/Physics-88-SP22', 'size': '9528'}