Hw2 Pr4: Numeric Integration

Lab Problem

(50 points, individual or pair)

This problem will ask you to write functions and a few text answers in a single file named hw2pr4.py.

Introduction

The goal of this lab is to build a program in Python that can integrate functions numerically. The lab is designed to exercise the following skills:

• using Python to write functions
• composing functions, i.e., extending or adapting old ones to meet new requirements
• applying functions to lists of inputs, i.e., list comprehensions or "map"
• breaking down a large problem (numerical integration) into several smaller problems

Before you begin, you need to set up your new file in the IDLE editor called hw2pr4.py. At the top of that file insert a header of comments that includes your name, the date, and the name of the file. For example:

# September 20-something, 2010
# _your name here_
#
# hw2pr4.py 

Overview

Integration is often described as finding the area between a mathematical function and the horizontal axis. This area is a single-value summary of the magnitude of the function across its independent variable. First, we present below a simple example of estimating an integral. Then, you will write an "integrator" function that works in general.

Recall the dbl function that we discussed last week:

def dbl(x):
    """ input: a number x (int or float)
        output: twice the input
    """
    return 2*x

Suppose you wanted to estimate the integral of the dbl function, between 0.0 and 10.0. You could do the following.

1. Divide the interval into 4 parts. Then the list [0, 2.5, 5, 7.5] represents the x value of the left-hand endpoint of each part of the interval.
2. Find the y-values from dbl(x) for each possible x in the list of left-hand endpoints above. These will be Y = [0, 5, 10, 15].
3. Add up the areas of the rectangles whose upper-left-hand corner is at these Y values. Each rectangle extends down to the x-axis.
4. Since there are four equal-width rectangles spanning a total of 10 units of width, each rectangle's individual width is 2.5. We find the rectangles' areas in the usual way. Their heights are contained in the list Y and their widths are 2.5, so the total area is
       0*2.5 + 5*2.5 + 10*2.5 + 15*2.5
   or
       (0 + 5 + 10 + 15)*2.5
   which is (30)*2.5, or 75.

Although this is quite a rough approximation for the integral, as the width of the rectangles gets smaller and smaller, it approximates the true integral more and more closely. Furthermore, the simple example above suggests a general way we can divide the problem of general numeric integration into pieces that we know how to write programs to solve. That is:

1. Divide the interval in question into N equal parts and create a list with the corresponding x values.
2. Find the y values of the function for each value of x calculated in step 1
3. Calculate the area of each rectangle under the curve (bounded by the two x values on either side and the y value of the leftmost x value)
4. Sum the areas and return the value

The rest of this lab involves writing functions to compute each of the four steps above (and optionally, to make the integral approximation slightly better) and then answering some questions about the results your integration approximation produces.

Step 1: Calculating the x values to be evaluated

At first you might think that a good way to approach this problem is to use the range function directly (recall that range(low, hi, step) returns a list of integers starting at low, up to but not including hi, step apart from each other.) Unfortunately, range returns a list of integers and we need a list of floats, so a single call to range will not suffice. Instead, we will again break up this problem into two pieces:

1. generate a list of N evenly spaced floating-point numbers starting at 0.0 and going up to but not including 1.0
2. adjust this list so that it spans the correct low, hi range.

Here are some examples:

>>> fracSteps(4) 
[0.0, 0.25, 0.5, 0.75]

>>> fracSteps(8)
[0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875]

>>> fracSteps(256)
(lots of floating-point numbers!)
        

NOTE: You should NOT use recursion to implement this function. A good strategy here is to use a list comprehension and the range command we looked at in class. For example:

[ float(x)/2 for x in range(5) ]

will produce the list

[0.0, 0.5, 1.0, 1.5, 2.0]

Recall that what's happening here is that range(5) is producing the list [0, 1, 2, 3, 4] and each of those values gets passed in as x to the function float(x) / 2, resulting in the list [0.0, 0.5, 1.0, 1.5, 2.0]. If we change either the list that gets produced, or the function that operates on the elements of the list, we will get a different result.

WARNING: Compare the above result to:

[ x / 2 for x in range(5) ]

which produces

[0, 0, 1, 1, 2]

Make sure you understand why.

Your list comprehension for this problem will look something like:

[                for x in range(N)]

If you want a step-by-step introduction to list comprehensions, check out the ListComprehension page.

[                  for x in fracSteps(N) ]
        

>>> steps(0,10,4)
[0.0, 2.5, 5.0, 7.5]

>>> steps(41,43,8)
[41.0, 41.25, 41.5, 41.75, 42.0, 42.25, 42.5, 42.75]
        

Step 2: Calculating the y values

>>> dbl_steps(0,1,4)
[0.0, 0.5, 1.0, 1.5]
        

>>> fsteps(dbl,0,1,4)
[0.0, 0.5, 1.0, 1.5]

>>> import math
>>> fsteps(math.sin,0,math.pi,6)
[0.0, 0.5, 0.87, 1.0, 0.87, 0.5]
(note: the above values are rounded versions of
          what actually gets displayed!)

>>> fsteps(math.exp,0,3,3)
[1.0, 2.7182818284590451, 7.3890560989306504]

# run dir(math) for the math functions (after "import math")
# run help(math.sin) for a bit more description of any f'n

>>> fsteps(dbl,-10,10,500)
[-20.0, ... 499 more list elements ...]
        

Steps 3+4: Calculate the area of the rectangles under the curve and putting it all together

Now that you can calculate the x values and the corresponding y values for functions of 1 input, you are ready to do the last two steps--use these values to calculate the area of the rectangles they form, and then sum these areas to compute the integral approximation.

def finteg(f,low,hi,N):
    """ finteg returns an estimate of the definite integral
        of the function f (the first input)
        with lower limit low (the second input)
        and upper limit hi (the third input)
        where N steps are taken (the fourth input)

	finteg simply returns the sum of the areas of rectangles
        under f, drawn at the left endpoint of the N steps taken
        from low to hi
    """

Hints: This function will be VERY short, but it can be a little tricky to write. Go very carefully through the example at the start of this lab, matching the concrete numbers in that example with the inputs to finteg . What corresponds to what?

• In particular, remember that the output of fsteps is a big list of y-values!
• You will want to multiply those y-values, which are the heights of the rectangles, by the width of those rectangles... what is that width?

Also, don't rewrite your fsteps function from the previous part of this lab - simply use it. If you use fsteps , the finteg function does not need list comprehensions or map or recursion at all...but sum will be useful.

Here is how sum works:

>>> sum( [10, 4, 1] )
15

>>> sum( range(0,101) )
5050

Here are some examples to check. You will need to add the two small functions

def dbl(x):
    """ doubler!  input: x, a number """
    return 2*x

def sq(x):
    """ squarer!  input: x, a number """
    return x**2

>>> finteg(dbl,0,10,4)
75.0

Help - I get 60 instead of 75!! You may be dividing an int by an int with a resulting rounding error... Multiply low by 1.0 to create a float value, for example.

>>> import math
>>> finteg(math.sin, 0, math.pi, 1000)
1.9999983550656628          # pretty good - should be 2.0

>>> finteg(sq,0,3,1000000)  # a million steps will give it pause
8.9999865000044963          # close! - should be 9.0

Questions to answer in your hw2pr4.py file

The last required part of the lab (and HW problem 1) involves using your finteg function to answer a number of questions. You should put your answers either within comments or within triple-quoted strings in your hw2pr4.py file. Strings are easier because you don't need to preface multiple lines with the comment character, #.

You don't have to answer this zeroth question, but we have placed the answer here in a string, as an example:

""" Q0. The "true value" of finteg(sq,0,3,1000000) -- that is,
     the integral of sq(x) from a low limit of 0 to a high
     limit of 3 -- equals

     x**3 / 3.0, evaluated at 3.0   MINUS
     x**3 / 3.0, evaluated at 0.0

     These evaluate to        27.0/3.0 - 0.0, which equals 9.0
"""

def c(x):
    return (4 - x**2)**0.5
    

>>> finteg(c,0,2,2)
3.7320508075688772

>>> finteg(c,0,2,20)
3.2284648797549815