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.,1.), b=(-1.,1.), c=(-10.,10.), d=(-10.,10.)):
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([-.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.
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.; p2=0.
for j in range(1, npts + 1):
p3 = p2; p2 = p1
p1 = ((2.*float(j)-1)*t*p2 - (float(j)-1.)*p3)/(float(j))
pp = npts*(t*p1 - p2)/(t*t - 1.)
t1=t
t=t1 - p1/pp
x[i - 1] = - t
x[npts - i] = t
w[i - 1] = 2./( (1. - 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. + (b + a)/2.
w[i] = w[i]*(b - a)/2.
if (job == 1):
for i in range(0, npts):
xi = x[i]
x[i] = a*b*(1. + xi) / (b + a - (b - a)*xi)
w[i] = w[i]*2.*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.-xi)
w[i] = w[i]*2.*(a + b)/( (1. - xi)*(1. - xi) )
def gaussint (no, min, max):
quadra = 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.
p2=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.
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!