Tukukkk

In the following code I have implemented composite Simpsons Rule in python. I have tested it by integrating f = sinx over [0,pi/2] and plotting the resulting absolute error as a function of n for a suitable range of integer values for n. As shown below. Now I am trying to verify that my method is of order 4. To do this instead of plotting the error vs n I need to plot n vs n^(-4) to show that both have the same slope.

from math import pi, cos, sin

from matplotlib import pyplot as plt





def simpson(f, a, b, n):

    """Approximates the definite integral of f from a to b by the composite Simpson's rule, using 2n subintervals """



    h = (b - a) / (2*n)

    s = f(a) + f(b)





    for i in range(1, 2*n, 2):

        s += 4 * f(a + i * h)

    for i in range(2, 2*n-1, 2):

        s += 2 * f(a + i * h)

    return s * h / 3





diffs = {}

exact = 1 - cos(pi/2)

for n in range(1, 100):

    result = simpson(lambda x: sin(x), 0.0, pi/2, n)

    diffs[2*n] = abs(exact - result)   # use 2*n or n here, your choice.





ordered = sorted(diffs.items())

x,y = zip(*ordered)

plt.autoscale()

plt.loglog(x,y)

plt.xlabel("Intervals")

plt.ylabel("Error")

plt.show()

this results in my error graph:

You can add another line to your plot by making another call to plt.loglog(). So just before your plt.show() you can add:

n = 

n_inv_4 = 

for i in range(1,100):

    n.append(2*i)

    n_inv_4.append(1.0 / ((2*i)**4))

plt.loglog(n, n_inv_4)

Note that the above calculations use (2*i) to match the (2*n) used in the simpson method.