# 汪群超 Chun-Chao Wang

Dept. of Statistics, National Taipei University, Taiwan

# Lesson 4: Numerical and Symbolic Calculus

### Objective

1. Learn the differentiation and integration in python.
2. Learn when and how to employ numerical and symbolic approaches.
3. Find applications in probability and statistics.

### Prerequisite:

For the function $f(x) = e^{2x}$, compute $f'(x), f''(x), f'''(x)$ at $x=a$.

Note:

1. The scipy module is used for numerical differentiation.
2. The numerical derivative is defined by $f'(x) = \frac{f(x+dx)-f(x)}{dx}$
3. The precision of $f'(a)$ is dependent on “dx”. Try on various “dx” by getting dx smaller and smaller.
import numpy as np
import matplotlib.pyplot as plt
from scipy.misc import derivative

f = lambda x: np.exp(2*x)
a = 0 # evaluate at a

fp = derivative(f, a, dx = 0.001, n = 1) # nth derivative
fpp = derivative(f, a, dx = 0.001, n = 2)
fppp = derivative(f, a, dx = 0.001, n = 3, order = 5) # use 5 points to calculate

print(fp, fpp, fppp)

Note:

1. Print out $f'(a)$ for each “dx” as shown below on the left.
2. Which “dx” gives the most accurate answer?
3. Try to make a line plot of $f'(1)$ vs. “dx” as shown below on the right where matplotlib.pyplot.semilogx() is employed.

Note:

1. It is necessary to separate two lines on a plot by “line style”.
2. Add legend with proper text on the plot.

Let $f(x) = e^{-x} \sin(x)$.

1. Derive the first and the second derivatives $f'(x)$ and $f''(x)$.
2. Evaluate $f'(x)$ at $x=a$ from a symbolic function $f'(x)$.
3. Evaluate symbolic function at more than one points and make a function plot as usual.

Note:
1. The sympy module is employed for symbolic differentiation.
2. The symbolic expression has its own function command such as sympy.exp(), sympy.sin() instead of np.exp(), np.sin().
3. In many situations, the symbolic expression has to be transferred to numerical expression for further computations.
import sympy as sym
import matplotlib.pyplot as plt

x = sym.Symbol('x')
f = sym.exp(-x) * sym.sin(x)
fp = f.diff(x) # the first derivative
fpp = f.diff(x, 2) # the second derivative
print(fp) # print out to see the symboolic function

# substitution; need to convert to float if necessary
print(float(fp.subs(x, 3))) # evaluate f'(3)
# evaluate symbolic function at more than one points
feval = sym.lambdify(x, fp) # transfer from symbolic function to lambda function
a = np.linspace(-3, 3, 1000)

plt.plot(a, feval(a)) # use symbolic computation to draw a figure
plt.show()
print(feval(3)) # compare with the previous one

Draw $f(x) = e^{-x} \sin(x)$ in symbolic way.

Note:

1. Sympy has its own plotting that is different from that of matplotlib.
2. Refer to Sympy document for more attributes in a figure.
import sympy as sym
from sympy.plotting import plot

x = sym.Symbol('x')
f = sym.exp(-x) * sym.sin(x)
fp = f.diff(x)
# draw two separate figures
# f1 = plot(f, xlim = [-4, 4], ylim = [-10, 5])
# f2 = plot(fp, xlim = [-4, 4], ylim = [-10, 5])
# draw two lines in a single figure
f1 = plot(f, xlim = [-4, 4], ylim = [-10, 6], \
line_color = 'b', label = 'f(x)', show = False)
f2 = plot(fp, xlim = [-4, 4], ylim = [-10, 6], \
line_color = 'g', label = "f'(x)", show = False)
f1.append(f2[0]) # append f2 to f1
f1.legend = True
f1.show()

Note:

1. Carefully select x and y range (i.e. xlim and ylim) in order to see the global extrema of both functions.

Compute the definite integral
$\int_0^1 \sin x \;dx$

Note:

1. Use the module Scipy.integrate and the quadratic rule (Simpson’s Rule) to compute the definite integral.
2. Practice to draw the area that the definite integral is defined.
The result from the definite integral by Scipy.integrate contains two numbers. You need to refer to the document for their meanings.
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integral

# define anonymous function
f = lambda x : np.sin(x**2)
lb, ub = 0, 1
print("The integration is {:.4f}".format(result[0]))

# draw the function line -------------------
x = np.linspace(0, np.pi/2, 100)
y = f(x)
fig, ax = plt.subplots(1, 1)
ax.plot(x, y)
# draw the area ----------------------------
x = np.linspace(lb, ub, 100)
y = f(x)
# fill color between y and 0
ax.fill_between(x, y, 0, alpha=0.3, color='b')
ax.set_xlabel('X')
ax.set_title("The area = {:.4f}".format(result[0]))
plt.grid(True, linestyle='--', which='major')
plt.show()

Note :

1. Try to put the formula of definite integral in the middle of the area.

Draw the function $F(x) = \int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}} e^{-t^2/2}\; dt$ for $-5 \leq x \leq 5$.

Note:

1. $F(x)$ is CDF of the standard normal distribution.
2. You already know what the function looks like.
3. The first thought come into mind to compute $F(x)$ for $-5 \leq x \leq 5$ is “for loop”. So do it this way by yourself.
4. This example demonstrates a vectorized function to compute $F(x)$ in a parallel fashion.
5. Learn to make a step plot.
6. To have a smoother line, increase the points evaluated, i.e. decrease the step size.
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integral

n = 100
f = lambda x : np.exp(-x**2/2) / np.sqrt(2*np.pi)
xlim = [-5, 5]
x = np.linspace(xlim[0], xlim[1], n)

def cdf(x):

vec_F = np.vectorize(cdf) # define a vectorized function
F = vec_F(x) # evalaute cdf for all x
plt.plot(x, F, drawstyle='steps-pre') # stairs plot
plt.title('A stair plot of CDF function by integration')
plt.show()

Express the following definite integration as a Riemann sum.

$\int_1^8 x^2 + 3x + 5 \; dx \approx \sum_{k=1}^{N} f(x_k)\Delta_k$

Note:

1. In this example, let the N intervals for the Riemann sum be equally spaced, i.e. $\Delta_k = (8-1)/N$.
2. Learn how to patch a rectangle.
3. Learn how to add animation for the visual effect.
4. Increase the number of intervals and see what happens (how close to the correct answer).
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integral

n = 100
lb, ub = 1, 8
f = lambda x : x**2 + 3*x +5
xlim = [0, 10]
x = np.linspace(xlim[0], xlim[1], n)
plt.plot(x, f(x))
N = 10 # how many intervals
rec = np.linspace(lb, ub, N+1)
for i in range(N) :
plt.fill_between([rec[i], rec[i+1]],
[f(rec[i+1]), f(rec[i+1])], 0, alpha=0.3, color='b')
# plt.pause(0.5) # animation not work in Jupyter notbook
plt.title('Integration by Riemann Sum')
plt.show()

1. In this previous example, the right endpoints are used for $f(x_k)$ as seen in the graph. Try left end points and midpoints and compare the results based on the same number of intervals.
2. Replace the rectangles by trapezoids and calculate the sum.
3. Scipy.integrate also has trapezoidal rule for definite integral. Try to use it.

For the function $f(x) = \int \sqrt{x^2 + 4x + 12}$

1. Compute $\int f(x) \; dx$
2. Evaluate $\int_{-2}^6 f(x) \; dx$.

Note:
1. In the following codes, the integration is called as a method of a symbolic function $f$. The other way is to use the command sym.integrate.
2. When the function is complicated, the symbolic result of integral is not seen as the case in this example. Try simpler functions such as $f(x) = \sin x$ or $f(x) = \sqrt{1+x}$ to see the explicit symbolic result of integration.
import sympy as sym

x = sym.Symbol('x')
f = sym.sqrt(x**2 + 4*x + 12)
# f = sym.sin(x) # check this simple function
# F = sym.integrate(f, x)
F = f.integrate(x) # indefinite integral
print(F)
# Definite integral with a symbolic result
I = sym.integrate(f, (x, -2, 6))

print("The definite integral is {:.4f}".format(float(I)))

Evaluate the following definite integrals:

1. $\int_0^5 100 e^x\; dx$
2. $\int_0^7 \frac{x^2}{2} e^{-x}\;dx$ (This is Gamma(3,1))
3. $\int_{-1.96}^{1.96} \frac{1}{\sqrt{2\pi}} e^{-x^2/2}\;dx$
4. $E[\theta] =\frac{\int_0^1 (2+\theta)^{125}(1-\theta)^{38}\theta^{35}\;d\theta}{\int_0^1 (2+\theta)^{125}(1-\theta)^{38}\theta^{34}\;d\theta}$

Note:
1. You may want to try both numerical and symbolic approaches.
2. Draw the function and the patch area for each definite integral.

1. Calculate the area enclosed by the functions $f(x) = 2x - x^2$ and $g(x) = x^2$.
2. Plot the functions and patch the enclosed area.