CODEWITHSUNDEEP
pythonnumerical-integration
user2763502
user2763502

Reputation: 11

Trapezoid Rule in Python

I am trying to write a program using Python v. 2.7.5 that will compute the area under the curve y=sin(x) between x = 0 and x = pi. Perform this calculation varying the n divisions of the range of x between 1 and 10 inclusive and print the approximate value, the true value, and the percent error (in other words, increase the accuracy by increasing the number of trapezoids). Print all the values to three decimal places.

I am not sure what the code should look like. I was told that I should only have about 12 lines of code for these calculations to be done.

I am using Wing IDE.

This is what I have so far

# base_n = (b-a)/n
# h1 = a + ((n-1)/n)(b-a)
# h2 = a + (n/n)(b-a)
# Trap Area = (1/2)*base*(h1+h2)
# a = 0, b = pi

from math import pi, sin

def TrapArea(n):
for i in range(1, n):
    deltax = (pi-0)/n
    sum += (1.0/2.0)(((pi-0)/n)(sin((i-1)/n(pi-0))) + sin((i/n)(pi-0)))*deltax
    return sum

for i in range(1, 11):
    print TrapArea(i)

I am not sure if I am on the right track. I am getting an error that says "local variable 'sum' referenced before assignment. Any suggestions on how to improve my code?

Upvotes: 1

Views: 8746

Answers (3)

Mike Housky
Mike Housky

Reputation: 4069

A "nicer" way to do the trapezoid rule with equally-spaced points...

Let dx = pi/n be the width of the interval. Also, let f(i) be sin(i*dx) to shorten some expressions below. Then interval i (in range(1,n)) contributes:

dA = 0.5*dx*( f(i) + f(i-1) )

...to the sum (which is an area, so I'm using dA for "delta area"). Factoring out the 0.5*dx, makes the whole some look like:

A = 0.5*dx * ( (f(0) + f(1)) + (f(1) + f(2)) + .... + (f(n-1) + f(n)) )

Notice that there are two f(1) terms, two f(2) terms, on up to two f(n-1) terms. Combine those to get:

A = 0.5*dx * ( f(0) + 2*f(1) + 2*f(2) + ... + 2*f(n-1) + f(n) )

The 0.5 and 2 factors cancel except in the first and last terms:

A = 0.5*dx(f(0) + f(n)) + dx*(f(1) + f(2) + ... + f(n-1))

Finally, you can factor dx out entirely to do just one multiplication at the end. Converting back to sin() calls, then:

def TrapArea(n):
    dx = pi/n
    asum = 0.5*(sin(0) + sin(pi))   # this is 0 for this problem, but not others
    for i in range(1, n-1):
        asum += sin(i*dx)
    return sum*dx

That changed "sum" to "asum", or maybe "area" would be better. That's mostly because sum() is a built-in function, which I'll use below the line.

Extra credit: The loop part of the sum can be done in one step with a generator expression and the sum builtin function:

def TrapArea2(n):
    dx = pi/n
    asum = 0.5*(sin(0) + sin(pi))
    asum += sum(sin(i*dx) for i in range(1,n-1))
    return asum*dx

Testing both of those:

>>> for n in [1, 10, 100, 1000, 10000]:
        print n, TrapArea(n), TrapArea2(n)

1 1.92367069372e-16 1.92367069372e-16
10 1.88644298557 1.88644298557
100 1.99884870579 1.99884870579
1000 1.99998848548 1.99998848548
10000 1.99999988485 1.99999988485

That first line is a "numerical zero", since math.sin(math.pi) evaluates to about 1.2e-16 instead of exactly zero. Draw the single interval from 0 to pi and the endpoints are indeed both 0 (or nearly so.)

Upvotes: 1

Shashank
Shashank

Reputation: 13869

You have some indentation issues with your code but that could just be because of copy paste. Anyways adding a line sum = 0 at the beginning of your TrapArea function should solve your current error. But as @Blender pointed out in the comments, you have another issue, which is the lack of a multiplication operator (*) after your floating point division expression (1.0/2.0).

Remember that in Python expressions are not always evaluated as you would expect mathematically. Thus (a op b)(c) will not automatically multiply the result of a op b by c like you would expect with a mathematical expression. Instead this is the function call notation in Python.

Also remember that you must initialize all variables before using their values for assignment. Python has no default value for unnamed variables so when you reference the value of sum with sum += expr which is equivalent to sum = sum + expr you are trying to reference a name (sum) that is not binded to any object at all.

The following revision to your function should do the trick. Notice how I place multiplication operators (*) between every expression that you intend to multiply.

def TrapArea(n):
    sum = 0
    for i in range(1, n):
        i = float(i)
        deltax = (pi-0)/n
        sum += (1.0/2.0)*(((pi-0)/n)*(sin((i-1)/n*(pi-0))) + sin((i/n)*(pi-0)))*deltax
    return sum

EDIT: I also dealt with the float division issue by converting i to float(i) within every iteration of the loop. In Python 2.x, if you divide one integer type object with another integer type object, the expression evaluates to an integer regardless of the actual value.

Upvotes: 1

Veedrac
Veedrac

Reputation: 60117

Your original problem and problem with Shashank Gupta's answer was /n does integer division. You need to convert n to float first:

from math import pi, sin

def TrapArea(n):
    sum = 0
    for i in range(1, n):
        deltax = (pi-0)/n
        sum += (1.0/2.0)*(((pi-0)/float(n))*(sin((i-1)/float(n)*(pi-0))) + sin((i/float(n))*(pi-0)))*deltax
    return sum

for i in range(1, 11):
    print TrapArea(i)

Output:

0
0.785398163397
1.38175124526
1.47457409274
1.45836902046
1.42009115659
1.38070223089
1.34524797198
1.31450259385
1.28808354

Note that you can heavily simplify the sum += ... part.

First change all (pi-0) to pi:

sum += (1.0/2.0)*((pi/float(n))*(sin((i-1)/float(n)*pi)) + sin((i/float(n))*pi))*deltax

Then do pi/n wherever possible, which avoids needing to call float as pi is already a float:

sum += (1.0/2.0)*(pi/n * (sin((i-1) * pi/n)) + sin(i * pi/n))*deltax

Then change the (1.0/2.0) to 0.5 and remove some brackets:

sum += 0.5 * (pi/n * sin((i-1) * pi/n) + sin(i * pi/n)) * deltax

Much nicer, eh?

Upvotes: 3

Related Questions