﻿ Programming for Physicists: Exercise 3

Exercise 2: Introduction to Fortran

Introduction:

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:

 a_real variable No Yes Credits: 2/-1
 a_real_variable No Yes Credits: 2/-1
 a-real_variable No Yes Credits: 2/-1
 1:variable No Yes Credits: 2/-1
 one:variable No Yes Credits: 2/-1
 one_very_long_variable_name No Yes Credits: 2/-1
 one_very_very_long_variable_name No Yes Credits: 2/-1

Intrinsic data types:

Classify the following constants according to their intrinsic data type.

 1.5 Integer Real Complex Character Logical Non Credits: 2/-1
 99 Integer Real Complex Character Logical Non Credits: 2/-1
 "(1.5, 2.4)" Integer Real Complex Character Logical Non Credits: 2/-1
 e4 Integer Real Complex Character Logical Non Credits: 2/-1
 3-1.5 Integer Real Complex Character Logical Non Credits: 2/-1
 4_6 Integer Real Complex Character Logical Non Credits: 2/-1
 (2.5, 1) Integer Real Complex Character Logical Non Credits: 2/-1

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

 real, dimension(10) :: a Credits: 2/-1
 real, dimension(-3:6) :: a Credits: 2/-1
 real, dimension(11:20) :: a Credits: 2/-1
 real, dimension(10,10) :: a Credits: 2/-1
 real, dimension(3,0:1,4) :: a Credits: 2/-1

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

 What does the ??? above have to be for an integer inside the range -1025 - 1025 Credits: 2/-1
 and for a real with 15 digits precision inside the range 10-120 - 10120 Credits: 2/-1

Value of expressions:

In the following you have to determine the value of the expressions as a Fortran code will evaluate them:

 5/4+2 Credits: 2/-1
 3.**4/2 Credits: 2/-1
 6**3/5 Credits: 2/-1
 2/4. - 1/4 Credits: 2/-1

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?

 a = b OK No Credits: 2/-1
 c = a(0,:) OK No Credits: 2/-1
 a(3,:) = c - 1.5 OK No Credits: 2/-1
 c = a(2,:) + 4*b(:,4) OK No Credits: 2/-1
 a(:,4) = c(2:7) + a(:,7) OK No Credits: 2/-1

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

 When does it go wrong? (value of eps) 0.001 0.00001 0.0000001 0.000000001 Credits: 2/-1
 Could you have predicted this? Yes, looking at the number of zeros written in the other answers. No, you are wrong! Credits: 1/-1
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.

 Fulfilling this criteria, what is the comparable exponential range (10^n) that the variables are defined in?  (give the value of n) Credits: 50/-1

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.

 How far can you go this time before the algorithm fails? 1e10 1e20 1e30 1e38 Credits: 50/-1

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:

dfi+0.5 = (f(xi+1)- f(xi))/(1.125 h)  +  (f(xi+2)- f(xi-1))/ (24 h)

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

 I have obtained the same results as the one given above Credits: 50/-1

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

f(x) = sin(2x) cos(3x)    (click to see a graphical representation of the function)

with a 50 points representation of the 2Pi x range.

Cleanup...

During the exercise you have created a number of files that
• take up disk space
• are easy to recreate
Therefore make it a good habit to cleanup redundant files when you are finished with a project.

\$Id: index.php,v 1.8 2009/09/11 05:47:15 kg Exp \$ kg@astro.ku.dk