#! /usr/bin/env python
"""
 Simpson.py is a program  which calculate the integral 
 of a function f(x) = x**4 - 2x +1 utilizing the Simpson's rule.

 Paul Eugenio
 PHZ4151C
 13 Feb 2019
"""
from __future__ import division, print_function
import numpy as np


 

def  simIntegrate(f, a, b, N):
    """
    Intgrate the function f(x) where a,b are the 
    lower,higher limits using the Simpson's rule for N steps

    Note this implimentation uses a 3 points over one step implimentation
    which is valid for odd & even number of steps where as the 
    3 points over two steps implimentations requires an even number of steps.
    
    This approach is like half stepping making the effective number of steps N
    equiliviant to 2*N.
    """
    h = (b-a)/N
    s1, s2 = 0.0, 0.0
    for k in range(0, N):
        if not k == 0:
            s1 += f(a+k*h)
        s2 += f(a+(k+0.5)*h)
    
    return (f(a) + f(b) + 2*s1 + 4*s2)*h/6

def  simIntegrate2(f, a, b, N):
    """
    This is the Simpson's rule implimentation following the 
    prescription given in Newman's textbook.  The results from
    this integration with N steps will give identical results
    to the function simIntegrate above but with N/2 steps.
    """
    
    h = (b-a)/N
    s1 = 0.0
    for k in range(1, N, 2):
        s1 += f(a+k*h)

    s2 = 0.0
    for k in range(2, N, 2):
        s2 += f(a+k*h)

    return (f(a) + f(b) + 4*s1 + 2*s2)*h/3



#
# main 
   
def fx(x):
    """ integrand function """
    return x**4 - 2*x + 1

N = 100        # number of steps
a = 0.0        # lower x limit
b = 2.0        # upper x limit
h = (b-a)/N    # step size
actualValue = 4.4

print("Nsteps\t I_Simp1\t\t I_Simp2\t\t Fractional Error")
for N in [10, 100, 1000]:
    I = simIntegrate(fx, a, b, N)
    I2 = simIntegrate2(fx, a, b, 2*N)    
    fractionalError = 100*abs(actualValue - I)/actualValue
 
    print("{0}\t {1:2.10e}\t {2:2.10e}\t {3:e}%".format(N, I, I2, fractionalError), sep="")

# Results
#
# physics 251% simpson.py
# Nsteps   I_Simp1         I_Simp2         Fractional Error
# 10   4.4000266667e+00    4.4000266667e+00    6.060606e-04%
# 100  4.4000000027e+00    4.4000000027e+00    6.060607e-08%
# 1000 4.4000000000e+00    4.4000000000e+00    6.075948e-12%
#
# from example 5.1:
#
# 10      4.450656     2%
# 100     4.40107      0.02%
# 1000    4.40001      0.0002%
#
# This shows that the error for the trapezodial rule is proportional 
# to h**2 or 1/N**2 (i.e. h=(b-a)/N) where as for the Simpson's rule 
# the error goes like h**4 or 1/N**4
#