Simpson's rule in Python

Question

For a numerical methods class, I need to write a program to evaluate a definite integral with Simpson's composite rule. I already got this far (see below), but my answer is not correct. I am testing the program with f(x)=x, integrated over 0 to 1, for which the outcome should be 0.5. I get 0.78746... etc. I know there is a Simpson's rule available in Scipy, but I really need to write it myself.

I suspect there is something wrong with the two loops. I tried "for i in range(1, n, 2)" and "for i in range(2, n-1, 2)" before, and this gave me a result of 0.41668333... etc. I also tried "x += h" and I tried "x += i*h". The first gave me 0.3954, and the second option 7.9218.

# Write a program to evaluate a definite integral using Simpson's rule with
# n subdivisions

from math import *
from pylab import *

def simpson(f, a, b, n):
    h=(b-a)/n
    k=0.0
    x=a
    for i in range(1,n/2):
        x += 2*h
        k += 4*f(x)
    for i in range(2,(n/2)-1):
        x += 2*h
        k += 2*f(x)
    return (h/3)*(f(a)+f(b)+k)

def function(x): return x

print simpson(function, 0.0, 1.0, 100)

Answer 1

You probably forget to initialize x before the second loop, also, starting conditions and number of iterations are off. Here is the correct way:

def simpson(f, a, b, n):
    h=(b-a)/n
    k=0.0
    x=a + h
    for i in range(1,n/2 + 1):
        k += 4*f(x)
        x += 2*h

    x = a + 2*h
    for i in range(1,n/2):
        k += 2*f(x)
        x += 2*h
    return (h/3)*(f(a)+f(b)+k)

Your mistakes are connected with the notion of a loop invariant. Not to get into details too much, it's generally easier to understand and debug cycles which advance at the end of a cycle, not at the beginning, here I moved the x += 2 * h line to the end, which made it easy to verify where the summation starts. In your implementation it would be necessary to assign a weird x = a - h for the first loop only to add 2 * h to it as the first line in the loop.

Answer 2

All you need to do to make this code work is add a variable for a and b in the function bounds() and add a function in f(x) that uses the variable x. You could also implement the function and bounds directly into the simpsonsRule function if desired... Also, these are functions to be implimented into a program, not a program itself.

def simpsonsRule(n):

    """
    simpsonsRule: (int) -> float
    Parameters:
        n: integer representing the number of segments being used to 
           approximate the integral
    Pre conditions:
        Function bounds() declared that returns lower and upper bounds of integral.
        Function f(x) declared that returns the evaluated equation at point x.
        Parameters passed.
    Post conditions:
        Returns float equal to the approximate integral of f(x) from a to b
        using Simpson's rule.
    Description:
        Returns the approximation of an integral. Works as of python 3.3.2
        REQUIRES NO MODULES to be imported, especially not non standard ones.
        -Code by TechnicalFox
    """

    a,b = bounds()
    sum = float()
    sum += f(a) #evaluating first point
    sum += f(b) #evaluating last point
    width=(b-a)/(2*n) #width of segments
    oddSum = float()
    evenSum = float()
    for i in range(1,n): #evaluating all odd values of n (not first and last)
        oddSum += f(2*width*i+a)
    sum += oddSum * 2
    for i in range(1,n+1): #evaluating all even values of n (not first and last)
        evenSum += f(width*(-1+2*i)+a)
    sum += evenSum * 4
    return sum * width/3

def bounds():
    """
    Description:
        Function that returns both the upper and lower bounds of an integral.
    """
    a = #>>>INTEGER REPRESENTING LOWER BOUND OF INTEGRAL<<<
    b = #>>>INTEGER REPRESENTING UPPER BOUND OF INTEGRAL<<<
    return a,b

def f(x):
    """
    Description:
        Function that takes an x value and returns the equation being evaluated,
        with said x value.
    """
    return #>>>EQUATION USING VARIABLE X<<<

Answer 3

def simps(f, a, b, N):   # N must be an odd integer
    """ define simpson method, a and b are the lower and upper limits of
    the interval, N is number of points, dx is the slice
        """
    integ = 0
    dx = float((b - a) / N)

    for i in range(1,N-1,2):
        integ += f((a+(i-1)*dx)) + 4*f((a+i*dx)) + f((a+(i+1)*dx))
        integral = dx/3.0 * integ
    
    # if number of points is even, then error araise
    if (N % 2) == 0:
      raise ValueError("N must be an odd integer.")
        

    return integral


def f(x):
    return x**2


integrate = simps(f,0,1,99)
print(integrate)

Answer 4

You can use this program for calculating definite integrals by using Simpson's 1/3 rule. You can increase your accuracy by increasing the value of the variable panels.

import numpy as np

def integration(integrand,lower,upper,*args):    
    panels = 100000
    limits = [lower, upper]
    h = ( limits[1] - limits[0] ) / (2 * panels)
    n = (2 * panels) + 1
    x = np.linspace(limits[0],limits[1],n)
    y = integrand(x,*args)
    #Simpson 1/3
    I = 0
    start = -2
    for looper in range(0,panels):
        start += 2
        counter = 0
        for looper in range(start, start+3):
            counter += 1
            if (counter ==1 or counter == 3):
                I += ((h/3) * y[looper])
            else:
                I += ((h/3) * 4 * y[looper])
    return I

For example:

def f(x,a,b):
    return a * np.log(x/b)
I = integration(f,3,4,2,5)
print(I)

will integrate 2ln(x/5) within the interval 3 and 4

Answer 5

There is my code (i think that is the most easy method). I done this in jupyter notebook. The easiest and most accurate code for Simpson method is 1/3.

Explanation

For standard method (a=0, h=4, b=12) and f=100-(x^2)/2

We got: n= 3.0, y0 = 100.0, y1 = 92.0, y2 = 68.0, y3 = 28.0,

So simpson method = h/3*(y0+4*y1+2*y2+y3) = 842,7 (this is not true). Using 1/3 rule we got:

h = h/2= 4/2= 2 and then:

n= 3.0, y0 = 100.0, y1 = 98.0, y2 = 92.0, y3 = 82.0, y4 = 68.0, y5 = 50.0, y6 = 28.0,

Now we calculate the integral for each step (n=3 = 3 steps):

Integral of the first step: h/3*(y0+4*y1+y2) = 389.3333333333333

Integral of the second step: h/3*(y2+4*y3+y4) = 325.3333333333333

Integral of the third step: h/3*(y4+4*y5+y6) = 197.33333333333331

Sum all, and we get 912.0 AND THIS IS TRUE

x=0
b=12
h=4
x=float(x)
h=float(h)
b=float(b)
a=float(x)
def fun(x): 
    return 100-(x**2)/2
h=h/2
l=0  #just numeration
print('n=',(b-x)/(h*2))
n=int((b-a)/h+1)
y = []   #tablica/lista wszystkich y / list of all "y"
yf = []
for i in range(n):
    f=fun(x)
    print('y%s' %(l),'=',f)
    y.append(f)
    l+=1
    x+=h
print(y,'\n')
n1=int(((b-a)/h)/2)  
l=1
for i in range(n1):
    nf=(h/3*(y[+0]+4*y[+1]+y[+2]))
    y=y[2:]  # with every step, we deleting 2 first "y" and we move 2 spaces to the right, i.e. y2+4*y3+y4
    print('Całka dla kroku/Integral for a step',l,'=',nf)
    yf.append(nf)
    l+=1
print('\nWynik całki/Result of the integral =', sum(yf) )

At the beginning I added simple data entry:

d=None
while(True):
    print("Enter your own data or enter the word "test" for ready data.\n")
    x=input ('Enter the beginning of the interval (a): ') 
    if x == 'test':
        x=0
        h=4  #krok (Δx)
        b=12 #granica np. 0>12  
        #w=(20*x)-(x**2)  lub   (1+x**3)**(1/2)
        break
    h=input ('Enter the size of the integration step (h): ')
    if h == 'test':
        x=0
        h=4 
        b=12 
        break
    b=input ('Enter the end of the range (b): ')
    if b == 'test':
        x=0
        h=4  
        b=12 
        break 
    d=input ('Give the function pattern: ')
    if d == 'test':
        x=0
        h=4  
        b=12
        break
    elif d != -9999.9:
        break

x=float(x)
h=float(h)
b=float(b)
a=float(x)

if d == None or d == 'test':
    def fun(x): 
        return 100-(x**2)/2 #(20*x)-(x**2)
else:
    def fun(x): 
        w = eval(d)
        return  w
h=h/2
l=0  #just numeration
print('n=',(b-x)/(h*2))
n=int((b-a)/h+1)
y = []   #tablica/lista wszystkich y / list of all "y"
yf = []
for i in range(n):
    f=fun(x)
    print('y%s' %(l),'=',f)
    y.append(f)
    l+=1
    x+=h
print(y,'\n')
n1=int(((b-a)/h)/2)  
l=1
for i in range(n1):
    nf=(h/3*(y[+0]+4*y[+1]+y[+2]))
    y=y[2:]
    print('Całka dla kroku/Integral for a step',l,'=',nf)
    yf.append(nf)
    l+=1
print('\nWynik całki/Result of the integral =', sum(yf) )

Answer 6

Example of implementing Simpson's rule for integral sinX with a = 0 and b = pi/4. And use 10 panels for the integration

def simpson(f, a, b, n):
  x = np.linspace(a, b, n+1)
  w = 2*np.ones(n+1); w[0] = 1.0; w[-1] = 1;
  for i in range(len(w)):
      if i % 2 == 1:
          w[i] = 4
  width = x[1] - x[0]
  area = 0.333 * width * np.sum( w*f(x))
  return area

f = lambda x: np.sin(x)
a = 0.0; b = np.pi/4

areaSim = simpson(f, a, b, 10)
print(areaSim)