Exercise 3: Conditional
In the previous
exercise you were introduced to the various data types and how to compose simple
arithmetic expressions using these. Restricted to these building blocks only
very simple programs can be written. The next important step is the ability to
control the execution flow of the program. This ability makes it possible to
solve much more complicated scenarios. The purpose of this exercise is to teach
you exactly this. There are a number of ways to control the flow in a program.
When writing a code, there is never one unique solution, but some ways of doing
it can be more elegant and efficient than others. Depending on your individual
choices while structuring the program, the result may appear very different from
you friends. If you both obtain the correct answers, then both program are
correct, but which solution is the best in general depends on other
circumstances. Such criteria can relate to the readability of the code, how user
friendly it is, cost for continued maintenance, the cost in development time and
also the expenses in CPU time when running the program.
The commands you
encounter below are the do, if, case,
where and forall statements. The
following exercises train you in using these, while you write a number of short
Getting the required files:
To get the files needed for this weeks exercise, you have to get
the following tar file.
the file in your Dat_F directory, and unpack it.
> cd ~/Dat-F
> tar -zxvf fortran3.tgz
> rm fortran3.tgz
creates a new directory ~/Dat_F/Exercise_3 that contains the files needed
for this exercise.
Second order equation: --- (working
with complex variables, do and if statements)
In the previous exercise
you wrote a program to solve the 2ed order equation. We only dealt with the
simple case having real solutions and saw that a stable solution could be
obtained by coding up the calculations in a particular way. The disadvantage
doing it that way was you need to know something about the solution to the
problem. Here you have to extent the program to handle a number of situations
automatically. In the Exercise_3 directory you have the same empty file
you started with last time, second.f90. This allows you to implement a
number of new requirements (these are not given in order, so read all of them
and decide in which order it is most easy to implement them!):
One point in writing a new code has to be mentioned here,
namely, that it is not cheating to copy parts of an old solution into a new
code. One only has to be careful doing it and remember that the pasted pieces of
code have to work together with the rest of the code!
- Deal with both real and complex solutions.
- Automatically take advantage of the trick introduced last time to get both
solutions with as high numerical precision as possible.
- Continue asking for new parameters until a default stop number is given
(say A=0, B=0, C=0)
- Test that it works for different parameters (testing all possible cases)
- Make it easy to change between different precisions - (need to recompile
To test that the
your program works correct for the different case, use the following
combinations for a, b and c -- using single precision for the variables:
lets leave the 2ed order equation for this time and go onto something
Driven oscillator ---
(working with do and if or case, scalars)
gives the angular
frequency and the period. This case is to idealised, it lacks the effect of
damping. In the simplest approximation equation (1) can be changed to
(2)where gamma represents a damping coefficient. When gamma is
positive the amplitude of the oscillation will decrease with time and eventually
the oscillations will stop. In this case the equation has been linearized (sin()
is represented by it's first order term in the Taylor series) and it is
therefore only valid for small angular amplitudes.
A typical way of
representing the solution, is to plot the trajectory in a phase diagram
representing (x,y) = (coordinate, velocity) of the pendulum. In such a diagram
the solution of (1) is represented by a circle, while the solution to (2)
becomes a spiral -- towards origin for positive gamma and infinity for negative
Finding the solution to (2) in phase space we rearrange (2) such
that the equation
and (4) represent two coupled ordinary differential equations that are easier to
solve than equation (2), while we at the same time get the solution in the
coordinates needed for plotting a phase diagram for the pendulum
In oscillation.f90 you have the skeleton for writing the
code that solves these equations. You have to
time integration you should try two different approaches. The simple
Euler integration where the solution is forwarded by simply doing
- Decide which input the code needs (start position, gamma, l, .....)
- Define all required variables
- Define how many integration time steps you need to take and a finish time
- Make a loop to do the time integration
- Save the updated (x,y) for later plotting with gnuplot
here the index n represents the time
When this works, use these values for (x,y)=(1,0), g=10, l=10,
gamma=0, t=10 and n=400 steps.
does the result changes as n is increased? One simple way to see this is to let
the program only integrate the phase space curve for a single resolution.
the integration is not all that accurate it tends to overshoot a little in the
orbit, you will measure the distance only in the x-direction after a a full
(5) does not take into account the curvature of the trajectory. To improve the
integration accuracy one can include information from the previous time steps. A
simple second order algorithm integration formula is given here:
n and n-1 on the time derivatives represents the corresponding
values at the appropriate times.
This formula has a problem with the
first time step. Here is how to solve that:
When it is working, you should use the same parameters as
for the first test case above [(x,y)=(1,0), g=10, l=10, gamma=0, t=10 and n=400
steps.] to answer this question.
- Use (5) as the first time step and then adopt (6) below for the following
- As equation (5) is not all that good it is custom to take a smaller
timestep with this formula than with the higher order to initiate the
integration. Here you should take 0.5 times the timestep used in the following
prove this is a much better choice for the integration repeat the convergence
test from above with the new version of the code:
different input values - especially change the value of gamma and plot the
result using gnuplot.
In nature you can find situations where an
oscillating system is driven by some external force. This can be represented in
a number of different ways. Here we take the simple representation where the
damping coefficient has been rewritten
Depending on the angular amplitude, the "damping"
coefficient will either amplify or damp the amplitude. One could directly
implement this into the formula you used before, but a more general approach is
to rewrite the equation by non dimensionalising the equation by taking
as the unit of the amplitude and
as the unit of time. The equation then
This make a small change to equation (4) such that it becomes
- Implement this in the program such that both versions work depending on
the given input
Try using different values for the initial
conditions to see how the phase trajectory changes depending on where in the
phase diagram the initial values are relative the closed trajectory.
- Make a gnuplot combining the solutions of
- the damped oscillator with parameters:
- l=10, g=10, x=1, y=0, gamma=0.1
- the forced oscillator with
- l=10, g=10, x=1, y=0, epsilon=1.
- and the free oscillator with
- l=10, g=10, x=1, y=0, gamma=0., epsilon=0
link shows the mathematical background for the Pendulum.
Game of Life -- (working
with do, where, cshift and arrays)
The game of life is not a game, it is a set of
simple rules that simulate the evolution of a population. It is an example of a
cellular automaton devised invented by the British mathematician John Horton
Conway in 1970. For the "game" you need a "large" 2D array on which you have a
distribution of initial cells containing "life". The distribution is change by
applying a set of rules by which the population number and distribution is
changed in time --- either growing, dying or oscillating.
The rules that
apply for our game are:
- Survival: Every grid point with two
or three neighbours survive for the next generation.
- Each grid point with 4 or more neighbours dies from over population.
- Each grid point with 1 or less neighbours dies from isolation.
- Birth: Each empty cell with exactly
three neighbours is a birth cell.
For each life circle these three
rules are applied to the data grid and the new generation on the grid is found.
The 2D setup implies that each cell has eight neighbours and we need to
test the above rules for each cell and keep track of the dead or birth of new
cells. A simple way of doing this is to use a number of 2D integer arrays and
manipulate these to contain the information need. The program life.f90
contains a skeleton for developing the program. You can use this, or if you want
start writing the program from scratch. Before starting you need a little more
Write a program
- random_number(number): is a
fortran standard routine that returns a random number between 0 and 1 in each
array entry. This can be used to initialise the "life" array by choosing the
locations above a given threshold to hold initial life.
- cshift(array,1,2): is a function
that takes an array and shift the values in different directions. The 1
represents a shift of one in positive index direction for the the second, 2,
index of the array. The shift assumes the array to be periodic in the given
index direction such that the last array value in the second index direction
is moved to the first index... If these are the index numbers then the value
associated with them are moved in the following manner:
1, 2, 3, 4, 5,
6, ...., n-1, n
n, 1, 2, 3, 4, 5, 6, ...., n-1
By changing the two
numbers the array can be shifted in different direction and in different array
indexes. A suitable combinations of cshift operations on the life array can be
used to give a variable that contains the number of neighbours for each grid
point. (You may also do this using the array syntax used in the derivative
routine from last exercise, but then it is a little more difficult to include
the assumption of a periodic domain. It requires two lines of code for each
cshift operation making the code less easy to read).
Plot the results using idl. In
principle it is possible to plot the data with gnuplot, but it requires a number
of items to work together, and as it is not a course in gnuplot, I am here take
the simplest approach....
- Dependent on a threshold value to pick out a number of cells to initially
- Evolve the grid for n life circles.
- Save each life circle to disk.
To plot the data in idl you have to do the
Start idl in the Exercise_2 directory by typing
inside idl you have to type the
following few lines:
a=lonarr(50,50) & number=128
IDL> close,10 &
IDL> for i=0,number-1 do begin readf,10,a &
The first line creates a 2D array of 50x50 elements and
defines the number of data snapshots to read. The second line first make sure
that the file pointer number 10 is free and then opens the data file and
associate the file pointer number 10 with the file. In the third line, we loop
over the number of snapshots in the file. For each loop the associated snapshot
is read and then plotted on the graphical screen. Each array getting it's own
50x50 pixel square on the display.
If you experience problems of the type
where it reads past the end of the file you can do
to bring you back again
to a workable setup. To save the image to a png file you do the
with your usrid',align=.5,/normal
get out of idl type
Choosing the value of the random threshold to be 0.8 for
identifying cells containing life. Then do 128 iterations to get a time
evolution like this, where each 50x50 pixels represent one instance in time. The
white dots represent the cells containing life. The way the random number
generator is used right now, it will each time reproduce the same numbers. It is
possible to change this to get different patterns by using the function
random_seed(). But we will not do that here.
the following two images are only obtained using ifort. Using any other compiler
and the initial random number sequence will be a different one, and therefore
also the evolution of the life game.
If you instead uses a threshold of
0.96 and 16 iterations you should get this:
Notice that the cyclic pattern at the
end of the series continues forever.
that you have reproduced the image above, we are going to change the initial
conditions and make a few specific tests of the code. There exist a number of
special structures that either oscillate or move systematically in a particular
direction. We are going to test one of these here, the Jelly fish. To do
this one you have to do the following:
- Change the grid resolution to 25x25.
- Allow for different initial conditions
- Use the following 5x5 image to reconstruct the initial condition where the
black pattern represents life cells.
- Place the initial 5x5 in the lower left corner of the data array.
- Include a change in the rules for the game of life:
- Survival: Every grid point with
two, four or five neighbours survive for the next generation.
- Each grid point with 6 or more neighbours dies from over population.
- Each grid point with 1 or less neighbours dies from isolation.
- Each populated grid point with 3 neighbours dies.
- Birth: Each empty cell with
exactly 3 neighbours is a birth cell.
- Run the experiment for 20 time steps and look at the image using idl.
Change the lonarr(50,50) to lonarr(25,25) and number to
20 in the lines above. Also change tvscl,a,i to
tvscl,rebin(a,n*25,n*25,/sample),i where n is an integer number that
expands the size of the image n times.
The Jelly fish is a
moving feature that has a cyclic behaviour. Use the idl plot to look at its
motion and how it changes its appearance.
are two, of many, links to more detailed information about the Game of Life
$Id: index.php,v 1.10 2009/09/15 08:24:43 kg Exp $ email@example.com