<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 ★ Numerical Differentiaion and Integration ★
Step1:  5.1 Numerical Differentiation
Two-point forward-difference formula
$f'(x) = \frac{f(x+h) - f(x)}{h} - \frac{h}{2}f''(c)$
where $c$ is between $x$ and $x+h$
Example
Use the two-point forward-difference formula with $h = 0.1$ to approximate the derivative of $f(x) = 1/x$ at $x = 2$
Step2:  Three-point centered-difference formula
$f'(x) =  \frac{f(x+h) - f(x-h)}{2h} - \frac{h^2}{6}f'''(c)$
where $x-h < c < x+h$
Example
Use the three-point centered-difference formula with $h = 0.1$ to approximate the derivative of $f(x) = 1 / x$ at $x = 2$
Step3:  Three-point centered-difference formula for second derivative
$f''(x) = \frac{f(x - h) - 2f(x) + f(x + h)}{h^2} - \frac{h^2}{12}f^{(iv)}(c)$
for some $c$ between $x - h$ and $x + h$
Rounding error
Example
Approximate the derivative of $f(x) = e^x$ at $x = 0$
Step4:  Extrapolation for order n formula
$ Q \approx \frac{2^nF(h/2) - F(h)}{2^n - 1} $
Step5:  5.2 Newton-Cotes Formulas For Numerical Integration
Trapezoid Rule
$\int_{x_0}^{x_1} f(x) dx = \frac{h}{2}(y_0 + y_1) - \frac{h^3}{12}f''(c)$
where $h = x_1 - x_0$ and $c$ is between $x_0$ and $x_1$
Simpson's Rule
$\int_{x_0}^{x_2} f(x) dx = \frac{h}{3}(y_0 + 4y_1 + y_2) - \frac{h^5}{90}f^{(iv)}(c)$
where $h = x_2 - x_1 = x_1 - x_0$ and $c$ is between $x_0$ and $x_2$
Example
Apply the Trapezoid Rule and Simpson's Rule to approximate $\int_{1}^{2} \ln(x) dx$ and find an upper bound for the error in your approximations
Step6:  Composite Trapezoid Rule
$\int_{a}^{b} f(x) dx = \frac{h}{2} \left ( y_0 + y_m + 2\sum_{i=1}^{m-1}y_i \right ) - \frac{(b-a)h^2}{12}f''(c)$
where $h = (b - a) / m $ and $c$ is between $a$ and $b$
Composite Simpson's Rule
$ \int_{a}^{b}f(x)dx = \frac{h}{3}\left [ y_0 + y_{2m} + 4\sum_{i=1}^{m}y_{2i-1} + 2\sum_{i=1}^{m - 1}y_{2i} \right ] - \frac{(b-a)h^4}{180}f^{(iv)}(c) $
where $c$ is between $a$ and $b$
Example
Carry out four-panel approximations of $\int_{1}^{2} \ln{x} dx$ using the composite Trapezoid Rule and composite Simpson's Rule
Step7:  Midpoint Rule
$ \int_{x_0}^{x_1} f(x)dx = hf(\omega) + \frac{h^3}{24}f''(c) $
where $ h = (x_1 - x_0) $, $\omega$ is the midpoint $ x_0 + h / 2 $, and $c$ is between $x_0$ and $x_1$
Composite Midpoint Rule
$ \int_{a}^{b} f(x) dx = h \sum_{i=1}^{m}f(\omega_{i}) + \frac{(b - a)h^2}{24} f''(c) $
where $h = (b - a) / m$ and $c$ is between $a$ and $b$. The $\omega_{i}$ are the midpoints of the $m$ equal subintervals of $[a,b]$
Example
Approximate $\int_{0}^{1} \frac{\sin x}{x} dx$ by using the composite Midpoint Rule with $m = 10$ panels
Step8:  5.3 Romberg Integration
Step9:  Example
Apply Romberg Integration to approximate $\int_{1}^{2} \ln{x}dx$
Step10:  5.4 Adaptive Quadrature
Step11:  Example
Use Adaptive Quadrature to approximate the integral $ \int_{-1}^{1} (1 + \sin{e^{3x}}) dx $
Step12:  5.5 Gaussian Quadrature
Step13:  Example
Approximate $\int_{-1}^{1} e^{-\frac{x^2}{2}}dx$ using Gaussian Quadrature
Step14:  Example
Approximate the integral $\int_{1}^{2} \ln{x} dx$ using Gaussian Quadrature
<ASSISTANT_TASK:>
Python Code:
# Import modules
import math
import sympy as sym
import numpy as np
import scipy 
import matplotlib.pyplot as plt
import plotly
import plotly.plotly as ply
import plotly.figure_factory as ply_ff
from IPython.display import Math
from IPython.display import display
# Startup plotly
plotly.offline.init_notebook_mode(connected=True)
''' Fix MathJax issue '''
# The polling here is to ensure that plotly.js has already been loaded before
# setting display alignment in order to avoid a race condition.
from IPython.core.display import display, HTML
display(HTML(
    '<script>'
        'var waitForPlotly = setInterval( function() {'
            'if( typeof(window.Plotly) !== "undefined" ){'
                'MathJax.Hub.Config({ SVG: { font: "STIX-Web" }, displayAlign: "center" });'
                'MathJax.Hub.Queue(["setRenderer", MathJax.Hub, "SVG"]);'
                'clearInterval(waitForPlotly);'
            '}}, 250 );'
    '</script>'
))
Explanation: ★ Numerical Differentiaion and Integration ★
End of explanation
# Parameters
x = 2
h = 0.1
# Symbolic computation
sym_x = sym.Symbol('x')
sym_deri_x1 = sym.diff(1 / sym_x, sym_x)
sym_deri_x1_num = sym_deri_x1.subs(sym_x, x).evalf()
# Approximation
f = lambda x : 1 / x
deri_x1 = (f(x + h) - f(x)) / h
# Comparison
print('approximate = %f, real value = %f, backward error = %f' %(deri_x1, sym_deri_x1_num, abs(deri_x1 - sym_deri_x1_num)) )
Explanation: 5.1 Numerical Differentiation
Two-point forward-difference formula
$f'(x) = \frac{f(x+h) - f(x)}{h} - \frac{h}{2}f''(c)$
where $c$ is between $x$ and $x+h$
Example
Use the two-point forward-difference formula with $h = 0.1$ to approximate the derivative of $f(x) = 1/x$ at $x = 2$
End of explanation
# Parameters
x = 2
h = 0.1
f = lambda x : 1 / x
# Symbolic computation
sym_x = sym.Symbol('x')
sym_deri_x1 = sym.diff(1 / sym_x, sym_x)
sym_deri_x1_num = sym_deri_x1.subs(sym_x, x).evalf()
# Approximation
deri_x1 = (f(x + h) - f(x - h)) / (2 * h)
# Comparison
print('approximate = %f, real value = %f, backward error = %f' %(deri_x1, sym_deri_x1_num, abs(deri_x1 - sym_deri_x1_num)) )
Explanation: Three-point centered-difference formula
$f'(x) =  \frac{f(x+h) - f(x-h)}{2h} - \frac{h^2}{6}f'''(c)$
where $x-h < c < x+h$
Example
Use the three-point centered-difference formula with $h = 0.1$ to approximate the derivative of $f(x) = 1 / x$ at $x = 2$
End of explanation
# Parameters
f = lambda x : math.exp(x)
real_value = 1
h_msg = "$10^{-%d}$"
twp_deri_x1 = lambda x, h : ( f(x + h) - f(x) ) / h
thp_deri_x1 = lambda x, h : ( f(x + h) - f(x - h) ) / (2 * h)
data = [
    ["h", 
     "$f'(x) \\approx \\frac{e^{x+h} - e^x}{h}$", 
     "error", 
     "$f'(x) \\approx \\frac{e^{x+h} - e^{x-h}}{2h}$", 
     "error"],
]
for i in range(1,10):
    h = pow(10, -i)
    twp_deri_x1_value = twp_deri_x1(0, h) 
    thp_deri_x1_value = thp_deri_x1(0, h)
    row = ["", "", "", "", ""]
    row[0] = h_msg %i
    row[1] = '%.14f' %twp_deri_x1_value
    row[2] = '%.14f' %abs(twp_deri_x1_value - real_value)
    row[3] = '%.14f' %thp_deri_x1_value
    row[4] = '%.14f' %abs(thp_deri_x1_value - real_value)
    data.append(row)
table = ply_ff.create_table(data)
plotly.offline.iplot(table, show_link=False)
Explanation: Three-point centered-difference formula for second derivative
$f''(x) = \frac{f(x - h) - 2f(x) + f(x + h)}{h^2} - \frac{h^2}{12}f^{(iv)}(c)$
for some $c$ between $x - h$ and $x + h$
Rounding error
Example
Approximate the derivative of $f(x) = e^x$ at $x = 0$
End of explanation
sym.init_printing(use_latex=True)
x = sym.Symbol('x')
dx = sym.diff(sym.exp(sym.sin(x)), x)
Math('Derivative : %s' %sym.latex(dx) )
Explanation: Extrapolation for order n formula
$ Q \approx \frac{2^nF(h/2) - F(h)}{2^n - 1} $
End of explanation
# Apply Trapezoid Rule
trapz = scipy.integrate.trapz([np.log(1), np.log(2)], [1, 2])
# Evaluate the error term of Trapezoid Rule
sym_x = sym.Symbol('x')
expr = sym.diff(sym.log(sym_x), sym_x, 2)
trapz_err = abs(expr.subs(sym_x, 1).evalf() / 12)
# Print out results
print('Trapezoid rule : %f and upper bound error : %f' %(trapz, trapz_err) )
# Apply Simpson's Rule
area = scipy.integrate.simps([np.log(1), np.log(1.5), np.log(2)], [1, 1.5, 2])
# Evaluate the error term
sym_x = sym.Symbol('x')
expr = sym.diff(sym.log(sym_x), sym_x, 4)
simps_err = abs( pow(0.5, 5) / 90 * expr.subs(sym_x, 1).evalf() )
# Print out results
print('Simpson\'s rule : %f and upper bound error : %f' %(area, simps_err) )
Explanation: 5.2 Newton-Cotes Formulas For Numerical Integration
Trapezoid Rule
$\int_{x_0}^{x_1} f(x) dx = \frac{h}{2}(y_0 + y_1) - \frac{h^3}{12}f''(c)$
where $h = x_1 - x_0$ and $c$ is between $x_0$ and $x_1$
Simpson's Rule
$\int_{x_0}^{x_2} f(x) dx = \frac{h}{3}(y_0 + 4y_1 + y_2) - \frac{h^5}{90}f^{(iv)}(c)$
where $h = x_2 - x_1 = x_1 - x_0$ and $c$ is between $x_0$ and $x_2$
Example
Apply the Trapezoid Rule and Simpson's Rule to approximate $\int_{1}^{2} \ln(x) dx$ and find an upper bound for the error in your approximations
End of explanation
# Apply composite Trapezoid Rule
x = np.linspace(1, 2, 5)
y = np.log(x)
trapz = scipy.integrate.trapz(y, x)
# Error term
sym_x = sym.Symbol('x')
expr = sym.diff(sym.log(sym_x), sym_x, 2)
trapz_err = abs( (2 - 1) * pow(0.25, 2) / 12 * expr.subs(sym_x, 1).evalf() )
print('Trapezoid Rule : %f, error = %f' %(trapz, trapz_err) )
# Apply composite Trapezoid Rule
x = np.linspace(1, 2, 9)
y = np.log(x)
area = scipy.integrate.simps(y, x)
# Error term
sym_x = sym.Symbol('x')
expr = sym.diff(sym.log(sym_x), sym_x, 4)
simps_err = abs( (2 - 1) * pow(0.125, 4) / 180 * expr.subs(sym_x, 1).evalf() )
print('Simpson\'s Rule : %f, error = %f' %(area, simps_err) )
Explanation: Composite Trapezoid Rule
$\int_{a}^{b} f(x) dx = \frac{h}{2} \left ( y_0 + y_m + 2\sum_{i=1}^{m-1}y_i \right ) - \frac{(b-a)h^2}{12}f''(c)$
where $h = (b - a) / m $ and $c$ is between $a$ and $b$
Composite Simpson's Rule
$ \int_{a}^{b}f(x)dx = \frac{h}{3}\left [ y_0 + y_{2m} + 4\sum_{i=1}^{m}y_{2i-1} + 2\sum_{i=1}^{m - 1}y_{2i} \right ] - \frac{(b-a)h^4}{180}f^{(iv)}(c) $
where $c$ is between $a$ and $b$
Example
Carry out four-panel approximations of $\int_{1}^{2} \ln{x} dx$ using the composite Trapezoid Rule and composite Simpson's Rule
End of explanation
# Parameters
m = 10
h = (1 - 0) / m
f = lambda x : np.sin(x) / x
mids = np.arange(0 + h/2, 1, h)
# Apply composite midpoint rule
area = h * np.sum(f(mids))
# Error term
sym_x = sym.Symbol('x')
expr = sym.diff(sym.sin(sym_x) / sym_x, sym_x, 2)
mid_err = abs( (1 - 0) * pow(h, 2) / 24 * expr.subs(sym_x, 1).evalf() )
# Print out
print('Composite Midpoint Rule : %.8f, error = %.8f' %(area, mid_err) )
Explanation: Midpoint Rule
$ \int_{x_0}^{x_1} f(x)dx = hf(\omega) + \frac{h^3}{24}f''(c) $
where $ h = (x_1 - x_0) $, $\omega$ is the midpoint $ x_0 + h / 2 $, and $c$ is between $x_0$ and $x_1$
Composite Midpoint Rule
$ \int_{a}^{b} f(x) dx = h \sum_{i=1}^{m}f(\omega_{i}) + \frac{(b - a)h^2}{24} f''(c) $
where $h = (b - a) / m$ and $c$ is between $a$ and $b$. The $\omega_{i}$ are the midpoints of the $m$ equal subintervals of $[a,b]$
Example
Approximate $\int_{0}^{1} \frac{\sin x}{x} dx$ by using the composite Midpoint Rule with $m = 10$ panels
End of explanation
def romberg(f, a, b, step):
    R = np.zeros(step * step).reshape(step, step)
    R[0][0] = (b - a) * (f(a) + f(b)) / 2
    for j in range(1, step):
        h = (b - a) / pow(2, j)
        summ = 0
        for i in range(1, pow(2, j - 1) + 1):
            summ += h * f(a + (2 * i - 1) * h)
        R[j][0] = 0.5 * R[j - 1][0] + summ
        
        for k in range(1, j + 1):
            R[j][k] = ( pow(4, k) * R[j][k - 1] - R[j - 1][k - 1] ) / ( pow(4, k) - 1 )
        
    return R[step - 1][step - 1]
Explanation: 5.3 Romberg Integration
End of explanation
f = lambda x : np.log(x)
result = romberg(f, 1, 2, 4)
print('Romberg Integration : %f' %(result) )
f = lambda x : np.log(x)
result = scipy.integrate.romberg(f, 1, 2, show=True)
print('Romberg Integration : %f' %(result) )
Explanation: Example
Apply Romberg Integration to approximate $\int_{1}^{2} \ln{x}dx$
End of explanation
''' Use Trapezoid Rule '''
def adaptive_quadrature(f, a, b, tol):
    return adaptive_quadrature_recursively(f, a, b, tol, a, b, 0)
    
def adaptive_quadrature_recursively(f, a, b, tol, orig_a, orig_b, deep):
    c = (a + b) / 2
    S = lambda x, y : (y - x) * (f(x) + f(y)) / 2
    if abs( S(a, b) - S(a, c) - S(c, b) ) < 3 * tol * (b - a) / (orig_b - orig_a) or deep > 20 :
        return S(a, c) + S(c, b)
    else:
        return adaptive_quadrature_recursively(f, a, c, tol / 2, orig_a, orig_b, deep + 1) + adaptive_quadrature_recursively(f, c, b, tol / 2, orig_a, orig_b, deep + 1)
''' Use Simpon's Rule '''
def adaptive_quadrature(f, a, b, tol):
    return adaptive_quadrature_recursively(f, a, b, tol, a, b, 0)
    
def adaptive_quadrature_recursively(f, a, b, tol, orig_a, orig_b, deep):
    c = (a + b) / 2
    S = lambda x, y : (y - x) * ( f(x) + 4 * f((x + y) / 2) + f(y) ) / 6
    if abs( S(a, b) - S(a, c) - S(c, b) ) < 15 * tol  or deep > 20 :
        return S(a, c) + S(c, b)
    else:
        return adaptive_quadrature_recursively(f, a, c, tol / 2, orig_a, orig_b, deep + 1) + adaptive_quadrature_recursively(f, c, b, tol / 2, orig_a, orig_b, deep + 1)
Explanation: 5.4 Adaptive Quadrature
End of explanation
f = lambda x : 1 + np.sin(np.exp(3 * x))
val = adaptive_quadrature(f, -1, 1, tol=1e-12)
print(val)
Explanation: Example
Use Adaptive Quadrature to approximate the integral $ \int_{-1}^{1} (1 + \sin{e^{3x}}) dx $
End of explanation
poly = scipy.special.legendre(2)
# Find roots of polynomials
comp = scipy.linalg.companion(poly)
roots = scipy.linalg.eig(comp)[0]
Explanation: 5.5 Gaussian Quadrature
End of explanation
f = lambda x : np.exp(-np.power(x, 2) / 2)
quad = scipy.integrate.quadrature(f, -1, 1)
print(quad[0])
# Parametes
a = -1
b = 1
deg = 3
f = lambda x : np.exp( -np.power(x, 2) / 2 )
x, w = scipy.special.p_roots(deg) # Or use numpy.polynomial.legendre.leggauss
quad = np.sum(w * f(x))
    
print(quad)
Explanation: Example
Approximate $\int_{-1}^{1} e^{-\frac{x^2}{2}}dx$ using Gaussian Quadrature
End of explanation
# Parametes
a = 1
b = 2
deg = 4
f = lambda t : np.log( ((b - a) * t + b + a) / 2) * (b - a) / 2
x, w = scipy.special.p_roots(deg)
np.sum(w * f(x))
Explanation: Example
Approximate the integral $\int_{1}^{2} \ln{x} dx$ using Gaussian Quadrature
End of explanation
<END_TASK>