Exercises

The exercise numbers refer to the numbers used in the book Python Scripting for Computational Science, by H. P. Langtangen. Exercise numbers of the form X.y, where y is an integer, are extra exercises not found in the book.

• 0.0 Make sure your set-up is ok

  In a terminal window, do the following steps to install the software associated with this course:

  cd $HOME
  mkdir scripting
  cd scripting
  export scripting=`pwd`
  wget http://folk.uio.no/hpl/scripting/TCSE3-3rd-examples.tar.gz
  wget http://folk.uio.no/hpl/scripting/scitools.tar.gz
  gunzip TCSE3-3rd-examples.tar.gz scitools.tar.gz
  tar xvf TCSE3-3rd-examples.tar
  rm TCSE3-3rd-examples.tar
  tar xvf scitools.tar
  rm scitools.tar
  

  Insert the following commands in your .bashrc start-up file:

  export scripting=$HOME/scripting
  export PYTHONPATH=$scripting/src/tools:$scripting/scitools/lib
  PATH=$PATH:$scripting/src/tools:$scripting/scitools/bin
  

  Start a new terminal window, or log in and out again, and check that you can do the following session: following:

  Unix> python $scripting/src/py/intro/hw.py 3.14159
  Hello, World! sin(3.14159)=2.65358979335e-06
  
  Unix> python -c 'from scitools.all import *'
  scitools.easyviz backend is gnuplot
  
  Unix> python -c 'import Tkinter, Pmw'
  

  To get access to the oscillator code, used in the simviz*.py family of introductory scripts, run the commands

  Unix> cat > $scripting/src/tools/oscillator
  #!/bin/bash
  python $scripting/src/app/oscillator/Python/oscillator.py
  Ctrl-D
  Unix> chmod a+rx $scripting/src/tools/oscillator
  

  Now, try to run the code:

  Unix> oscillator
  

  The program just hangs (waiting for input). Kill it with Ctrl-C. If you encounter problems with installing the oscillator code, it is not that important - it means that you cannot run the simviz1.py script, but you can work with all the other examples.

• X.1 Use Python to do some serious calculations

  Write a Python script that computes 1+1 and writes out the result.

• 2.1 Become familiar with electronic Python documentation

  Write a script that prints a uniformly distributed random number between -1 and 1 on the screen. The number should be written with four decimals as implied by the %.4f format.

  To create the script file, you can use a standard editor such as Emacs or Vim on Unix-like systems. On Windows you must use an editor for pure text files - Notepad is a possibility, but I prefer to use Emacs or the ``IDLE'' editor that comes with Python (you usually find IDLE on the start menu, choose File--New Window to open up the editor). IDLE supports standard key bindings from Unix, Windows, or Mac (choose Options--Configure IDLE... and Keys to get a menu where you can choose between the three classes of key bindings).

  The standard Python module for generation of uniform random numbers is called random. To figure out how to use this module, you can look up the description of the module in the Python Library Reference. Load the file $scripting/doc.html into a Web browser and click on the link Python Library Reference: Index. You will then see the index of Python functions, modules, data structures, etc. Find the item random (standard module) in the index and follow the link. This will bring you to the manual page for the random module. In the bottom part of this page you will find information about functions for drawing random numbers from various distributions (do not use the classes in the module, use plain functions). Also apply pydoc to look up documentation of the random module: just write pydoc random on the command line.

• 2.2 Extend Exercise 2.1 with a loop

  Extend the script from Exercise 2.1 such that you draw n random uniformly distributed numbers, where n is given on the command line, and compute the average of these numbers.

• 2.3 Find five errors in a script

  The file src/misc/averagerandom2.py contains the following Python code:

  #!/usr/bin/ env python
  
  import sys, random
  def compute(n):
      i = 0; s = 0
      while i <= n:
          s += random.random()
           i += 1
      return s/n
  
  n = sys.argv[1]
  print 'the average of %d random numbers is %g" % (n, compute(n))
  

  There are five errors in this file - find them!

• 2.4 Basic use of control structures

  To get some hands-on experience with writing basic control structures in Python, we consider an extension of the Scientific Hello World script hw.py from Chapter 2.1. The script is now supposed to read an arbitrary number of command-line arguments and write the sine of each number to the screen. Let the name of the new script be hw2a.py. As an example, we can write

  python hw2a.py 1.4 -0.1 4 99
  

  and the program writes out

  Hello, World!
  sin(1.4)=0.98545
  sin(-0.1)=-0.0998334
  sin(4)=-0.756802
  sin(99)=-0.999207
  

  Traverse the command-line arguments using a for loop. The complete list of the command-line arguments can be written sys.argv[1:] (i.e., the entries in sys.argv, starting with index 1 and ending with the last valid index). The for loop can then be written as for r in sys.argv[1:].

  Make an alternative script, hw2b.py, where a while loop construction is used for handling each number on the command line.

  In a third version of the script, hw2c.py, you should take the natural logarithm of the numbers on the command line. Look up the documentation of the math module in the Python Library Reference to see how to compute the natural logarithm of a number. Include an if test to ensure that you only take the logarithm of positive numbers. Running, for instance,

  python hw2c.py 1.4 -0.1 4 99
  

  should give this output:

  Hello, World!
  ln(1.4)=0.336472
  ln(-0.1) is illegal
  ln(4)=1.38629
  ln(99)=4.59512
  

• 2.5 Use standard input/output instead of files

  Modify the datatrans1.py script such that it reads its numbers from standard input, sys.stdin, and writes the results to standard output, sys.stdout. You can work with sys.stdin and sys.stdout as the ordinary file objects you already have in datatrans1.py, except that you do not need to open and close them.

  You can feed data into the script directly from the terminal window (after you have started the script, of course) and terminate input with Ctrl-D. Alternatively, you can on Unix systems send a file into the script using a pipe, and if desired, redirect output to a file:

  cat inputfile | datatrans1stdio.py > res
  

  (datatrans1stdio.py is the name of the modified script.) A suitable input file for testing the script is src/py/intro/.datatrans_infile.

• 3.15 Annotate a filename with the current date

  Write a function that adds the current date to a filename. For example, calling the function with the text myfile as argument results in the string myfile_Aug22_2010 being returned if the current date is August 22, 2010. Read about the time module in the Python Library Reference to see how information about the date can be obtained. Exercise 3.16 has a useful application of the function from the present exercise, namely a script that takes backup of files (on your laptop or home computer) and annotes backup directories with the date.

• X.2 Use Fortran to do some serious calculations

  In Exercise X.1, let the computation of 1+1 be carried out in a Fortran function. The Python program should look like this:

  from f77module import calc
  print calc()
  

• X.3 Use Fortran to do some more serious calculations

  Write a Fortran subroutine that computes the sum and difference of two numbers, given as arguments to the subroutine. Call the subroutine from Python as follows:

  from f77module2 import sumdiff
  a = 1; b = -10
  sum, diff = sumdiff(a, b)
  print 'The sum of %g and %g is %g and the difference is %g' % (a, b, sum, diff)
  

• 5.1 Implement a numerical integration rule in F77

  Implement the Trapezoidal rule (see equation (4.1) in Exercise 4.5 in the course book) in Fortran 77 along with a function to integrate and a main program. Verify that the program works (check, e.g., that a linear function is integrated exactly, i.e., the error is zero to machine precision). Thereafter, interface this code from Python and write a new main program in Python calling the integration rule in F77 (the function to be integrated is still implemented in F77).

• 5.2 Implement a numerical integration rule in C

  As Exercise 5.1, but implement the numerical integration rule and the function to be integrated in C.

• 6.1 Modify the Scientific Hello World GUI

  

  Create a GUI as shown in figure above, where the user can write a number and click on sin, cos, tan, or sqrt to take the sine, cosine, etc. of the number. The GUI consists of an entry field, four buttons, and a label (with sunken relief). After the GUI is functioning, adjust the layout such that the computed number appears to the right in the label field. (Hint: Look up the man page for the Label widget. The 'Introduction to Tkinter' link in doc.html is a starting point.)

  Name of scriptfile: hwGUI1b.py (exercise class: GUI)

• 6.2 Change the layout of the GUI in Exercise 6.1

  Change the GUI in Exercise 6.1 on page \pageref{intro:pythontk:exer1} such that the layout becomes similar to the one in the figure below. Now there is only one input/output field (and you can work with only one StringVar or DoubleVar variable), just like a calculator. A Courier 18pt bold font is used in the text entry field.

  

  Name of scriptfile: hwGUI1c.py (exercise class: GUI)

• 4.1 Matrix-vector multiply with NumPy arrays

  Define a matrix and a vector, e.g.,

  A = array([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
  b = array([-3, -2, -1])
  

  Use the NumPy manual to find a function that computes the standard matrix-vector product A times b.

• 4.2 Work with slicing and matrix multiplication

  Extract the 2x2 matrix in the lower right corner of the matrix A in Exercise 4.1 as a slice. Add this slice to another 2x2, multiply the result by a 2x2 matrix, and insert this final result in the upper left corner of the original matrix A. Control the result by hand calculations.

• X.4 Compare two formulas in a plot

  The exact relation between Celsius and Fahrenheit degrees reads C = (F-32)*5/9. An approximate formula used by many people is C = (F-30)/2. Plot the two C(F) curves for F between -20 and 120. (Hint: watch out for integer division in the conversion formulas).

• X.5 Plot a wave packet

  The function f(x, t) = exp(-(x-3t)**2)*sin(3*pi*(x-t)) describes for a fixed value of t a wave localized in space. Make a program that visualizes this function as a function of x on the interval [-4,4] when t=0.

• X.6 Animate a wave packet

  Display an animation of the function f(x,t) in Exercise X.5 by plotting f as a function of x on [-6,6] for a set of t values in [-1,1]. Also make an animated GIF file. A suitable resolution can be 1000 intervals (1001 points) along the x axis, 60 intervals (61 points) in time, and 6 frames per second in the animated GIF file.

  (You will see that f(x,t) models waves that are moving to the right (when x is a space coordinate and t is time). The velocity of the individual waves and the packet is different, demonstrating a case in physics when the phase velocity of waves is different from the group velocity.)

• X.7 Judge a plot

  Assume you have the following program for plotting a parabola:

  from scitools.all import *
  x = linspace(0, 2, 20)
  y = x*(2 - x)
  plot(x, y)
  

  Then you alter the second statement to y = cos(18*pi*x). Judge the resulting plot. Is it correct? Display the cos(18*x) function with 1000 points in the same plot.

• X.11 Implement a very simple class

  Make a class Simple with one attribute i, one method double, which replaces the value of i by i+i, and a constructor that initializes the attribute. Try out the following code for testing the class:

  s1 = Simple(4)
  for i in range(4):
      s1.double()
  print s1.i
  
  s2 = Simple('Hello')
  s2.double();  s2.double()
  print s2.i
  s2.i = 100
  print s2.i
  

  Before you run this code, convince yourself what the output of the print statements will be.

• 8.26 Implement a class for vectors in 3D

  The purpose of this exercise is to program with classes and special methods. Create a class Vec3D with support for the inner product, cross product, norm, addition, subtraction, etc. The following application script demonstrates the required functionality:

  >>> from Vec3D import Vec3D
  >>> u = Vec3D(1, 0, 0)  # (1,0,0) vector
  >>> v = Vec3D(0, 1, 0)
  >>> str(u)        # pretty print
  '(1, 0, 0)'
  >>> repr(u)       # u = eval(repr(u))
  'Vec3D(1, 0, 0)'
  >>> u.len()       # Eucledian norm
  1.0
  >>> u[1]          # subscripting
  0.0
  >>> v[2]=2.5      # subscripting w/assignment
  >>> print u**v    # cross product
  (0, -2.5, 1)      # (output applies __str__)
  >>> u+v           # vector addition
  Vec3D(1, 1, 2.5)  # (output applies __repr__)
  >>> u-v           # vector subtraction
  Vec3D(1, -1, -2.5)
  >>> u*v           # inner (scalar, dot) product
  0.0
  

  We remark that class Vec3D is just aimed at being an illustrating exercise. Serious computations with a class for 3D vectors should utilize either a NumPy array (see Chapter 4), or better, the Vector class in the Scientific.Geometry.Vector module, which is a part of ScientificPython (see Chapter 4.4.1).

  Name of scriptfile: vec3d.py

• 8.27 Extend the class from Exericse 8.26

  Extend and modify the Vec3D class from Exericse 8.26 such that operators like + also work with scalars:

  u = Vec3D(1, 0, 0)
  v = Vec3D(0, -0.2, 8)
  a = 1.2
  u+v  # vector addition
  a+v  # scalar plus vector, yields (1.2, 1, 9.2)
  v+a  # vector plus scalar, yields (1.2, 1, 9.2)
  

  In the same way we should be able to do a-v, v-a, a*v, v*a, and v/a (a/v is not defined).

  Name of scriptfile: vec3d_ext.py

• 4.4 Vectorize a constant function

  The function

  def initial_condition(x):
      return 3.0
  

  does not work properly when x is a NumPy array. In that case the function should return a NumPy array with the same shape as x and with all entries equal to 3.0. Perform the necessary modifications such that the function works for both scalar types and NumPy arrays.

• X.8 Solve a linear system

  Solve the linear system A*x=b with respect to x, where the matrix A and the vector b are given as in Exercise 4.1.

• 2.9 Estimate the chance of an event in a dice game

  What is the probability of getting at least one 6 when throwing two dice? This question can be analyzed theoretically by methods from probability theory (the result is 1/6 + 1/6 - (1/6)*(1/6) = 11/36). However, a much simpler and much more general alternative is to let a computer program 'throw' two dice a large number of times and count how many times a 6 shows up. Such type of computer experiments, involving uncertain events, is often called Monte Carlo simulation (see also Exercise 4.14).

  Create a script that in a loop from 1 to n draws two uniform random numbers between 1 and 6 and counts how many times p a 6 shows up. Write out the estimated probability p/float(n) together with the exact result 11/36. Run the script a few times with different n values (preferably read from the command line) and determine from the experiments how large n must be to get the first three decimals (0.306) of the probability correct.

  Hint: Check out the random module to see how to draw random uniformly distributed integers in a specified interval.

  Remark. Division of integers (p/n) yields in this case zero, since p is always smaller than n. Integer p divided by integer n implies in Python, and in most other languages, integer division, i.e., p/n is the largest integer that when multiplied by n becomes less than or equal to p. Converting at least one of the integers to float (p/float(n)) ensures floating-point division, which is what we need.

• 4.8 Implement Exercise 2.9 using NumPy arrays

  Solve the same problem as in Exercise 2.9, but use Numerical Python and a vectorized algorithm. That is, generate a (long) random vector e of 2n uniform integer numbers ranging from 1 to 6, find the entries that are 6 by using where(e == 6, 1, 0), reshape the vector to a two-dimensional 2xn array, add the two rows of this array to a new array e2, count how many of the elements in e2 that are greater than zero (these are the events where at least one die shows a 6) by sum(where(e2 > 0, 1, 0)). Estimate the probability from this count. Insert CPU-time measurements in the scripts (see Chapter 4.1.4 or 8.10.1) and compare the plain Python loop and the standard random module with the vectorized version utilizing random, where, and sum from numpy.

• 2.10 Determine if you win or loose a hazard game

  Somebody suggests the following game. You pay 1 unit of money and are allowed to throw four dice. If the sum of the eyes on the dice is less than 9, you win 10 units of money, otherwise you loose your investment. Should you play this game? Hint: Use the simulation method from Exercise 2.9.

• 9.1 Extend Exercise 5.1 with a callback to Python

  Modify the solution of Exercise 5.1 such that the function to be integrated is implemented in Python (i.e., perform a callback to Python) and transferred to the Fortran code as a subroutine or function argument. Test different types of callable Python objects: a plain function, a lambda function, and a class instance with a __call__ method.

• 10.2 Investigate the efficiency of various DAXPY implementations

  A DAXPY operation performs the calculation u = a*x+y, where u, x, and y are vectors and a is a scalar. Implement the DAXPY operation in various ways:
  1. a plain Python loop over the vector indices
  2. a NumPy vector expression u = a*x + y
  3. a Fortran 77 or 90 program with a subroutine for computing u = a*x + y
  Run m DAXPY operations with vector length n, such that n=2**(2k), k=1,...,11, and mn=const (i.e., the total CPU time is ideally the same for each test). Show for each implementation the logarithm of the (normalized) CPU time versus the logarithm of n (in a plot or table). (The smallest arrays will probably lead to a blow-up of the CPU time of the Python implementations, and that is why it might be convenient to use the logarithm of the CPU time.)

  (The name DAXPY originates from the name of the subroutine in the standard BLAS library offering this computation.)

• X.9: Python-Fortran extension of Exercise 10.2

  Implement a Python interface to the Fortran 77 or Fortran 90 subroutines for computing u = a*x+y. Allocate the arrays in Python and send them to Fortran (let u be in/out and not only out, to avoid repeated allocation of a new u array every time the subroutine is called). Add experiments with the Python-Fortran version to the presentation you made in Exercise 10.2.