Exercise 3: Conditional controls




Introduction:

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



Getting the required files:

To get the files needed for this weeks exercise, you have to get the following tar file.
Save the file in your Dat_F directory, and unpack it.

  > cd ~/Dat-F
  > tar -zxvf fortran3.tgz
  > rm fortran3.tgz

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

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:

a=1, b=2, c=1 Credits: 2/-1
a=1, b=10000, c=1 Credits: 2/-1
a=1, b=1, c=1 Credits: 2/-1
a=0, b=1, c=1 Credits: 2/-1

OK, lets leave the 2ed order equation for this time and go onto something new.



Driven oscillator --- (working with do and if or case, scalars)

     The simple pendulum is a well known free oscillator. It consist of a mass m in a gravitational field of acceleration g, suspended from a point O by a rigid wire of length l. The mass oscillate in the vertical plane. If theta(t) measures the angle between the vertical line through O, then using Newton's 2ed law, a 2ed order ordinary differential equation for theta can be written:
    (1)
This equation has the solution,
where
and
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 gamma.

Finding the solution to (2) in phase space we rearrange (2) such that the equation becomes:

                                                              (3)
            (4)
Equations (3) 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 motion.

In oscillation.f90 you have the skeleton for writing the code that solves these equations. You have to
For the time integration you should try two different approaches. The simple Euler integration where the solution is forwarded by simply doing
  (5)

here the index n represents the time step.
When this works, use these values for (x,y)=(1,0), g=10, l=10, gamma=0, t=10 and n=400 steps.

For the case with gamma=0, what does the curve looks like? Credits: 10/-5

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

How long does it take the pendulum to make a full resolution? Credits: 2/-1

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

For the 1st order algorithm, how many integration steps does it take to reach a distance in x below 1e-3 after one period? Credits: 10/-5

Equation (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:

(6)
The indexes 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.

What does the curve looks like now? Credits: 10/-5

To prove this is a much better choice for the integration repeat the convergence test from above with the new version of the code:

For the 2nd order algorithm, how many integration steps does it take to get the distance in x below 1e-3 after one period? Credits: 10/-5


Try 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

(7)


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 becomes,
(8)

This make a small change to equation (4) such that it becomes
(9)


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.



I have created the oscillator.png image. Credits: 20/-10

How long time does it take the pendulum to do 5 full oscillations in the decaying case? Credits: 5/-3

How many periods does the forced pendulum in equilibrium do in a quarter of the previous questions time? Credits: 5/-3

This 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:
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 information.
Write a program that:
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....
To plot the data in idl you have to do the following:
Start idl in the Exercise_2 directory by typing

$ idl

inside idl you have to type the following few lines:

IDL> a=lonarr(50,50) & number=128
IDL> close,10 & openr,10,'life.dat'
IDL> for i=0,number-1 do begin readf,10,a & tvscl,a,i

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

IDL> retall

to bring you back again to a workable setup. To save the image to a png file you do the following:

IDL> xyouts,.5,.01,'replace with your usrid',align=.5,/normal
IDL> write_png,'life.png',tvrd()

To get out of idl type

IDL> exit

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.

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

I have created the life.png image. Credits: 50/-25


Now 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:


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.

What is the length of the life circle for the jelly fish? Credits: 50/-25


Below are two, of many, links to more detailed information about the Game of Life
http://en.wikipedia.org/wiki/Conway%27s_Game_of_Life
http://www.math.com/students/wonders/life/life.html


$Id: index.php,v 1.10 2009/09/15 08:24:43 kg Exp $ kg@astro.ku.dk