#!/usr/bin/env python
# coding: utf-8

# In[1]:


#simple trapezoid rule
trpz = lambda f, a, b: (b-a) * (f(a) + f(b))/2


# In[2]:


#example integral
from numpy import exp
f = lambda x: exp(x)
a = 0
b = 1


# In[3]:


trpz(f, a, b)


# In[4]:


#composite trapezoid rule with n equal segments
def trpz_c(f, a, b, n):
    from numpy import arange, sum
    return (b - a) * (f(a) + 2*sum(f(a + (b-a)*arange(1, n)/n)) + f(b)) / (2*n)


# In[5]:


trpz_c(f, a, b, 1)


# In[6]:


trpz_c(f, a, b, 10)


# In[7]:


trpz_c(f, a, b, 100)


# In[11]:


#composite trapezoid rule with 2n equal segments
#if the one with n segments was already calculated (Tn)
def trpz_c2(f, a, b, n, Tn):
    from numpy import arange, sum
    return (Tn + (b - a) * (sum(f(a + (b-a)*arange(1, 2*n, 2)/(2*n)))) / n) / 2


# In[14]:


T5 = trpz_c(f, a, b, 5)
T10 = trpz_c(f, a, b, 10)
T10_2 = trpz_c2(f, a, b, 5, T5)


# In[15]:


T10, T10_2 #should be the same except for roundoff error


# In[17]:


#Romberg integration
from numpy import zeros, nan
def romberg(f, a, b, n):
    R = zeros((n, n))
    R[:] = nan
    for i in range(n):
        if i > 0:
           R[i, 0] = trpz_c2(f, a, b, 2**(i-1), R[i-1, 0])
        else: #for i = 0 (simple trapezoid rule)
           R[0, 0] = (b-a) * (f(a) + f(b))/2 
        for j in range(1, n):
            R[i-j, j] = (4**j/(4**j-1))*R[i-j+1, j-1] - (1/(4**j-1))*R[i-j, j-1]
    return R  


# In[19]:


R = romberg(f, a, b, 4)


# In[20]:


R - (exp(1) - 1)


# In[28]:


#2-point Gauss quadrature
from numpy import sqrt
w = [1, 1]
x = [-1, 1] / sqrt(3)
((b-a)/2) * sum (w * f(a + (x+1)*(b-a)/2))


# In[32]:


#3-point Gauss quadrature
from numpy import array
w = array([5, 8, 5]) / 9.
x = array([-0.774596669, 0, 0.774596669])
((b-a)/2) * sum (w * f(a + (x+1)*(b-a)/2))


# In[35]:


#using the SciPy function
from scipy.integrate import fixed_quad
fixed_quad(f, a, b, n=3)[0]


# In[40]:


#composite trapezoid rule
from scipy.integrate import trapezoid
from numpy import linspace
x = linspace(a, b, 11)
y = f(x)
trapezoid(y, x)


# In[42]:


from scipy.integrate import quad
I, err = quad(f, a, b)


# In[45]:


#examples of double integration from SciPy documentation
from scipy.integrate import dblquad
f = lambda y, x: x*y**2
dblquad(f, 0, 2, 0, 1)


# In[46]:


from numpy import pi, sin, cos
f = lambda y, x: 1
dblquad(f, 0, pi/4, sin, cos)


# In[ ]: