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')

Week 3: Second lecture

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

%matplotlib inline
%config InlineBackend.figure_format='retina'
In [2]:
def func(x):
    return x**2 +2*x-1
In [3]:
func(2)
Out[3]:
7
In [4]:
X = np.arange(-2,2,0.01)
In [5]:
plt.plot(X,func(X))
plt.show()
In [6]:
F = lambda c: c**2 +2*c -1
plt.plot(X,F(X))
Out[6]:
[<matplotlib.lines.Line2D at 0x7fde604a5fd0>]
In [7]:
XX = np.linspace(-2,2,5000)
In [8]:
plt.plot(XX,F(XX))
Out[8]:
[<matplotlib.lines.Line2D at 0x7fdeb0926f90>]
In [9]:
import Pythonic_files as PF
In [10]:
PF.func(2)
Out[10]:
7
In [15]:
PF.Generate_Trap_int(-2,4,F,50)

The area using trapezoidal integration is : 30.014399999999995
In [16]:
PF.Generate_MidP_int(-2,4,F,50)

The area using Mid-point integration is : 29.992799999999995

This is a bit weird,

In [22]:
answer = 1+2
In [23]:
expected = 3
In [24]:
answer == expected
Out[24]:
True
In [25]:
answer = 0.1 + 0.2
In [26]:
expected = 0.3
In [27]:
answer == expected
Out[27]:
False
In [28]:
answer
Out[28]:
0.30000000000000004
In [43]:
%%time
PF.Midpoint_int(-2,2,F,500000)

CPU times: user 232 ms, sys: 18.1 ms, total: 250 ms
Wall time: 303 ms
Out[43]:
1.3333333333119477
In [44]:
%%time
PF.Trapezoid_int(-2,2,F,500000)

CPU times: user 429 ms, sys: 18.2 ms, total: 447 ms
Wall time: 476 ms
Out[44]:
1.3333333333759765
In [45]:
def mid_point_Vectorization(a,b,f,n):
    Deltax= (b-a)/n
    X = np.arange(a+Deltax/2,b+Deltax/2,Deltax)
    return Deltax*sum(f(X))
In [46]:
%%time
mid_point_Vectorization(-2,2,F,500000)

CPU times: user 98 ms, sys: 9.55 ms, total: 108 ms
Wall time: 113 ms
Out[46]:
1.3333333333386423

Lets look at some 3-d

In [47]:
from mpl_toolkits import mplot3d

import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

%config InlineBackend.figure_format='retina'
In [53]:
fig = plt.figure() # this generates something that you will add to, its like plt.show()

ax = plt.axes(projection ="3d") # makes sure that its a 3-d

ax.view_init(50, 35)
In [71]:
fig = plt.figure() 
ax = plt.axes(projection ="3d") 

X=np.linspace(-15,15,200)
Y=np.linspace(-15,15,200)

X,Y = np.meshgrid(X,Y)

R = np.sqrt(X**2 + Y**2)

Z=np.sin(R)*np.exp(-R/10)

ax.plot_surface(X,Y,Z,linewidth=0, antialiased=False)

ax.view_init(70, 35)
In [77]:
fig = plt.figure() 
ax = plt.axes(projection ="3d") 
ax.view_init(70, 35)

ax.contour3D(X,Y,Z,50,cmap='binary')
ax.set_xlabel('This is x')
Out[77]:
Text(0.5, 0, 'This is x')
In [81]:
fig = plt.figure() 
ax = plt.axes(projection ="3d") 
ax.view_init(70, 35)

ax.contour3D(X,Y,Z,50,cmap='inferno')
Out[81]:
<matplotlib.contour.QuadContourSet at 0x7fde81626c10>
In [96]:
fig = plt.figure() 
ax = plt.axes(projection ="3d") 

zline = np.linspace(0,45,1000)
xline = np.sin(zline)*np.exp(zline/20)
yline= np.cos(zline)*np.exp(zline/20)

ax.plot3D(xline,yline,zline)

zpoints = 45*np.random.random(300)
xpoints = np.sin(zpoints)*np.exp(zpoints/20) + 1*np.random.random()
ypoints = np.cos(zpoints)*np.exp(zpoints/20) + 1*np.random.random()

ax.view_init(30, 35)

ax.scatter3D(xpoints,ypoints,zpoints,c=(max(zpoints)-zpoints),cmap='inferno')
Out[96]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fdeb2059710>
In [89]:
np.random.random(5)
Out[89]:
array([0.55222841, 0.24222912, 0.78834489, 0.26655526, 0.6982878 ])
In [97]:
def func3D(x,y):
    return np.exp(-(x**2+y**2)/100)*(np.cos(x**2/16+y**2/4))
In [101]:
fig = plt.figure() 
ax = plt.axes(projection ="3d") 

X = np.linspace(-10,10,1000)
Y = np.linspace(-10,10,1000)

X,Y = np.meshgrid(X,Y)

Z = func3D(X,Y)

ax.contour3D(X,Y,Z,500,cmap='inferno')

ax.view_init(60, 55)
In [102]:
fig = plt.figure() 
ax = plt.axes(projection ="3d") 

X = np.linspace(-10,10,100)
Y = np.linspace(-10,10,100)

X,Y = np.meshgrid(X,Y)

Z = func3D(X,Y)

ax.contour3D(X,Y,Z,500,cmap='plasma')

ax.view_init(60, 55)
In [109]:
fig = plt.figure() 
ax = plt.axes(projection ="3d") 
ax.view_init(60, 55)

X = np.linspace(-10,10,100)
Y = np.linspace(-10,10,100)

X,Y = np.meshgrid(X,Y)

Z = func3D(X,Y)


ax.plot_surface(X,Y,Z,cmap='viridis',edgecolor='none',rstride=1, cstride=1)
Out[109]:
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7fdd903b5990>
In [113]:
fig = plt.figure() 
ax = plt.axes(projection ="3d")
ax.view_init(60, 55)

X = np.linspace(-10,10,100)
Y = np.linspace(-10,10,100)

X,Y = np.meshgrid(X,Y)

Z = func3D(X,Y)

ax.plot_wireframe(X,Y,Z,color='black',alpha=0.1)

ax.set_title('wireframe')
Out[113]:
Text(0.5, 0.92, 'wireframe')
In [116]:
fig = plt.figure() 
ax = plt.axes(projection ="3d")
ax.view_init(60, 55)

X = np.linspace(-10,10,100)
Y = np.linspace(-10,10,100)

X,Y = np.meshgrid(X,Y)

Z = func3D(X,Y)

ax.scatter(X,Y,Z,alpha=0.1)
Out[116]:
<mpl_toolkits.mplot3d.art3d.Path3DCollection at 0x7fddd110b450>
In [128]:
fig = plt.figure() 
ax = plt.axes(projection ="3d")
ax.view_init(60, 55)
X = np.linspace(-10,10,100)
Y = np.linspace(-10,10,100)

X,Y = np.meshgrid(X,Y)

X=np.ravel(X)
Y=np.ravel(Y)

Z = func3D(X,Y)

ax.plot_trisurf(X,Y,Z)
Out[128]:
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7fddd0af1b90>
In [121]:
X = np.linspace(-2,2,5)
Y = np.linspace(-2,2,5)
In [122]:
Y
Out[122]:
array([-2., -1.,  0.,  1.,  2.])
In [123]:
X,Y = np.meshgrid(X,Y)
In [125]:
Y
Out[125]:
array([[-2., -2., -2., -2., -2.],
       [-1., -1., -1., -1., -1.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 1.,  1.,  1.,  1.,  1.],
       [ 2.,  2.,  2.,  2.,  2.]])
In [126]:
X=np.ravel(X)
In [127]:
X
Out[127]:
array([-2., -1.,  0.,  1.,  2., -2., -1.,  0.,  1.,  2., -2., -1.,  0.,
        1.,  2., -2., -1.,  0.,  1.,  2., -2., -1.,  0.,  1.,  2.])
In [130]:
u = np.linspace(0, 2.0 * np.pi, endpoint=True, num=500) 
v = np.linspace(-0.5, 0.5, endpoint=True, num=100)

u, v = np.meshgrid(u, v)
u, v = u.flatten(), v.flatten()

x = (1 + 0.5 * v * np.cos(u / 2.0)) * np.cos(u) 
y = (1 + 0.5 * v * np.cos(u / 2.0)) * np.sin(u) 
z = 0.5 * v * np.sin(u / 2.0)

ax = plt.axes(projection='3d')
fig = plt.figure()

ax.scatter(x, y, z, c=z, cmap='viridis', linewidth=0.1)

ax.view_init(60, 55)

<Figure size 432x288 with 0 Axes>
In [ ]: