This chapter is taken from the book A Primer on Scientific Programming with Python by H. P. Langtangen, 5th edition, Springer, 2016.
Make a class F that implements the function
$$
\begin{equation*} f(x; a, w)=e^{-ax}\sin(wx) \thinspace . \end{equation*}
$$
A value(x) method computes values of \( f \), while \( a \) and \( w \) are data
attributes.
Test the class in an interactive session:
>>> from F import F
>>> f = F(a=1.0, w=0.1)
>>> from math import pi
>>> print f.value(x=pi)
0.013353835137
>>> f.a = 2
>>> print f.value(pi)
0.00057707154012
Filename: F.
Add a data attribute transactions to the Account class from
the section Bank accounts. The new attribute counts the number
of transactions done in the deposit and withdraw methods.
Print the total number of transactions in the dump
method. Write a test function test_Account() for testing that the
implementation of the extended class Account is correct.
Filename: Account2.
In class AccountP from the section Bank accounts,
introduce a list self._transactions, where each element holds a dictionary
with the amount of a transaction and the point of time the transaction
took place. Remove the _balance attribute and use instead the
_transactions list to compute the balance in the method
get_balance. Print out a nicely formatted table of all
transactions, their amounts, and their time in a method
print_transactions.
Use the time or datetime module to get the date and local time.
Filename: Account3.
Observe that the computation of the balance is implemented
in a different way in the present version of class
AccountP compared
to the version in the section Bank accounts, but the usage
of the class, especially the get_balance method,
remains the same. This is one of the great advantages
of class programming: users are supposed to use the methods only, and
the implementation of data structures and computational techniques
inside methods can be changed without affecting existing
programs that just call the methods.
The purpose of this exercise is to create classes like
class Circle from the section A circle for
representing other geometric figures: a rectangle with width
\( W \), height \( H \), and lower left corner \( (x_0, y_0) \); and
a general triangle specified by its
three vertices \( (x_0,y_0) \), \( (x_1,y_1) \), and
\( (x_2,y_2) \).
Provide three methods: __init__ (to initialize the geometric data),
area, and perimeter. Write test functions test_Rectangle()
and test_Triangle() for
checking that the results produced by area and perimeter
coincide with exact values within a small tolerance.
Filename: geometric_shapes.
Consider a quadratic function \( f(x;a,b,c)=ax^2 + bx + c \).
Make a class Quadratic
for representing \( f \), where \( a \), \( b \), and \( c \) are data attributes, and
the methods are
__init__ for storing the attributes \( a \), \( b \), and \( c \),value for computing a value of \( f \) at a point \( x \),table for writing out a table of \( x \) and \( f \) values for \( n \) $x$ values in the interval \( [L,R] \),roots for computing the two roots.Quadratic and corresponding demonstrations and/or
tests should be organized as a module such that other programs can
do a from Quadratic import Quadratic to use the class.
Also equip the file with a test function for verifying the implementation
of value and roots.
Filename: Quadratic.
Make a class Line whose constructor takes two points
p1 and p2 (2-tuples or 2-lists) as input.
The line goes through these two
points.
A value(x) method computes a value on the line
at the point x. Also make
a function test_Line() for verifying the implementation.
Here is a demo in an interactive session:
>>> from Line import Line, test_Line
>>> line = Line((0,-1), (2,4))
>>> print line.value(0.5), line.value(0), line.value(1)
0.25 -1.0 1.5
>>> test_Line()
Filename: Line.
The constructor in class Line in Exercise 6: Make a class for straight lines takes
two points as arguments. Now we want to have more flexibility in the
way we specify a straight line: we can give two points, a point and a
slope, or a slope and the line's interception with the \( y \) axis.
Write this extended class and a test function for checking that the
increased flexibility does work.
Let the constructor take two arguments p1 and p2
as before, and test with
isinstance whether the
arguments are float versus tuple or list to determine
what kind of data the user supplies:
if isinstance(p1, (tuple,list)) and isinstance(p2, (float,int)):
# p1 is a point and p2 is slope
self.a = p2
self.b = p1[1] - p2*p1[0]
elif ...
Filename: Line2.
The purpose of this exercise is to make a class interface to an
already existing set of functions implementing Lagrange's
interpolation method from
the exercise named ``Implement Lagrange's interpolation formula'' in the document Arrays and plotting
[4].
We want to
construct a class LagrangeInterpolation with a typical usage like:
import numpy as np
# Compute some interpolation points along y=sin(x)
xp = np.linspace(0, np.pi, 5)
yp = np.sin(xp)
# Lagrange's interpolation polynomial
p_L = LagrangeInterpolation(xp, yp)
x = 1.2
print 'p_L(%g)=%g' % (x, p_L(x)),
print 'sin(%g)=%g' % (x, np.sin(x))
p_L.plot() # show graph of p_L
The plot method visualizes \( p_L(x) \) for \( x \) between the first
and last interpolation point (xp[0] and xp[-1]).
In addition to writing the class itself, you should write code to verify
the implementation.
The class does not need much code as it can call
the functions p_L from
ref{sec:class:ex27a} and
graph from ref{sec:class:ex27b}, available in
the Lagrange_poly2 module made in the latter exercise.
Filename: Lagrange_poly3.
Instead of manually computing the interpolation points,
as demonstrated in Exercise 8: Wrap functions in a class, we now want
the constructor in class LagrangeInterpolation to also
accept some Python function f(x) for computing
the interpolation points.
Typically, we would like to write
this code:
from numpy import exp, sin, pi
def myfunction(x):
return exp(-x/2.0)*sin(x)
p_L = LagrangeInterpolation(myfunction, x=[0, pi], n=11)
With such a code, \( n=11 \) uniformly distributed \( x \) points between
\( 0 \) and \( \pi \) are computed, and the corresponding \( y \) values are
obtained by calling myfunction.
The Lagrange interpolation polynomial is then constructed from
these points. Note that the previous types of calls,
LangrangeInterpolation(xp, yp), must still be valid.
The constructor in class LagrangeInterpolation must now
accept two different sets of arguments: xp, yp vs. f, x, n.
You can use the isinstance(a, t) function to test if object a is
of type t. Declare the constructor with three
arguments arg1, arg2, and arg3=None. Test
if arg1 and arg2 are arrays
(isinstance(arg1, numpy.ndarray)), and
in that case, set xp=arg1 and yp=arg2. On the other hand,
if arg1 is a function (callable(arg1) is True),
arg2 is a list or tuple (isinstance(arg2, (list,tuple))),
and arg3 is an integer, set
f=arg1, x=arg2, and n=arg3.
Filename: Lagrange_poly4.
Write a class Hello that behaves as illustrated in the following
session:
>>> a = Hello()
>>> print a('students')
Hello, students!
>>> print a
Hello, World!
Filename: Hello.
Modify the class from Exercise 1: Make a function class such that the following interactive session can be run:
>>> from F import F
>>> f = F(a=1.0, w=0.1)
>>> from math import pi
>>> print f(x=pi)
0.013353835137
>>> f.a = 2
>>> print f(pi)
0.00057707154012
>>> print f
exp(-a*x)*sin(w*x)
Filename: F2.
The task in this exercise is to calculate a sum \( S(x)=\sum_{k=M}^N f_k(x) \),
where \( f_k(x) \) is some user-given formula for the terms in the sum.
The following snippet demonstrates the typical use and functionality
of a class Sum for computing \( S(x)= \sum_{k=0}^N (-x)^k \):
def term(k, x):
return (-x)**k
S = Sum(term, M=0, N=3)
x = 0.5
print S(x)
print S.term(k=4, x=x) # (-0.5)**4
a)
Implement class Sum such that the code snippet above works.
b)
Implement a test function test_Sum() for verifying the results of
the various methods in class Sum for a specific choice of \( f_k(x) \).
c)
Apply class Sum to compute the Taylor polynomial approximation to
\( \sin x \) for \( x=\pi \) and some chosen \( x \) and \( N \).
Filename: Sum.
Isolate class Derivative from the section Example: Automagic differentiation
in a module file. Also isolate class Y from
the section Representing a function as a class in a module file.
Make a program that imports class Derivative and class
Y and applies the former to differentiate the
function \( y(t)=v_0t - \frac{1}{2}gt^2 \) represented by class Y.
Compare the computed derivative with the exact value for
\( t=0, \frac{1}{2}v_0/g, v_0/g \).
Filenames: dYdt.py, Derivative.py, Y.py.
An anthropologist was asking a primitive tribesman about arithmetic. When the anthropologist asked, What does two and two make? the tribesman replied, Five. Asked to explain, the tribesman said, If I have a rope with two knots, and another rope with two knots, and I join the ropes together, then I have five knots.
a)
Make a class Rope for representing a rope with a given number of knots.
Implement the addition operator in this class such that we can join two
ropes together in the way the tribesman described:
>>> from Rope import Rope
>>> rope1 = Rope(2)
>>> rope2 = Rope(2)
>>> rope3 = rope1 + rope2
>>> print rope3
5
As seen, the class also features a __str__ method for returning
the number of knots on the rope.
b) Equip the module file with a test function for verifying the implementation of the addition operator.
Filename: Rope.py.
+= and -= operators
As alternatives to the deposit and withdraw methods
in class Account class from
the section Bank accounts, we could use the operation += for
deposit and -= for withdraw.
Implement the += and -= operators, a
__str__ method, and preferably a
__repr__ method in class Account. Write a test_Account()
function to verify the implementation of all functionality
in class Account.
The special methods __iadd__ and __isub__
implement the += and -= operators, respectively.
For instance, a -= p implies a call to a.__isub__(p).
One important feature of __iadd__ and __isub__
is that they must return self to work properly,
see the documentation of these
methods in Chapter 3 of the Python Language Reference.
Filename: Account4.
A widely used formula for numerical differentiation of a function \( f(x) \) takes the form $$ \begin{align} f'(x) & \approx {f(x+h) - f(x-h)\over 2h} \tp \tag{8} \end{align} $$ This formula usually gives more accurate derivatives than (1) because it applies a centered, rather than a one-sided, difference.
The goal of this exercise is to use the formula (8)
to automatically differentiate a mathematical function \( f(x) \) implemented
as a Python function f(x). More precisely, the following
code should work:
def f(x):
return 0.25*x**4
df = Central(f) # make function-like object df
# df(x) computes the derivative of f(x) approximately
x = 2
print 'df(%g)=%g' % (x, df(x))
print 'exact:', x**3
a)
Implement class Central and test that the code above works.
Include an optional argument h to the constructor in class
Central so that \( h \) in the
approximation (8) can be specified.
b)
Write a test function test_Central() to verify the implementation.
Utilize the fact that the formula (8) is
exact for quadratic polynomials (provided \( h \) is not too small, then
rounding errors in (8) require use of a
(much) larger tolerance than the expected machine precision).
c)
Write a function table(f, x, h=1E-5) that prints a
table of errors in the numerical derivative (8)
applied to a function f at some points x.
The argument f is a sympy expression for
a function. This f object can be transformed to a Python function and fed
to the constructor of class Central, and f can be used to compute
the exact derivative symbolically.
The argument x is a list or array of points \( x \), and h is the \( h \)
in (8).
The following session demonstrates how sympy can differentiate a
mathematical expression and turn the result into a Python
function:
>>> import sympy
>>> x = sympy.Symbol('x')
>>> f_expr = 'x*sin(2*x)'
>>> df_expr = sympy.diff(f_expr)
>>> df_expr
2*x*cos(2*x) + sin(2*x)
>>> df = sympy.lambdify([x], df_expr) # make Python function
>>> df(0)
0.0
d) Organize the file with the class and functions such that it can be used a module.
Filename: Central.
Consider this program file for computing a backward difference approximation
to the derivative of a function f(x):
from math import *
class Backward(object):
def __init__(self, f, h=e-9):
self.f, self.h = f, h
def __call__(self, x):
h, f = self.h, self.f
return (f(x) - f(x-h))/h # finite difference
dsin = Backward(sin)
e = dsin(0) - cos(0); print 'error:', e
dexp = Backward(exp, h=e-7)
e = dexp(0) - exp(0); print 'error:', e
The output becomes
error: -1.00023355634 error: 371.570909212
Is the approximation that bad, or are there bugs in the program?
Filename: find_errors_class.
Make the two data attributes h and f of class Derivative from
the section Example: Automagic differentiation protected as explained in the section Bank accounts. That is, prefix h and f with an underscore
to tell users that these attributes should not be accessed
directly. Add two methods get_precision() and set_precision(h) for
reading and changing h. Make a separate test function for checking
that the new class works as intended.
Filename: Derivative_protected.
a) Use a class to implement the discontinuous Heaviside function and smoothed continuous version, as defined in the exercise "Implement a smoothed Heaviside function" in the document Functions and branching, such that the following code works:
H = Heaviside() # original discontinous Heaviside function
print H(0.1)
H = Heaviside(eps=0.8) # smoothed continuous Heaviside function
print H(0.1)
b)
Extend class Heaviside such that
array arguments are allowed:
H = Heaviside() # original discontinous Heaviside function
x = numpy.linspace(-1, 1, 11)
print H(x)
H = Heaviside(eps=0.8) # smoothed Heaviside function
print H(x)
Use ideas from the section ref{sec:vec:Heaviside}.
c)
Extend class Heaviside such that it supports plotting:
H = Heaviside()
x, y = H.plot(xmin=-4, xmax=4) # x in [-4, 4]
from matplotlib.pyplot import plot
plot(x, y)
H = Heaviside(eps=1)
x, y = H.plot(xmin=-4, xmax=4)
plot(x, y)
Techniques from the section ref{sec:plot:pwisefunc} must in the first case
be used to return arrays x and y such that the discontinuity is
exactly reproduced. In the continuous (smoothed) case, one needs to
compute a sufficiently fine resolution (x) based on the eps
parameter, e.g., 201/$\epsilon$ points in the interval \( [-\epsilon,
\epsilon] \), with a coarser set of coordinates outside this interval
where the smoothed Heaviside function is almost constant, 0 or 1.
d)
Write a test function test_Heaviside() for verifying the result
of the various methods in class Heaviside.
Filename: Heaviside_class.
The purpose of this exercise is the make a class implementation of the indicator function from ref[ref{sec:basic:exH3}[ in [5]][the exercise named "Implement an indicator function" in the document Functions and branching [5]]. Let the implementation be based on expressing the indicator function in terms of Heaviside functions. Allow for an \( \epsilon \) parameter in the calls to the Heaviside function, such that we can easily choose between a discontinuous and a smoothed, continuous version of the indicator function:
I = Indicator(a, b) # indicator function on [a,b]
print I(b+0.1), I((a+b)/2.0)
I = Indicator(0, 2, eps=1) # smoothed indicator function on [0,2]
print I(0), I(1), I(1.9)
Note that if you build on the version of class Heaviside
in Exercise 19: Make a class for the Heaviside functionb, any Indicator instance
will accept array arguments too.
Filename: Indicator.
The purpose of this exercise is to make a class implementation of a piecewise constant function, as defined in ref[ref{sec:basic:exH4}[ in [5]][the exercise named "Implement a piecewise constant function" in the document Functions and branching [5]].
a) Implement the minimum functionality such that the following code works:
f = PiecewiseConstant([(0.4, 1), (0.2, 1.5), (0.1, 3)], xmax=4)
print f(1.5), f(1.75), f(4)
x = np.linspace(0, 4, 21)
print f(x)
b)
Add a plot method to class PiecewiseConstant
such that we
can easily plot the graph of the function:
x, y = f.plot()
from matplotlib.pyplot import plot
plot(x, y)
Filename: PiecewiseConstant.
The observant reader may have noticed that our Integral class
from the section Example: Automagic integration
is very inefficient if we want to tabulate or plot a
function \( F(x)=\int_a^xf(x) \) for several consecutive values of \( x \):
\( x_0 < x_1 < \cdots < x_m \). Requesting \( F(x_k) \) will recompute
the integral computed for \( F(x_{k-1}) \), and this is of course
waste of computer work.
Modify
the __call__ method such that if x is an array,
assumed to contain coordinates of increasing value:
\( x_0 < x_1 < \cdots < x_m \), the method returns an array with
\( F(x_0), F(x_1),\ldots,F(x_m) \) with the minimum computational work.
Also write a test function to verify that the implementation is correct.
The \( n \) (n) parameter in the constructor of the Integral class can be
taken as the total number of trapezoids (intervals) that are to be used to
compute the final \( F(x_m) \) value. The integral over an interval
\( [x_k, x_{k+1}] \) can then be computed by the trapezoidal function
(or an Integral object) using an appropriate fraction of the \( n \)
total trapezoids. This fraction can be \( (x_{k+1}-x_k)/(x_m-a) \)
(i.e., \( n_k = n(x_{k+1}-x_k)/(x_m-a) \)) or
one may simply use a constant \( n_k=n/m \) number of trapezoids for
all the integrals over \( [x_k, x_{k+1}] \), \( k=0,\ldots,m-1 \).
Filename: Integral_eff.
The Taylor polynomial of degree \( N \) for the exponential function \( e^x \)
is given by
$$
\begin{equation*} p(x) = \sum_{k=0}^N {x^k\over k!} \thinspace . \end{equation*}
$$
Make a program that (i) imports class
Polynomial from the section Example: Class for polynomials, (ii) reads \( x \)
and a series of \( N \) values
from the command line, (iii) creates a
Polynomial object for each
\( N \) value for computing with the given Taylor polynomial, and (iv)
prints the values of \( p(x) \) for all the given \( N \) values as well as the
exact value \( e^x \).
Try the program out with
\( x=0.5, 3, 10 \) and \( N=2,5,10, 15, 25 \).
Filename: Polynomial_exp.
Go through this alternative implementation of class
Polynomial from the section Example: Class for polynomials
and explain each line in detail:
class Polynomial(object):
def __init__(self, coefficients):
self.coeff = coefficients
def __call__(self, x):
return sum([c*x**i for i, c in enumerate(self.coeff)])
def __add__(self, other):
maxlength = max(len(self), len(other))
# Extend both lists with zeros to this maxlength
self.coeff += [0]*(maxlength - len(self.coeff))
other.coeff += [0]*(maxlength - len(other.coeff))
result_coeff = self.coeff
for i in range(maxlength):
result_coeff[i] += other.coeff[i]
return Polynomial(result_coeff)
The enumerate function, used in the __call__ method,
enables us to iterate over a list somelist with
both list indices and list elements: for index, element in enumerate(somelist).
Write the code above in a file, and demonstrate that adding two polynomials
does not work. Find the bug and correct it.
Filename: Polynomial_error.
Implement the special method __sub__ in class
Polynomial from the section Example: Class for polynomials.
Add a test for this functionality in function test_Polynomial.
Study the __add__ method in class Polynomial and treat the
two cases, where the lengths of the lists in the polynomials differs,
separately.
Filename: Polynomial_sub.
Verify the functionality of the __str__ method in class Polynomial from
the section Example: Class for polynomials by writing a new test function
test_Polynomial_str().
Filename: Polynomial_test_str.
Introducing an array instead of a list in class Polynomial does not
enhance the efficiency of the implementation unless the mathematical
computations are also vectorized. That is, all explicit Python loops
must be substituted by vectorized expressions.
a)
Go through class Polynomial.py and make sure the coeff attribute
is always a numpy array with float elements.
b)
Update the test function test_Polynomial to make use of the fact
that the coeff attribute is always a
numpy array with float elements. Run test_Polynomial
to check that the new implementation is correct.
c)
Vectorize the __add__ method by adding the
common parts of the coefficients arrays and then appending
the rest of the longest array to the result.
Appending an array
a to an array b can be done by concatenate(a, b).
d)
Vectorize the __call__ method by
observing that evaluation of a polynomial, \( \sum_{i=0}^{n-1}c_ix^i \),
can be computed as the inner product of two arrays:
\( (c_0,\ldots,c_{n-1}) \) and \( (x^0,x^1,\ldots,x^{n-1}) \).
The latter array can be computed by x**p, where
p is an array with powers \( 0,1,\ldots,n-1 \), and x
is a scalar.
e)
The differentiate method can be vectorized by
the statements
n = len(self.coeff)
self.coeff[:-1] = linspace(1, n-1, n-1)*self.coeff[1:]
self.coeff = self.coeff[:-1]
Show by hand calculations in a case where n is 3 that
the vectorized statements produce the same result as the
original differentiate method.
Filename: Polynomial_vec.
The __mul__ method is more challenging to vectorize so you may leave
this unaltered. Check that the
vectorized versions of __add__,
__call__, and differentiate work as intended by calling
the test_Polynomial function.
Use a dictionary (instead of a list)
for the coeff attribute in class Polynomial from
the section Example: Class for polynomials such that self.coeff[k] holds the
coefficient of the \( x^k \) term. The advantage with a dictionary is
that only the nonzero coefficients in a polynomial need to be stored.
a)
Implement a constructor and the __call__ method for evaluating the
polynomial. The following demonstration code should work:
from Polynomial_dict import Polynomial
p1_dict = {4: 1, 2: -2, 0: 3} # polynomial x^4 - 2*x^2 + 3
p1 = Polynomial(p1_dict)
print p1(2) # prints 11 (16-8+3)
b)
Implement the __add__ method. The following demonstration code should work:
p1 = Polynomial({4: 1, 2: -2, 0: 3}) # x^4 - 2*x^2 + 3
p2 = Polynomial({0: 4, 1: 3} # 4 + 3*x
p3 = p1 + p2 # x^4 - 2*x^2 + 3*x + 7
print p3.coeff # prints {0: 7, 1: 3, 2: -2, 4: 1}
The structure of __add__ may be
class Polynomial(object):
...
def __add__(self, other):
"""Return self + other as a Polynomial object."""
result = self.coeff.copy()
for exponent in result:
if exponent in other.coeff:
# add other's term to result's term
else:
result[exponent] = other[exponent]
# return Polynomial object based on result dict
c)
Implement the __sub__ method. The following demonstration code should work:
p1 = Polynomial({4: 1, 2: -2, 0: 3}) # x^4 - 2*x^2 + 3
p2 = Polynomial({0: 4, 1: 3} # 4 + 3*x
p3 = p1 - p2 # x^4 - 2*x^2 - 3*x - 1
print p3.coeff # prints {0: -1, 1: -3, 2: -2, 4: 1}
d)
Implement the __mul__ method. The following demonstration code should work:
p1 = Polynomial({0: 1, 3: 1}) # 1 + x^3
p2 = Polynomial({1: -2, 2: 3}) # -2*x + 3*x^2
p3 = p1*p3
print p3.coeff # prints {1: -2, 2: 3, 4: -2, 5: 3}
Study the __mul__ method in class Polynomial based on a list
representation of the data in the polynomial and adapt to
a dictionary representation.
e)
Write a test function for each of the methods __call__, __add__,
and __mul__.
Filename: Polynomial_dict.
The Vec2D class from the section Example: Class for vectors in the plane supports
addition and subtraction, but only addition and subtraction of two
Vec2D objects. Sometimes we would like to add or subtract a
point that is represented by a list or a tuple:
u = Vec2D(-2, 4)
v = u + (1,1.5)
w = [-3, 2] - v
That is, a list or a tuple must be allowed in the right or left operand.
Implement such an extension of class Vec2D.
Ideas are found in the sections Mixing complex and real numbers and Special methods for "right" operands.
Filename: Vec2D_lists.
Extend the implementation of class Vec2D
from the section Example: Class for vectors in the plane to a class Vec3D for vectors
in three-dimensional space. Add a method cross for computing
the cross product of two 3D vectors.
Filename: Vec3D.
The internal code in class Vec2D from the section Example: Class for vectors in the plane
can be valid for vectors in any space dimension if we represent the
vector as a NumPy array in the class instead of separate variables x
and y for the vector components. Make a new class Vec where you
apply NumPy functionality in the methods. The constructor should be
able to treat all the following ways of initializing a vector:
a = array([1, -1, 4], float) # numpy array
v = Vec(a)
v = Vec([1, -1, 4]) # list
v = Vec((1, -1, 4)) # tuple
v = Vec(1, -1) # coordinates
In the constructor, use variable number of arguments as
described in the document Variable number of
function arguments in Python
[2]. All arguments are then available as
a tuple, and if there is only one element in the tuple, it should be
an array, list, or tuple you can send through asarray to get a NumPy
array. If there are many arguments, these are coordinates, and the
tuple of arguments can be transformed by array to a NumPy array.
Assume in all operations that the involved vectors have equal
dimension (typically that other has the same dimension as
self). Recall to return Vec objects from all arithmetic
operations, not NumPy arrays, because the next operation with the
vector will then not take place in Vec but in NumPy. If self.v is
the attribute holding the vector as a NumPy array, the addition
operator will typically be implemented as
class Vec(object):
...
def __add__(self, other):
return Vec(selv.v + other.v)
Filename: Vec.
Consider the function \( f(x)=x/(1+x) \) on \( [1,2] \). Find the variation of \( f \)
over \( [1,2] \). Use interval arithmetics from
the section Example: interval arithmetic to compute the variation of
\( f \) when \( x\in [1,2] \).
Filename: interval_arithmetics.
In this case, interval arithmetics overestimates the variation in \( f \). The reason is that \( x \) occurs more than once in the formula for \( f \) (the so-called dependency problem).
Use classes to reimplement the summarizing problem in
the section "Example: A file database"
in the document Dictionaries and strings
[4].
More precisely, introduce a class Student and a class
Course. Find appropriate attributes. The classes
should have a __str__ method for pretty-printing of the contents.
Filename: Student_Course.
Extreme points of a function \( f(x) \) are normally found by solving \( f'(x)=0 \). A much simpler method is to evaluate \( f(x) \) for a set of discrete points in the interval \( [a, b] \) and look for local minima and maxima among these points. We work with \( n+1 \) equally spaced points \( a=x_0 < x_1 < \cdots < x_{n}=b \), \( x_i=a+ih \), \( h=(b-a)/n \).
First we find all local extreme points in the interior of the domain. Local minima are recognized by $$ \begin{equation*} f(x_{i-1}) > f(x_i) < f(x_{i+1}),\quad i=1,\ldots,n-1 \thinspace . \end{equation*} $$ Similarly, at a local maximum point \( x_i \) we have $$ \begin{equation*} f(x_{i-1}) < f(x_i) > f(x_{i+1}),\quad i=1,\ldots,n-1 \thinspace . \end{equation*} $$ Let \( P_{\min} \) be the set of \( x \) values for local minima and \( F_{\min} \) the set of the corresponding \( f(x) \) values at these minima. Two sets \( P_{\max} \) and \( F_{\max} \) are defined correspondingly for the maxima.
The boundary points \( x=a \) and \( x=b \) are for algorithmic simplicity also defined as local extreme points: \( x=a \) is a local minimum if \( f(a) < f(x_1) \), and a local maximum otherwise. Similarly, \( x=b \) is a local minimum if \( f(b) < f(x_{n-1}) \), and a local maximum otherwise. The end points \( a \) and \( b \) and the corresponding function values must be added to the sets \( P_{\min},P_{\max},F_{\min},F_{\max} \).
The global maximum point is defined as the \( x \) value corresponding to the maximum value in \( F_{\max} \). The global minimum point is the \( x \) value corresponding to the minimum value in \( F_{\min} \).
a)
Make a class MinMax with the following functionality:
__init__ takes \( f(x) \), \( a \), \( b \), and \( n \) as arguments, and calls a method _find_extrema to compute the local and global extreme points._find_extrema implements the algorithm above for finding local and global extreme points, and stores the sets \( P_{\min},P_{\max},F_{\min},F_{\max} \) as list attributes in the (self) instance.get_global_minimum returns the global minimum point as a pair
\( (x, f(x)) \).get_global_maximum returns the global maximum point as a pair
\( (x, f(x)) \).get_all_minima returns a list or array of all \( (x,f(x)) \) minima.get_all_maxima returns a list or array of all \( (x,f(x)) \) maxima.__str__ returns a string where a nicely formatted table
of all the min/max points are listed, plus the global extreme points.MinMax:
def f(x):
return x**2*exp(-0.2*x)*sin(2*pi*x)
m = MinMax(f, 0, 4, 5001)
print m
The output becomes
All minima: 0.8056, 1.7736, 2.7632, 3.7584, 0 All maxima: 0.3616, 1.284, 2.2672, 3.2608, 4 Global minimum: 3.7584 Global maximum: 3.2608
Make sure that the program also works for functions without local extrema, e.g., linear functions \( f(x)=ax+b \).
b)
The algorithm sketched above finds local extreme
points \( x_i \), but all we know is that the true extreme point is in the
interval \( (x_{i-1},x_{i+1}) \). A more accurate algorithm may take this
interval as a starting point and run a Bisection method (see
the final part of the document User
input and error handling
[6]) to find the extreme point \( \bar x \) such
that \( f'(\bar x)=0 \). Add a method
_refine_extrema in class MinMax,
which goes through all the interior local minima
and maxima and solves \( f'(\bar x)=0 \). Compute \( f'(x) \) using the
Derivative class (the section Example: Automagic differentiation with \( h \ll
x_{i+1}-x_{i-1} \).
Filename: minmaxf.
The company PROD produces two different products, P$_1$ and P$_2$, based on three different raw materials, \( \mbox{M}_1 \), \( \mbox{M}_2 \) and \( \mbox{M}_3 \). The following table shows how much of each raw material \( \mbox{M}_i \) that is required to produce a single unit of each product P$_j$:
| P$_1$ | P$_2$ | |
| \( \mbox{M}_1 \) | 2 | 1 |
| \( \mbox{M}_2 \) | 5 | 3 |
| \( \mbox{M}_3 \) | 0 | 4 |
For instance, to produce one unit of \( \mbox{P}_2 \) one needs 1 unit of \( \mbox{M}_1 \), 3 units of \( \mbox{M}_2 \) and 4 units of \( \mbox{M}_3 \). Furthermore, PROD has available 100, 80 and 150 units of material \( \mbox{M}_1 \), \( \mbox{M}_2 \) and \( \mbox{M}_3 \) respectively (for the time period considered). The revenue per produced unit of product \( \mbox{P}_1 \) is 150 NOK, and for one unit of \( \mbox{P}_2 \) it is 175 NOK. On the other hand the raw materials \( \mbox{M}_1 \), \( \mbox{M}_2 \) and \( \mbox{M}_3 \) cost 10, 17 and 25 NOK per unit, respectively. The question is: how much should PROD produce of each product? We here assume that PROD wants to maximize its net revenue (which is revenue minus costs).
a) Let \( x \) and \( y \) be the number of units produced of product \( \mbox{P}_1 \) and \( \mbox{P}_2 \), respectively. Explain why the total revenue \( f(x,y) \) is given by $$ \begin{equation*} f(x,y)=150x-(10\cdot 2+17\cdot 5)x+ 175y-(10\cdot 1+17\cdot 3 + 25\cdot 4)y \end{equation*} $$ and simplify this expression. The function \( f(x,y) \) is linear in \( x \) and \( y \) (make sure you know what linearity means).
b) Explain why PROD's problem may be stated mathematically as follows: $$ \begin{equation} \tag{9} \begin{array}{lrrrrr} \mbox{\rm maximize} &\multicolumn{3}{c}{f(x,y)} \\ \mbox{\rm subject to} & \\ &2x &+ &y &\le &100 \\ &5x &+ &3y &\le &80 \\ & & &4y &\le &150 \\ & &\multicolumn{4}{l}{x\ge 0, y \ge 0.} \end{array} \end{equation} $$ This is an example of a linear optimization problem.
c) The production \( (x,y) \) may be considered as a point in the plane. Illustrate geometrically the set \( T \) of all such points that satisfy the constraints in model (9). Every point in this set is called a feasible point.
For every inequality determine first the straight line obtained by replacing the inequality by equality. Then, find the points satisfying the inequality (a half-plane), and finally, intersect these half-planes.
d) Make a program for drawing the straight lines defined by the inequalities. Each line can be written as \( ax + by = c \). Let the program read each line from the command line as a list of the \( a \), \( b \), and \( c \) values. In the present case the command-line arguments will be
'[2,1,100]' '[5,3,80]' '[0,4,150]' '[1,0,0]' '[0,1,0]'
Perform an eval on the elements of sys.argv[1:] to
get \( a \), \( b \), and \( c \) for each line as a list in the program.
e) Let \( \alpha \) be a positive number and consider the level set of the function \( f \), defined as the set $$ \begin{equation*} L_{\alpha}= \{(x,y) \in T: f(x,y)=\alpha\}. \end{equation*} $$ This set consists of all feasible points having the same net revenue \( \alpha \). Extend the program with two new command-line arguments holding \( p \) and \( q \) for a function \( f(x,y)=px + qy \). Use this information to compute the level set lines \( y=\alpha /q - px/q \), and plot the level set lines for some different values of \( \alpha \) (use the \( \alpha \) value in the legend for each line).
f) Use what you saw in e) to solve the problem (9) geometrically. This solution is called an optimal solution.
How large can you choose \( \alpha \) such that \( L_{\alpha} \) is nonempty?
g) Assume that we have other values on the revenues and costs than the actual numbers in a). Explain why (9), with these new parameter values, still has an optimal solution lying in a corner point of \( T \). Extend the program to calculate all the corner points of a region \( T \) in the plane determined by the linear inequalities like those listed above. Moreover, the program shall compute the maximum of a given linear function \( f(x,y)=ax+by \) over \( T \) by calculating the function values in the corner points and finding the smallest function value.
Filename: optimization.