$$ \newcommand{\tp}{\thinspace .} $$




This chapter is taken from the book A Primer on Scientific Programming with Python by H. P. Langtangen, 5th edition, Springer, 2016.


Exercise 1: Make a function class

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)
>>> f.a = 2
>>> print f.value(pi)

Filename: F.

Exercise 2: Add a data attribute to a class

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.

Exercise 3: Add functionality to a class

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.

Exercise 4: Make classes for a rectangle and a triangle

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.

Exercise 5: Make a class for quadratic functions

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

The file with class 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.

Exercise 6: Make a class for straight lines

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.

Exercise 7: Flexible handling of function arguments

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.

Exercise 8: Wrap functions in a class

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.

Exercise 9: Flexible handling of function arguments

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.

Exercise 10: Deduce a class implementation

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.

Exercise 11: Implement special methods in a class

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)
>>> f.a = 2
>>> print f(pi)
>>> print f

Filename: F2.

Exercise 12: Make a class for summation of series

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.

Exercise 13: Apply a numerical differentiation class

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.

Exercise 14: Implement an addition operator

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

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.

Exercise 15: Implement in-place += 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.

Exercise 16: Implement a class for numerical differentiation

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)

d) Organize the file with the class and functions such that it can be used a module.

Filename: Central.

Exercise 17: Examine a program

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.

Exercise 18: Modify a class for numerical differentiation

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.

Exercise 19: Make a class for the Heaviside function

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.

Exercise 20: Make a class for the indicator function

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.

Exercise 21: Make a class for piecewise constant functions

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.

Exercise 22: Speed up repeated integral calculations

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.

Exercise 23: Apply a class for polynomials

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.

Exercise 24: Find a bug in a class for polynomials

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.

Exercise 25: Implement subtraction of polynomials

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.

Exercise 26: Test the functionality of pretty print of polynomials

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.

Exercise 27: Vectorize a class for polynomials

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.

Exercise 28: Use a dict to hold polynomial coefficients

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
                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.

Exercise 29: Extend class Vec2D to work with lists/tuples

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.


Filename: Vec2D_lists.

Exercise 30: Extend class Vec2D to 3D vectors

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.

Exercise 31: Use NumPy arrays in class Vec2D

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.

Exercise 32: Impreciseness of interval arithmetics

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).

Exercise 33: Make classes for students and courses

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.

Exercise 34: Find local and global extrema of a function

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:

Here is a sample code using class 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.

Exercise 35: Find the optimal production for a company

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.


  1. H. P. Langtangen. Object-oriented programming, \emphhttp://hplgit.github.io/primer.html/doc/pub/oo, http://hplgit.github.io/primer.html/doc/pub/oo.
  2. H. P. Langtangen. Variable number of function arguments in Python, \emphhttp://hplgit.github.io/primer.html/doc/pub/varargs, http://hplgit.github.io/primer.html/doc/pub/varargs.
  3. H. P. Langtangen. Debugging in Python, \emphhttp://hplgit.github.io/primer.html/doc/pub/debug, http://hplgit.github.io/primer.html/doc/pub/debug.
  4. H. P. Langtangen. Array computing and curve plotting, \emphhttp://hplgit.github.io/primer.html/doc/pub/plot, http://hplgit.github.io/primer.html/doc/pub/plot.
  5. H. P. Langtangen. Functions and branching, \emphhttp://hplgit.github.io/primer.html/doc/pub/funcif, http://hplgit.github.io/primer.html/doc/pub/funcif.
  6. H. P. Langtangen. User input and error handling, \emphhttp://hplgit.github.io/primer.html/doc/pub/input, http://hplgit.github.io/primer.html/doc/pub/input.