Week 3: Numerical Intergration

Consider a function lets say $f(x) = xe^{-x^2} $. How can we determine the area under the curve from say 1 to 10

In mathematical terms this is just $ \int_1^{100} x e^{-x^2} \,dx $. We simply can do this by U-substitution and solve.

For us that would be $ \frac{-e^{-x^2}}{2} |_1^{10} $ which is approx 0.18393972

In [4]:
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
%config InlineBackend.figure_format='retina'

Lets plot the graph

In [6]:
X = np.arange(0,10,0.1)
plt.plot(X,X)
plt.show()
In [7]:
def func(x):
    return np.exp(-x**2)*x
In [9]:
X = np.arange(-2,10,0.1)
plt.plot(X,func(X))
plt.show()
In [12]:
X = np.arange(-2,10,0.1)
Y = np.arange(1,10,0.1)
plt.plot(X,func(X))
plt.fill_between(Y,func(Y),color='r')
plt.xlim(-2,4)
plt.show()

Trapezoidal Rule

In [20]:
X = np.arange(-2,10,0.1)
a=2
b=3
Y=np.arange(a,b,0.1)
plt.plot(X,X)
plt.fill_between(Y,Y,color='r')
plt.plot([a,b],[a,0],c='black')
plt.show()

$$ \text{Area} = \frac{a}{2}\Delta x+\frac{b}{2}\Delta x $$

Or another way to say this with functions

$$ \text{Area} = \frac{f(x_n)+f(x_{n+1})}{2}\Delta x $$

In [21]:
def Trapezoid_int(a,b,f,n):
    Deltax = (b-a)/n
    Area = 0
    for c in range(0,n):
        Area += (f(a+c*Deltax)+f(a+(c+1)*Deltax))/2*Deltax
    return Area
In [22]:
Trapezoid_int(1,10,func,3)
Out[22]:
0.5518205121792601
In [23]:
X = np.arange(-2,10,0.1)
Y = np.arange(1,10,0.1)
plt.plot(X,func(X))
plt.fill_between(Y,func(Y),color='r')
plt.xlim(-2,4)
plt.show()
In [24]:
def Trap_top_line(x1,x2,f,x):
    y1 = f(x1)
    y2 = f(x2)
    return (y2-y1)/(x2-x1)*(x-x1)+y1
In [32]:
a=1
b=10
n=10
Deltax = (b-a)/n

X = np.arange(a-3,b+1,0.1)

plt.plot(X,func(X))
for c in range(0,n):
    Y = np.arange(a+c*Deltax,a+(c+1)*Deltax,0.01)
    plt.fill_between(Y,Trap_top_line(a+c*Deltax,a+(c+1)*Deltax,func,Y),color='r')
    plt.vlines(a+c*Deltax,0,func(a+c*Deltax),alpha=0.3)


plt.show()
In [31]:
Trapezoid_int(1,10,func,10)
Out[31]:
0.2128002306378246

Lets look at $ \sin (x) $

In [38]:
def func2(x):
    return np.sin(x)**2
In [40]:
Trapezoid_int(-np.pi,np.pi,func2,20)
Out[40]:
3.141592653589793
In [44]:
a=-np.pi
b=np.pi
n=100
Deltax = (b-a)/n

X = np.arange(a-1,b+1,0.1)

plt.plot(X,func2(X))
for c in range(0,n):
    Y = np.arange(a+c*Deltax,a+(c+1)*Deltax,0.01)
    plt.fill_between(Y,Trap_top_line(a+c*Deltax,a+(c+1)*Deltax,func2,Y),color='r')
    plt.vlines(a+c*Deltax,0,func2(a+c*Deltax),alpha=0.3)


plt.show()
In [49]:
Trapezoid_int(-np.pi,np.pi,func2,3)
Out[49]:
3.1415926535897922
In [53]:
for c in range(2,50):
    plt.plot(c,Trapezoid_int(1,10,func,c),'.',c='b')
    plt.ylim(0.15,0.3)
In [56]:
Trapezoid_int(1,10,func,7)
Out[56]:
0.25232692530495565

Midpoint Rule

In [57]:
def Midpoint_int(a,b,f,n):
    Deltax = (b-a)/n
    Area = 0
    for c in range(0,n):
        Area += (f(a+(c+0.5)*Deltax)*Deltax)
        
    return Area
In [58]:
Midpoint_int(1,10,func,7)
Out[58]:
0.1428092081588567
In [59]:
def Box_line(x1,x2,f,x):
    y1 = f((x1+x2)/2)
    return y1
In [61]:
a=1
b=10
n=100
Deltax = (b-a)/n

X = np.arange(a-3,b+1,0.1)
plt.plot(X,func(X))

for c in range(0,n):
    Y = np.arange(a+c*Deltax,a+(c+1)*Deltax,0.01)
    plt.fill_between(Y,Box_line(a+c*Deltax,a+(c+1)*Deltax,func,Y),color='r')
    plt.vlines(a+c*Deltax,0,func(a+c*Deltax),alpha=0.3)


plt.show()
In [63]:
Midpoint_int(1,10,func,1000)
Out[63]:
0.1839384789632739
In [64]:
Trapezoid_int(1,10,func,1000)
Out[64]:
0.18394220380547277
In [66]:
for c in range(2,50):
    plt.plot(c,Trapezoid_int(1,10,func,c),'.',c='b')
    plt.plot(c,Midpoint_int(1,10,func,c),'.',c='r')
plt.ylim(0.15,0.3)
Out[66]:
(0.15, 0.3)
In [67]:
for c in range(2,50):
    plt.plot(c,Trapezoid_int(1,10,func2,c),'.',c='b')
    plt.plot(c,Midpoint_int(1,10,func2,c),'.',c='r')
In [ ]: