import numpy as np
import matplotlib.pyplot as plt

# Given a tuple quad = (xs, coeff) that describes a quadrature formula for the
# interval [0,1] uses that chained to calculate the integral of the given
# function f over the interval [a,b], iterating quad a total of n times.
def chainedQuad(quad, a, b, n, f):
    xs, coeff = quad
    bma = (b - a) / n

    points = np.linspace(a, b, n+1)[:-1] + bma * np.array(xs, ndmin=2).T
    return np.sum(bma * (coeff * f(points.T)))

def chainedTrapezoid(a, b, n, f):
    xs, coeff = [0, 1], [1/2, 1/2]
    return chainedQuad((xs, coeff), a, b, n, f)

def chainedSimpson(a, b, n, f):
    xs, coeff = [0, 1/2, 1], [1/6, 2/3, 1/6]
    return chainedQuad((xs, coeff), a, b, n, f)

def chainedMilne(a, b, n, f):
    xs, coeff = [0, 1/4, 1/2, 3/4, 1], [7/90, 32/90, 12/90, 32/90, 7/90]
    return chainedQuad((xs, coeff), a, b, n, f)

if __name__ == "__main__":
    f = lambda x: np.cos(x - 1)**2 + np.sin(0.7 * x) + x
    a, b = -1, 2 * np.pi

    # Integrate[f[x], {x, a, b}] ~ 23.998
    exactResult = (-10 + 10*np.sqrt(5) + 28*np.pi + 56*np.pi**2 + 7*(np.sin(4) - np.sin(2)) + 40*np.cos(7/10))/28

    err = lambda n: np.abs(exactResult - np.array([chainedTrapezoid(a, b, n, f),
                                                   chainedSimpson(a, b, n, f),
                                                   chainedMilne(a, b, n, f)]))

    ns = 5 * (np.arange(30) + 1)
    plt.semilogy(ns, [err(n) for n in ns], label=["Trapez-Regel", "Simpson-Regel", "Milne-Regel"])
    plt.ylabel("Approximationsfehler")
    plt.xlabel("Anzahl der Iterationen der Quadraturformel")
    plt.legend()
    plt.savefig("image.png")