How to Integrate Using Simpsons Rule in Python?

Question

I am trying to integrate this function: x^4 - 2x + 1 from 0 to 2

I wrote this program:

def f(x):
    return (x**4)-(2*x)+1

N=10 
a=0.0
b=2.0
h=(b-a)/N

s=f(a)+f(b)

for k in range(1,N/2):
    s+=4*f(a+(2*k-1)*h)

for k in range(1,N/(2-1)):
    s1 +=f(a+(2*k*h)

M=(s)+(2*s1)
print((1/3.0)*h)*(3)

But I got this error:

File "<ipython-input-29-6107592420b6>", line 17
   M=(s)+(2*s1):
   ^
SyntaxError: invalid syntax

I tried writing it in different forms but I always get an error in M

BTW, your Simpson's rule formula isn't quite right. My old answer here may be helpful http://stackoverflow.com/a/33715116/4014959 — PM 2Ring, Apr 17 '16 at 16:51
You might want to declare `s1=0` before using it in the loop. Or use `s += 2*f(...)`. — Lutz Lehmann, Apr 17 '16 at 17:06

score 4 · Answer 1 · answered Apr 17 '16 at 16:22

4

You forgot a closing parentheses in your second for loop here: s1 += f(a+(2*k*h). It should be:

s1 += f(a + (2 * k * h)) # <<< here it is

answered Apr 17 '16 at 16:22

ForceBru

43,482
10
63
98

2

Just to add - Typically when you get an invalid syntax error, it is *usually* the line above that is the issue. – idjaw Apr 17 '16 at 16:39

Joe Gavin · Answer 2 · 2016-04-17T17:49:51.377

For future reference you might also think about using scipy.integrate. Look here for some methods which might have better accuracy depending on the nature and resolution of your data set.

A code might look like this:

import scipy.integrate as int
x = [ ii/10. for ii in range(21)]
y = [ xi**4 - 2*xi + 1 for xi in x]
tahdah = int.simps(y,x,even='avg')
print(tahdah)

Which yields and answer of 4.4 that you can confirm with pencil and paper.

score 1 · Answer 3 · answered Jul 29 '16 at 12:47

Have you seen the code example of Simpsons Rule on Wikipedia (written in python)? I will repost it here for the benefit of future readers.

#!/usr/bin/env python3
from __future__ import division  # Python 2 compatibility

def simpson(f, a, b, n):
    """Approximates the definite integral of f from a to b by the
    composite Simpson's rule, using n subintervals (with n even)"""

    if n % 2:
        raise ValueError("n must be even (received n=%d)" % n)

    h = (b - a) / n
    s = f(a) + f(b)

    for i in range(1, n, 2):
        s += 4 * f(a + i * h)
    for i in range(2, n-1, 2):
        s += 2 * f(a + i * h)

    return s * h / 3

# Demonstrate that the method is exact for polynomials up to 3rd order
print(simpson(lambda x:x**3, 0.0, 10.0, 2))       # 2500.0
print(simpson(lambda x:x**3, 0.0, 10.0, 100000))  # 2500.0

print(simpson(lambda x:x**4, 0.0, 10.0, 2))       # 20833.3333333
print(simpson(lambda x:x**4, 0.0, 10.0, 100000))  # 20000.0

How to Integrate Using Simpsons Rule in Python?

3 Answers3