Simpson Integration to Python : Can't get the output as expected

125 Views Asked by At

I have already implement this even() & simpson(), but can't get the desired output. Secondly getting error ([lambda x:1/x, 1, 11, 6]]) ZeroDivisionError: division by zero, can't understand this error

import math
import matplotlib
import matplotlib.pyplot as plt
import numpy as np

def even(x):
    if x % 2 == 0:
        return 2
    else:
        return 4

def simpson(fn,a,b,n):
    ans = 0

    for x in range(1, n):
        ans = ans + even(x) * fn(x)
    return ((b - a) / (3 * n)) * (ans + fn(0) + fn(n))

    

if __name__ == "__main__":
    """
    The code in "__main__" is not being graded, but a tool for you to test 
    your code outside of the `test_a8.py`. Feel free to add print statements. 
    """

    data = [[lambda x:3*(x**2)+1, 0,6,2],
            [lambda x:x**2,0,5,6],
            [lambda x:math.sin(x), 0,math.pi, 4],
            [lambda x:1/x, 1, 11, 6]]

    for d in data:
        f,a,b,n = d
        print(simpson(f,a,b,n))

    t = np.arange(0.0, 10.0,.1)
    fig,ax = plt.subplots()
    s = np.arange(0,6.1,.1)
    ax.plot(t, (lambda t: 3*(t**2) + 1)(t),'g')
    plt.fill_between(s,(lambda t: 3*(t**2) + 1)(s))
    ax.grid()
    ax.set(xlabel ="x", ylabel=r"$f(x)=3x^2 + 1$",
        title = r"Area under the curve $\int_0^6\,f(x)$")

    plt.show()

Expected Output I am trying for like this

222.0
41.66666666666667
2.0045597549844207
2.4491973405016885
1

There are 1 best solutions below

0
STerliakov On

Do I properly interpret third number as total points count, including ends? If yes - the following works. You made an even function too heavy (no if's are requred, only integer evaluation). Main problem of your code is that you loop x over integer points only and start from zero (parameters a and b should be bounds, of course).

import math
import numpy as np

def simpson(fn,a,b,n):
    # n = n + 2 # uncomment this line if only inner points are included in n
    ans = sum((2-i%2)*fn(x) for i,x in enumerate(np.linspace(a,b,n)[1:-1]))
    return ((b - a) / (3 * n)) * (fn(a) + 2*ans + fn(b))

data = [[lambda x:3*(x**2)+1, 0,6,2],
        [lambda x:x**2,0,5,6],
        [lambda x:math.sin(x), 0,math.pi, 4],
        [lambda x:1/x, 1, 11, 6]]

for d in data:
    f,a,b,n = d
    print(simpson(f,a,b,n))

Don't expect the output you provided. You should specify at least 30-50 points per range to get result close to accurate, 2-6 is too small. If you use 100's of points, result will be close to expected.