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

Exercises

Exercise 1: Demonstrate the magic of inheritance

Consider class Line from the section A class for straight lines and a subclass Parabola0 defined as

class Parabola0(Line):
    pass
That is, class Parabola0 does not have any own code, but it inherits from class Line. Demonstrate in a program or interactive session, using dir and looking at the __dict__ object, that an instance of class Parabola0 contains everything (i.e., all attributes) that an instance of class Line contains. Filename: dir_subclass.

Exercise 2: Make polynomial subclasses of parabolas

The task in this exercise is to make a class Cubic for cubic functions $$ \begin{equation*} c_3x^3 + c_2x^2 + c_1x + c_0\end{equation*} $$ with a call operator and a table method as in classes Line and Parabola from the section Inheritance and class hierarchies. Implement class Cubic by inheriting from class Parabola, and call up functionality in class Parabola in the same way as class Parabola calls up functionality in class Line.

Make a similar class Poly4 for 4-th degree polynomials $$ \begin{equation*} c_4x^4 + c_3x^3 + c_2x^2 + c_1x + c_0\end{equation*} $$ by inheriting from class Cubic. Insert print statements in all the __call__ methods such that you can easily watch the program flow and see when __call__ in the different classes is called.

Evaluate cubic and a 4-th degree polynomial at a point, and observe the printouts from all the superclasses. Filename: Cubic_Poly4.

Remarks

This exercise follows the idea from the section Inheritance and class hierarchies where more complex polynomials are subclasses of simpler ones. Conceptually, a cubic polynomial is not a parabola, so many programmers will not accept class Cubic as a subclass of Parabola; it should be the other way around, and Exercise 2: Make polynomial subclasses of parabolas follows that approach. Nevertheless, one can use inheritance solely for sharing code and not for expressing that a subclass is a kind of the superclass. For code sharing it is natural to start with the simplest polynomial as superclass and add terms to the inherited data structure as we make subclasses for higher degree polynomials.

Exercise 3: Implement a class for a function as a subclass

Implement a class for the function \( f(x)=A\sin (wx) + ax^2 + bx + c \). The class should have a call operator for evaluating the function for some argument \( x \), and a constructor that takes the function parameters A, w, a, b, and c as arguments. Also a table method as in classes Line and Parabola should be present. Implement the class by deriving it from class Parabola and call up functionality already implemented in class Parabola whenever possible. Filename: sin_plus_quadratic.

Exercise 4: Create an alternative class hierarchy for polynomials

Let class Polynomial from the document Introduction to classes in Python [2] be a superclass and implement class Parabola as a subclass. The constructor in class Parabola should take the three coefficients in the parabola as separate arguments. Try to reuse as much code as possible from the superclass in the subclass. Implement class Line as a subclass specialization of class Parabola.

Which class design do you prefer, class Line as a subclass of Parabola and Polynomial, or Line as a superclass with extensions in subclasses? (See also remark in Exercise 2: Make polynomial subclasses of parabolas.) Filename: Polynomial_hier.

Exercise 5: Make circle a subclass of an ellipse

The document Introduction to classes in Python [2] presents class Circle. Make a similar class Ellipse for representing an ellipse. Then create a new class Circle that is a subclass of Ellipse. Filename: Ellipse_Circle.

Exercise 6: Make super- and subclass for a point

A point \( (x,y) \) in the plane can be represented by a class:

class Point(object):
    def __init__(self, x, y):
        self.x, self.y = x, y

    def __str__(self):
        return '(%g, %g)' % (self.x, self.y)
We can extend the Point class to also contain the representation of the point in polar coordinates. To this end, create a subclass PolarPoint whose constructor takes the polar representation of a point, \( (r,\theta) \), as arguments. Store \( r \) and \( \theta \) as data attributes and call the superclass constructor with the corresponding \( x \) and \( y \) values (recall the relations \( x=r\cos\theta \) and \( y=r\sin\theta \) between Cartesian and polar coordinates). Add a __str__ method in class PolarPoint which prints out \( r \), \( \theta \), \( x \), and \( y \). Write a test function that creates two PolarPoint instances and compares the four data attributes x, y, r, and theta with the expected values. Filename: PolarPoint.

Exercise 7: Modify a function class by subclassing

Consider a class F implementing the function \( f(t;a,b) = e^{-at}\sin (bt) \):

class F(object):
    def __init__(self, a, b):
        self.a, self.b = a, b
    def __call__(self, t):
        return exp(-self.a*t)*sin(self.b*t)
We now want to study how the function \( f(t;a,b) \) varies with the parameter \( b \), given \( t \) and \( a \). Mathematically, this means that we want to compute \( g(b;t,a)=f(t;a,b) \). Write a subclass Fb of F with a new __call__ method for evaluating \( g(b;t,a) \). Do not reimplement the formula, but call the __call__ method in the superclass to evaluate \( f(t;a,b) \). The Fs should work as follows:

f = Fs(t=2, a=4.5)
print f(3)  # b=3

Hint. Before calling __call__ in the superclass, the data attribute b in the superclass must be set to the right value.

Filename: Fb.

Exercise 8: Explore the accuracy of difference formulas

The purpose of this exercise is to investigate the accuracy of the Backward1, Forward1, Forward3, Central2, Central4, Central6 methods for the function $$ \begin{equation*} v(x)={1-e^{x/\mu}\over 1-e^{1/\mu}}\tp\end{equation*} $$ Compute the errors in the approximations for \( x=0, 0.9 \) and \( \mu =1, 0.01 \). Illustrate in a plot how the \( v(x) \) function looks like for these two \( \mu \) values.

Hint. Modify the src/oo/Diff2_examples.py program which produces tables of errors of difference approximations as discussed at the end of the section Extensions.

Filename: boundary_layer_derivative.

Exercise 9: Implement a subclass

Make a subclass Sine1 of class FuncWithDerivatives from the section Superclass for defining an interface for the \( \sin x \) function. Implement the function only, and rely on the inherited df and ddf methods for computing the derivatives. Make another subclass Sine2 for \( \sin x \) where you also implement the df and ddf methods using analytical expressions for the derivatives. Compare Sine1 and Sine2 for computing the first- and second-order derivatives of \( \sin x \) at two \( x \) points. Filename: Sine12.

Exercise 10: Make classes for numerical differentiation

Carry out the exercise called Implement a class for numerical differentiation in the document Introduction to classes in Python [2]. Find the common code in the classes Derivative, Backward, and Central. Move this code to a superclass, and let the three mentioned classes be subclasses of this superclass. Compare the resulting code with the hierarchy shown in the section Classes for differentiation. Filename: numdiff_classes.

Exercise 11: Implement a new subclass for differentiation

A one-sided, three-point, second-order accurate formula for differentiating a function \( f(x) \) has the form $$ \begin{equation} f'(x)\approx {f(x-2h) -4f(x-h) + 3f(x)\over 2h}\tp \tag{17} \end{equation} $$ Implement this formula in a subclass Backward2 of class Diff from the section Class hierarchy for numerical differentiation. Compare Backward2 with Backward1 for \( g(t)=e^{-t} \) for \( t=0 \) and \( h=2^{-k} \) for \( k=0,1,\ldots, 14 \) (write out the errors in \( g'(t) \)). Filename: Backward2.

Exercise 12: Understand if a class can be used recursively

Suppose you want to compute \( f''(x) \) of some mathematical function \( f(x) \), and that you apply some class from the section Class hierarchy for numerical differentiation twice, e.g.,

ddf = Central2(Central2(f))
Will this work?

Hint. Follow the program flow, and find out what the resulting formula will be. Then see if this formula coincides with a formula you know for approximating \( f''(x) \) (actually, to recover the well-known formula with an \( h \) parameter, you would use \( h/2 \) in the nested calls to Central2).

Exercise 13: Represent people by a class hierarchy

Classes are often used to model objects in the real world. We may represent the data about a person in a program by a class Person, containing the person's name, address, phone number, date of birth, and nationality. A method __str__ may print the person's data. Implement such a class Person.

A worker is a person with a job. In a program, a worker is naturally represented as class Worker derived from class Person, because a worker is a person, i.e., we have an is-a relationship. Class Worker extends class Person with additional data, say name of company, company address, and job phone number. The print functionality must be modified accordingly. Implement this Worker class.

A scientist is a special kind of a worker. Class Scientist may therefore be derived from class Worker. Add data about the scientific discipline (physics, chemistry, mathematics, computer science, ...). One may also add the type of scientist: theoretical, experimental, or computational. The value of such a type attribute should not be restricted to just one category, since a scientist may be classified as, e.g., both experimental and computational (i.e., you can represent the value as a list or tuple). Implement class Scientist.

Researcher, postdoc, and professor are special cases of a scientist. One can either create classes for these job positions, or one may add an attribute (position) for this information in class Scientist. We adopt the former strategy. When, e.g., a researcher is represented by a class Researcher, no extra data or methods are needed. In Python we can create such an empty class by writing pass (the empty statement) as the class body:

class Researcher(Scientist):
    pass
Finally, make a demo program where you create and print instances of classes Person, Worker, Scientist, Researcher, Postdoc, and Professor. Print out the attribute contents of each instance (use the dir function).

Remark

An alternative design is to introduce a class Teacher as a special case of Worker and let Professor be both a Teacher and Scientist, which is natural. This implies that class Professor has two superclasses, Teacher and Scientist, or equivalently, class Professor inherits from two superclasses. This is known as multiple inheritance and technically achieved as follows in Python:

class Professor(Teacher, Scientist):
    pass
It is a continuous debate in computer science whether multiple inheritance is a good idea or not. One obvious problem in the present example is that class Professor inherits two names, one via Teacher and one via Scientist (both these classes inherit from Person). Filename: Person.

Exercise 14: Add a new class in a class hierarchy

a) Add the Monte Carlo integration method from the document Random numbers and simple games [6] as a subclass MCint in the Integrator hierarchy explained in the section Class hierarchy for numerical integration. Import the superclass Integrator from the integrate module in the file with the new integration class.

b) Make a test function for class MCint where you fix the seed of the random number generator, use three function evaluations only, and compare the result of this Monte Carlo integration with results calculated by hand using the same three random numbers.

c) Run the Monte Carlo integration class in a case with known analytical solution and see how the error in the integral changes with \( n=10^k \) function evaluations, \( k=3,4,5,6 \).

Filename: MCint_class.

Exercise 15: Compute convergence rates of numerical integration methods

Numerical integration methods can compute "any" integral \( \int_a^b f(x)dx \), but the result is not exact. The methods have a parameter \( n \), closely related to the number of evaluations of the function \( f \), that can be increased to achieve more accurate results. In this exercise we want to explore the relation between the error \( E \) in the numerical approximation to the integral and \( n \). Different numerical methods have different relations.

The relations are of the form $$ \begin{equation*} E = Cn^r,\end{equation*} $$ where and \( C \) and \( r < 0 \) are constants to be determined. That is, \( r \) is the most important of these parameters, because if Simpson's method has a more negative \( r \) than the Trapezoidal method, it means that increasing \( n \) in Simpson's method reduces the error more effectively than increasing \( n \) in the Trapezoidal method.

One can estimate \( r \) from numerical experiments. For a chosen \( f(x) \), where the exact value of \( \int_a^bf(x)dx \) is available, one computes the numerical approximation for \( N+1 \) values of \( n \): \( n_0 < n_1 < \cdots < n_N \) and finds the corresponding errors \( E_0,E_1,\ldots,E_N \) (the difference between the exact value and the value produced by the numerical method).

One way to estimate \( r \) goes as follows. For two successive experiments we have $$ \begin{equation*} E_i = Cn_i^r\tp\end{equation*} $$ and $$ \begin{equation*} E_{i+1} = Cn_{i+1}^r\end{equation*} $$ Divide the first equation by the second to eliminate \( C \), and then take the logarithm to solve for \( r \): $$ \begin{equation*} r = {\ln (E_{i}/E_{i+1})\over\ln (n_{i}/n_{i+1})}\tp \end{equation*} $$ We can compute \( r \) for all pairs of two successive experiments. Say \( r_i \) is the \( r \) value found from experiment \( i \) and \( i+1 \), $$ \begin{equation*} r_i = {\ln (E_{i}/E_{i+1})\over\ln (n_{i}/n_{i+1})},\quad i=0,1,\ldots,N-1\tp \end{equation*} $$ Usually, the last value, \( r_{N-1} \), is the best approximation to the true \( r \) value. Knowing \( r \), we can compute \( C \) as \( E_in_i^{-r} \) for any \( i \).

Use the method above to estimate \( r \) and \( C \) for the Midpoint method, the Trapezoidal method, and Simpson's method. Make your own choice of integral problem: \( f(x) \), \( a \), and \( b \). Let the parameter \( n \) be the number of function evaluations in each method, and run the experiments with \( n=2^{k}+1 \) for \( k=2,\ldots,11 \). The Integrator hierarchy from the section Class hierarchy for numerical integration has all the requested methods implemented. Filename: integrators_convergence.

Exercise 16: Add common functionality in a class hierarchy

Suppose you want to use classes in the Integrator hierarchy from the section Class hierarchy for numerical integration. to calculate integrals of the form $$ \begin{equation*} F(x) = \int_a^x f(t)dt\tp\end{equation*} $$ Such functions \( F(x) \) can be efficiently computed by the method from the exercise "Speed up repeated integral calculations" in the document Introduction to classes [2]. Implement this computation of \( F(x) \) in an additional method in the superclass Integrator. Test that the implementation is correct for \( f(x)=2x-3 \) for all the implemented integration methods (the Midpoint, Trapezoidal and Gauss-Legendre methods, as well as Simpson's rule, integrate a linear function exactly). Filename: integrate_efficient.

Exercise 17: Make a class hierarchy for root finding

Given a general nonlinear equation \( f(x)=0 \), we want to implement classes for solving such an equation, and organize the classes in a class hierarchy. Make classes for three methods: Newton's method, the Bisection method, and the Secant method.

It is not obvious how such a hierarchy should be organized. One idea is to let the superclass store the \( f(x) \) function and its derivative \( f'(x) \) (if provided - if not, use a finite difference approximation for \( f'(x) \)). A method

def solve(start_values=[0], max_iter=100, tolerance=1E-6):
    ...
in the superclass can implement a general iteration loop. The start_values argument is a list of starting values for the algorithm in question: one point for Newton, two for Secant, and an interval \( [a,b] \) containing a root for Bisection. Let solve define a list self.x holding all the computed approximations. The initial value of self.x is simply start_values. For the Bisection method, one can use the convention \( a, b, c \) = self.x[-3:], where \( [a,b] \) represents the most recently computed interval and \( c \) is its midpoint. The solve method can return an approximate root \( x \), the corresponding \( f(x) \) value, a boolean indicator that is True if \( |f(x)| \) is less than the tolerance parameter, and a list of all the approximations and their \( f \) values (i.e., a list of \( (x, f(x)) \) tuples). Filename: Rootfinders.

Exercise 18: Make a calculus calculator class

Given a function \( f(x) \) defined on a domain \( [a,b] \), the purpose of many mathematical exercises is to sketch the function curve \( y=f(x) \), compute the derivative \( f'(x) \), find local and global extreme points, and compute the integral \( \int_a^b f(x)dx \). Make a class CalculusCalculator which can perform all these actions for any function \( f(x) \) using numerical differentiation and integration, and a simple search method over a large number of function values for finding extrema.

Here is an interactive session with the class where we analyze \( f(x)=x^2e^{-0.2x}\sin (2\pi x) \) on \( [0,6] \) with a grid (set of \( x \) coordinates) of 700 points:

>>> from CalculusCalculator import *
>>> def f(x):
...     return x**2*exp(-0.2*x)*sin(2*pi*x)
...
>>> c = CalculusCalculator(f, 0, 6, resolution=700)
>>> c.plot()                # plot f
>>> c.plot_derivative()     # plot f'
>>> c.extreme_points()

All minima: 0.8052, 1.7736, 2.7636, 3.7584, 4.7556, 5.754, 0
All maxima: 0.3624, 1.284, 2.2668, 3.2604, 4.2564, 5.2548, 6
Global minimum: 5.754
Global maximum: 5.2548

>>> c.integral
-1.7353776102348935
>>> c.df(2.51)    # c.df(x) is the derivative of f
-24.056988888465636
>>> c.set_differentiation_method(Central4)
>>> c.df(2.51)
-24.056988832723189
>>> c.set_integration_method(Simpson)  # more accurate integration
>>> c.integral
-1.7353857856973565
Design the class such that the above session can be carried out.

Hint. Use classes from the Diff and Integrator hierarchies (the sections Class hierarchy for numerical differentiation and Class hierarchy for numerical integration) for numerical differentiation and integration (with, e.g., Central2 and Trapezoidal as default methods for differentiation and integration). The method set_differentiation_method takes a subclass name in the Diff hierarchy as argument, and makes a data attribute df that holds a subclass instance for computing derivatives. With set_integration_method we can similarly set the integration method as a subclass name in the Integrator hierarchy, and then compute the integral \( \int_a^bf(x)dx \) and store the value in the attribute integral. The extreme_points method performs a print on a MinMax instance, which is stored as an attribute in the calculator class.

Filename: CalculusCalculator.

Exercise 19: Compute inverse functions

Extend class CalculusCalculator from Exercise 18: Make a calculus calculator class to offer computations of inverse functions.

Hint. A numerical way of computing inverse functions is explained in the document Sequences and difference equations [7]. Other, perhaps more attractive methods are described in the exercises on inverse functions in the document Programming of ordinary differential equations [4].

Filename: CalculusCalculator2.

Exercises

Exercise 20: Make line drawing of a person; program

A very simple sketch of a human being can be made of a circle for the head, two lines for the arms, one vertical line, a triangle, or a rectangle for the torso, and two lines for the legs. Make such a drawing in a program, utilizing appropriate classes in the Shape hierarchy. Filename: draw_person.

Exercise 21: Make line drawing of a person; class

Use the code from Exercise 20: Make line drawing of a person; program to make a subclass of Shape that draws a person. Supply the following arguments to the constructor: the center point of the head and the radius \( R \) of the head. Let the arms and the torso be of length \( 4R \), and the legs of length \( 6R \). The angle between the legs can be fixed (say 30 degrees), while the angle of the arms relative to the torso can be an argument to the constructor with a suitable default value. Filename: Person.

Exercise 22: Animate a person with waving hands

Make a subclass of the class from Exercise 21: Make line drawing of a person; class where the constructor can take an argument describing the angle between the arms and the torso. Use this new class to animate a person who waves her/his hands. Filename: waving_person.

References

  1. H. P. Langtangen. Debugging in Python, \emphhttp://hplgit.github.io/primer.html/doc/pub/debug, http://hplgit.github.io/primer.html/doc/pub/debug.
  2. H. P. Langtangen. Introduction to classes in Python, \emphhttp://hplgit.github.io/primer.html/doc/pub/class, http://hplgit.github.io/primer.html/doc/pub/class.
  3. H. P. Langtangen. Unit testing with nose, \emphhttp://hplgit.github.io/primer.html/doc/pub/nose, http://hplgit.github.io/primer.html/doc/pub/nose.
  4. H. P. Langtangen. Programming of ordinary differential equations, \emphhttp://hplgit.github.io/primer.html/doc/pub/ode2, http://hplgit.github.io/primer.html/doc/pub/ode2.
  5. 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.
  6. H. P. Langtangen. Random numbers and simple games, \emphhttp://hplgit.github.io/primer.html/doc/pub/random, http://hplgit.github.io/primer.html/doc/pub/random.
  7. H. P. Langtangen. Sequences and difference equations, \emphhttp://hplgit.github.io/primer.html/doc/pub/diffeq, http://hplgit.github.io/primer.html/doc/pub/diffeq.