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

# In[8]:


#apply the composite trapezoid rule over a set of points x and function values f(x) = y
def trapezoid(x, y):
    from numpy import sum
    w = x[1:] - x[:-1]
    a = (y[1:] + y[:-1])/2
    return sum(w * a)
    


# In[9]:


from numpy import array
x = array([2, 3, 5, 7, 9])
y = array([1, 2, 4, 5, 4])
trapezoid(x, y)


# In[5]:


from


# In[6]:


x[0:-1]


# In[11]:


from numpy import sqrt
f = lambda x: sqrt(x)
a = 1; b = 2


# In[14]:


x = array([a, b])
y = f(x)
T1 = trapezoid(x, y)
T1


# In[15]:


x = array([a, (a+b)/2, b])
y = f(x)
T2 = trapezoid(x, y)
T2


# In[16]:


M = (b-a) * f((a+b)/2) #simple midpoint rule
M


# In[17]:


S = (b-a) * (f(a) + 4*f((a+b)/2) + f(b))/6 #simple Simpson rule
S


# In[18]:


(4*T2 - T1)/3 #same as Simpson rule


# In[19]:


from numpy import linspace
n = 100
x = linspace(a, b, n+1)
y = f(x)
Tn = trapezoid(x, y)


# In[21]:


#Romberg integration
def romberg(f, a, b, jmax):
    from numpy import array, zeros, linspace
    A = zeros((jmax+1, jmax+1))
    for i in range(jmax+1):
        n = 2**i
        x = linspace(a, b, n+1)
        y = f(x)
        A[i, 0] = trapezoid(x, y)
        for j in range(1, i+1):
            A[i-j, j] = ((4**j) * A[i-j+1, j-1] - A[i-j, j-1]) / (4**j - 1) 
    print(A)
    return A[0, jmax]


# In[22]:


romberg(f, a, b, 3)


# In[23]:


#analytic integral for the example function
(2/3) * (sqrt(8) - 1)


# In[25]:


#Gauss quadrature with n = 2
w = array([1, 1])
x = array([-1, 1]) / sqrt(3)
G = ((b-a)/2) * sum (w * f(a + ((b-a)/2)*(x + 1)))


# In[26]:


G


# In[27]:


#Gauss quadrature with n = 3
w = array([5, 8, 5]) / 9.
x = array([-1, 0, 1]) * sqrt(0.6)
G = ((b-a)/2) * sum (w * f(a + ((b-a)/2)*(x + 1)))


# In[28]:


G


# In[30]:


#same thing, with a built-in function
from scipy.integrate import fixed_quad
fixed_quad(f, a, b, n=3)[0]


# In[31]:


#general-purpose adaptive quadrature
from scipy.integrate import quad
quad(f, a, b)


# In[ ]: