Trapezoidal Rule
Passing functions around
Goal:
Making use of functions as objects – functions that act on functions.
Trapezoidal rule
The “trapezoidal rule”:
https://en.wikipedia.org/wiki/Trapezoidal_rule
Is one of the easiest “quadrature” methods.
Otherwise known as computing a definite integral, or, simply:
Computing the area under a curve.
The task
Your task is to write a trapz() function that will compute the area under an arbitrary function, using the trapezoidal rule.
The function will take another function as an argument, as well as the start and end points to compute, and return the area under the curve.
Example:
def line(x):
'''a very simple straight horizontal line at y = 5'''
return 5
area = trapz(line, 0, 10)
area
50
About the simplest “curve” you can have is a horizontal straight line, in this case, at y = 5. The area under that line from 0 to 10 is a rectangle that is 10 wide and 5 high, so with an area of 50.
Of course in this case, it’s easiest to simply multiply the height times the width, but we want a function that will work for Any curve.
HINT: this simple example could be a good test case!
The Solution:
Your function definition should look like:
def trapz(fun, a, b):
"""
Compute the area under the curve defined by
y = fun(x), for x between a and b
:param fun: the function to evaluate
:type fun: a function that takes a single parameter
:param a: the start point for the integration
:type a: a numeric value
:param b: the end point for the integration
:type b: a numeric value
"""
pass
In the function, you want to compute the following equation:
So you will need to:
create a list of x values from a to b (maybe 100 or so values to start)
compute the function for each of those values and double them
add them all up
multiply by the half of the difference between a and b divided by the number of steps.
Note that the first and last values are not doubled, so it may be more efficient to rearrange it like this:
NOTE: for those of you confused by that weird big greek letter (no it’s a name for a fraternity…) – see: A bit about math symbology
Can you use comprehensions for this?
NOTE: range() only works for integers – how can you deal with that?
Once you have that, it should work for any function that can be evaluated between a and b.
Try it for some built-in math functions, like math.sin
tests
Do this using test-drive development.
A few examples of analytical solutions you can use for tests:
A simple horizontal line – see above.
A sloped straight line:
The quadratic:
The sine function:
Computational Accuracy
In the case of the linear functions, the result should theoretically be exact. But with the vagaries of floating point math may not be.
And for non-linear functions, the result will certainly not be exact.
So you want to check if the answer is close to what you expect.
In py3.5+ – there is an isclose() function (PEP485)
Stage Two:
Some functions need extra parameters to do their thing. But the above will only handle a single parameter. For example, a quadratic function:
Requires values for A, B, and C in order to compute y from an given x.
You could write a specialized version of this function for each A, B, and C:
def quad1(x):
return 3 * x**2 + 2*x + 4
But then you need to write a new function for any value of these parameters you might need.
Instead, you can pass in A, B and C each time:
def quadratic(x, A=0, B=0, C=0):
return A * x**2 + B * x + C
Nice and general purpose.
But how would we compute the area under this function?
The function we wrote above only passes x in to the function it is integrating.
Passing arguments through:
Update your trapz() function so that you can give it a function that takes arbitrary extra arguments, either positional or keyword, after the x.
So you can do:
trapz(quadratic, 2, 20, A=1, B=3, C=2)
or
trapz(quadratic, 2, 20, 1, 3, C=2)
or
coef = {'A':1, 'B':3, 'C': 2}
trapz(quadratic, 2, 20, **coef)
NOTE: Make sure this will work with ANY function, with ANY additional positional or keyword arguments – not just this particular function.
This is pretty conceptually challenging – but it’s very little code!
If you are totally lost – look at the lecture notes from previous classes – how can you both accept and pass arbitrary arguments to/from a function?
You want your trapz function to take ANY function that can take ANY arbitrary extra arguments – not just the quadratic function, and not just A,B, and C. So good to test with another example.
The generalized sine function is:
where \(A\) is the amplitude, and \(\omega\) is the frequency of the function. In this case, the area under the curve from a to b is:
The test code has a test for this one, too.
Currying
Another way to solve the above problem is to use the original trapz, and create a custom version of the quadratic() function instead.
Write a function that takes A, B, and C as arguments, and returns a function that evaluates the quadratic for those particular coefficients.
Try passing the results of this into your trapz() and see if you get the same answer.
partial
Do the above with functools.partial as well.
Extra credit
This isn’t really the point of the exercise, but see if you can make it dynamically accurate.
How accurate it is depends on how small the chunks are that you break the function up into.
See if you can think of a way to dynamically determine how small a step you should use.
This is one for the math and computational programming geeks!
A bit about math symbology
Those of you without a lot of math background may be confused by the symbols. So here’s a quick intro to the “Summation Symbol” (Greek Capital sigma)
Is shorthand for “add up a bunch of values, with varying i from a to b”. where each x is a different value each time. Translating this into code you get:
x = a_list_of_numbers
total = 0
for i in range(a, b+1):
total += x[i]
or, in more compact python:
x = an iterable_of_numbers
total = sum(x[a:b+1])
So the full expression used above:
Can be written as:
sum(fun(x) for x in list_of_x[1:-1])