<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 A primer on numerical differentiation
In order to numerically evaluate a derivative $y'(x)=dy/dx$ at point $x_0$, we approximate is by using finite differences
Step1:  Why is it that the sequence does not converge? This is due to the round-off errors in the representation of the floating point numbers. To see this, we can simply type
Step2:  Let's try using powers of 1/2
Step3:  In addition, one could consider the midpoint difference, defined as
Step4:  A more in-depth discussion about round-off errors in numerical differentiation can be found <a href="http
Step5:  Notice above that gradient() uses forward and backward differences at the two ends.
Step6:  More discussion about numerical differenciation, including higher order methods with error extrapolation can be found <a href="http
Step7:  One way to improve the roundoff errors is by simply using the decimal package
<ASSISTANT_TASK:>
Python Code:
dx = 1.
x = 1.
while(dx > 1.e-10):
    dy = (x+dx)*(x+dx)-x*x
    d = dy / dx
    print("%6.0e %20.16f %20.16f" % (dx, d, d-2.))
    dx = dx / 10.
    
Explanation: A primer on numerical differentiation
In order to numerically evaluate a derivative $y'(x)=dy/dx$ at point $x_0$, we approximate is by using finite differences:
Therefore we find: $$\begin{eqnarray}
&& dx \approx \Delta x &=&x_1-x_0, \
&& dy \approx \Delta y &=&y_1-y_0 = y(x_1)-y(x_0) = y(x_0+\Delta_x)-y(x_0),\end{eqnarray}$$
Then we re-write the derivative in terms of discrete differences as:
$$\frac{dy}{dx} \approx \frac{\Delta y}{\Delta x}$$
Example
Let's look at the accuracy of this approximation in terms of the interval $\Delta x$. In our first example we will evaluate the derivative of $y=x^2$ at $x=1$.
End of explanation
((1.+0.0001)*(1+0.0001)-1)
Explanation: Why is it that the sequence does not converge? This is due to the round-off errors in the representation of the floating point numbers. To see this, we can simply type:
End of explanation
dx = 1.
x = 1.
while(dx > 1.e-10):
    dy = (x+dx)*(x+dx)-x*x
    d = dy / dx
    print("%8.5e %20.16f %20.16f" % (dx, d, d-2.))
    dx = dx / 2.
Explanation: Let's try using powers of 1/2
End of explanation
from math import sin, sqrt, pi
dx = 1.
while(dx > 1.e-10):
    x = pi/4.
    d1 = sin(x+dx) - sin(x); #forward
    d2 = sin(x+dx*0.5) - sin(x-dx*0.5); # midpoint
    d1 = d1 / dx;
    d2 = d2 / dx;
    print("%8.5e %20.16f %20.16f %20.16f %20.16f" % (dx, d1, d1-sqrt(2.)/2., d2, d2-sqrt(2.)/2.) )
    dx = dx / 2.
Explanation: In addition, one could consider the midpoint difference, defined as:
$$ dy \approx \Delta y = y(x_0+\frac{\Delta_x}{2})-y(x_0-\frac{\Delta_x}{2}).$$
For a more complex function we need to import it from the math module. For instance, let's calculate the derivative of $sin(x)$ at $x=\pi/4$, including both the forward and midpoint differences.
End of explanation
%matplotlib inline
import numpy as np
from matplotlib import pyplot
y = lambda x: x*x
x1 = np.arange(0,10,1)
x2 = np.arange(0,10,0.1)
y1 = np.gradient(y(x1), 1.)
print (y1)
pyplot.plot(x1,np.gradient(y(x1),1.),'r--o');
pyplot.plot(x1[:x1.size-1],np.diff(y(x1))/np.diff(x1),'b--x');
Explanation: A more in-depth discussion about round-off errors in numerical differentiation can be found <a href="http://www.uio.no/studier/emner/matnat/math/MAT-INF1100/h10/kompendiet/kap11.pdf">here</a>
Special functions in numpy
numpy provides a simple method diff() to calculate the numerical derivatives of a dataset stored in an array by forward differences. The function gradient() will calculate the derivatives by midpoint (or central) difference, that provides a more accurate result.
End of explanation
pyplot.plot(x2,np.gradient(y(x2),0.1),'b--o');
Explanation: Notice above that gradient() uses forward and backward differences at the two ends.
End of explanation
from scipy.misc import derivative
y = lambda x: x**2
dx = 1.
x = 1.
while(dx > 1.e-10):
    d = derivative(y, x, dx, n=1, order=3)
    print("%6.0e %20.16f %20.16f" % (dx, d, d-2.))
    dx = dx / 10.
Explanation: More discussion about numerical differenciation, including higher order methods with error extrapolation can be found <a href="http://young.physics.ucsc.edu/115/diff.pdf">here</a>. 
The module scipy also includes methods to accurately calculate derivatives:
End of explanation
from decimal import Decimal
dx = Decimal("1.")
while(dx >= Decimal("1.e-10")):
    x = Decimal("1.")
    dy = (x+dx)*(x+dx)-x*x
    d = dy / dx
    print("%6.0e %20.16f %20.16f" % (dx, d, d-Decimal("2.")))
    dx = dx / Decimal("10.")
Explanation: One way to improve the roundoff errors is by simply using the decimal package
End of explanation
<END_TASK>