- A 11.0 Introduction
- A 11.1 Exercise 1 (Runge-Kutta)
- A 11.2 Exercise 2 (Horner's scheme and files)
- A 11.3 Exercise 3 (Factorial and Bessel)
- A 11.4 Exercise 4 (Factorial and Runge-Kutta)
- A 11.5 Exercise 5 (Linear systems of equations and files)

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 of the complete old name
`/undv/mai/numt/fort/kursbib ` or the new name
`/mailocal/lab/numt/TANA70/`

Locally is both
Fortran 77 (`f77`) and Fortran 90
(`f90`) available on our workstations.
Both the workstations and the
file servers are from Digital Equipment Corporation, they are DEC Alpha stations
running Digital UNIX.
Fortran 77 is the DEC compiler, Fortran 90 is also a DEC compiler.
Previously we used DEC MIPS stations, running DEC Fortran 77 and NAG Fortran 90.

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 have been replaced during 1996, from DEC Ultrix on MIPS machines to DEC UNIX on Alpha machines. The new computers also mean a change in location of the course directory. You define the course directory by writing

setenv KURSBIB /mailocal/lab/numt/TANA70/A suitable environment is obtained with the

Further information on the DEC Fortran 90, and its system parameters, hints on compilation, and some examples are available.

The following material shall be handed to the teacher:

- Listing of the translated program.
- 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.125For 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.If you wish to run this program you have to replace the first line

program RK1;with the following extended line

program RK1(INPUT, OUTPUT);when running on a DEC UNIX system, but not on a Sun Solaris system.

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:

- Listing of the corrected and modified program.
- Listing (at UNIX from script) of a run with the three test cases (x*x + 2 = 0),
`lab21.dat`and`lab22.dat`.

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

- Read the program and correct the errors you find.
- Run the program to see if it works satisfactorily.
- 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. - 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

- Correct the additional errors found now, and return to item 2.
- When the program seems to work properly, start to modify the program to handle also input data from a file.
- The commands to open and close a file are given in Appendix 2.
- 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)`. - 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)

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

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 -lF77and on DEC ULTRIX with

f90 prog.f -lnag -lfor -lutil -li -lots

On DEC UNIX ` -lfor -lutil -lots ` are available.

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.

Write a program in Fortran 77 for tabulating the Bosse function Bo(x). Write the main program in such a way that it asks for an interval for which the function values are to be evaluated.

Make a test run and print a table with *x* and *Bo(x)* for
*x* = 0.0, 0.1, ... 1.0.
Make the output to look nice! Investigate also what happens with the
arguments +800 and -900!

The Bosse function `BO ` is in the library `libbosse.a `
on the course directory and has the arguments
`X ` and `IFAIL ` and in this
order, where `X ` is the argument for the Bosse 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. If it
is 0 at exit, everything is OK, if it is 1 at exit, the argument
is too large, if it is 2 at exit, the argument is too small. If the
error parameter is zero on entry and an error occurs the program will stop with an
error message from the library error handling subroutine, and not return to
the user's program. The error handling is similar to the one in the
NAG Fortran 77 library.

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

f77 prog.f -lbosse -L$KURSBIBwhere

If you use Fortran 90 for your own program but you still have the Bosse library available only in Fortran 77 special care is necessary, see chapter 15, Method 2, see also the discussion above in exercise 3b.

To make life simpler we have however made also a Fortran 90 version of the Bosse library, which is used with

f90 prog.f90 -lbosse90 -L$KURSBIBIt is very important to note in which precision the Bosse library is made available on your system.

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

Compare with Exercise 3 a.

y'(t) = z(t) y(0) = 1 z'(t) = 2*y*(1+y^2) + t*sin(z) z(0) = 2Calculate an approximation to

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

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
`SAVE`d 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 -359A =-24 -10 -57b =281 -8 -4 -17 85

**Suggestion 1.**Test also on the trivial case with dimension equals 1 for both the matrix and the vector.**Suggestion 2.**The provided routine may only be changed from single to double precision. The user therefore has to use an`INTERFACE`.**Suggestion 3.**I recommend you to perform the first compilation one routine (program unit) at a time, in order not to be flooded with error messages. When all the program units compile four different ways of combining them into a program exist.I assume that the main program is on the file

`huvud.f90`, the subroutine`LASMAT`on the file`lasmat.f90`, and so on, and that a possible complete program is on the file`labb5.f90`.- Compile all the routines separately but at the same time
and with the same command,
f90 huvud.f90 lasmat.f90 lasvek.f90 ... los.f90 solve.f90

Execute witha.out

- Compile all the routines separately and use the object modules for linking,
f90 -c huvud.f90 f90 -c lasmat.f90 f90 -c lasvek.f90 ... f90 -c los.f90 f90 -c solve.f90 f90 huvud.o lasmat.o lasvek.o ... los.o solve.o

Execute witha.out

When one routine has been changed, now only that one has to be recompiled, followed by linking of all the routines. This method is much faster! - Combine all the files into one large and compile,
f90 labb5.f90

Execute witha.out

This method is very slow! - Compile all the routines separately but using the
`make`command. Execute with`labb5`. This is advanced UNIX! A suitable`makefile`is below and in the course directory.**This is the best method!**

- Compile all the routines separately but at the same time
and with the same command,

labb5: huvud.o lasmat.o lasvek.o matin.o vekin.o \ matut.o vekut.o los.o solve.o f90 -o labb5 huvud.o lasmat.o lasvek.o matin.o vekin.o \ matut.o vekut.o los.o solve.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 solve.o: solve.f90 f90 -c solve.f90The above is used by moving the file

If you have used other file names you have to make the appropriate modifications
to the file `makefile`.

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

At the the specification of the array to store the matrix, you may chose either the method in a) or b) above, or perhaps invent a third method.

If the `makefile`
is used, the reference to `solve.f90` must be changed.

Last modified: 27 June 1997