[NOTES] Numerical Integration
Fellas, sometimes you just can't integrate analytically. That's when you gotta get a bit freaky.
In early calculus classes I was taught that area under the curve can be thought of as the sum of several tiny rectangles (Reimann's Rule). But the Big-Education kept lying to me and forced me to evaluate integrals analytically. This is an exposé on Big-Educations dark secrets (and a code dump).
Trapezoid Rule
Trapezoid's are more flexible than rectangles so they calculate the area under the curve a bit more accurately. This is how it's defined:
Simpson's Rule
This one is better because it was invented by Bart Simpson from the Simpsons. Here is how it looks:
Code Example
Let's start with integrals of Sin and Cos to establish a benchmark. We already know the integrals so code would look something like this:
import numpy as np
import math
def int_sin(a, b):
return -np.cos(b) + np.cos(a)
def int_cos(a, b):
return np.sin(b) - np.sin(a)
int_sin(0, 1), int_cos(0, 1)
# Output: (0.45969769413186023, 0.8414709848078965)
Here is how we apply both rules:
def num_integration(f, a, b, n = 10, h = 0.1, method = "simpsons"):
if method == "simpsons":
acc = 0
for i in range(1, n//2 + 1):
acc += f(a + (2*i - 2)*h) + 4*f(a + (2*i - 1)*h) + f(a + (2*i)*h)
return h*acc/3
else:
# use trapezoidal method
t1 = h/2 * f(a)
t2 = h/2 * f(b)
acc = 0
for i in range(1, n):
acc += h * f(a + i * h)
return t1 + acc + t2
int_sin(0, 1), num_integration(np.sin, 0, 1), num_integration(np.sin, 0, 1, method="trapezoidal")
# (0.45969769413186023, 0.4596979498238207, 0.4593145488579765)
int_cos(0, 1), num_integration(np.cos, 0, 1), num_integration(np.cos, 0, 1, method="trapezoidal")
# (0.8414709848078965, 0.8414714528488904, 0.8407696420884198)
Too Bad
Too bad that numerical integration is not scalable to high-dimensional functions. You gotta use Monte Carlo Integration for this. But that's an exposé for another time.