<jupyter_start><jupyter_text>First exploration of finite-difference approximations to derivativesNote: Prof. Brown's [differentiation notebook](https://github.com/cu-numcomp/numcomp-class/blob/master/Differentiation.ipynb) for CSCI-3656 is quite good, and covers [algorithmic differentiation](https://en.wikipedia.org/wiki/Automatic_differentiation) using JAX.  [![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/cu-numcomp/numcomp-class/blob/master/Differentiation.ipynb)- Task 1: try the simple forward differences formula to approximate the derivative of $f(x)=\sin(x)$. The formula is:$$f'(x) \approx \frac{ f(x+h)-f(x) }{h}$$Try this for a range of stepsizes $h$. Can we make the approximation (aka "truncation") error arbitrarily small?<jupyter_code>import numpy as np
import matplotlib.pyplot as plt

# Approximate the derivative of sin(x)
f   = lambda x : np.sin(x)
trueDerivative = lambda x : np.cos(x)

def twoPointForwardDifferences( f, x, h ):
  """ The simplest, plain vanilla finite difference approx"""
  return ( f(x+h) - f(x) )/h

# Pick some point x at which we want the derivative
x   = 2.
df      = trueDerivative(x)
twoPointF= lambda h : twoPointForwardDifferences( f, x, h)

# Pick stepsizes h
hList   = np.logspace(-1,-12,30)

# Plot
plt.loglog( hList, np.abs( df - twoPointF(hList) ), 'o-', label="2 Pt Forward Diff" );
plt.xlabel("h");
plt.ylabel("Error");<jupyter_output><empty_output><jupyter_text>Getting fancier: adding in other methods- Task 2: try other finite difference methods (e.g, 3 point forward differences; 3 point centered differences; 5 point centered differences)<jupyter_code>import numpy as np
import matplotlib.pyplot as plt

# Approximate the derivative of sin(x)
f   = lambda x : np.sin(x)
trueDerivative = lambda x : np.cos(x)

def twoPointForwardDifferences( f, x, h ):
  """ The simplest, plain vanilla finite difference approx"""
  return ( f(x+h) - f(x) )/h
def threePointForwardDifferences( f, x, h):
  """Eq. (4.4) in the Burden and Faires book"""
  return (-3*f(x)+4*f(x+h)-f(x+2*h))/(2*h)
def threePointCenteredDifferences( f, x, h ):
  """ Book calls it 3 points, since it sort of includes f(x) """
  return ( f(x+h) - f(x-h) )/(2*h)
def fivePointCenteredDifferences( f, x, h ):
  """ Eq. (4.6) in Burden and Faires book """
  return ( f(x-2*h) - 8*f(x-h) + 8*f(x+h) - f(x+2*h) )/(12*h)

# Pick some point x at which we want the derivative
x   = 2.
df      = trueDerivative(x)
twoPointF= lambda h : twoPointForwardDifferences( f, x, h)
threePointC= lambda h : threePointCenteredDifferences( f, x, h)
threePointF= lambda h : threePointForwardDifferences( f, x, h)
fivePointC= lambda h : fivePointCenteredDifferences( f, x, h)

# Pick stepsizes h
hList   = np.logspace(-1,-12,20)

# Plot
plt.figure(figsize=(10,8)) 
plt.loglog( hList, np.abs( df - twoPointF(hList) ), 'o-', label="2 Pt Forward Diff" );
plt.loglog( hList, np.abs( df - threePointC(hList) ), 'o-', label="3 Pt Centered Diff" );
plt.loglog( hList, np.abs( df - threePointF(hList) ), 'o-', label="3 Pt Forward Diff" );
plt.loglog( hList, np.abs( df - fivePointC(hList) ), 'o-', label="5 Pt Centered Diff" )
plt.legend();
plt.xlabel("h");
plt.ylabel("Error");<jupyter_output><empty_output><jupyter_text>Can we estimate how low the error can be?It's a trade off between truncation error (low $h$ is good), and roundoff error (low $h$ is bad)- Task 3: can we predict how low each method gets? Assume error looks like$$\epsilon/h + h^k$$ where $k$ is the order, and use `sympy` to minimize this (e.g., differentiate it and set it equal to 0).  You may want to declare $\epsilon$ as a symbol like:  `e = sym.symbols('\epsilon',positive=True)`<jupyter_code>import sympy as sym
from sympy import init_printing
init_printing()

h = sym.symbols('h')
e = sym.symbols('\epsilon',positive=True)
g = e/h + h
d = sym.diff( g,  h)
d
h0 = list( sym.solveset( d, h, domain=sym.S.Reals) )[0]
print("Order 1, (optimal h, value-at-optimal)")
h0, g.subs(h,h0)
g = e/h + h**2
d = sym.diff( g,  h)
h0 = list( sym.solveset( d, h, domain=sym.S.Reals) )[0]
print("Order 2")
h0, g.subs(h,h0)
g = e/h + h**3
d = sym.diff( g,  h)
h0 = list( sym.solveset( d, h, domain=sym.S.Reals) )[0]
print("Order 3")
h0, g.subs(h,h0)
g = e/h + h**4
d = sym.diff( g,  h)
h0 = list( sym.solveset( d, h, domain=sym.S.Reals) )[0]
print("Order 4")
h0, g.subs(h,h0)<jupyter_output>Order 4<jupyter_text>And put these in the plot- 2 pt forward diff is $O(h)$- 3 pt forward diff is $O(h^2)$- 3 pt centered diff is $O(h^2)$- 5 pt centered diff is $O(h^4)$<jupyter_code># So add to the plot
plt.figure(figsize=(10,8)) 
plt.loglog( hList, np.abs( df - twoPointF(hList) ), 'o-', label="2 Pt Forward Diff" );
plt.loglog( hList, np.abs( df - threePointC(hList) ), 'o-', label="3 Pt Centered Diff" );
plt.loglog( hList, np.abs( df - threePointF(hList) ), 'o-', label="3 Pt Forward Diff" );
plt.loglog( hList, np.abs( df - fivePointC(hList) ), 'o-', label="5 Pt Centered Diff" )
plt.legend();
plt.xlabel("h");
plt.ylabel("Error");

# And can we predict the values of h which gives us the smallest amount?
eps = np.finfo(float).eps
eps**(1/2)
eps**(1/4)
plt.vlines( eps**(1/2), *plt.gca().get_ylim() ,linestyles='--',colors="C0");
plt.vlines( eps**(1/3), *plt.gca().get_ylim() ,linestyles='--',colors="C1");
plt.vlines( eps**(1/5), *plt.gca().get_ylim() ,linestyles='--',colors="C3");

plt.hlines( eps**(1/2), *plt.gca().get_xlim() ,linestyles='--',colors="C0");
plt.hlines( eps**(2/3), *plt.gca().get_xlim() ,linestyles='--',colors="C1");
plt.hlines( eps**(4/5), *plt.gca().get_xlim() ,linestyles='--',colors="C3");<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1779541446208113, 'fraction_numerical': 0.018871252204585538, 'mean_word_length': 4.676676676676677, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 16, 'https://': 4, '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': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '420370', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2302, 'n_tokens_neox': 2179, 'n_words': 625, 'path': 'Demos/Ch4_FiniteDifferences.ipynb', 'repo_name': 'skhadem/numerical-analysis-class', 'size': '5670'}