#-----------------------------------------------------------
# Lecture 14 -- Example -- Trapezoidal Rule
#-----------------------------------------------------------

import numpy as np

# Use the composite trapezoidal rule to evaluate the integral:
#   I = Int ( sin(pi*x)*cos(pi*x)^2 ), from x = 0 to x = 1
# with N = 10 points.
#
# The exact result is 2/(3 pi). How accurate is your result?
# How does it change when you use N = 20 points?

# Steps to using the trapezoidal rule
# (1) Define a function for the integrand
# (2) Make an array for the x points and calculate dx
# (3) Make a loop to calculate the f(x[i])
# (4) Sum the integral using the trapezoidal rule formula:
#       I = dx/2(f0 + 2*f1 + 2*f2 + ... + 2*fn-1 + fn)

def f(x):
    return np.sin(np.pi*x)*np.cos(np.pi*x)**2

N = 10 # number of points
x = np.linspace(0, 1, N)
dx = x[1]-x[0]

I = 0
for i in range(N):
    if (i > 0 and i < N-1):
        I += 2*f(x[i])
    else:
        I += f(x[i])
        
I *= dx/2

print('Exact value = ', 2./(3*np.pi))
print('Trapezoidal rule = ', I)
print('Error = ', np.abs(2./(3*np.pi) - I))