Standard imports
%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt
Aside: I wanted to setup a fun function to integrate. Here's how I played with a polynomial till I got one I liked. First, make a (mostly empty) plot to put our interactive output:
from ipywidgets import interact
x = np.linspace(0, 10)
fig, ax = plt.subplots()
(lines,) = plt.plot([], [])
ax.set_xlim(0, 10)
ax.set_ylim(-10, 10)
plt.show()
def f(a=(-1.0, 1.0), b=(-1.0, 1.0), c=(-10.0, 10.0), d=(-10.0, 10.0)):
lines.set_data(x, np.poly1d([a, b, c, d, 0])(x))
interact(f)
Okay, I have something I like. Let's compute:
$$ Ax^4 + Bx^3 + Cx^2 + Dx + E $$Where $A=-0.1, B=1.0, C=-3.1, D=3.6, E=0.5$. I've just picked nice numbers to give an interesting plot and integral.
x = np.linspace(0, 5.5, 500)
f = np.poly1d([-0.1, 1.0, -3.1, 3.6, 0.5])
y = f(x)
x_i = np.linspace(0, 5.5, 4)
y_i = f(x_i)
fig, ax = plt.subplots()
plt.plot(x, y)
plt.plot(x_i, y_i, "o-r")
plt.plot(x_i, y_i * 0, "o-r")
plt.vlines(x_i, y_i * 0, y_i, "r")
plt.show()
Let's try integrating with several different values of N (picking odd values for later use - we could have even numbers if we wanted). Let's also compare with np.trapz.
for i in 3, 5, 7, 11, 21, 101, 1001:
x_i = np.linspace(0, 5.5, i)
y_i = f(x_i)
h = np.diff(x_i)[0]
my_intg = 1 / 2 * h * np.sum(y_i[:-1] + y_i[1:])
of_intg = np.trapz(y_i, x_i)
print(f"{i: 5} {my_intg: 16.12f} {of_intg: 16.12f}")
We can also calculate the exact value, since it's just a polynomial. In fact, since we built it with a specialized class in numpy, we can just call .integ to get the analytic integral (and then we call that with our points):
f.integ()(5.5)
from scipy import integrate
Let's make sure we are selecting the times we want with our indexing expressions:
index = list(range(7))
index[:-1:2], index[1::2], index[2::2]
Looks good! Now let's try the rule, comparing with integrate.simps:
for i in 3, 5, 7, 11, 21, 101, 1001:
x_i = np.linspace(0, 5.5, i)
y_i = f(x_i)
h = np.diff(x_i)[0]
my_intg = 1 / 3 * h * np.sum(y_i[:-1:2] + 4 * y_i[1::2] + y_i[2::2])
of_intg = integrate.simps(y_i, x_i)
print(f"{i: 5} {my_intg: 16.12f} {of_intg: 16.12f}")
p = np.poly1d([7, -8, -3, 3])
x = np.linspace(-1, 1)
fig, ax = plt.subplots()
ax.plot(x, p(x))
plt.show()
p.integ()(1) - p.integ()(-1)
p(-np.sqrt(1 / 3)) + p(np.sqrt(1 / 3))
from sympy.abc import x
import sympy
sympy.init_printing()
p(x).expand()
p(x).integrate(x).expand()
Slightly improved version:
max_in = 11
w = np.zeros(2001, dtype=np.float64)
x = np.zeros(2001, dtype=np.float64)
def func_e(x):
return np.exp(-x)
def gauss(npts, job, a, b, x, w, vmin=0.0, vmax=0.0):
m = i = j = t = t1 = pp = p1 = p2 = p3 = 0.0
eps = 3e-14 # Accuracy: ******ADJUST THIS*******!
m = (npts + 1) // 2
for i in range(1, m + 1):
t = np.cos(np.pi * (float(i) - 0.25) / (float(npts) + 0.5))
t1 = 1
while abs(t - t1) >= eps:
p1 = 1.0
p2 = 0.0
for j in range(1, npts + 1):
p3 = p2
p2 = p1
p1 = ((2.0 * float(j) - 1) * t * p2 - (float(j) - 1.0) * p3) / (
float(j)
)
pp = npts * (t * p1 - p2) / (t * t - 1.0)
t1 = t
t = t1 - p1 / pp
x[i - 1] = -t
x[npts - i] = t
w[i - 1] = 2.0 / ((1.0 - t * t) * pp * pp)
w[npts - i] = w[i - 1]
if job == 0:
for i in range(0, npts):
x[i] = x[i] * (b - a) / 2.0 + (b + a) / 2.0
w[i] = w[i] * (b - a) / 2.0
if job == 1:
for i in range(0, npts):
xi = x[i]
x[i] = a * b * (1.0 + xi) / (b + a - (b - a) * xi)
w[i] = (
w[i]
* 2.0
* a
* b
* b
/ ((b + a - (b - a) * xi) * (b + a - (b - a) * xi))
)
if job == 2:
for i in range(0, npts):
xi = x[i]
x[i] = (b * xi + b + a + a) / (1.0 - xi)
w[i] = w[i] * 2.0 * (a + b) / ((1.0 - xi) * (1.0 - xi))
def gaussint(no, min, max):
quadra = 0.0
gauss(no, 0, min, max, x, w)
for n in range(0, no):
quadra += func_e(x[n]) * w[n]
return quadra
for i in range(3, max_in + 1, 2):
result = gaussint(i, 0.0, 1.0)
print(" i ", i, " err ", abs(result - 1 + 1 / np.exp(1)))
Much better version:
def gauss_new(npts, job, a, b, vmin=0.0, vmax=0.0, eps=3e-14):
x = np.zeros(npts, dtype=np.float64)
w = np.zeros(npts, dtype=np.float64)
for i in range((npts + 1) // 2):
t = np.cos(np.pi * (i + 0.75) / (npts + 0.5))
t1 = 1
while abs(t - t1) >= eps:
p1 = 1.0
p2 = 0.0
for j in range(npts):
p3 = p2
p2 = p1
p1 = ((2 * j + 1) * t * p2 - j * p3) / (j + 1)
pp = npts * (t * p1 - p2) / (t ** 2 - 1)
t1 = t
t = t1 - p1 / pp
x[i] = -t
x[~i] = t
w[i] = 2 / ((1 - t ** 2) * pp ** 2)
w[~i] = w[i]
if job == 0:
w = w * (b - a) / 2
x = x * (b - a) / 2 + (b + a) / 2
elif job == 1:
# Use old x for w and x calcs
w = w * 2 * a * b ** 2 / (b + a - (b - a) * x) ** 2
x = a * b * (1 + x) / (b + a - (b - a) * x)
elif job == 2:
# Use old x for w and x calcs
w = w * 2 * (a + b) / (1 - x) ** 2
x = (b * x + b + 2 * a) / (1 - x)
else:
raise RuntimeError(f"Job ({job}) must be 0, 1, or 2")
return x, w
def gaussint_new(f, no, min, max):
quadra = 0.0
x, w = gauss_new(no, 0, min, max)
for n in range(no):
quadra += f(x[n]) * w[n]
return quadra
for i in range(3, max_in + 1, 2):
result = gaussint_new(func_e, i, 0, 1)
print("i", i, "err", abs(result - 1 + 1 / np.exp(1)))
Using the classic QUADPACK from the Fortran days...
result, *_ = integrate.quad(func_e, 0, 1)
print(abs(result - 1 + 1 / np.exp(1)))
Using fixed_quad, which should be fairly similar to our algorithm:
for i in range(3, max_in + 1, 2):
result, *_ = integrate.fixed_quad(func_e, 0, 1, n=i)
print("i", i, "err", abs(result - 1 + 1 / np.exp(1)))
Look at the "See Also" section of the docstrings for the integration functions!