Week 2, Problem 1: Numeric
Integration
[35
points; individual or pair]
Note: There is a lot of explanation and text in this
lab. All the functions that you have to write in your hw2pr1.py
file are labeled with PN.
(where N is the problem number, e.g., 1, 2, 3, etc) and set off with a
horizontal line like this:
P0. This indicates a
function to write in your hw2pr1.py file.
The goal of this lab is to build a program in Python
that can do numerical integration of functions. The lab is designed to exercise
the following skills:
- using Python
to write functions
- composing
functions, i.e., extending or adapting old one 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.
Please create a new file in the IDLE editor called hw2pr1.py
and 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:
# April 14,2009
# John Smith
# hw2pr1.py -- done in lab
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.
- 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.
- 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].
- 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.
- 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:
- Divide the
interval in question into
Nequal parts and create a list with the correspondingxvalues. - Find the
yvalues of the function for each value ofxcalculated in step 1 - Calculate the
area of each rectangle under the curve (bounded by the two
xvalues on either side and theyvalue of the leftmostxvalue) - 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
Our goal here is to write the function steps(low,hi,N),
which takes in two numbers low and hi and a nonnegative integer N and
returns a regularly-spaced list of N floats starting at low and going up to, but not including, hi itself. This list will give us our x values (an number we want) within a given
range (whatever range we want).
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:
- generate a
list of
Nevenly spaced floating-point numbers starting at0.0and going up to but not including1.0 - adjust this
list so that it spans the correct
low,hirange.
P1. To solve step 1 above, in your hw2pr1.py file, write fracSteps(N),
which
should output a list of N evenly spaced floating-point values starting at 0.0 and going up to, but not including 1.0. N will be a positive
integer.
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:
[
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
comprehesions, check out the ListComprehension page.
P2. Now write steps(low,hi,N), which should take in two numbers low and hi and a nonnegative integer N. Then, steps should return a regularly-spaced list of N floats,
starting at low and going up to, but
not including, hi itself.
Hint:
You will want to incorporate low into TWO places in
your code. What is the importance of hi-low ?
We
would suggest getting rid of range at this point. Try writing
your return statement to include a list comprehension that looks like the
following (you'll need to fill in the missing portion (marked by ___), of
course!)
[_______________for x in fracSteps(N) ]
In
essence, steps is basically a
floating-point version of the built-in function, range. Here
are some examples of how your steps function should work:
>>>
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
Now that we have our x values, we need calculate the y
value of the function at each of these x positions. We'll again use list
comprehensions to make this process simple and fast! Although our goal is to
handle arbitrary functions, we'll start with a concrete function and build up.
P3. Write a function dbl_steps(low,hi,N)
that will double each of the values in the list that your steps
function creates. Again, try to use a list comprehension that takes advantage
of your steps
function. That is, consider how you could use [ ... for x in steps(low,hi,N)]
in writing dbl_steps.
This may turn out to be less difficult than it seems at first glance.
Here is an example:
>>>
dbl_steps(0,1,4)
[0.0, 0.5, 1.0, 1.5]
P4. Write fsteps(f,low,hi,N),
which should do the same thing as dbl_steps, but with an arbitrary
function f
that is passed in as the first input. You might be best off copying dbl_steps
and changing it. Here are some examples:
>>>
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 ...]
You might note that if you wrote fsteps
first, then you wouldn't need to write dbl_steps -- you could just use fsteps
with dbl
as the initial input.
This kind of generality is sometimes a good thing... . However, some people
find that they like to write a couple of specific examples before generalizing.
Now that you can generate x
and y
values for your functions, it would be nice to see them plotted graphically.
Please follow the link below that describes how to plot your functions in
python. There's nothing to submit in this part, but throughout the rest of the
lab it will be nice to have the ability to visualize your functions and your
integration approximations. When you understand how to plot functions, please
return to this page and continue the lab.
Click here
to go the EECS 110 Plot Page
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.
P5. Write the function
finteg (described in more detail
below). Its behavior just generalizes the step-by-step integration example
depicted and described above. You will want to use built-in functions and
functions that you've already written in previous parts of this lab! Feel
free to simply copy this function signature and docstring as-is and then fill
in the return statement:
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?
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.
>>>
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 hw2pr1.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 hw2pr1.py
file. Strings are easier because you don't need to preface multiple lines with
the comment character, #. The first problem is done in a string as an example:
Q0. Explain why the
true value of the last example should be 9.
The answer in the hw2pr1.py
file would look like:
""" 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
"""
Q1. The following
function, c,
traces the upper arc of a semicircle from -2 to 2 :
def c(x):
return (4 - x**2)**0.5
Place this function into your hw2pr1.py
file and confirm the values of
>>> finteg(c,0,2,2)
3.7320508075688772
>>> finteg(c,0,2,20)
3.2284648797549815
Next, find the values of finteg(c,0,2,200)
and finteg(c,0,2,2000) .
As N
goes to infinity, what does the value of this integral become?
Q2. Sometimes,
integrals without closed-form solutions are so important that we simply define
a function to be their solution! The natural log, sometimes written ln(x)
is one such example. It is defined as
This function is a built-in part of the math library under the name math.log(x).
Write a Python function ln(x,N) that uses your finteg
to compute the natural log of x using N steps. How accurate are ln(4.2,100)
and ln(4.2,10000)
in comparison with the real value of math.log(4.2)?
Note: You should include ln(x, N)
in your hw2pr1.py
file, but what we are really interested in is the answer to the above question.
(You just need to write this function to answer the question.)
Q3. Using finteg
as a guide, write a new function, finteg_tr, that uses a trapezoid
for each small interval, instead of rectangles.
Submission
You
should submit your hw2pr1.py file at Canvas .
Next