### Review of Numerical Differentiation

#### Q. What orders are the errors?



In [ ]:

def forward_difference(f, x, h):
return (f(x + h) - f(x)) / h

def central_difference(f, x, h):
return (f(x + h) - f(x - h)) / (2. * h)



# Integration: Chalk Board

## Numerical Integration

Trapezoidal integration based on our notes:


In [2]:

def trapezoidal(f, x_L, x_U, n):
"""
Integrate function f from x_L to x_U using n points.
"""

# Step size (or grid spacing):
h = (x_U - x_L) / float(n)  # h is dx in our notes

# Include start and end points in summation
tot = f(x_L) + f(x_U)

# Loop over all the trapezoids:
# We don't include k = 0 (x = x_L) or k = n-1 (x = x_U) because
# they're included in the sum before the loop.
for k in range(1, n, 1):
x = x_L + (k * h)
tot += 2. * f(x)

tot *= h/2.




In [3]:

%matplotlib inline




In [4]:

import numpy as np
import matplotlib.pyplot as pl



Test case:

$$g(x) = \tau A \cos(\tau x)$$

with $\tau=2\pi$



In [7]:

tau= 2*np.pi
g = lambda x, A=1: tau*A*np.cos(tau*x)




In [8]:

x = np.linspace(0, 2, 100)

pl.plot(x, g(x), color='k')
pl.xlabel('x')
pl.ylabel('g(x)')




Out[8]:

<matplotlib.text.Text at 0x1187143c8>



#### Q. What should we get if we integrate from 0 to 1?

We can verify this analytically:

$\int_{x_i}^{x_f} g(x) dx = \tau A \int_{x_i}^{x_f} \cos(\tau x) dx$

Substitute $u = \tau x$, in which case $du = \tau dx$,

$$\int_{x_i}^{x_f} g(x) dx = A \int_{\tau x_i}^{\tau x_f} \cos(u) du$$

Finally,

$$\int_{x_i}^{x_f} g(x) dx = A \left[\sin(u) \right]_{\tau x_i}^{\tau x_f} du = A \left[\sin(\tau x_f) - sin(\tau x_i) \right]$$

For $x_i = 0$, $x_f = 1$,

$$\int_0^1 g(x) dx = A \left[\sin(\tau) - sin(0) \right] = 0$$
Now, in Python:


In [9]:

from math import sin, pi

int_analytic = lambda A, x_i, x_f:  A * (sin(tau*x_f) - sin(tau*x_i))




In [13]:

print('Integral of g(x) from x_i=%g to x_f=%g is %g ' \
% (0., 0.2, int_analytic(1., 0., 0.2)))




Integral of g(x) from x_i=0 to x_f=0.2 is 0.951057


Now, let's calculate integral using the trapezoidal rule:


In [14]:

x_L = 0.
x_U = 0.2
n = 5
print('Approximate integral is %g' % trapezoidal(g, x_L, x_U, n))




Approximate integral is 0.946045



How about from $x_i=0$ to $x_f=0.2$?

#### Q. Before we do this, do you expect the trapezoidal sum version of the integral to be an overestimate or an underestimate? Hint: see figure below.



In [15]:

A = 1.0
pl.figure()
pl.plot(x, g(x, A=A), color='k')
pl.xlabel('x')
pl.ylabel('g(x)')
pl.xlim(0, 0.2)




Out[15]:

(0, 0.2)




In [16]:

print('Analytic integral is %g' % int_analytic(A, 0.0, 0.2))
print('Approximate integral is %g' \
% trapezoidal(lambda x: g(x, A=A), 0.0, 0.2, 5))




Analytic integral is 0.951057
Approximate integral is 0.946045



#### Q. How can we improve the accuracy of this numerical integral?

Plot error as a function of n up to n = 100


In [17]:

n = 100        # Max number of grid points

x_L = 0.0      # limits of integration
x_U = 0.2

# Grid point and integral approximation arrays
n_grid = np.arange(1, n)  # Starts with 1:  need at least one trapezoid.
approx = np.zeros(n-1)

# Loop over all n starting with n = 1 (n_grid[0] = 1):
for i in range(n-1):
integral = trapezoidal(g, x_L, x_U, n_grid[i])
approx[i] = integral

# Compute relative error in percent
error = 100 * (approx - int_analytic(1., x_L, x_U)) / int_analytic(1., x_L, x_U)

# Plot results
pl.figure()
pl.plot(n_grid, error)
pl.xlabel('Number of Grid Points')
pl.ylabel('Percent Error')




Out[17]:

<matplotlib.text.Text at 0x118a0c9e8>



#### Expanding the y-axis using logarithm



In [18]:

error_abs = np.abs(error)  # absolute value of the error

pl.semilogy(n_grid, error_abs)
pl.xlabel('Number of Grid Points')
pl.ylabel('Percent Relative Error')

# Plot straight lines corresponding to 1% and 10% error
pl.axhline(y=1, color='k', ls='--')
pl.axhline(y=1e-1, color='k', ls=':')
pl.axhline(y=1e-2, color='k', ls='-.')




Out[18]:

<matplotlib.lines.Line2D at 0x118c05c88>



#### Q. How could we find the number of grid points required to achieve 0.01% accuracy more precisely than the visual approach above?



In [19]:

# location of point closest to 0.01% error
i = np.abs(error_abs - 0.01).argmin() # This yields the index of the point
print(np.abs(error_abs - 0.01).min()) # This prints the point closest to 0.01%.
pl.semilogy(np.abs(error_abs - 0.01))




0.000154120202821

Out[19]:

[<matplotlib.lines.Line2D at 0x1193be320>]




In [20]:

print(i, n_grid[i], error_abs[i])




35 36 0.0101541202028



#### Q. Are there any downsides to letting n be very large?



In [ ]: