Appendix 11

Computer Exercises

This Appendix with Computer Exercises is also available in Swedish

A 11.0 Introduction

The source of this appendix is my complete Fortran 90 tutorial in Swedish, Lärobok i Fortran 90.

These computer exercises are small variants of those used at the teaching of Fortran at Linköping University, both for the students of technology and for the students in the mathematics and science program.

When I write course directory below it is intended as an advice to the local teacher to make the required files available to the students. Of course, it is also possible to type them from the textbook, extract them from this HTML-file, or get them electronically from our course directory. Using UNIX it is easy to introduce symbolic names, so I am calling the course directory $KURSBIB instead for the complete name /undv/mai/numt/fort/kursbib

Locally is only Fortran 77 (f77) available on our workstations , and Fortran 90 (f90) is only available on the file server limes. The local user therefore has to make a remote login, using the command rlogin limes. A good alternative is to only make a temporary visit to limes with the command rsh limes command, where command may be f90 lab4a.f90. All user files are automatically available also on the the file server. Both the workstations and the file server are from Digital Equipment Corporation, they are DEC MIPS stations. Fortran 77 is the DEC compiler, Fortran 90 is the NAG compiler.

Students of technology are supposed to solve the exercises 2 to 5, while mathematicians are supposed to solve the exercises 1, 2, 3, 4a and 5.

Further information on the course (including experiences from earlier students) is available (partly in Swedish, partly in English) on URL=http://www.nsc.liu.se/~boein/edu/edu.html

The computers used for teaching at the Department of Mathematics will be replaced during 1996. This may mean a change in location of the course directory.

A 11.1 EXERCISE 1

Runge-Kutta

The Pascal program below for solving an ordinary differential equation with the Runge-Kutta method is on the course directory, as the file $KURSBIB/lab1_e.p. The assignment is to translate it from Pascal to Fortran. The output is required to include not only the numeric values but also suitable explanatory text. It is also valuable if the program can use other input values than those given in the Pascal code, 1 and 2, respectively.

The following material shall be handed to the teacher:

  1. Listing of the translated program.
  2. Listing (for example using the UNIX script) from a complete run of the four test cases
        Number of steps      Step length

             1               1
             2               0.5
             4               0.25
             8               0.125
For those who do not know Pascal the recommended alternative is to get an elementary textbook on numerical methods and code Runge-Kutta directly in Fortran.
program RK1;
    (* A simple program in Pascal for Runge-Kutta's method for a
       first order differential equation.
                dy/dx = x^2 + sin(xy)
                y(1) = 2                                 *)
var     number, i : integer;
        h, k1, k2, k3, k4, x, y : real;
        function f(x,y : real) : real;
            begin
                        f := x*x + sin(x*y)
            end;
begin
        number := 1;
        while number > 0 do
        begin
            x := 1.0;
            y := 2.0;
            writeln(' Give the number of steps ');
            read(number);
            if number >= 1 then
                begin
                        writeln(' Give the step length ');
                        read(h);
                        writeln('       x              y');
                        writeln(x, y);
                        for i := 1 to number do
                        begin
                                k1 := h*f(x,y);
                                k2 := h*f(x+0.5*h,y+0.5*k1);
                                k3 := h*f(x+0.5*h,y+0.5*k2);
                                k4 := h*f(x+h,y+k3);
                                x  := x + h;
                                y  := y + (k1+2*k2+2*k3+k4)/6;
                                writeln(x, y);
                        end;
                end;
        end;
end.

A 11.2 EXERCISE 2

Horner's scheme and files

The program below is in Fortran 77 (or fix form Fortran 90) and is in the course directory as the file $KURSBIB/lab2_e.f. The assignment is first to correct the errors in the program. There are no errors in the numerical algorithm (Horner's scheme) or in the comments in the program. But some programming errors are included in the presented program, several of these are based on a mixture of Fortran and Pascal.

When the program is believed to be correct (test with the equation x*x + 2 = 0), the second assignment is to modify the program so that it can accept data from either the keyboard or from a file $KURSBIB/lab21.dat. The modifications are required to be done in such a way that the user can select if data are supposed to be input directly from the keyboard or from a file, in which case the user also will be required to give the name of the file.

Finally, check that the program also works for the file $KURSBIB/lab22.dat.

Note that the program uses complex numbers, and that such values are input as pairs within parenthesis.

The following material shall be handed to the teacher:

  1. Listing of the corrected and modified program.
  2. Listing (at UNIX from script) of a run with the three test cases (x*x + 2 = 0), lab21.dat and lab22.dat.
The file lab2_e.f for computer exercise 2. Find the errors in the program!
******************************************************************************
**      This program calculates all roots (real and/or complex) to an       **
**      N:th-degree polynomial with complex coefficients. (N <= 10)         **
**                                                                          **
**                 n        n-1                                             **
**       P(z) = a z  + a   z    + ... + a z + a                             **
**               n      n-1              1     0                            **
**                                                                          **
**      The execution terminates if                                         **
**                                                                          **
**                 1) Abs (Z1-Z0) < EPS       ==>   Root found = Z1         **
**                 2) ITER > MAX              ==>   Slow convergence        **
**                                                                          **
**      The program sets EPS to 1.0E-7 and MAX to 30                        **
**                                                                          **
**      The NEWTON-RAPHSON method is used:   z1 = z0 - P(z0)/P'(z0)         **
**      The values of P(z0) and P'(z0) are calculated with HORNER'S SCHEME. **
**                                                                          **
**      The array A(0:10) contains the complex coefficients of the          ** 
**      polynomial P(z).                                                    **
**      The array B(1:10) contains the complex coefficients of the          **
**      polynomial Q(z),                                                    **
**                  where P(Z) = (z-z0)*Q(z) + P(z0)                        **
**      The coefficients to Q(z) are calculated with HORNER'S SCHEME.       **
**                                                                          **
**      When the first root has been obtained with the NEWTON-RAPHSON 	    **
**      method, it is divided away (deflation).                             **
**      The quotient polynomial is Q(z).                                    **
**      The process is repeated, now using the coefficients of Q(z) as      **
**      input data.                                                         **
**      As a starting approximation the value STARTV = 1+i is used          **
**      in all cases.                                                       **
**      Z0 is the previous approximation to the root.                       **
**      Z1 is the latest approximation to the root.                         **
**      F0 = P(Z0)                                                          **
**      FPRIM0 = P'(Z0)                                                     **
******************************************************************************
      COMPLEX      A(0:10), B(1:10), Z0, Z1, STARTV
      INTEGER      N, I, ITER, MAX
      REAL         EPS
      DATA  EPS/1E-7/, MAX /30/, STARTV /(1,1)/
******************************************************************************
20    WRITE(6,*)  'Give the degree of the polynomial'
      READ (5,*) N
      IF (N .GT. 10) THEN
         WRITE(6,*) 'The polynomial degree is maximum 10'
      GOTO 20
      WRITE (6,*) 'Give the coefficients of the polynomial,',
     , ' as complex constants'
      WRITE (6,*) 'Highest degree coefficient first'
      DO 30 I = N, 0, -1
          WRITE (6,*) 'A(' , I, ') = '
          READ  (5,*)  A(I)
30    CONTINUE
      WRITE (5,*) '    The roots are','        Number of iterations'
******************************************************************************
40    IF (N GT 0) THEN
C     ******     Find the next root                                     ******
          Z0 = (0,0)
          ITER = 0
          Z1 = STARTV
50        IF (ABS(Z1-Z0) .GE. EPS) THEN 
C         ++++++ Continue the iteration until two estimates             ++++++
C         ++++++ are sufficiently close.                                ++++++
              ITER = ITER + 1
              IF (ITER .GT. MAX) THEN
C             ------ Too many iterations  ==> TERMINATE                 ------
                     WRITE (6,*) 'Too many iterations.'
                     WRITE (6,*) 'The latest approximation to the root is ',Z1
              GOTO 200
              ENDIF
              Z0 = Z1
              HORNER (N, A, B, Z0, F0, FPRIM0)
C             ++++++   NEWTON-RAPHSON'S METHOD                          ++++++
              Z1 = Z0 - F0/FPRIM0
          GOTO 50
          ENDIF
C         ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
100       WRITE (6,*) Z1, '         ',Iter
C         ****** The root is found. Divide it away and seek the next one *****
          N = N - 1
          FOR I = 0 TO N DO
              A(I) = B(I+1)
      GOTO 40
      ENDIF
200   END

      SUBROUTINE  HORNER(N, A, B, Z, F, FPRIM)
******    For the meaning of the parameters - see the main program      ******
******    BI and CI are temporary variable.                             ******

      INTEGER N, I
      COMPLEX A(1:10), B(0:10), Z, F, FPRIM, BI, CI

      BI = A(N)
      B(N) = BI
      CI = BI

      DO 60  I = N-1, 1, -1
             BI = A(I) + Z*BI
             CI = BI + Z*CI
             B(I) = BI
C            ++++++ BI = B(I) for the calculation of P(Z)               ++++++
C            ++++++ CI  for the calculation of P'(Z)                    ++++++
60    CONTINUE

      F = A(0) + Z*BI
      FPRIM = CI

      RETURN
      END
******   END HORNER'S SCHEME                                            ******

******   This program is composed by Ulla Ouchterlony 1984              ******
******   The program is translated by Bo Einarsson 1995                 ******
Comments. I recommend the following process at the solution of this exercise.
  1. Read the program and correct the errors you find.
  2. Run the program to see if it works satisfactorily.
  3. If it does not work properly, correct any new errors and return to item 2. If you did not find any new errors, recompile the program using a check of not declared variables. This can on our DEC system be done with
                    f77 -u lab2_e.f
                    a.out
    
    Using a modern system, for example Fortran 90, you put the statement IMPLICIT NONE first in booth the main program and the subroutine.
  4. Correct the additional errors found now, and return to item 2. If you did not discover any new errors, recompile using checking of array indices. On some systems, as on our DEC system and with the NAG Fortran 90, this can be done with
                    f77 -C lab2_e.f
                    a.out
    
  5. Correct the additional errors found now, and return to item 2.
  6. When the program seems to work properly, start to modify the program to handle also input data from a file.
  7. The commands to open and close a file are given in Appendix 2.
  8. At list-directed input of text the text characters should preferably be given within apostrophes (is required by some systems). This is not user-friendly! I therefore recommend formatted input in this case. If only one character is to be given, use FORMAT(A1), if thirteen characters are to be given, use FORMAT(A13).
  9. Remember that each read or write statement starts with a new record (line). This is also the case within an explicit DO-loop. If several data values are on the same line, they have to be read with the same statement. This can be achieved by using the implicit DO-loop.

    The explicit DO-loop below will write seven lines, each one with one value of X

    	DO I = 1, 13, 2
    	   WRITE(*,*) X(I)
    	END DO
    

    The implicit DO-loop below will write only one line, but with seven values of X

     	WRITE(*,*) ( X(I), I = 1, 13, 2)
    

A 11.3 EXERCISE 3

Factorial and Bessel

This is computer exercise number three, and it consists of writing two small programs in the programming language Fortran 77.

Exercise 3 a)

Write a function in Fortran 77 with a main program for evaluation of the factorial function. Use only integers. Write the main program in such a way that it asks the user for an integer, for which the factorial is to be calculated. Make a test run and evaluate 10!

Exercise 3 b)

Write a program in Fortran 77 for tabulating the Bessel function J0(x). Use for example the NAG library (or a similar library) for obtaining values of this function. Write the main program in such a way that it asks for an interval for which the functions is to be evaluated.

Make a test run and print a table with x and J0(x) for x = 0.0, 0.1, ... 1.0. Make the output to look nice!

At the use of the NAG library the NAG User Notes are very useful. They describe the system specific parts of the library. Some of these are available here as

The required Bessel function in the NAG library is S17AEF and has the arguments X and IFAIL and in this order, where X is the argument for the Bessel function in single or double precision, and where IFAIL is an error parameter (integer). The error parameter is given the value 1 before entry, and is checked at exit. It it is 0 at exit, everything is OK, if it is 1 at exit, the magnitude of the argument is too large.

If the NAG library is installed in the normal way it is linked with the command

        f77 prog.f -lnag
where prog.f is your program.

If you use Fortran 90 for your own program but has the NAG library available only in Fortran 77 special care is necessary, see chapter 15, Method 2.

On Sun you compile and link with

        f90 prog.f -lnag -lF77
and on DEC with
        f90 prog.f -lnag -lfor -lutil -li -lots
It is very important to note in which precision the NAG library is made available on your system. Using the wrong precision will usually give completely wrong results!

Local information: The NAG library is in double precision on the DEC system.

A 11.4 EXERCISE 4

Factorial and Runge-Kutta

This is computer exercise number four, and means full transition from Fortran 77 to Fortran 90. Free form of the source code is required!

Exercise 4 a)

Write a recursive function in Fortran 90 with a main program for evaluation of the factorial function. Use only integers. Write the main program in such a way that it asks the user for an integer, for which the factorial is to be calculated. Make a test run and evaluate 10!

Compare with Exercise 3 a.

Exercise 4 b)

Write a program in Fortran 90 for the numerical solution of a system of ordinary differential equations using the classical Runge-Kutta method. The following system is given
        y'(t) = z(t)                        y(0) = 1
        z'(t) = 2*y*(1+y^2) + t*sin(z)      z(0) = 2
Calculate an approximation to y(0.2) using the step size 0.04. Repeat the calculation with the step sizes 0.02 and 0.01. The functions representing the differential functions are to be given as internal functions. It is therefore not permitted (in this exercise) to use the usual external functions.

Starting values of t, y and z are to be given interactively, as well as the step size h.

Note that the requirement not to use external functions makes it a little more difficult to use a subprogram for the Runge-Kutta steps (in this case also the subprogram has to be internal). The aim of this requirement is to give the student an opportunity to use the new internal functions of Fortran 90.

The formulas for Runge-Kutta at systems are available here as an dvi-file and here as an HTML-file (which requires indices)!

A 11.5 EXERCISE 5

Linear systems of equations and files

This is computer exercise number five. It consists of writing in Fortran 90 one main program and a number of subprograms, all in double precision, for the solution of a common numerical problem. The problem is to solve a linear system of equations Ax = b, using an available routine for the numerical task. The available routine is however given in single precision, and therefore requires manual conversion into double. At least four hours are suggested for this exercise.

The first subprogram LASMAT is used to interactively input the floating point values of a quadratic matrix. The subprogram is required to ask for the dimension and give the user the choice of providing from the keyboard either all the values of the matrix, or only those which are different from zero, The matrix shall then be stored in a file, using maximal accuracy and minimal storage space. The possible sparsity of the matrix is however not permitted to be exploited at this storage process. The subprogram shall ask for the name of the file, but the file type (extension) will be .mat. The user is not permitted to give the file type.

The second subprogram LASVEK is used to interactively input the floating point values of a vector. The subprogram is required to ask for the dimension and give the user the choice of providing from the keyboard either all the values of the vector, or only those which are different from zero, The vector shall then be stored in a file, using maximal accuracy and minimal storage space. The possible sparsity of the vector is however not permitted to be exploited at this storage process. The subprogram shall ask for the name of the file, but the file type (extension) will be .vek. The user is not permitted to give the file type.

The third subprogram MATIN reads a matrix from the file with the given name and stores it as an array (matrix) in the subprogram. The user is not permitted to give the file type.

The fourth subprogram VEKIN reads a vector from the file with the given name and stores it as an array (vector) in the subprogram. The user is not permitted to give the file type.

The fifth subprogram MATUT prints a matrix on paper, with rows and columns in the normal way if the dimension is at most 10, and in an easily understandable way if the dimension is larger than 10. The available 132 positions of a typical line printer are to be utilized as well as possible. Output on paper is under UNIX not done directly, so you will have to first write to a file, which is then moved to paper with the UNIX print command lpr or preferably fpr or asa, see the end of Appendix 2. In order to get 132 positions you can usually add -w132 to the print command.

The sixth subprogram VEKUT prints a vector on paper, preferably in the transposed form (line-wise).

The seventh subprogram LOS solves the system of equations through a call to the provided routine SOLVE_LINEAR_SYSTEM, which is given in double precision. No changes are permitted in that routine! If the dimension of the matrix and the vector do not agree, then the subprogram LOS or the main program is required to give an error message. An error message is also required if the routine SOLVE_LINEAR_SYSTEM detects that the matrix is singular.

The main program HUVUD utilizes LASMAT, LASVEK, MATIN and VEKIN in order to get the matrix A and the vector b, calls the solver subprogram LOS and prints the matrix A , the right hand side b and the solution x with the subprograms MATUT and VEKUT.

The dimension in all the programs shall use the dynamic possibilities of Fortran 90. One alternative is to use pointers at the specification, see chapter 12. Another possibility, suggested by Arie ten Cate is to introduce SAVEd arrays in a MODULE. A simple example of this is available.

The program is required to return to the starting point after each calculation, and ask for a new matrix and/or vector. It has to be possible to use a stored set of matrix and vector, without manual input of a lot of values.

In the exercise it is included the essential task of performing any additional specifications, test the reasonableness of the given values, and to construct suitable test matrices and vectors.

A compulsory test example follows

		 33	 16	 72			-359
	A =	-24	-10	-57		b =	 281
		 -8	 -4	-17			  85
The students are requried to hand in a complete program listing, and results from actual runs with matrices of the order 3, 10, and 12, including input of a sparse matrix.
labb5:  huvud.o lasmat.o lasvek.o matin.o vekin.o \
                matut.o vekut.o los.o loes.o
        f90 -o labb5 huvud.o lasmat.o lasvek.o matin.o vekin.o \
                matut.o vekut.o los.o loes.o
huvud.o: huvud.f90
        f90 -c huvud.f90
lasmat.o: lasmat.f90
        f90 -c lasmat.f90
lasvek.o: lasvek.f90
        f90 -c lasvek.f90
matin.o: matin.f90
        f90 -c matin.f90
vekin.o: vekin.f90
        f90 -c vekin.f90
matut.o: matut.f90
        f90 -c matut.f90
vekut.o: vekut.f90
        f90 -c vekut.f90
los.o:  los.f90
        f90 -c los.f90
loes.o: loes.f90
        f90 -c loes.f90
The above is used by moving the file makefile to your directory and when you wish to compile you write only make in the terminal window. Then only those files that have been changed as compared with the previous make command are recompiled (or all files if it is the first time), then all the object modules are linked to a program, ready to run with the command labb5.

If you have used other file names you have to make the appropriate modifications to the file makefile. Especially note that the provided routine is loes.f90 in the Swedish case but solve.f90 in the English case.

In principle make works in such a way that if an item which is after a colon : has a later time than that before the colon, the command on the next line is performed.

Remark. The backslash \ indicates a continuation line in UNIX.

The routine SOLVE_LINEAR_SYSTEM is available in single precision in the course directory with the name solve1.f90 and in double precision with the name solve.f90, and looks as follows. The single precision version is discussed at length in the Solution to Exercise (11.1).

SUBROUTINE SOLVE_LINEAR_SYSTEM(A, X, B, ERROR)
        IMPLICIT NONE
        ! Array specifications
        DOUBLE PRECISION, DIMENSION (:, :), INTENT (IN)  :: A
        DOUBLE PRECISION, DIMENSION (:),    INTENT (OUT) :: X
        DOUBLE PRECISION, DIMENSION (:),    INTENT (IN)  :: B
        LOGICAL, INTENT (OUT)                            :: ERROR

        ! The work area M is A extended with B
        DOUBLE PRECISION, DIMENSION (SIZE (B), SIZE (B) + 1) :: M
        INTEGER, DIMENSION (1)                               :: MAX_LOC
        DOUBLE PRECISION, DIMENSION (SIZE (B) + 1)           :: TEMP_ROW 
        INTEGER                                              :: N, K

        ! Initiate M
        N = SIZE (B)
        M (1:N, 1:N) = A
        M (1:N, N+1) = B

        ! Triangularization phase
        ERROR = .FALSE.
        TRIANG_LOOP: DO K = 1, N - 1
                ! Pivotering
                MAX_LOC = MAXLOC (ABS (M (K:N, K)))
                IF ( MAX_LOC(1) /= 1 ) THEN
                    TEMP_ROW (K:N+1 ) = M (K, K:N+1)
                    M (K, K:N+1) = M (K-1+MAX_LOC(1), K:N+1)
                    M (K-1+MAX_LOC(1), K:N+1) = TEMP_ROW (K:N+1)
                END IF

                IF (M (K, K) == 0.0D0) THEN
                        ERROR  = .TRUE.     ! Singular matrix A
                        EXIT TRIANG_LOOP
                ELSE
                        TEMP_ROW (K+1:N) = M (K+1:N, K) / M (K, K)
                        M (K+1:N, K+1:N+1) = M (K+1:N, K+1:N+1)  &
                           - SPREAD( TEMP_ROW (K+1:N), 2, N-K+1) &
                           * SPREAD( M (K, K+1:N+1), 1, N-K)
                        M (K+1:N, K) = 0 ! These values are not used
                END IF
        END DO TRIANG_LOOP

        IF (M (N, N) == 0.0D0) ERROR  = .TRUE.  ! Singular matrix A

        ! Back substitution
        IF (ERROR) THEN
           X = 0.0D0
        ELSE
           DO K = N, 1, -1
              X (K) = M (K, N+1) - SUM (M (K, K+1:N) * X (K+1:N))
              X (K) = X (K) / M (K, K)
           END DO
        END IF

END SUBROUTINE SOLVE_LINEAR_SYSTEM


Last modified: 27 January 1996
boein@nsc.liu.se