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

# # Composite corrected Trapezoid rule

# In[1]:


import numpy as np

# n = number of intervals
# There are n+1 points
def trapz(a,b,n,f,df):
    h = (b-a)/n
    x = np.linspace(a,b,n+1)
    y = f(x)
    res = np.sum(y[1:n]) + 0.5*(y[0] + y[n])
    return h*res, h*res - (h**2/12)*(df(b) - df(a))


# ## Example
# $$
# f(x) = \exp(x)\cos(x), \qquad x \in [0,\pi]
# $$
# The exact integral is $-\frac{1}{2}(1+\exp(\pi))$.

# In[2]:


f  = lambda x: np.exp(x)*np.cos(x)
df = lambda x: np.exp(x)*(np.cos(x) - np.sin(x))
qe = -0.5*(1.0 + np.exp(np.pi)) # Exact integral

n,N = 4,10
e1,e2 = np.zeros(N),np.zeros(N)
for i in range(N):
    e1[i],e2[i] = trapz(0.0,np.pi,n,f,df) - qe
    if i > 0:
        print('%6d %24.14e %10.5f %24.14e %10.5f'%(n,e1[i],e1[i-1]/e1[i],e2[i],e2[i-1]/e2[i]))
    else:
        print('%6d %24.14e %10.5f %24.14e %10.5f'%(n,e1[i],0,e2[i],0))
    n = 2*n