P6.1.3Gaussian functions and their derivatives
Using NumPy, it is possible to do this exercise without using a single (Python) loop.
The normalized Gaussian function with mean \(\mu\) and standard deviation \(\sigma\) is
\[g(x) = \frac{1}{\sigma\sqrt{2\pi}}\exp\left( -\frac{(x-\mu)^2}{2\sigma^2}\right).\]
Write a program to calculate and plot the Gaussian functions with \(\mu=0\) and the three values \(\sigma = 0.5, 1, 1.5\). Use a grid of 1000 points in the interval \(-10 \le x \le 10\).
Verify (by direct summation) that the functions are normalized with area 1.
Finally, calculate the first derivative of these functions on the same grid using the first-order central difference approximation:
\[g'(x) \approx \frac{g(x+h) - g(x-h)}{2h}\]
for some suitably-chosen, small \(h\).
Show solution
The following solution avoids Python loops by storing the three Gaussian functions in a single array, y, with shape (1000,3).
gaussians
import numpy as np
import matplotlib.pyplot as plt
def gaussian(x, s, m):
return (
1.0 / (np.sqrt(2.0 * np.pi) * s) * np.exp(-0.5 * ((x - m) / s) ** 2)
)
m = 0
s = np.array([0.5, 1, 1.5])
x, dx = np.linspace(-5, 5, 1000, retstep=True)
x = x[:, np.newaxis]
y = gaussian(x, s, m)
# Integrate approximately by direct summation along the x-axis
int_y = np.sum(gaussian(x, s, m), axis=0) * dx
print("If these numbers are 1 (or close to it), then the g(x) is normalized:")
print(int_y)
# Numerical differentiation by the first-order central difference formula
h = 1.0e-6
dydx = (gaussian(x + h, s, m) - gaussian(x - h, s, m)) / 2 / h
plt.plot(x, y)
plt.plot(x, dydx)
plt.show()
