Programming for Physicists: Exercise 4

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

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

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

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

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:

Use (5) as the first time step and then adopt (6) below for the following integration steps
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 integration.

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)

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

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:

Survival: Every grid point with two or three neighbours survive for the next generation.
Deaths:
- 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 information.

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

Write a program that:

Dependent on a threshold value to pick out a number of cells to initially host life.
Evolve the grid for n life circles.
Save each life circle to disk.

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:

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

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