# coding: utf-8

# # Accuracy of Simpson's Rule

# In[1]:

import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as pt


# In[2]:

simpson_weights = np.array([1/6, 4/6, 1/6])
simpson_weights


# Here is a function and its definite integral:
# 
# $$\text{int_f}(x)=\int_0^x f(\xi)d\xi$$

# In[1]:

alpha = 4

def f(x):
    return np.cos(alpha*x)
def int_f(x):
    return 1/alpha*(np.sin(alpha*x)-np.sin(alpha*0))


# Plotted:

# In[53]:

plot_x = np.linspace(0, 1, 200)

pt.plot(plot_x, f(plot_x), label="f")
pt.fill_between(plot_x, 0*plot_x, f(plot_x),alpha=0.3)
pt.plot(plot_x, int_f(plot_x), label="$\int f$")
pt.grid()
pt.legend(loc="best")


# This here plots the function, the interpolant, and the area under the interpolant:

# In[61]:

# fix nodes
h = 1
x = np.linspace(0, h, 3)

# find interpolant
V = np.array([
  1+0*x,
  x,
  x**2
]).T
a, b, c = la.solve(V, f(x))
interpolant = a + b*plot_x + c*plot_x**2

# plot
pt.plot(plot_x, f(plot_x), label="f")
pt.plot(plot_x, interpolant, label="Interpolant")
pt.fill_between(plot_x, 0*plot_x, interpolant, alpha=0.3, color="green")
pt.plot(x, f(x), "og")
pt.grid()
pt.legend(loc="best")


# In[55]:

true_val = int_f(h)
error = (f(x).dot(simpson_weights) * h - true_val)/true_val
print(error)


# * Compare the error for $h=1,0.5,0.25$. What order of accuracy do you observe?