This
exercise is the first introduction into the Fortran language. Here we are going
to deal with the basic concepts of variables, their different classes, the usage
of scalar and arrays and assignments. These issues provide one of the basic
building blocks of a numerical code. Initially there are a number of questions
related to understanding the various properties related to variables. Following
this the emphasis changes and you will start writing your first simple Fortran
programs.
Notice that you have online access to an older version of the
Metcalf f90 book, which contains all that is required for this
exercise. Names of variables:
In this part
of the exercise we investigate the syntax of simple expressions. Knowing the
correct syntax is important for writing the code correctly right from the start,
and will save you lots of time in the compiling phase.
Which of the
following names are allowed for Fortran variables:
Intrinsic data types:
Classify
the following constants according to their intrinsic data type.
Arrays:
Assume the following
definitions, what are the second, the fifth and the last array
element? Write the answer separated by a "," and without blank
spaces: ( as here a(#),a(#),a(#), where # is replaced with the appropriate
numbers).
Precision - Range:
Write the
specification that defines a scalar parameter (long) that can be used to define
the precision and range for other variables:
integer, parameter :: long=??????? real (kind=long) :: a,b,c.....
Value of expressions:
In the
following you have to determine the value of the expressions as a Fortran code
will evaluate them:
Array statements:
Using the
following definitions of the three arrays, a, b and c
real, dimension(6,9) :: a, b real, dimension(9) ::
c
Which of the following expressions are valid?
Download the exercise
files:
Before progressing to the next part of the exercises, you
have to download a file and unpack it content. To do this execute the following
commands in a command tool (xterm) window:
$ cd ~/Dat_F $ mkdir Exercise_2 $ cd Exercise_2
Then
download this file. Save it
in the new directory, and unpack the file:
$ tar -zxvf
fortran2.tgz $ rm
fortran2.tgz
For the rest of this exercise you will be work in the
"Exercise_2" directory.
First program:
The first program
you are going to write, or rather finalise, will determine the roots of a second
order equation:
a x² + b x + c = 0.
It is well known
that the mathematical solution to this equation can be written
as:
x = (-b +/- sqrt(b² - 4 a c))/ 2a
Open
second.f90 with your favorite editor. This file already contains the
skeleton of the program, but you will have to add the missing parts:
Define the needed variables correctly
Include missing read and print statement required to reading the
coefficients (a,b,c) and print the result x1, x2
Implement the solution given above
When this is done and you think
you have a code that can do the required task, compile the program by
typing:
$ ifort second.f90 -o
second.x
Depending on your success with implementing the above
conditions, you either get a smooth compilation of the code and the file second.x is created. If you have made one or
more errors, you will be presented with a list of error messages. These tells
you where and what may be wrong. Start to look at the errors from the top of the
list and correct them one by one. In many cases an errors in the code tricker
more errors further down in the program!
When it compiles without errors,
then try the following values as input for the three coefficients in the
equation:
a=1, b=2+eps, c=1
where eps represents a small deviation from the
value 2. What happens as the value of eps approaches 0? (try to change eps, starting with 0.1 and continue to
decrease it with a factor of 10 until the program fails.). To run the program,
you simply type:
$ ./second.x
and give the required input
(you can use the arrow keys to get back to the old commands and reuse
it).
If
we want to improve the precision of the calculations there is a simple way to do
this. You should know how from one of the previous questions. Implement this
change to the code such that the precision of the result improves to 15 digits.
There
are two routines you can add to the program that tells you what the actual
representation of the floating point numbers is (Check chapter 2.6 or 8.7. in
Metcalf or page 114 in Birtz.). Include them in the code to see what the
numerical representation of the real variables actually is.
Numerical precision:
In the
above example all coefficients were well behaved and the problem arises as the
determinant approaches zero when eps approaches zero. The problem is that
the numerical representation of eps eventually becomes to small compared to the
value of b in order to make a numerical change to the solution of the equation.
Another problem arises when the two contributions inside the determinant spans a
to large numerical range.
Try to make the same analysis as above, this
time setting a=1, c=-1 and b=10. This gives one positive and one negative root.
This time increase the value of b with a factor of 10 (using the low numerical
representation as you initially did in the first experiment). This time it goes
wrong when b=1e4 (e.i. one of the roots become zero.). The reason this time is
that the value of b² becomes so large that subtracting 4ac does not make a
difference to the number. Mathematically this is wrong, but remember that the
computer has a limited precision of numbers, and therefore when one adds a very
small number to a large number, their respective domains of precision may not
overlap and the operation will have no effect on the calculation. In this
particular case there is a way around the problem. We can write the second order
polynomial as:
(x - x{1})(x-x{2})
where x{1}
and x{2} represents the two roots. Evaluate this expression we
get
x² + (x_{1}+x_{2}) x + x_{1} x_{2}.
From
this we can identify
x_{1} x_{2}= c/a
From the
previous expression we can always get one numerically valid root. Knowing which
one (the numerically largest one), we can use the identity given above to
evaluate the second root. Now implement this change to deriving the solutions in
the code. Still using the default (low) definition of precision in the
calculations.
Let
this be a warning! There are situations where the most strait forward solution
makes mathematically sense, but when applied in a numerical approach may turnout
to give an unreliable solution. In many cases one can predict this and take
appropriate action to prevent wrong results.
Simple array handling:
As we
have only started the introduction to Fortran, it is limited what we can do with
arrays at present. One issue that can be addressed is the calculation of a
spacial derivatives of a function known in discreet points. Mathematically a
functional derivative is defined as the limit of the following expression as h
approaches zero:
df(x)/dx =
(f(x+h)-f(x-h))/(2 h)
From a numerical point of view there are a
number of issues that have importance for calculating a derivative of a
function, but this is not the time or place for discussing these. If we assume
the function f to be defined on a
uniform grid with a grid spacing h, and
the individual grid points are indexed by an i going from 1. to N. Then
we can define a numerical derivative using this
formula:
dfi = (f(xi+1) -
f(xi-1))/(2 h),
where the index i referees to
the array index in the grid. As you can see, then dfi depends
only on data points on either side of the point itself. This is a general
disadvantage, as the two sets of points (1+2 j) and (2+2 j) where j=0,N/2-1 do
not couple, and one can get a situation where the function looks like a sawtooth
function (alternating values between neighboring points) where the actual
derivative is large, but this approach gives a near zero value. This problem can
be eliminated if we instead calculate the derivative on a data grid that is
centered at half distance between the original grid.
The light green dots represent the
position of grid points, with the top line representing the centered grid, while
the lower line represents a half staggered grid. The formula calculating the
derivative using points on the top grid, returning the result on the lower grid
is given here:
dfi+0.5 = (f(xi+1)-
f(xi))/ h
This is a fist order approximation to the
derivative. One can also make higher order derivatives by combining more
gridpoints. A third order Taylor approximation to a first order derivative is
given by:
The file derivative.f90 contains the skeleton for
writing a code that calculates the derivative of the data given in the data file
function.dat. The data files contain the
grid and function values as 2 columns in ASCII format (the format gnuplot uses
for plotting). Expand the skeleton code, derivative.f90, to use the data to calculate
the three derivative approximations given above and the x values of the half
centered grid. Write the three datasets (x,df) to different files. Notice that
it is not possible to define the derivative for all gridpoints. There are a
number of points close to the ends of the data array that can't be defined,
unless periodicity of the data will be assumed - which we will ignore
here.
Use gnuplot to plot the three derivatives together with the
analytical data given in derivative.dat. Make a script file for the plot.
You can check the result against this image:
From this it can be seen that the
centered first order gives the worsted representation of the three, especially
in the regions where the derivative is changing fast. When you have got the
correct result, save the plot in png format with the filename derivatives.png
The
plot will be checked before you are credited for the 50 points.
The
underlying function for deriving the derivative shown in the plot above is given
by
