Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.mrao.cam.ac.uk/~rachael/compphys/CompPhys1.pdf Äàòà èçìåíåíèÿ: Fri Sep 21 15:20:19 2007 Äàòà èíäåêñèðîâàíèÿ: Tue Oct 2 05:34:02 2012 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: http astrokuban.info astrokuban

University of Cambridge Department of Physics

Computational Physics
Handout 1 Introduction to Computational Physics

Dr Rachael Padman Michaelmas 2007


1. Introduction
This course aims to expose you to the types of computational problems that occur in real scientific work. We will use the language Fortran 90/95 as a vehicle for learning how to compute effectively and will encourage the use of program packages such as the NAG library. It is also an opportunity for you to develop skills at working reasonably independently on a project. This course will assume no prior knowledge of programming or Unix, and no knowledge of computing beyond that already taught in Part I. These skills will be acquired by working through a set of self-study guides during the afternoon sessions in the Physics PWF when there will be demonstrators available to provide help and advice. 1.1 What is Computational Physics? We can split computational physics (roughly) into three types of problem: A. Some problems can be solved numerically which are intractable analytically. This is a major area of computational physics. In particular problems involving the solution of differential equations and/or large linear algebra problems: · a differential equation which can only be solved numerically · an eigenvalue problem involving very large matrices Examples · SchrÆdinger's equation · fluid mechanics · diffusion equation B. Some problems are difficult to reduce to a set of differential equations. We may have some idea of the physics, but it may be difficult to formulate (let alone solve) the problem analytically. This is often the case with highly non-linear systems. In these cases one uses the technique of simulation. Here we approximate the system to a computational model. The difference from A is that each run of the simulation will have to be treated much like the results of an experiment. We will need to average the results of these computer experiments (or make them large enough to be representative) to model the physics. Examples · N-body simulations of galaxies and structure formation in the Universe · Molecular Dynamics · Cellular Automata (problem is highly abstracted) C. The analysis of data using fitting a complicated function minimisation. Here we seek to analytically or algorithmically. (e.g. local or false minima). sophisticated numerical techniques. For example this might involve to real (or perhaps simulated) data. This introduces us to the idea of find a (usually global) minimum in a function which may be described Many techniques are available to us to avoid finding incorrect answers

Examples · Monte Carlo techniques · Simulated Annealing

1


1.2 Relationship to Numerical Analysis In solving these problems we will often call upon techniques to perform numerically well-defined (but complicated) tasks familiar to us from mathematics. · · · · · matrix manipulation (Linear Algebra) Fourier transforms minimisation algorithms special mathematical functions solution of equations (root finding etc.)

This will introduce us to the concept of "Numerical Analysis". There is considerable overlap between Computational Physics and Numerical Analysis: often we can use tried and tested algorithms to help solve our physical problem. Many of you will have had some exposure to numerical analysis in last year's IB maths course. For the professional Computational Physicist it is often essential to have a good understanding of such numerical techniques so as to be aware of their limitations. Going into this area in such depth is beyond the scope of the present course. Instead we will use algorithms in a pre-packed form, usually from the NAG library which we will have much to say about later. This is NOT a course on numerical analysis. 1.3 Our Approach As physicists we use a variety of computer-based/computational tools in addition to standard tools such as word processors etc.. Broadly speaking they fall into the following types: Computer algebra packages. [Maple, Mathcad, Mathematica] These have become vital tools and can be used for complex problems. They are not only the domain of the theorist, but are useful to all physicists. They are also suitable for small-scale numerical computations. Data analysis and computational packages. [Excel, IDL, Matlab] There are many packages available to help perform data analysis tasks, image processing and medium-scale numerical computations (useful, for example, when model-fitting data). Data visualisation. [Excel, gnuplot] There are a huge number of packages available to help visualise data; indeed data visualisation is fast becoming a specialist discipline of its own. These range from basic graph plotting to complex 3D representations requiring special hardware. Even virtual reality hardware/software is now being used. Programming languages [Fortran, C, C++] Despite all the available packages, as scientists we still often need to resort to programming in a traditional programming language. This is especially true for a computational physicist wishing to solve a complex physical problem such as we will be discussing later in the course. However, all physicists are likely to come across the need to write a simple program because the supplied application doesn't quite do what we want. This may not be in a traditional language, but may be in a language associated with a particular application (many "programs" are written now in Matlab, Maple, Mathematica etc.). Traditional programming languages are still very useful for the more computationally demanding problems. These are NOT alternatives, but complementary tools.

2


1.3.1 A Typical Project A typical project in computational physics will involve many steps: · · · · Use a computer algebra package to help when developing a physical model Develop a simulation or numerical solution to the physical model, say in a programming language (e.g. Fortran) Produce results Use a package to analyse and display the results (e.g. Excel)

In IA you learnt to use Excel and Mathcad; still useful tools. This course will teach Fortran 90/95 (unlike Fortran 77 this is a modern programming language) which is well-suited to computational physics. We will also introduce other tools (the text editor emacs and the graphics package gnuplot). 1.3.2 Why Unix? Note we shall be using a variant of Unix, Linux, for this course running on the PWF. But why Unix? · · Unix is a widely used Operating System (OS) (c.f. Windows Operating System) on large machines as well as PCs. Linux is a variant of Unix. Trained scientists should be familiar with Unix (also used in business on servers).

You are able to access the software and other course material, as well as run programs, on the Physics PWF, the Physics PWF Server, or perhaps your College or Personal PC (Linux is now a very common and widespread operating system). We shall assume that you have no experience whatsoever of programming languages, Unix/Linux and only a very limited experience of using a Windows PC. Almost all of you will in fact have more experience than we shall assume. 1.4 Structure of the Course The structure of this course will be somewhat different from most you have so far attended. The lectures and practical classes are very closely linked, and the practical sessions are extremely important. The practical part of the course will use Linux running on the Physics PWF. · · · Practical sessions will run each weekday in weeks 3 to 8, from 2.00 ­ 5.00pm. You should sign up for one practical session per week. The Physics PWF is located between the Part IA practical classes.

At each practical session there will be a head of class and a demonstrator to help you with any problems whatsoever you are having. To start with many of these may quite straightforward, for example how to log on! Do not be afraid to ask for any help you need: this way of working will be very new to almost all of you. You should go at your own pace in the practical sessions. To help you do this there are two self-study guides. One for PWF Linux and one for programming in Fortran 95. You choose how to use this material. The lectures will include mention of Linux and Fortran 95, BUT:

3


You will only learn Unix/Linux and programming by DOING. It is essential that you work through the self study guides both in the classes and your own time to familiarise yourself with the basics of Linux and Fortran 95.

As a guide, we expect that you will work through these self-study guides during the first four practical sessions. The course is assessed by asking you to complete a programming exercise (a small project). Again as a guide we envisage that during the second four practical sessions (plus thinking time outside the class) you will be working on the problem. During this period the head of class and demonstrator will still be on hand to answer your problems; make use of this help as much as you want to. Also, as in other practicals, you are encouraged to talk to other students and discuss problems. The final program that you produce and write up must, as always, be your own work, but you may make use of any help you have obtained in the practical sessions and any insight gained from discussions among yourselves. We do not expect you to become experts in the use of Unix or programming during this course and a good mark for the assessment can be achieved with a working program and a well constructed writeup. Lecture 1 will serve as an introduction to the course, provide a brief overview of Fortran, and will provide general advice on how to approach the course work. Lectures 2-7 will provide an introduction to computational physics, including: · · · · · Working with Data Ordinary Differential Equations Linear Algebra N-body simulations Abstract Simulation

Lecture 8 (time permitting) will be of a more general nature, and will discuss computational physics in the wider context.

4


1.5 Timetable Lectures: Practical sessions: Exercise: Tuesdays and Thursdays 9.00am (first eight lectures of Term) 14.00-17.00 starting Thursday 18th October 2007 Deadline for handing in exercise is 17.00 Thursday 17th January 2008

Further details are on the course web site: http://www.mrao.cam.ac.uk/~rachael/compphys 1.6 Reference Material The written material accompanying the course will be: Two self-study guides · PWF Linux · Programming in Fortran 95 Two handouts · An Introduction to Computational Physics (this one!) · The Physics of Computational Physics Other material · Problem sheet and administrative notes Recommended books on Fortran programming and Unix: · · · · · Fortran 90 For Scientists and Engineers, B.D. Hahn, Arnold Programming in Fortran 90, I.M. Smith, Wiley Fortran 90/95 explained, Metcalf and Reid, pub. Oxford Unix for Fortran Programmers, M. Loukides, pub. O'Reilly and Associates Unix in a Nutshell by Gilly et al., pub. O'Reilly and Associates

Recommended books on computational physics and numerical analysis: · · · · · Numerical methods for scientists and engineers, Hamming, McGraw Hill (A classic: probably out of print, but there may be copies in college libraries). Computational Physics, Thijssen, Cambridge Computational techniques in physics, MacKeown & Newman, Adam Hilger Computer simulation methods, Heerman, Springer-Verlag Numerical recipes in FORTRAN, 2nd Edition. Press et al., Cambridge (An advanced text useful for reference)

For on-line material an important starting place is the web site for this course.

5


2. A Brief History of Digital Computers
Charles Babbage is generally regarded as the father of modern computers. He was educated at Trinity, and became the Lucasian Professor of Mathematics. More importantly, in 1821 he invented the "Difference Engine", a mechanical computer for calculating mathematical tables. Although Babbage's version never progressed beyond a prototype, in 1854 George Scheutz constructed a version which saw successful and useful service.

Charles Babbage (1791 ­ 1871) The first programmer is also English: Ada Countess Lovelace. Ada, a daughter of Lord Byron, met Babbage in 1834, and wrote an article in 1843 which described how to use computers such as Babbage had invented for solving mathematical problems.

Ada Countess Lovelace (1815-52) Despite these early successes it was the advent of the electronic computer which really heralds the start of computing as we know it. This occurred at Bletchley Park in 1944 when Colossus was built. Again, many exCambridge people were involved (Turing, Newman, Tutte and others). Colossus used thermionic valves, not transistors, so was rather large. It was also not programmable: it was built for one specific task (decrypting a certain class of German codes). The Manchester Small Scale Experimental Machine (1948) was the first genuinely programmable, fully electronic digital computer. A year later Cambridge followed with EDSAC. Other early machines include ENIAC (US, 1945), an electronic an nearly-programmable computer, and Zuse Z3 (Germany, 1944), an electromechanical programmable computer.

EDSAC-1

6


EDSAC ran its first program in May 1949. It was run by the forerunner of the University's Computing Service, and in 1950 it started doing real scientific work. Early fields of study included X-ray crystallography and radio astronomy. In 1958 EDSAC was replaced by EDSAC-2, another Cambridge-built valve-based machine, which lasted until 1965. EDSAC-2's successor, Titan, was programmable in Fortran. The 1950s saw the introduction of transistors to computing, and the 1960s the introduction of integrated circuits. Companies such as IBM and DEC were very active at this time. The 1970s saw the introduction of the VLSI (very large scale integration) chips that all modern computers use, as well as the rise of a company called Intel.

3. Inside the Computer
It is helpful to have some general picture of how a computer works in order to appreciate some aspects of program development and design. One part of a computer, the CPU, is responsible for actually performing operations. Another, loosely called "memory", is responsible for storing the data being processed, any intermediate results, and the instructions which the CPU is executing. Memory has taken many forms over the years: paper tape, punched cards, magnetised ferrite beads, steel wire, magnetic tape, magnetic disks, and semiconductor-based memory. The last two are relevant for the computers used in this course! Standard semiconductor memory (RAM) costs about $150 for a gigabyte (8x109 bits), whereas disk drives are about 20 times cheaper. RAM is much faster than disk drives, but is volatile (its contents are lost when the computer is turned off). The term "memory" usually refers exclusively to RAM. The CPU is the brains of the computer and everything else is subordinate to it. A typical modern CPU understands two main classes of data: integer and floating point. Within those classes it may understand some additional subclasses, such as different precisions. It can perform basic arithmetic operations, governed by a sequence of instructions, or a program. It can also perform comparisons, the result of which can change the execution path through the program. To perform calculations the CPU must work with data: both input to it and returned results. The CPU itself holds very little data and the data must be transferred between the CPU and memory; this occurs via the "bus". In many calculations we find that it is not the speed at which the calculations themselves are done, but rather the speed with which this data transfer occurs which limits the overall performance of the calculation. Aside: there are many indications of how fast a particular computer is at performing calculations. Some of these include · · · MIPS: Millions of Instructions Per Second. Theoretical peak speed of basic operations. MTOPS: Millions of Theoretical Operations Per Second. Current favourite of the US Government. FLOPS: Floating Point Operations Per Second. Theoretical peak speed for floating point instructions. (M, G and T FLOPS: 106, 109, 1012 FLOPS). 7


3.1 Representing Numbers Computers store data as bits, each of which can represent either a 0 or 1; i.e. numbers are stored as binary. Clearly any integer can be represented exactly using a binary coding. Bits are processed together in groups of (usually) eight, called bytes. Most CPUs can process integers of different sizes, typically of 1, 2, 4 and possibly 8 bytes long. The larger the integer the more bits/bytes are needed to represent it. All CPUs will have a maximum number of bytes/bits they can process as an integer which means there is a maximum size to the integer which can be represented. The finite number of bytes used to represent an integer means that it is possible that the result of a calculation is larger than the integer that the computer can represent ­ this is called "overflow" and is something that must be avoided. Both positive and negative integers can be represented using one of a number of techniques (for example one bit, say the first bit can be used to represent the sign and remaining bits to represent the magnitude of the integer). To represent decimal, or "floating point" numbers the computer stores the leading digits (mantissa) together with an exponent and some indication of the sign of the number and the sign of the exponent. The maximum range of a floating point number is determined by the number of bits used to represent the exponent, while the precision will be determined mainly by the number of bits representing the significant digits. We may think of this in the following way: ± 0.XXXXXXX â 10
±YY

where we have assumed 6 digits for the mantissa and 2 for the exponent (of course the limits will actually be set by a binary representation). Overflow is still possible if the calculated number cannot be represented by the range set by the exponent. Similarly numbers of very small magnitude cannot be represented leading to a condition called "under flow". Even more important is the fact that the computer only stores a number to finite precision; this has important consequences which we have to carefully consider when solving problems computationally. The precision is, for ordinary floating point variables actually very low and for many calculations we have to use higher precision to store values than the default. 3.2 Representing Text Text can be represented in a similar way to integers. A pattern of bits can particular character. English only needs 7 bits to store any character (including all symbols on the keyboard). For accented characters in European languages still for textually rich languages such as Chinese. The most common system of of bits to characters is the ASCII mapping. be used to represent a digits, punctuation and 8 are needed and more assigning such patterns

8


4. The elements of programming
The CPU at the heart of a computer although able to perform calculations very efficiently is not in any way intelligent ­ it only carries out instructions given to it in the form of a program. A program is just a set of instructions to the computer to perform a series of operations. Those operations will often be mathematical calculations, decisions based on equalities and inequalities, or special instructions to say write output to the screen. The program consists of "source code" which is "stored" in a text file. This code contains the instructions in a highly structured form ­ a "Computer Language". Each computer language has a different set of rules (or syntax) for specifying these operations, although the logical operations we ask the computer to do are independent of the language. Here we will only consider the FORTRAN 90/95 (F95 for short) programming language and syntax. In order to appreciate what a program is and how to design a program we will consider a simply physical problem which we are going to solve using a computer. 4.1 The Problem We will simulate an ideal air table. Imagine we have a 2D rectangular region in which a frictionless circular object, a "puck" is able to move. The 2D region (the "table") is bounded on all sides by "walls". All collisions with the walls are elastic. Given some initial conditions we want to calculate the motion of the puck. Clearly this is very simple and we know the answer. However we will design a program to do the job in such a way as to allow for the possibility of adding additional pucks to the table which may collide with one another.

y x

9


4.2 Analysis of the problem We shall now consider how we might go about designing a program to compute the motion of the puck. Firstly, however, we need to consider exactly what the computer does (we have touched on this slightly already). Basically we can think of the computer as a very advanced calculator. Imagine how we might break down the problem if we were going to solve it without a computer, but instead (tediously) repeating calculations on a hand-held calculator. We need to consider not only the calculations we would perform, but all the decisions we would make along the way. We will then be in a position to convert this analysis into a "program" in which we will encode the calculations that the computer should make and rules for decisions etc. in a completely logical and completely prescribed fashion. The results will be a specific set of instructions to the computer for our particular problem. To fix our ideas we will decide to proceed by calculating the new position of the puck after a short, but finite time. Clearly we could thing of calculating the position after quite a long time interval since we know that the motion will be a perfect straight line until there is a collision with the wall, however our approach will enable us to extend the problem in the future. The starting point is then easy to see. If the initial position is r0 and the velocity is v then the new position after a time t can be calculated trivially r = r0 + v t So with our imaginary calculator we could do the sum and write down the new position, r. Clearly the next step must just be to repeat this calculation; a better notation than the equation we have above is perhaps to think of updating the position, r: r r + v t where the can read for the right arrow "becomes" or "is updated to". This simple calculation can be repeated again and again (each time recording the answer for future reference) until the puck reaches the wall. Let's consider what happens when the puck hits the wall in more detail. First we need to define where the walls are; we need four coordinates wy1, wy2, wx1, wx2. The puck will hit the wall, say the upper wall at wy2 when the y-component of the puck's position after the next update is greater than wy2 or within a puck radius of it. Defining the puck radius as R, and with r = (x,y), the condition can be written as y > (wy2 - R) similar conditions can be determined for collisions with the other walls.
y x

10


If the puck hits the wall during the next update of position as determined by the above criterion, then we must (a) account for the elastic collision by simply reversing the sign of the appropriate component of the puck's velocity and (b) determine the new position of the puck allowing for the reflection happening during the interval t. At this point we can decide to approximate the physical model. Reversing the sign of the velocity is easy, but calculating the new position correctly is a bit more tricky. A simpler approximation is to ignore the change in position of the puck during the time interval in which the collision occurs. We can ensure that the error is not too large by requiring that |v t| < R We now know what calculations and checks we need to do. The list of actions includes: · · · · · updating positions and velocities mathematical calculations recording positions and velocities repeating the calculations checking for certain conditions

In order to get the computer to do these calculations for us we need a way of specifying these same sorts of actions; it is precisely this form of description which a programming language gives us so that we can specify the actions the computer is to perform.

4.3 Development of the program In this section we will consider how each of the above actions is translated into F95. By the end of this section we will have covered the main elements of a program. More details and practice of using these ideas is in the self-study guides. 4.3.1 The overall program We firstly need a framework in which we can specify actions and commands that the computer is to perform. This is how we do it in F95: program table ! ! My first program ! write(*,*) 'Start of the Table Program' end program table

11


This is a complete F95 program. We have introduced here: · the overall structure of the program · lines beginning with a `!' are not instructions to the computer but comments for us so that we can document what we are doing · the write statement is a command telling the computer to output something to the screen, in this case the text `Start of the Table Program'. It is important to realise that F95 (unlike Unix) is not case sensitive so within the source code write and WRITE are identical for example. The symbolic description of what we want the computer to do, the "source code", must now be turned into an "executable" program. To do this we must "compile" the program. This will determine if there are any syntax errors in our program. Suppose the code is in a file called table.f90 then we compile the program as follows: f95 -o table table.f90 We run the program by typing: ./table 4.3.2 Variables and expressions We need some way of storing the current values of quantities such as position and velocity in our example and of performing mathematical operations on them. F95 provides this via variables and expressions which may include constants and variables. Variables may be of different types, the most important being integer, real (holding floating point variables), character and complex. In our example we shall only need real and integer. The next program performs a single step of our calculation assuming we do not encounter a wall. program table implicit none ! define variables used to hold position and velocity real :: x, y, vx, vy, dt ! set the initial values for the position and velocity x = 0.0; y = 0.0 vx = 1.5; vy = 2.0 ! and the time step dt = 0.1 x = x + vx * dt y = y + vy * dt write(*,*) x, y end program table

12


Our program now has two types of statements. The line near the top beginning real is called a declaration statement and is used to define the variables in our program. The remaining lines are called executable statements and are instructions to the computer to perform certain operations, in this case a few simple calculations. There are few important things to realise: · · · · statements end at the end of the line or at the comment character `!' two or more statements can be separated on a single line by semicolons `;' all variables must have a value before we can use them, our computer is just like the calculator of our problem it must work on data, we are not doing symbolic manipulations lines such as: x = x + vx * dt can be thought of us updating the value of a variable; the `=' is not an equality sign but rather an assignment operator the statement implicit none says that we require all variables we use to be declared

·

A full table of numeric operators can be found in section 1.3 of self-study guide 2. 4.3.3 Repetition: loops The next step in our program is to get the program to repeat the calculation many times. We do this using a loop. Loops are one of the key ideas in programming enabling us to repeat operations many times. The statement to implement a loop is called do and is terminated by end do. The statements between these two are repeated for a specified number of times. Here is a loop which will execute ten times: do i=1,10 write(*,*) i end do An integer variable, i, is used to keep count of how many times we have gone round the loop. It is assigned an initial value of 1 and the loop continues until i exceeds 10. The following would also loop for ten iterations, but i would run from 11 to 20: do i=11,20 write(*,*) i end do

13


Now let's modify our program: program table implicit none ! define variables used to hold position and velocity real :: x, y, vx, vy, dt integer :: i ! set the initial values for the position and velocity x = 0.0; y = 0.0 vx = 1.5; vy = 2.0 ! and the time step dt = 0.1 ! loop updating the position and outputting the value do i=1,30 x = x + vx * dt y = y + vy * dt write(*,*) x, y end do end program table Note · · ·

how all the declarations must come at the top of the program we put the output statement, i.e. the write inside the loop to keep track of the answer we have indented the lines between do and end do. This is good practice but not essential, however it provides a clear visual clue that these statements are within the "do loop" and will be repeated a number of times (30 in this case).

4.3.4 More about variables: arrays So far we have used an explicit component representation of the position and velocity in our programme. However, F95 provides us with the concept of an array variable which for our example will provide a very natural way of defining a vector. The syntax for defining an array variable is easy: real :: r(2), v(2) defines two array variables each with two elements of type real; in effect each defines two variables which we can refer to as r(1) and r(2) and similarly for v.

14


Let's see how the program looks using these array variables to represent out vectors: program table implicit none ! define variables used to hold position and velocity real :: r(2), v(2), dt integer :: i, j ! set the initial values for the position and velocity r(1) = 0.0; r(2) = 0.0 v(1) = 1.5; v(2) = 2.0 ! and the time step dt = 0.1 ! loop updating the position and outputting the value do i=1,30 do j=1,2 r(j) = r(j) + v(j) * dt write(*,*) r(1), r(2) end do end do end program table Note how we now use a second do loop to iterate over the elements of the vectors. In this example using arrays doesn't help a great deal but it is essential when we use arrays to hold large numbers of variables. Arrays may have more than one dimension. For example to define a matrix we use something like: real :: mat(3,3) 4.3.5 Conditional code: if The final missing element we need to implement our complete problem as a program is to do the check for a collision against one of the walls. F95 provides us with a way of conditionally executing some statements depending on the outcome of a numerical test, the if statement which in its most general form looks like: if (logical expression) then xxx else if (logical expression) then xxx else xxx end if

15


We can omit parts we do not need. For example to set a variable with a negative value to zero: if (x < 0) x = 0 ! which replaces a negative value in x by zero Returning to our example we introduce some variables to represent the positions of the walls and the radius of the puck real :: wx1, wx2, wy1, wy2, radius The tests that we need to perform can then be written as follows: if ((r(2)+radius) >= wy2) then v(2) = -v(2) ! reverse y-component at upper wall end if We have now used two comparison (or logical) operators, `<' for less than and `>=' for greater than or equal to. The full set of such operators is in section 1.6 of self-study guide 2. 4.3.6 Functions So far the calculations we have done have involved simple numerical expressions. We often want to perform more complicated calculations which use for example trigonometric functions, square roots etc. F95 provides us with the idea of a function which returns the result of such a calculation given some parameters. We can write our own functions, but F95 also provides some built-in or intrinsic functions. They are easy to use as we shall see. A list of the most common is in section 1.5 of selfstudy guide 2. Using them in a program is easy as this example of finding the sine of the variable x or the cosine of the angle 0.98 radians shows: real :: x, y x = 0.47 y = sin(x) y = cos(0.98) Importantly we can also define our own functions as you will learn later.

16


4.3.7 The final program We can now put everything together: program table implicit none ! define variables with fixed values real, parameter :: wx1=-1.5, wx2=1.5, wy1=-1.0, wy2=1.0 real, parameter :: radius=0.1 ! define variables which are updated in the program real :: r(2), rnew(2), v(2), dt integer :: i, imax, j, hit_wall ! set the initial values for the position and velocity r(1) = 0.0; r(2) = 0.0 v(1) = 1.5; v(2) = 2.0 ! calculate a suitable time step dt = 0.25 * radius / sqrt(v(1)*v(1) + v(2)*v(2)) ! a maximum number of steps to give a path of 10 units imax = 10 / (0.25 * radius) ! loop updating the position and outputting the value do i=1,imax do j=1,2 rnew(j) = r(j) + v(j) * dt end do ! check for wall collisions hit_wall = 0 if (((rnew(1)+radius) >= wx2) .or. & ((rnew(1)-radius) <= wx1)) then v(1) = -v(1); hit_wall = 1 end if if (((rnew(2)+radius) >= wy2) .or. & ((rnew(2)-radius) <= wy1)) then v(2) = -v(2); hit_wall = 1 end if ! output results and reset r if we have not hit a wall if (hit_wall == 0) then do j=1,2 r(j) = rnew(j) end do write(*,*) r(1), r(2) end if end do end program table 17


There are lots of things to note here: · · · We have used some new syntax on the declaration of some of the variables: real, parameter :: wx1=-1.5, wx2=1.5, wy1=-1.0, wy2=1.0 The parameter key word ensures that the variables wx1 etc. cannot have their values changed in the rest of the program; this stops us from making possible mistakes later on Some of the statements would have been very long and not easy to read (the ones associated with the logical tests); we have split these statements across multiple lines for ease of reading by ending the line with a `&' which means that the end of the line is not the end of the statement. We have used the intrinsic function sqrt to find a square root of a number or in this case an expression which is passed to it. Note how we now self-consistently find a value for the time step and then the number of iterations we want the program to run for; note how the limits of the do loop may themselves be variables. We use indentation to show the logical flow of the program. We have introduced a variable hit_wall which we use as an indicator to show if the results of one of our tests have found a collision with a wall; we only report results if there has been no collision.

· ·

· ·

More details on all of the topics covered in this section can be found in the self-study guides.

18