AUTHORED BY KONSTANTINOS N. ANAGNOSTOPOULOS
Physics Department, National Technical University of Athens, Zografou Campus, 15780 Zografou, Greece
konstant@mail.ntua.gr, www.physics.ntua.gr/~konstant/

PUBLISHED BY KONSTANTINOS N. ANAGNOSTOPOULOS
and the
NATIONAL TECHNICAL UNIVERSITY OF ATHENS

Book Website:
www.physics.ntua.gr/~konstant/ComputationalPhysics

©Konstantinos N. Anagnostopoulos 2014, 2016

First Published 2014
Second Edition 2016
Version¹ 2.0.20161206202800

Cover: Design by K.N. Anagnostopoulos. The front cover picture is a snapshot taken during Monte Carlo simulations of hexatic membranes. Work done with Mark J. Bowick. Relevant video at youtu.be/Erc7Q6YXfLk

This book and its cover(s) are subject to copyright. They are licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit
creativecommons.org/licenses/by-sa/4.0/

The book is accompanied by software available at the book’s website. All the software, unless the copyright does not belong to the author, is open source, covered by the GNU public license, see www.gnu.org/licenses/. This is explicitly mentioned at the end of the respective source files.

ISBN 978-1-365-58322-3 (lulu.com, vol. I)
ISBN 978-1-365-58338-4 (lulu.com, vol. II)

Contents

Foreword to the Second Edition
Foreword to the First Edition
1 The Computer
 1.1 The Operating System
  1.1.1 Filesystem
  1.1.2 Commands
  1.1.3 Looking for Help
 1.2 Text Processing Tools – Filters
 1.3 Programming with Emacs
  1.3.1 Calling Emacs
  1.3.2 Interacting with Emacs
  1.3.3 Basic Editing
  1.3.4 Cut and Paste
  1.3.5 Windows
  1.3.6 Files and Buffers
  1.3.7 Modes
  1.3.8 Emacs Help
  1.3.9 Emacs Customization
 1.4 The C++ Programming Language
  1.4.1 The Foundation
 1.5 Gnuplot
 1.6 Shell Scripting
2 Kinematics
 2.1 Motion on the Plane
  2.1.1 Plotting Data
  2.1.2 More Examples
 2.2 Motion in Space
 2.3 Trapped in a Box
  2.3.1 The One Dimensional Box
  2.3.2 Errors
  2.3.3 The Two Dimensional Box
 2.4 Applications
 2.5 Problems
3 Logistic Map
 3.1 Introduction
 3.2 Fixed Points and Cycles
 3.3 Bifurcation Diagrams
 3.4 The Newton-Raphson Method
 3.5 Calculation of the Bifurcation Points
 3.6 Liapunov Exponents
 3.7 Problems
4 Motion of a Particle
 4.1 Numerical Integration of Newton’s Equations
 4.2 Prelude: Euler Methods
 4.3 Runge–Kutta Methods
  4.3.1 A Program for the 4th Order Runge–Kutta
 4.4 Comparison of the Methods
 4.5 The Forced Damped Oscillator
 4.6 The Forced Damped Pendulum
 4.7 Appendix: On the Euler–Verlet Method
 4.8 Appendix: 2nd order Runge–Kutta Method
 4.9 Problems
5 Planar Motion
 5.1 Runge–Kutta for Planar Motion
 5.2 Projectile Motion
 5.3 Planetary Motion
 5.4 Scattering
  5.4.1 Rutherford Scattering
  5.4.2 More Scattering Potentials
 5.5 More Particles
 5.6 Problems
6 Motion in Space
 6.1 Adaptive Stepsize Control for RK Methods
  6.1.1 The rksuite Suite of RK Codes
  6.1.2 Interfacing C++ Programs with Fortran
  6.1.3 The rksuite Driver
 6.2 Motion of a Particle in an EM Field
 6.3 Relativistic Motion
 6.4 Problems
7 Electrostatics
 7.1 Electrostatic Field of Point Charges
 7.2 The Program – Appetizer and ... Desert
 7.3 The Program – Main Dish
 7.4 The Program - Conclusion
 7.5 Electrostatic Field in the Vacuum
 7.6 Results
 7.7 Poisson Equation
 7.8 Problems
8 Diffusion Equation
 8.1 Introduction
 8.2 Heat Conduction in a Thin Rod
 8.3 Discretization
 8.4 The Program
 8.5 Results
 8.6 Diffusion on the Circle
 8.7 Analysis
 8.8 Problems
9 The Anharmonic Oscillator
 9.1 Introduction
 9.2 Calculation of the Eigenvalues of
 9.3 Results
 9.4 The Double Well Potential
 9.5 Problems
10 Time Independent Schrödinger Equation
 10.1 Introduction
 10.2 The Infinite Potential Well
 10.3 Bound States
 10.4 Measurements
 10.5 The Anharmonic Oscillator - Again...
 10.6 The Lennard–Jones Potential
 10.7 Problems
11 The Random Walker
 11.1 (Pseudo)Random Numbers
 11.2 Using Pseudorandom Number Generators
 11.3 The MIXMAX Random Number Generator
 11.4 Random Walks
 11.5 Problems
12 Monte Carlo Simulations
 12.1 Statistical Physics
 12.2 Entropy
 12.3 Fluctuations
 12.4 Correlation Functions
 12.5 Sampling
  12.5.1 Simple Sampling
  12.5.2 Importance Sampling
 12.6 Markov Processes
 12.7 Detailed Balance Condition
 12.8 Problems
13 Simulation of the Ising Model
 13.1 The Ising Model
 13.2 Metropolis
 13.3 Implementation
  13.3.1 The Program
  13.3.2 Towards a Convenient User Interface
 13.4 Thermalization
 13.5 Autocorrelations
 13.6 Statistical Errors
  13.6.1 Errors of Independent Measurements
  13.6.2 Jackknife
  13.6.3 Bootstrap
 13.7 Appendix: Autocorrelation Function
 13.8 Appendix: Error Analysis
  13.8.1 The Jackknife Method
  13.8.2 The Bootstrap Method
  13.8.3 Comparing the Methods
 13.9 Problems
14 Critical Exponents
 14.1 Critical Slowing Down
 14.2 Wolff Cluster Algorithm
 14.3 Implementation
  14.3.1 The Program
 14.4 Production
 14.5 Data Analysis
 14.6 Autocorrelation Times
 14.7 Temperature Scaling
 14.8 Finite Size Scaling
 14.9 Calculation of
 14.10 Studying Scaling with Collapse
 14.11 Binder Cumulant
 14.12 Appendix: Scaling
  14.12.1 Binder Cumulant
  14.12.2 Scaling
  14.12.3 Finite Size Scaling
 14.13 Appendix: Critical Exponents
  14.13.1 Definitions
  14.13.2 Hyperscaling Relations
 14.14 Problems
Bibliography
Index

The book’s website is at
http://www.physics.ntua.gr/~konstant/ComputationalPhysics/
From there, you can can download the accompanying software, which contains, among other things, all the programs presented in the book.

Some conventions: Text using the font shown below refers to commands given using a shell (the “command line”), input and output of programs, code written in Fortran (or any other programming language), as well as to names of files and programs:

When a line starts with the prompt

then the text that follows is a command, which can be given from the command line of a terminal. The second line, Hello World, is the output of the command.

The contents of a file with C++ code is listed below:

What you need in order to work on your PC:

• An operating system of the GNU/Linux family and its basic tools.
• A Fortran compiler. The gfortran compiler is freely available for all major operating systems under an open source license at http://www.gfortran.org.
• An advanced text editor, suitable for editing code in several programming languages, like Emacs² .
• A good plotting program, suitable for data analysis, like gnuplot³ .
• The shell tcsh⁴ .
• The programs awk⁵ , grep, sort, cat, head, tail, less. Make sure that they are available in your computer environment.

If you have installed a GNU/Linux distribution on your computer, all of the above can be installed easily. For example, in a Debian like distribution (Ubuntu, ...) the commands

install all the necessary tools.

If you don’t wish to install GNU/Linux on your computer, you can try the following:

• Boot your computer using a usb/DVD live GNU/Linux, like Ubuntu⁶ . This will not make any permanent changes in your hard drive but it will start and run slower. On the other hand, you may save all your computing environment and documents and use it on any computer you like.
• Install Cygwin⁷ in your Microsoft Windows. It is a very good solution for Microsoft-addicted users. If you choose the full installation, then you will find all the tools needed in this book.
• Mac OS X is based on Unix. It is possible to install all the software needed in this book and follow the material as presented. Search the internet for instructions, e.g. google “gfortran for Mac”, “emacs for Mac”, “tcsh for Mac”, etc.

Foreword to the Second Edition

This book has been out “in the wild” for more than two years. Since then, its pdf version has been downloaded 2-5000 times/month from the main server and has a few thousand hits from sites that offer science e-books for free. I have also received positive feedback from students and colleagues from all over the world and that gave me the encouragement to devote some time to create a C++ version of the book. As far as scientific programming is concerned, the material has not changed apart from some typo and error corrections⁸ .

I have to make it clear that by using this book you will not learn much on the advanced features of C++. Scientific computing is usually simple at its core and, since it must be made efficient and accurate, it needs to go down to the lowest levels of programming. This also partly the reason of why I chose to use Fortran for the core programming in the first edition of the book: It is a language designed for numerical programming and high performance computing in mind. It is simple and a scientist or engineer can go directly into programming her code. C++ is not designed for scientific applications⁹ in mind and this reflects on some trivial omissions in its standard. Still, many scientific groups are now using C++ for programming and the C++ compilers have improved quite a lot. There is still an advantage in performance using a Fortran compiler on a supercomputer, but this is not going to last for much longer.

Still, for a scientist, the programming language is a tool to solve her scientific problems. One should not bind herself to a specific language. The treasures of today are the garbage of tomorrow, and the time scale for this happening is small in today’s computing environments. What has really lasting value is the ability to solve problems using a computer and this is what needs to be emphasized. Consistent with this idea is that, in the course of reading this book, you will also learn how to make your C++ code interact with code written in Fortran, like in the case of the popular library Lapack. This will improve your “multilingual skills” and flexibility with interacting with legacy code.

The good news for us scientists is that numerical code usually needs simple data structures and programming is similar in any language. It was simple for me to “translate” my book from Fortran to C++. Unfortunately I will not touch on all this great stuff in true object oriented programming but you may be happy to know that you will most likely not need it¹⁰ .

So, I hope that you will enjoy using my book and I remind you that I love fan mail and I appreciate comments/corrections/suggestions sent to me. Now, if you want to learn about the structure and educational procedure in this book, read the foreword to the first edition, otherwise skip to the real fun of solving scientific problems numerically.

Athens, 2016.

Foreword to the First Edition

This book is the culmination of my ten years’ experience in teaching three introductory, undergraduate level, scientific computing/computational physics classes at the National Technical University of Athens. It is suitable mostly for junior or senior level science courses, but I am currently teaching its first chapters to sophomores without a problem. A two semester course can easily cover all the material in the book, including lab sessions for practicing.

Why another book in computational physics? Well, when I started teaching those classes there was no bibliography available in Greek, so I was compelled to write lecture notes for my students. Soon, I realized that my students, majoring in physics or applied mathematics, were having a hard time with the technical details of programming and computing, rather than with the physics concepts. I had to take them slowly by the hand through the “howto” of computing, something that is reflected in the philosophy of this book. Hoping that this could be useful to a wider audience, I decided to translate these notes in English and put them in an order and structure that would turn them into “a book”.

I also decided to make the book freely available on the web. I was partly motivated by my anger caused by the increase of academic (e)book prices to ridiculous levels during times of plummeting publishing costs. Publishers play a diminishing role in academic publishing. They get an almost ready-made manuscript in electronic form by the author. They need to take no serious investment risk on an edition, thanks to print-on-demand capabilities. They have virtually zero cost ebook publishing. Moreover, online bookstores have decreased costs quite a lot. Academic books need no advertisement budget, their success is due to their academic reputation. I don’t see all of these reflected on reduced book prices, quite the contrary, I’m afraid.

My main motivation, however, is the freedom that independent publishing would give me in improving, expanding and changing the book in the future. It is great to have no length restrictions for the presentation of the material, as well as not having to report to a publisher. The reader/instructor that finds the book long, can read/print the portion of the book that she finds useful for her.

This is not a reference book. It uses some interesting, I hope, physics problems in order to introduce the student to the fundamentals of solving a scientific problem numerically. At the same time, it keeps an eye in the direction of advanced and high performance scientific computing. The reader should follow the instructions given in each chapter, since the book teaches by example. Several skills are taught through the solution of a particular problem. My lectures take place in a (large) computer lab, where the students are simultaneously doing what I am doing (and more). The program that I am editing and the commands that I am executing are shown on a large screen, displaying my computer monitor and actions live. The book provides no systematic teaching of a programming language or a particular tool. A very basic introduction is given in the first chapter and then the reader learns whatever is necessary for the solution of her problem. There is more than one way to do it¹¹ and the problems can be solved by following a basic or a fancy way, depending on the student’s computational literacy. The book provides the necessary tools for both. A bibliography is provided at the end of the book, so that the missing pieces of a puzzle can be sought in the literature.

This is also not a computational physics playground. Of course I hope that the reader will have fun doing what is in the book, but my goal is to provide an experience that will set the solid foundation for her becoming a high performance computing, number crunching, heavy duty data analysis expert in the future. This is why the programming language of the core numerical algorithms has been chosen to be Fortran, a highly optimized, scientifically oriented, programming language. The computer environment is set in a Unix family operating system, enriched by all the powerful GNU tools provided by the FSF¹² . These tools are indispensable in the complicated data manipulation needed in scientific research, which requires flexibility and imagination. Of course, Fortran is not the best choice for heavy duty object oriented programming, and is not optimal for interacting with the operating system. The philosophy¹³ is to let Fortran do what is best for, number crunching, and leave data manipulation and file administration to external, powerful tools. Tools, like awk, shell scripting, gnuplot, Perl and others, are quite powerful and complement all the weaknesses of Fortran mentioned before. The plotting program is chosen to be gnuplot, which provides very powerful tools to manipulate the data and create massive and complicated plots. It can also create publication quality plots and contribute to the “fun part” of the learning experience by creating animations, interactive 3d plots etc. All the tools used in the book are open source software and they are accessible to everyone for free. They can be used in a Linux environment, but they can also be installed and used in Microsoft Windows and Mac OS X.

The other hard part in teaching computational physics to scientists and engineers is to explain that the approach of solving a problem numerically is quite different from solving it analytically. Usually, students of this level are coming with a background in analysis and fundamental physics. It is hard to put them into the mode of thinking about solving a problem using only additions, multiplications and some logical operations. The hardest part is to explain the discretization of a model defined analytically, which can be done in many ways, depending on the accuracy of the approximation. Then, one has to extrapolate the numerical solution, in order to obtain a good approximation of the analytic one. This is done step by step in the book, starting with problems in simple motion and ending with discussing finite size scaling in statistical physics models in the vicinity of a continuous phase transition.

The book comes together with additional material which can be found at the web page of the book¹⁴ . The accompanying software contains all the computer programs presented in the book, together with useful tools and programs solving some of the exercises of each chapter. Each chapter has problems complementing the material covered in the text. The student needs to solve them in order to obtain hands on experience in scientific computing. I hope that I have already stressed enough that, in order for this book to be useful, it is not enough to be read in a café or in a living room, but one needs to do what it says.

Hoping that this book will be useful to you as a student or as an instructor, I would like to ask you to take some time to send me feedback for improving and/or correcting it. I would also appreciate fan mail or, if you are an expert, a review of the book. If you use the book in a class, as a main textbook or as supplementary material, I would also be thrilled to know about it. Send me email at konstantmail.ntua.gr and let me know if I can publish, anonymously or not, (part of) what you say on the web page (otherwise I will only use it privately for my personal ego-boost). Well, nothing is given for free: As one of my friends says, some people are payed in dollars and some others in ego-dollars!

Have fun computing scientifically!

Athens, 2014.

Chapter 1
The Computer

The aim of this chapter is to lay the grounds for the development of the computational skills which are necessary in the following chapters. It is not an in depth exposition but a practical training by example. For a more systematic study of the topics discussed, we refer to the bibliography. Many of the references are freely available n the web.

The are many choices that one has to make when designing a computer project. These depend on the needs for numerical efficiency, on available programming hours, on the needs for extensibility and upgradability and so on. In this book we will get the flavor of a project that is mostly scientifically and number crunching oriented. One has to make the best of the available computing resources and have powerful tools available for a productive analysis of the data. Such an environment, found in most of today’s supercomputers, that offers flexibility, dependability, simplicity, powerful tools for data analysis and effective compilers is provided by the family of the Unix operating systems. The GNU/Linux operating system is a Unix variant that is freely available and most of its utilities are open source software. The voluntary work of millions of excellent programmers worldwide has built the most stable, fastest and highest quality software available for scientific computing today. Thanks to the idea of the open source software pioneered by Richard Stallman¹ this giant collaboration has been made possible.

Another choice that we have to make is the programming language. In this edition of the book we will be programming in C++. C++ is a language with very high level of abstraction designed for projects where modular programming and the use of complicated data structures is of very high priority. A large and complicated project should be divided into independent programming tasks (modules), where each task contains everything that it needs and does not interfere with the functionality of other modules. Although it has not been designed for high performance numerical applications, it is becoming more and more popular in the recent years.

C++, as well as other languages like C, Java and Fortran, is a language that needs to be compiled by a compiler. Other languages, like python, perl, awk, shell scripting, Macsyma, Mathematica, Octave, Matlab, , are interpreted line by line. These languages can be simple in their use, but they can be prohibitively slow when it comes to a numerically demanding program. A compiler is a tool that analyzes the whole program and optimizes the computer instructions executed by the computer. But if programming time is more valuable, then a simple, interpreted language can lead to faster results.

Another choice that we make in this book, and we mention it because it is not the default in most Linux distributions, is the choice of shell. The shell is a program that “connects” the user to the operating system. In this book, we will teach how to use a shell² to “send” commands to the operating system, which is the most effective way to perform complicated tasks. We will use the shell tcsh, although most of the commands can be interpreted by most popular shells. Shell scripting is simpler in this shell, although shells like bash provide more powerful tools, mostly needed for complicated system administration tasks. That may cause a small inconvenience to some readers, since tcsh is not preinstalled in Linux distributions³ .

1.1 The Operating System

The Unix family of operating systems offer an environment where complicated tasks can be accomplished by combining many different tools, each of which performs a distinct task. This way, one can use the power of each tool, so that trivial but complicated parts of a calculation don’t have to be programmed. This makes the life of a researcher much easier and much more productive, since research requires from us to try many things before we understand how to compute what we are looking for.

In the Unix operating system everything is a file, and files are organized in a unique and unified filesystem. Documents, pictures, music, movies, executable programs are files. But also directories or devices, like hard disks, monitors, mice, sound cards etc, are, from the point of view of the operating system, files. In order for a music file to be played by your computer, the music data needs to be written to a device file, connected by the operating system to the sound card. The characters you type in a terminal are read from a file “the keyboard”, and written to a file “the monitor” in order to be displayed. Therefore, the first thing that we need to understand is the structure of the Unix filesystem.

1.1.1 Filesystem

There is at least one path in the filesystem associated with each file. There are two types of paths, relative paths and absolute paths. These are two examples:

The paths shown above may refer to the same or a different file. This depends on “where we are”. If “we are” in the directory /home/george, then both paths refer to the same file. If on the other way “we are” in a directory /home/john or /home/george/CompPhys, then the paths refer⁴ to two different files. In the last two cases, the paths refer to the files

respectively. How can we tell the difference? An absolute path always begins with the / character, whereas a relative path does not. When we say that “we are in a directory”, we refer to a position in the filesystem called the current directory, or working directory. Every process in the operating system has a unique current directory associated with it.

The filesystem is built on its root and looks like a tree positioned upside down. The symbol of the root is the character / The root is a directory. Every directory is a file that contains a list of files, and it is connected to a unique directory, its parent directory . Its list of files contains other directories, called its subdirectories, which all have it as their parent directory. All these files are the contents of the directory. Therefore, the filesystem is a tree of directories with the root directory at its top which branch to its subdirectories, which in their turn branch into other subdirectories and so on. There is practically no limit to how large this tree can become, even for a quite demanding environment⁵ .

A path consists of a string of characters, with the characters / separating its components, and refers to a unique location in the filesystem. Every component refers to a file. All, but the last one, must be directories in a hierarchy, from parent directory to subdirectory. The only exception is a possible / in the beginning, which refers to the root directory. Such an example can be seen in figure 1.1.

In a Unix filesystem there is complete freedom in the choice of the location of the files⁶ . Fortunately, there are some universally accepted conventions respected by almost everyone. One expects to find home directories in the directory /home, configuration files in the directory /etc, application executables in directories with names such as /bin, /usr/bin, /usr/local/bin, software libraries in directories with names such as /lib, /usr/lib etc.

There are some important conventions in the naming of the paths. A single dot “.” refers to the current directory and a double dot “..” to the parent directory. Similarly, a tilde “~” refers to the home directory of the user. Assume, e.g., that we are the user george running a process with a current directory /home/george/Music/Rock (see figure 1.1). Then, the following paths refer to the same file /home/george/Doc/lyrics.doc:

Notice that ~ and ~george refer to the home directory of the user george (ourselves), whereas ~mary refer to the home directory of another user, mary.

We are now going to introduce the basic commands for filesystem navigation and manipulation⁷ . The command cd (=change directory) changes the current directory, whereas the command pwd (=print working directory) prints the current directory:

The argument of the command cd is an absolute or a relative path. If the path is correct and we have the necessary permissions, the command changes the current directory to this path. If no path is given, then the current directory changes to the home directory of the user. If the character - is given instead of a path, then the command changes the current directory to the previous current directory.

The command mkdir creates new directories, whereas the command rmdir removes empty directories. Try:

Note that the command mkdir cannot create directories more than one level down the filesystem, whereas the command mkdir -p can. The “switch” -p makes the behavior of the command different than the default one.

In order to list the contents of a directory, we use the command ls (=list):

The first command is given without an argument and it lists the contents of the current directory. The second one, lists the contents of the subdirectory of the current directory Programs. If the argument is a list of paths pointing to regular files, then the command prints the names of the paths. Another way of giving the command is

The switch -l makes ls to list the contents of the current directory together with useful information on the files in 9 columns. The first column lists the permissions of the files (see below). The second one lists the number of links of the files⁸ . The third one lists the user who is the owner of each file. The fourth one lists the group that is assigned to the files. The fifth one lists the size of the file in bytes (=8 bits). The next three ones list the modification time of the file and the last one the paths of the files.

File permissions⁹ are separated in three classes: owner permissions, group permissions and other permissions. Each class is given three specific permissions, r=read, w=write and x=execute. For regular files, read permission effectively means access to the file for reading/copying, write permission means permission to modify the contents of the file and execute permission means permission to execute the file as a command¹⁰ . For directories, read permission means that one is able to read the names of the files in the directory (but not make it as current directory with the cd command), write permission means to be able to modify its contents (i.e. create, delete, and rename files) and execute permission grants permission to access/modify the contents of the files (but not list the names of the files, this is granted by the read permission).

The command ls -l lists permissions in three groups. The owner (positions 2-4), the group (positions 5-7) and the rest of the world (others - positions 8-10). For example

In the first case, the owner has read and write but not execute permissions and the group+others have only read permissions. In the second case, the user has read, write and execute permissions, the group has read permissions and others have no permissions at all. In the last case, the user has read, write and execute permissions, whereas the group and the world have only execute permissions. The first character d indicates a special file, which in this case is a directory. All special files have this position set to a character, while regular files have it set to -.

File permissions can be modified by using the command chmod:

Using the first command, the owner (u user) obtains (+) permission to execute (x) the file named file. Using the second one, the rest of the world (o others) and the group (ggroup) loose (-) the write (w) permission to the files named file1 and file2. Using the third one, everyone (aall) obtain read (r) permission on the file named file.

We will close this section by discussing some commands which are used for administering files in the filesystem. The command cp (copy) copies the contents of files into other files:

If the file file2.cpp does not exist, the first command copies the contents of file1.cpp to a new file file2.cpp. If it already exists, it replaces its contents by the contents of the file file2.cpp. In order for the second command to be executed, Programs needs to be a directory. Then, the contents of the files file1.cpp, file2.cpp, file3.cpp are copied to indentical files in the directory Programs. Of course, we assume that the user has the appropriate privileges for the command to be executed successfully.

The command mv “moves”, or renames, files:

The first command renames the file file1.cpp to file2.cpp. The second one moves files file1.cpp, file2.cpp, file3.cpp into the directory Programs.

The command rm (remove) deletes files¹¹ . Beware, the command is unforgiving: after deletion, a file cannot be restored into the filesystem¹² . Therefore, after executing successfully the following commands

the files file1.cpp, file2.cpp, file3.cpp do not exist in the filesystem anymore. A more prudent use of the command demands the flag -i. Then, before deletion we are asked for confirmation:

When we type y, the file is deleted, when we type n, the file is not deleted.

We cannot remove directories the same way. It is possible to use the command rmdir in order to remove empty directories. In order to delete directories together with their contents (including subdirectories and their contents) use the command¹³ rm -r. For example, assume that the contents of the directories dir1 and dir1/dir2 are the files:

Then the results of the following commands are:

The last command removes all files (assuming that we have write permissions for all directories and subdirectories). Alternatively, we can empty the contents of all directories first, and then remove them with the command rmdir:

Note that by using a semicolon, we can execute two or more commands on the same line.

1.1.2 Commands

Commands in a Unix operating system are files with execute permission. When we write a sentence on the command line, like

the shell reads its and interprets it. The shell is a program that creates a interface between a user and the operating system. The first word (ls) of the sentence is interpreted as a command. The rest of the words are the arguments of the command and the program can use them (or not) at the discretion of its programmer. There is a special convention for arguments that begin with a - (e.g. -l, --help, --version, -O3). They are called options or switches, and they act as virtual switches that make the program act in a particular way. We have already seen that the program ls gives a different output with the switch -l.

In order for a command to be executed, the shell looks for a file that has the same name as the command (here a file named ls). In order to understand where the shell looks for such a file, we should digress a little bit and explain the use of shell variables and environment variables. These have a name, which is a string of permissible characters, and their values are obtained by preceding their name with the $ character. For example the variable PATH has value $PATH. The values of the environment variables can be set with the command¹⁴ setenv and of the shell variables with the command set:

Two special variables are the variables PATH and path:

The first one is an environment variable and the second one is a shell variable. Their values are set by the shell, and we don’t need to worry about them, unless we want to change them. Their value is a string of characters whose components should be valid paths to directories. In the first case, the components are separated by a :, while in the second case, by one or more spaces. In the example shown above, the shell searches each component of the path or PATH variables (in this order) until it finds a file ls in their contents. If it succeeds and the file has execute permissions, then the program in this file is executed. If it fails, then it prints an error message. Try the commands:

We see that the program that the ls command executes the program in the file /bin/ls.

The arguments of a command are passed on to the program that the command executes for possible interpretation. For example:

The argument -l is the switch that results in a long listing of the files. The arguments test.cpp and test.dat are interpreted by the program ls as paths that it will look up for file information.

You can use the * (wildcard) character as a shorthand notation for a group of files. For example, in the command shown below

the shell will expand *.cpp and *.dat to a list of all files whose names end with .cpp or .dat. Therefore, if the current directory contains the files test.cpp, test1.cpp, myprog.cpp, test.dat, hello.dat, the arguments that will be passed on to the command ls are

For each command there are three special files associated with it. The first one is the standard input (stdin), the second one is the standard output (stdout) and the third one the standard error (stderr). These are files where the program can print or read data from. By default, these files are the terminal that the user uses to execute the command. In this case, when the program reads data from the stdin, then it reads the data that we type to the terminal using the keyboard. When the program writes data to the stdout or to the stderr, then the data is written to the terminal.

The advantage of using these special files in order to read/write data is that the user can redirect the input/output to these files to any file she wants. Using the character > at the end of a command redirects the stdout to the file whose name is written after >. For example:

The first of the above commands, prints the contents of the current working directory to the terminal. The second command redirects data written to the stdout to the file results. After executing the command, the file results is created and its contents are the names of the files file1.cpp file2.cpp file3.cpp file4.csh. If the file results does not exist (as in the above example), the file is created. If it already exists, it is truncated and its contents replaced by the data written to the stdout of the command. If we want to append data without erasing the existing contents, then we should use the string of characters >>. Therefore, if we give the command

after executing the previous commands, then the contents of the file results will be

The redirection of the stdin is accomplished by the use of the character < while that of the stderr by the use of the string of characters¹⁵ >&. We will see more examples in section 1.2.

It is possible to redirect the stdout of a command to be the stdin of another command. This is very useful for creating filters. A filter is a command that creates a flow of data between two or more programs. This process is called piping. Pipes are creating by using the character |

Using the syntax shown above, the stdout of the command cmd1 is redirected to the stdin of the command cmd2, the stdout of the command cmd2 is redirected to the stdin of the command cmd3 etc. More examples will be presented in section 1.2.

1.1.3 Looking for Help

Unix got itself a reputation for not being user friendly. This is far from the truth. Although there is a steep learning curve, detailed documentation for almost everything is available online.

The key for a comfortable ride is to learn how to use the help system available on your computer and on the internet. Most of the commands are self documented. A simple test, like the one shown below, will help you with the basic usage of most of the commands:

For example, try the command ls --help. For a window application, start from the menu “Help”. You should not be afraid and/or lazy and you should proceed with careful searching and reading.

For example, let’s assume that you have heard about a command that sounds like printf, or something like that. The first level of online help is the man (=manual) command that searches the “man pages”. Read the output of the command

The command info usually provides more detailed and user friendly documentation. It has basic browsing capabilities like the browsers you use to read pages from the internet. Try the command

Furthermore, the commands

will inform you that there are other, possibly related, commands with names like fprintf, fwprintf, wprintf, sprintf...:

The second column printed by the whatis command is the “section” of the man pages. In order to gain access to the information in a particular section, you have to give it as an argument to the man command:

Section 1 of the man pages contains information of ordinary command line commands, section 3 contains information on functions in libraries of the C language. Section 2 contains information on commands used for system administration. You may browse the directory /usr/share/man, or read the man page of the man command (use the command man man for that!).

By using the command

we obtain plenty of memory refreshing information. The command

shows us many files related to the command printf. The commands

give information on the location of the executable(s) of the command printf.

Another useful feature of the shell is the command or it filename completion. This means that we can write only the first characters of the name of a command or filename and then press simultaneously the keys [Ctrl-d]¹⁶ (i.e. press the key Ctrl and the key of the letter d at the same time). Then the shell will complete the name of the command up to the point that is is unique with the given string of characters¹⁷ :

Try to type an x on the command line and then type [Ctrl-d]. You will learn all the commands that are available and whose name begins with an x: xterm, xeyes, xclock, xcalc, ...

Finally, the internet contains a wealth of information. Google your blues... and you will be rewarded!

1.2 Text Processing Tools – Filters

For doing data analysis, we will need powerful tools for manipulating data in text files. These are files that consist solely of printable characters. Some tools that can be used in order to construct complicated and powerful filters are the programs cat, less, head, tail, grep, sort and awk.

Suppose that we have data in a file named data¹⁸ which contains information on the contents of a food warehouse and their prices:

The command

prints the contents of the file data to the stdout. In general, this command prints the contents of all files given in its arguments or the stdin if none is given. Since the stdin and the stdout can be redirected, the command

takes the contents of the file data from the stdin and prints them to the stdout, which in this case is the file data1. This command has the same result as the command:

The command

prints the contents of the file data and then the contents of the file data1 to the stdout. Since the stdout is redirected to the file data2, data2 contains the data of both files.

By giving the command

you can browse the data contained in the file gfortran.txt one page at a time. Press [space] in order to “turn” a page, [b] to turn back a page. Press the up and down arrows to move one line backwards/forward. Press [g] in order to jump to the beginning of the file and press [G] in order to jump to the end. Press [h] in order to get a help message and press [q] in order to quit.

The commands

print the first line, the last two lines and the second to the last line of the file data to the stdout respectively. Note that, by piping the stdout of the command tail to the stdin of the command head, we are able to construct the filter “print the line before the last one”.

The command sort sorts the contents of a file by comparing each line of its text with all others. The sorting is alphabetical, unless otherwise set by using options. For example

For reverse sorting, try sort -r data. We can also sort by comparing specific fields of each line. By default, fields are words separated by one or more spaces. For example, in order to sort w.r.t. the second column of the file data, we can use the switch -k 2 (=second field). Furthermore, we can use the switch -n for numerical sorting:

If we omit the switch -n, the comparison of the lines is performed based on character sorting of the second field and the result is

The last column contains floating point numbers (not integers). In order to sort by the values of such numbers we should use the switch -g:

The command grep processes a text file line by line, searching for a given string of characters. When this string is found anywhere in a line, this line is printed to the stdout. The command

prints each line containing the string “kilos”. If we want to search for all line not containing the string “kilos”, then we add the switch -v:

We can use a regular expression for searching a whole family of strings of characters. These monsters need a full book for discussing them in detail! But it is not hard to learn how to use some simple forms of regular expressions. Here are some examples:

The first one, prints each line whose first character is a b. The second one, prints each line that ends with a 0. The third one, prints each line contaning the strings 32 or 34.

By far, the strongest tool in our toolbox is the awk program. By default, awk analyzes a text file line by line. Each word (or field in the awk jargon) of these lines is stored in a set of variables with names $1, $2, .... The variable $0 contains the full line currently processed, whereas the variable NF counts the number of fields in the current line. The variable NR counts the number of lines of the file processed so far by awk.

An awk program can be written in the command line. A set of commands within { ... } is executed for each line of input. The constructs BEGIN{ ... } and END{ ... } contain commands executed, only once, before and after the processing of the file respectively. For example, the command

prints the name of the product (1st column = $1) and the total value stored in the warehouse (2nd column = $2) (4th column = $4). More examples are given below:

The first one calculates the total value of all products: The processing of each line results in the increment (+=) of the variable value by the product of the second and fourth fields. In the end (END{ ... }), the string Total= is printed, together with the final value of the variable value. This is an easy way for computing the sum of the values calculated for each line. The second command, calculates and prints an average. The sum is calculated in each line and stored in the variable av. In the end, we print the quotient of the sum of all values by the number of lines that have been processed (NR). The last command shows a (crazy) mathematical expression based on numerical values found in each line of the file data: It computes the square of the second field times the sine of the fourth field plus the exponential of the fourth field.

There is much more potential in the commands presented above. Reading the documentation and getting experience by using them will provide you with very strong tools in order to accomplish complicated tasks.

1.3 Programming with Emacs

For a programmer that spends many hours programming every day, the environment and the tools available for editing the commands of a large and complicated program determine, to a large extent, the quality of her life! An editor edits the contents of a text file, that consists solely of printable characters. Such editors, available in most Linux environments, are the programs gedit, vim, pico, nano, zile... They provide basic functionality such as adding, removing or changing text within a file as well as more complicated functions, such as copying, pasting, searching and replacing text etc. There are many functions that are particularly useful to a programmer, such as detecting and formatting keywords of a particular programming language, pretty printing, closing scopes etc, which can be very useful for comfortable programming and for spotting errors. A very powerful and “knowledgeable” editor, offering many such functions for several programming languages, is the GNU Emacs editor¹⁹ . Emacs is open source software, it is available for free and can be used in most available operating systems. It is programmable²⁰ and the user can automate most of her everyday repeated tasks and configure it to her liking. There is a full interaction with the operating system, in fact Emacs has been built with the ambition of becoming an operating system. For example, a programmer can edit a C++ file, compile it, debug it and run it, everything done with Emacs commands.

1.3.1 Calling Emacs

In the command line type

Note the character & at the end of the line. This makes the particular command to run in the background. Without it, the shell waits until a command exits in order to return the prompt.

In a desktop environment, Emacs starts in its own window. For a quick and dirty editing session, or in the case that a windows environment is not available²¹ , we can run Emacs in a terminal mode. Then, we omit the & at the end of the line and we run the command

The switch -nw forces Emacs to run in terminal mode.

1.3.2 Interacting with Emacs

We can interact with Emacs in various ways. Newbies will prefer buttons and menus that offer a simple and intuitive interface. For advanced usage, however, we recommend that you make an effort to learn the keyboard shortcuts. There are also thousands of functions available to be used interactively. They are called from a “command line”, called the minibuffer in the Emacs jargon.

Keyboard shortcuts are usually combinations of keystrokes that consist of the simultaneous pressing of the Ctrl or Alt keys together with other keys. Our convention is that a key sequence starting with a C- means that the characters that follow are keys simultaneously pressed with the Ctrl key. A key sequance starting with a M- means that the characters that follow are keys simultaneously pressed with the Alt key²² . Some commands have shortcuts consisting of two or more composite keystrokes. For example by C-x C-c we mean that we have to press simultaneously the Ctrl key together with x and then press simultaneously the Ctrl key together with c. This sequence is a shortcut to the command that exits Emacs. Another example is C-x 2 which means to press the Ctrl key together with x and then press only the key 2. This is a shortcut to the command that splits a window horizontally to two equal parts.

The most useful shortcuts are M-x (press the Alt key siumutaneously with the x key) and C-g. The first command takes us to the minibuffer where we can give a command by typing its name. For example, type M-x and then type save-buffers-kill-emacs in the minibuffer (this will terminate Emacs). The second one is an “SOS button” that interrupts anything Emacs does and returns control to the working buffer. This can be pretty handy when a command hangs or destroys our work and we need to interrupt it.

Figure 1.4: The basic menus found in Emacs when run in a desktop environment. We can see the basic commands and the keyboard shortcut reminders in the parentheses. E.g. the command File Visit New File can be given by typing C-x C-f. Note the commands File Visit New File (open a file), FileSave (write contents of a buffer to a file), FileExit Emacs, File Split Window (split window in two), FileNew Frame (open a new Emacs desktop window) and of course the well known commands Cut, Copy, Paste, Undo from the Edit menu. We can choose different buffers from the menu Buffers, which contain the contents of other files that we have opened for editing. We recommend trying the Emacs Tutorial and Read Emacs Manual in the Help menu.

The conventions for the mouse events are as follows: With Mouse-1, Mouse-2 and Mouse-3 we denote a simple click with the left, middle and right buttons of the mouse respectively. With Drag-Mouse-1 we mean to press the left button of the mouse and at the same time drag the mouse.

We summarize the possible ways of giving a command in Emacs with the following examples that have the same effect: Open a file and put its contents in a buffer for editing.

• By pressing the toolbar button that looks like a white sheet of paper (see figure 1.2).
• By choosing the FileVisit New File menu entry.
• By typing the keyboard shortcut C-x C-f.
• By typing the name of the command in the minibuffer: M-x find-file

The number of available commands increases from the top to the bottom of the above list.

1.3.3 Basic Editing

In order to edit a file, Emacs places the contents of a file in a buffer. Such a buffer is a chunk of computer memory where the contents of the file are copied and it is not the file itself. When we make changes to the contents of a buffer, the file remains intact. For our changes to take effect and be written to the file, we have to save the buffer. Then, the contents of the buffer are written back to the file. It is important to understand the following cycle of events:

• Read a file’s contents to a buffer.
• Edit buffer contents.
• Write (save) buffer’s contents back into the file.

Emacs may have more than one buffers open for editing simultaneously. By default, the name of the buffer is the same as the name of the file that is edited, although this is not necessary²³ . The name of a buffer is written in the modeline of the window of the buffer, as can be seen in figure 1.3.

If Emacs crashes or exits before we save our edits, it is possible to recover (part of) them. There is a command M-x recover-file that will guide us through the necessary recovery steps, or we can look for a file that has the same name as the buffer we were editing surrounded by two #. For example, if we were editing the file file.cpp, the automatically saved changes can be found in the file #file.cpp#. Auto saving is done periodically by Emacs and its frequency can be controlled by the user.

The point where we insert text while editing is called “the point”. This is right before the blinking cursor²⁴ . Each buffer has another position marked by “the mark”. A point and the mark define a “region” in the buffer. This is a part of the text in the buffer where the functions of Emacs can act (e.g. copy, cut, change case, spelling etc.). We can set the region by setting a point and then press C-SPC²⁵ or give the command M-x set-mark-command. This defines the current point to be the mark. Then we can move the cursor to another point which will define a region together with the mark that we set. Alternatively we can use Drag-Mouse-1 (hold the left mouse button and drag the mouse) and mark a region. The mark can be set with Mouse-3, i.e. with a simple click of the right button of the mouse. Therefore by Mouse-1 at a point and then Mouse-3 at a different point will set the region between the two points.

We can open a file in a buffer with the command C-x C-f, and then by typing its path. If the file already exists, its contents are copied to a buffer, otherwise a new buffer is created. Then:

• We can browse the buffer’s contents with the Up/Down/Left/Right arrows. Alternatively, by using the commands C-n, C-p, C-f and C-b.
• If the buffer is large, we can browse its contents one page at a time by using the Page Up/Page Dn keys. Alternatively, by using the commands C-v, M-v.
• Enter text at the points simply by typing it.
• Delete characters before the point by using the Backspace key and after the point by using the Delete key. The command C-d deletes a forward character.
• Erase all the characters in a line that lie ahead of the point by using the command C-k.
• Open a new line by using Enter or C-o.
• Go to the first character of a line by using Home and the last one by using End. Alternatively, by using the commands C-a and C-e, respectively.
• Go to the first character of the buffer with the key C-Home and the last one with the key C-End. Alternatively, with M-x beginning -of -buffer and M-x end-of-buffer.
• Jump to any line we want: Type M-x goto-line and then the line number.
• Search for text after the point: Press C-s and then the text you are looking for. This is an incremental search and the point jumps immediately to the first string that matches the search. The same search can be repeated by pressing C-s repeatedely.

When we finish editing (or frequently enough so that we don’t loose our work due to an unfortunate event), we save the changes in the buffer, either by pressing the save icon on the toolbar, or by pressing the keys C-s, or by giving the command M-x save-buffer.

1.3.4 Cut and Paste

Use the instructions below for slightly more advanced editing:

• Undo! Some of the changes described below can be catastrophic. Emacs has a great Undo function that keeps in its memory many of the changes inflicted by our editing commands. By repeatedely pressing C-/, we undo the changes we made. Alternatively, we can use C-x u or the menu entry EditUndo. Remember that C-g interrupts any Emacs process currently running in the buffer.
• Cut text by using the mouse: Click with Mouse-1 at the point before the beginning of the text and then Mouse-3 at the point after the end. A second Mouse-3 and the region is ... gone (in fact it is written in the “kill ring” and it is available for pasting)!
• Cut text by using a keyboard shortcut: Set the mark by C-SPC at the point before the beginning of the text that you want to cut. Then move the cursor after the last character of the text that you want to cut and type C-w.
• Copy text by using the mouse: Drag the mouse Drag-Mouse-1 and mark the region that you want to copy. Alternatively, Mouse-1 at the point before the beginning of the text and then Mouse-3 at the point after the end.
• Copy text by using a keyboard shortcut: Set the mark at the beginning of the text with C-SPC and then move the cursor after the last character of the text. Then type M-w.
• Pasting text with the mouse: We click the middle button²⁶ Mouse-2 at the point that we want to insert the text from the kill ring (the copied text).
• Pasting text with a keyboard shortcut: We move the point to the desired insertion point and type C-y.
• Pasting text from previous copying: A fast choice is the menu entry EditPaste from kill manu and then select from the copied texts. The keyboard shortcut is to first type C-y and then M-y repeatedly, until the text that we want is yanked.
• Insert the contents of a file: Move the point to the desired place and type C-x i and the path of the file. Alternatively, give the command M-x insert-file.
• Insert the contents of a buffer: We can insert the contents of a whole buffer at a point by giving the command M-x insert-buffer.
• Replace text: We can replace text interactively with the command M-x query-replace, then type the string we want to replace, and then the replacement string. Then, we will be asked whether we want the change to be made and we can answer by typing y (yes), n (no), q (quit the replacements). A , (comma) makes only one replacement and quits (useful if we know that this is the last change that we want to make). If we are confident, we can change all string in a buffer, no questions asked, by giving the command M-x replace-string.
• Change case: We can change the case in the words of a region with the commands M-x upcase-region, M-x capitalize-region and M-x downcase-region. Try it.

We note that cutting and pasting can be made between different windows of the same or different buffers.

1.3.5 Windows

Sometimes it is very convenient to edit one or more different buffers in two or more windows. The term “windows” in Emacs refers to regions of the same Emacs desktop window. In fact, a desktop window running an Emacs session is referred to as a frame in the Emacs jargon. Emacs can split a frame in two or more windows, horizontally or/and vertically. Study figure 1.5 on page 81 for details. We can also open a new frame and edit several buffers simultaneously²⁷ . We can manipulate windows and frames as follows:

• Position the point at the center of the window and clear the screen from garbage: C-l (careful: l not 1).
• Split a window in two, horizontally: C-x 2.
• Split a window in two, vertically: C-x 3.
• Delete all other windows (remain only with the current one): C-x 1.
• Delete the current windows (the others remain): C-x 0.
• Move the cursor to the other window: Mouse-1 or C-x o.
• Change the size of window: Use Drag-Mouse-1 on the line separating two windows (the mode line). Use C-^, C-} for making a change of the horizontal/vertical size of a window respectively.
• Create a new frame: C-x 5 2.
• Delete a frame: C-x 5 0.
• Move the cursor to a different frame: With Mouse-1 or with C-x 5 o.

You can have many windows in a dumb terminal. This is a blessing when a dekstop environment is not available. Of course, in that case you cannot have many frames.

1.3.6 Files and Buffers

• Open a file: C-x C-f or M-x find-file.
• Save a buffer: C-x C-s or M-x save buffer. With C-x C-c or M-x save-buffers-kill-emacs we can also exit Emacs. From the menu: FileSave. From the toolbar: click on the save icon.
• Save buffer contents to a different file: C-x C-w or M-x write-file. From the menu: FileSave As. From the toolbar: click on the “save as” icon.
• Save all buffers: C-x s or M-x save-some-buffers.
• Connect a buffer to a different file: M-x set-visited-filename.
• Kill a buffer: C-x k.
• Change the buffer of the current window: C-x b. Also, use the menu Buffers, then choose the name of the buffer.
• Show the list of all buffers: C-x C-b. From the menu: Buffers List All Buffers. By typing Enter next to the name of the buffer, we make it appear in the window. There are several buffer administration commands. Learn about them by typing C-h m when the cursor is in the Bufer List window.
• Recover data from an edited buffer: If Emacs crashed, do not despair. Start a new Emacs and type M-x recover-file and follow the instructions. The command M-x recover-session recovers all unsaved buffers.
• Backup files: When you save a buffer, the previous contents of the file become a backup file. This is a file whose path is the same as the original’s file with a ~ appended in the end. For example a file test.cpp will have as a backup the file test.cpp~. Emacs has version control, and you can configure it to keep as many versions of your edits as you want.
• Directory browsing and directory administration commands: C-x d or M-x dired. You can act on the files of a directory (open, delete, rename, copy etc) by giving appropriate commands. When the cursor is in the dired window, type C-h m to read the relevant documentation.

1.3.7 Modes

Each buffer can be in different modes. Each mode may activate different commands or editing environment. For example each mode can color keywords relevant to the mode and/or bind keys to different commands. There exist major modes, and a buffer can be in only one of them. There are also minor modes, and a buffer can be in one or more of them. Emacs activates major and minor modes by default for each file. This usually depends on the filename but there are also other ways to control this. The user can change both major and minor modes at will, using appropriate commands.

Active modes are shown in a parenthesis on the mode line (see figures 1.3 and 1.5.

• M-x c++-mode: This mode is of special interest in this book since we will edit a lot of C++ code. We need it activated in buffers that contain a C++ program and its most useful characteristics are automatic code alignment by pressing the key TAB, the coloring of C++ statements, variables and other structural constructs (classes, if statements, for loops, variable declarations, comments etc). Another interesting function is the one that comments out a whole region of code, as well as the inverse function.
• M-x c-mode: For files containing programs written in the C language. Related modes are the java-mode, perl-mode, awk-mode, python-mode, makefile-mode, octave-mode, gnuplot-mode and others.
• latex-mode: For files containing LATEX  text formatting commands.
• text-mode: For editing simple text files (.txt).
• fundamental-mode: The basic mode, when one that fits better doesn’t exist...

Some interesting minor modes are:

• M-x auto-fill-mode: When a line becomes too long, it is wrapped automatically. A related command to do that for the whole region is M-x fill-region, and for a paragraph M-x fill-paragraph.
• M-x overwite-mode: Instead of inserting characters at the point, overwrite the existing ones. By giving the command several times, we toggle between activating and deactivating the mode.
• M-x read-only mode: When visiting a file with valuable data that we don’t want to change by mistake, we can activate this mode so that changes will not be allowed by Emacs. When we open a file with the command C-x C-r or M-x find-file-read-only this mode is activated. We can toggle this mode on and off with the command C-x C-q (M-x toggle-read-only). See the mode line of the buffer jack.c in figure 1.5 which contains a string %%. By clicking on the %% we can toggle the read-only mode on and off.
• flyspell-mode: Spell checking as we type.
• font-lock-mode: Colors the structural elements of the buffer which are defined by the major mode (e.g. the commands of a C++ program).

In a desktop environment, we can choose modes from the menu of the mode line. By clicking with Mouse-3 on the name of a mode we are offered options for (de)activating minor modes. With a Mouse-1 we can (de)activate the read-only mode with a click on :%% or :-- respectively. See figure 1.5.

1.3.8 Emacs Help

Emacs’ documentation is impressive. For newbies, we recommend to follow the mini course offered by the Emacs tutorial. You can start the tutorial by typing C-h t or select Help Emacs Tutorial from the menu. Enjoy... The Emacs man page (give the man emacs command in the command line) will give you a summary of the basic options when calling Emacs from the command line.

A quite detailed manual can be found in the Emacs info pages²⁸ . Using info needs some training, but using the Emacs interface is quite intuitive and similar to using a web browser. Type the command C-h r (or choose HelpEmacs Tutorial from the menu) and you will open the front page of the emacs manual in a new window. By using the keys SPC and Backspace we can read the documentation page by page. When you find a link (similar to web page hyperlinks), you can click on it in order to open to read the topic it refers to. Using the navigation icons on the toolbar, you can go to the previous or to the next pages, go up one level etc. There are commands that can be given by typing single characters. For example, type d in order to jump to the main info directory. There you can find all the available manuals in the info system installed on your computer. Type g (emacs) and go to the top page of the Emacs manual. Type g (info) and read the info manual.

Emacs is structured in an intuitive and user friendly way. You will learn a lot from the names of the commands: Almost all names of Emacs commands consist of whole words, separated by a hyphen “-”, which almost form a full sentence. These make them quite long sometimes, but by using auto completion of their names this does not pose a grave problem.

• auto completion: The names of the commands are auto completed by typing a TAB one or more times. E.g., type M-x in order to go to the minibuffer. Type capi[TAB] and the command autocompletes to capitalize-. By typing [TAB] for a second time, a new window opens and offers the options for completing to two possible commands: capitalize-region and capitalize-word. Type an extra r[TAB] and the command auto completes to the only possible choice capitalize-region. You can see all the commands that start with an s by typing M-x s[TAB][TAB]. Sure, there are many... Click on the *Completions* buffer and browse the possibilities. A lot will become clear just by reading the names of the commands. By typing M-x [TAB][TAB], all available commands will appear in your buffer!
• keyboard shortcuts: If you don’t remember what happens when you type C-s, no problem: Type C-h k and then the ... forgotten key sequence C-s. Conversely, have you forgotten what is the keyboard shortcut of the command save-buffer? Type C-h w and then the command.
• functions: Are you looking for a command, e.g. save-something -I-forgot? Type C-h f and then save-[TAB] in order to browse over different choices. Use Mouse-2 in order to select the command you are interested in, or type and complete the rest of its name (you may use [TAB] again). Read about the function in the *Help* buffer that opens.
• variables: Do the same after typing C-h v in order to see a variable’s value and documentation.
• command apropos: Have you forgotten the exact name of a command? No problem... Type C-h a and a keyword. All commands related to the keyword you typed will appear in a buffer. Use C-h d for even more information.
• modes: When in a buffer, type C-h m and read information about the active modes of the buffer.
• info: Type C-h i
• Have you forgotten everything mentioned above? Just type C-h ?

1.3.9 Emacs Customization

You can customize everything in Emacs. From key bindings to programming your own functions in the Elisp language. The most common way for a user to customize her Emacs sessions, is to put all her customization commands in the file /.emacs in her home directory. Emacs reads and executes all these commands just before starting a session. Such a .emacs file is given below:

Everything after a ; is a comment. Functions/commands are enclosed in parentheses. The first three ones bind the keys F1, C-c s and M-s to the commands save-buffer, save-some-buffers and isearch-forward respectively. The next one defines an alias of a command. This means that, when we give the command M-x is in the minibuffer, then the command isearch-forward will be executed. The last two commands are the definitions of the functions (fm) and (sign), which can be called interactively from the minibuffer.

For more complicated examples google “emacs .emacs file” and you will see other users’ .emacs files. You may also customize Emacs from the menu commands OptionsCustomize Emacs. For learning the Elisp language, you can read the manual “Emacs Lisp Reference Manual” found at the address
www.gnu.org/software/emacs/manual/elisp.html

1.4 The C++ Programming Language

In this section, we give a very basic introduction to the C++ programming language. This is not a systematic exposition and you are expected to learn what is needed in this book by example. So, please, if you have not already done it, get in front of a computer and do what you read. You can find many good tutorials and books introducing C++ in a more complete way in the bibliography.

1.4.1 The Foundation

The first program that one writes when learning a new programming language is the “Hello World!” program. This is the program that prints “Hello World!” on your screen:

Commands, or statements, in C++ are strings of characters separated by blanks (“words”) and end with a semicolon (;). We can put more than one command on each line by separating them with a semicolon.

Everything after two slashes (//) is a comment. Proliferation of comments is necessary for documenting our code. Good documentation of our code is an integral part of programming. If we plan to have our code read by others (or by us) at a later time, we have to make sure to explain in detail what each line is supposed to do. You and your collaborators will save a lot of time in the process of debugging, improving and extending your code.

The first line of the code shown above is a preprocessor directive. These lines start with a # and are interpreted by a separate program, the preprocessor. The #include directive, inserts the contents of a file replacing the line where the directive is. This acts like an editor! Actually, the code that will be compiled is not the one shown above, but the result of adding the contents of a file whose name is iostream²⁹ . iostream is an example of a header file that has many definitions of functions and symbols used by the program. The particular header has the necessary definitions in order to perform standard input and standard output operations.

The execution of a C++ program starts by calling a function whose name is main(). Therefore, the line int main(){ shows how to actually define a function in C++. Its name is the word before the parentheses () and the keyword int specifies that the function returns a value of integer type³⁰ . Within the parentheses placed after the name of the function, we put the arguments that we pass to the function. In our case the parentheses contain nothing, showing how to define a function without arguments.

The curly brackets { ... } define the scope or the body of the function and contain the statements to be executed when the function is called.

The line

is the only line that contains an executable statement that actually does something. Notice that it ends with a semicolon. This statement performs an output operation printing a string to the standard output. The sentence Hello World!∖n is a constant string and contains a sequence of printable characters enclosed in double quotes. The last character ∖n is a newline character, that prints a new line to the stdout.

cout identifies the standard character output device, which gives access to the stdout. The characters << indicate that we write to cout the expression to the right. In order to make cout accessible to our program, we need both the inclusion of the header file iostream and the statement using namespace std³¹ .

Statements in C++ end with a semicolon. Splitting them over a number of lines is only a matter of making the code legible. Therefore, the following parts of the code have equivalent effect as the one written above:

Finally notice that, for C++, uppercase and lowercase characters are different. Therefore main(), Main() and MAIN() are names of different functions.

In order to execute the commands in a program, it is necessary to compile it. This is a job done by a program called the compiler that translates the human language programming statements into binary commands that can be loaded to the computer memory for execution. There are many C++ compilers available, and you should learn which compilers are available for use in your computing environment. Typical names for C++ compilers are g++, c++, icc, .... You should find out which compiler is best suited for your program and spend time reading its documentation carefully. It is important to learn how to use a new compiler so that you can finely tune it to optimize the performance of your program.

We are going to use the open source and freely available compiler g++, which can be installed on most popular operating systems³² . The compilation command is:

The extension .cpp to the file name hello.cpp is important and instructs the compiler that the file contains source code in C++. Use your editor and edit a file with the name hello.cpp with the program shown above before executing the above command.

The switch -o defines the name of the executable file, which in our case is hello. If the compilation is successful, the program runs with the command:

The ./ is not a special symbol for running programs. The dot is the current working directory and ./hello is the full path to the file hello.

Now, we will try a simple calculation. Given the radius of a circle we will compute its length and area. The program can be found in the file area_01.cpp:

The first two statements in main() declare the values of the variables PI and R. These variables are of type double, which are floating point numbers³³ .

The following two commands have two effects: Computing the length and the area of the circle and printing the results. The expressions 2.0*PI*R and PI*R*R are evaluated before being printed to the stdout. Note the explicit decimal points at the constants 2.0 and 4.0. If we write 2 or 4 instead, then these are going to be constants of the int type and by using them the wrong way we may obtain surprising results³⁴ . We compile and run the program with the commands:

Now we will try a process that repeats itself for many times. We will calculate the length and area of 10 circles of different radii , . We will store the values of the radii in an array R[10] of the double type. The code can be found in the file area_02.cpp:

The declaration double R[10] defines an array of length 10. This way, the elements of the array are referred to by an index that takes values from 0 to 9. For example, R[0] is the first, R[3] is the fourth and R[9] is the last element of the array.

Between the lines

we can write commands that are repeatedly executed while the int variable i takes values from 1 to 9 with increasing step equal to 1. The way it works is the following: In the round brackets after the keyword for, there exist three statements separated by semicolons. The first, i=1, is the statement executed once before the loop starts. The second, i<10, is a statement that is evaluated each time before the loop repeats itself. If it is true, then the statements in the loop are executed. If it is false, the control of the program is transferred after the end of the loop. The last statement, i++, is evaluated each time after the last statement in the loop has been executed. The operator ++ is the increment operator, and its effect is equivalent to the statement:

The value of i is increased by one. The command:

defines the i-th radius from the value R[i-1]. For the loop to work correctly, we must define the initial value of R[0], otherwise³⁵ R[0]=0.0. The second loop uses the defined R-values in order to do the computation and print of the results.

Now we will write an interactive version of the program. Instead of hard coding the values of the radii, we will interact with the user asking her to give her own values. The program will read the 10 values of the radii from the standard input (stdin). We will also see how to write the results directly to a file instead of the standard output (stdout). The program can be found in the file area_03.cpp:

In the above program, the size of the array R is defined by a const int. A const declares a variable to be a parameter whose value does not change during the execution of the program and, if it is of int type, it can be used to declare the size of an array.

The array elements R[i] are read using the command:

cin is the standard input stream, the same way that cout is the standard output stream³⁶ . We read input using the >> operator, which indicates that input is written to the variable on the right.

In order to interact with ordinary files, we need to include the header

In this header, the C++ class ofstream is defined and it can be used in order to write to files (output stream). An object in this class, like myfile, is defined (“instantiated”) by the statement:

This object’s constructor is called by placing the parentheses ("AREA.DAT"), and then the output stream myfile is directed to the file AREA.DAT. Then we can write output to the file the same way we have already done with cout:

When we are done writing to the file, we can close the stream with the statement:

Reading from files is done in a similar way by using the class ifstream instead of ofstream.

The next step will be to learn how to define and use functions. The program below shows how to define a function area_of_circle(), which computes the length and area of a circle of given radius. The following program can be found in the file area_04.cpp:

The calculation of the length and the area of the circle is performed by the function

Calling a function, transfers the control of the program to the statements within the body of the function. The above function has arguments (R[i], perimeter, area). The argument R[i] is intended to be only an input variable whose value is not going to change during the calculation. The arguments perimeter and area are intended for output. Upon return of the function to the main program, they store the result of the computation. The user of a function must learn how to use its arguments in order to be able to call it in her program. These must be documented carefully by the programmer of the function.

In order to use a function, we need to declare it the same way we do with variables or, as we say, to provide its prototype. The prototype of a function can be declared without providing the function’s definition. We may provide just enough details that determine the types of its arguments and the value returned. In our program this is done on the line:

This is the same syntax used later in the definition of the function, but replacing the body of the function with a semicolon. The argument list does not need to include the argument names, only their types. We could have also used the following line in order to declare the function’s prototype:

We could also have used different names for the arguments, if we wished so. Including the names is a matter of style that improves legibility of the code.

The argument R is intended to be left unchanged during the function execution. This is why we used the keyword const in its declaration. The arguments L and R, however, will return a different value to the calling program. This is why const is not used for them.

The actual program executed by the function is between the lines:

The type of the value returned by a function is declared by the keyword before its name. In our case, this is void which declares that the function does not return a value.

The arguments (R,L,A) must be declared in the function and need not have the same names as the ones that we use when we call it. All arguments are declared to be of type double. The character & indicates that they are passed to the function by reference. This makes possible to change their values from within the function.

If & is omitted, then the arguments will be passed by value and a statement like L = 2.0*PI*R will not change the value of the variable passed by the calling program. This happens because, in this case, only the value of the variable L of the calling program is copied to a local variable which is used only within the function. This is important to understand and you are encouraged to run the program with and without the & and check the difference in the computed results.

The names of variables in a function are only valid within the scope of the function, i.e. between the curly brackets that contain the body of the function. Therefore the variable const int N is valid only within the scope of main(). You may use any name you like, even if it is already used outside the scope of the function. The names of arguments need not be the same as the ones used in the calling program. Only their types have to match.

Variables in the global scope are accessible by all functions in the same file³⁷ . An example of such a variable is PI, which is accessible by main(), as well as by area_of_circle().

We summarize all of the above in a program trionymo.cpp, which computes the roots of a second degree polynomial:

The program reads the coefficients of the polynomial . After a check whether , it computes the discriminant by calling the Discriminant(a,b,c).

The type of the value returned must be declared at the function’s prototype

and at the function’s definition

The value returned to the calling program is the value of the expression given as an argument to the return statement. return has also the effect of transferring the control of the program back to the calling statement.

1.5 Gnuplot

Plotting data is an indispensable tool for their qualitative, but also quantitative, analysis. Gnuplot is a high quality, open source, plotting program that can be used for generating publication quality plots, as well as for heavy duty analysis of a large amount of scientific data. Its great advantage is the possibility to use it from the command line, as well as from shell scripts and other programs. Gnuplot is programmable and it is possible to call external programs in order manipulate data and create complicated plots. There are many mathematical functions built in gnuplot and a fit command for non linear fitting of data. There exist interactive terminals where the user can transform a plot by using the mouse and keyboard commands.

This section is brief and only the features, necessary for the following chapters, are discussed. For more information visit the official page of gnuplot http://gnuplot.info. Try the rich demo gallery at http://gnuplot.info/screenshots/, where you can find the type of graph that you want to create and obtain an easy to use recipe for it. The book  [16] is an excellent place to look for many of gnuplot’s secrets³⁸ .

You can start a gnuplot session with the gnuplot command:

There is a welcome message and then a prompt gnuplot> is issued waiting for your command. Type a command an press [Enter]. Type quit in order to quit the program. In the following, when we show a prompt gnuplot>, it is assumed that the command after the prompt is executed from within gnuplot.

Plotting a function is extremely easy. Use the command plot and x as the independent variable of the function³⁹ . The command

plots the function which is a straight line with slope 1. In order to plot many functions simultaneously, you can write all of them in one line:

The above command plots the functions , , and . Within the square brackets [:], we set the limits of the and axes, respectively. The bracket [-5:5] sets and the bracket [-2:4] sets . You may leave the job of setting such limits to gnuplot, by omitting some, or all of them, from the respective positions in the brackets. For example, typing [1:][:5] changes the lower and upper limits of and and leaves the upper and lower limits unchanged⁴⁰ .

In order to plot data points , we can read their values from files. Assume that a file data has the following numbers recorded in it:

The first line is taken by gnuplot as a comment line, since it begins with a #. In fact, gnuplot ignores everything after a #. In order to plot the second column as a function of the first, type the command:

The name of the file is within double quotes. After the keyword using, we instruct gnuplot which columns to use as the and coordinates, respectively. The keywords with points instructs gnuplot to add each pair to the plot with points.

The command

plots the third column as a function of the first, and the keywords with lines instruct gnuplot to connect each pair with a straight line segment.

We can combine several plots together in one plot:

The first line plots the 1st and 3rd columns in the file data together with the function . The second line adds the plot of the 1st and 2nd columns in the file data and the third line adds the plot of the function .

There are many powerful ways to use the keyword using. Instead of column numbers, we can put mathematical expressions enclosed inside brackets, like using (...):(...). Gnuplot evaluates each expression within the brackets and plots the result. In these expressions, the values of each column in the file data are represented as in the awk language. $i are variables that expand to the number read from columns i=1,2,3,.... Here are some examples:

The first line plots the 1st column of the file data together with the value , where , and are the numbers in the 2nd, 1st and 3rd columns respectively. The second line adds the plot of the function .

The first line plots the logarithm of the 1st column together with the logarithm of the square of the 2nd column.

We can plot the data written to the standard output of any command. Assume that there is a program called area that prints the perimeter and area of a circle to the stdout in the form shown below:

The interesting data is at the second and fourth columns. These can be plotted directly with the gnuplot command:

All we have to do is to type the full command after the < within the double quotes. We can create complicated filters using pipes as in the following example:

The filter produces data to the stdout, by combining the action of the commands area, sort and awk. The data printed by the last program is in two columns and we plot the results using 1:2.

In order to save plots in files, we have to change the terminal that gnuplot outputs the plots. Gnuplot can produce plots in several languages (e.g. PDF, postscript, SVG, LATEX, jpeg, png, gif, etc), which can be interpreted and rendered by external programs. By redirecting the output to a file, we can save the plot to the hard disk. For example:

The first line makes the plot as usual. The second one sets the output to be in the JPEG format and the third one sets the name of the file to which the plot will be saved. The fourth lines repeats all the previous plotting commands and the fifth one closes the file data.jpg. The last line chooses the interactive terminal qt to be the output of the next plot. High quality images are usually saved in the PDF, encapsulated postcript or SVG format. Use set terminal pdf,postscript eps or svg, respectively.

And now a few words for 3-dimensional (3d) plotting. The next example uses the command splot in order to make a 3d plot of the function . After you make the plot, you can use the mouse in order to rotate it and view it from a different perspective:

If you have data in the form and you want to create a plot of , write the data in a file, like in the following example:

Note the empty line that follows the change of the value of the first column. If the name of the file is data3, then you can plot the data with the commands:

We close this section with a few words on parametric plots. A parametric plot on the plane (2-dimensions) is a curve , where is a parameter. A parametric plot in space (3-dimensions) is a surface , where are parameters. The following commands plot the circle and the sphere :

1.6 Shell Scripting

A typical GNU/Linux environment offers very powerful tools for complicated system administration tasks. They are much more simple to use than to incorporate them into your program. This way, the programmer can concentrate on the high performance and scientific computing part of the project and leave the administration and trivial data analysis tasks to other, external, programs.

One can avoid repeating the same sequence of commands by coding them in a file. An example can be found in the file script01.csh:

This is a very simple shell script. The first line instructs the operating system that the lines that follow are to be interpreted by the program /bin/tcsh⁴¹ . This can be any program in the system, which in our case is the tcsh shell. The following lines are valid commands for the shell, one in each line. They compile the C++ programs found in the files that we created in section 1.4 with g++, and then they run the executable ./area. In order to execute the commands in the file, we have to make sure that the file has the appropriate execute permissions. If not, we have to give the command:

Then we simply type the path to the file script01.csh

and the above commands are run the one after the other. Some of the versions of the programs that we wrote are asking for input from the stdin, which, normally, you have to type on the terminal. Instead of interacting directly with the program, we can write the input data to a file Input, and run the command

A more convenient solution is to use the, so called, “Here Document”. A “Here Document” is a section of the script that is treated as if it were a separate file. As such, it can be used as input to programs by sending its “contents” to the stdin of the command that runs the program⁴² . The “Here Document” does not appear in the filesystem and we don’t need to administer it as a regular file. An example of using a “Here Document” can be found in the file script02.csh:

The stdin of the command ./area is redirected to the contents between the lines

The string EOF marks the beginning and the end of the “Here Document”, and can be any string you like. The last EOF has to be placed exactly in the beginning of the line.

The power of shell scripting lies in its programming capabilities: Variables, arrays, loops and conditionals can be used in order to create a complicated program. Shell variables can be used as discussed in section 1.1.2: The value of a variable name is $name and it can be set with the command set name = value. An array is defined, for example, by the command

and its data can be accessed using the syntax $R[1] ... $R[10].

Lets take a look at the following script:

The first two lines of the script define the values of the arrays files (4 values) and R (10 values). The command echo echoes its argument to the stdin. $USER is the name of the user running the script. ‘date‘ is an example of command substitution: When a command is enclosed between backquotes and is part of a string, then the command is executed and its stdout is pasted back to the string. In the example shown above, ‘date‘ is replaced by the current date and time in the format produced by the date command.

The foreach loop

is executed once for each of the 4 values of the array files. Each time the value of the variable file is set equal to one of the values area_01.cpp, area_02.cpp, area_03.cpp, area_04.cpp. These values can be used by the commands in the loop. Therefore, the command g++ $file -o area compiles a different file each time that it is executed by the loop.

The last line in the loop

is a conditional. It executes the command cat AREA.DAT if the condition -f AREA.DAT is true. In this case, -f constructs a logical expression which is true when the file AREA.DAT exists.

We close this section by presenting a more complicated and advanced script. It only serves as a demonstration of the shell scripting capabilities. For more information, the reader is referred to the bibliography  [18, 19, 20, 21, 22]. Read carefully the commands, as well as the comments which follow the # mark. Then, write the commands to a file script04.csh⁴³ , make it an executable file with the command chmod u+x script04.csh and give the command

The script will run with the words “This is my first serious tcsh script” as its arguments. Notice how these arguments are manipulated by the script. Then, the script asks for the values of the radii of ten or more circles interactively, so that it will compute their perimeter and area. Type them on the terminal and then observe the script’s output, so that you understand the function of each command. You will not regret the time investment!

Chapter 2
Kinematics

2.1 Motion on the Plane

When a particle moves on the plane, its position can be given in Cartesian coordinates . These, as a function of time, describe the particle’s trajectory. The position vector is , where and are the unit vectors on the and axes respectively. The velocity vector is where

Figure 2.1: The trajectory of a particle moving in the plane. The figure shows its position vector , velocity and acceleration and their Cartesian components in the chosen coordinate system at a point of the trajectory.

In this section we study the kinematics of a particle trajectory, therefore we assume that the functions are known. By taking their derivatives, we can compute the velocity and the acceleration of the particle in motion. We will write simple programs that compute the values of these functions in a time interval , where is the initial and is the final time. The continuous functions are approximated by a discrete sequence of their values at the times such that .

We will start the design of our program by forming a generic template to be used in all of the problems of interest. Then we can study each problem of particle motion by programming only the equations of motion without worrying about the less important tasks, like input/output, user interface etc. Figure 2.2 shows a flowchart of the basic steps in the algorithm. The first part of the program declares variables and defines the values of the fixed parameters (like , , etc). The program starts by interacting with the user (“user interface”) and asks for the values of the variables , , , , . The program prints these values to the stdout so that the user can check them for correctness and store them in her data.

The main calculation is performed in a loop executed while . The values of the positions and the velocities are calculated and printed in a file together with the time . At this point we fix the format of the program output, something that is very important to do it in a consistent and convenient way for easing data analysis. We choose to print the values t, x, y, vx, vy in five columns in each line of the output file.

The specific problem that we are going to solve is the computation of the trajectory of the circular motion of a particle on a circle with center and radius with constant angular velocity . The position on the circle can be defined by the angle , as can be seen in figure 2.3. We define the initial position of the particle at time to be .

The equations giving the position of the particle at time are

The data structure is quite simple. The constant angular velocity is stored in the double variable omega. The center of the circle , the radius of the circle and the angle are stored in the double variables x0, y0, R, theta. The times at which we calculate the particle’s position and velocity are defined by the parameters and are stored in the double variables t0, tf, dt. The current position is calculated and stored in the double variables x, y and the velocity in the double variables vx, vy. The declarations of the variables are put in the beginning of the program:

The user interface of the program is the interaction of the program with the user and, in our case, it is the part of the program where the user enters the parameters omega, x0, y0, R, t0, tf, dt. The program issues a prompt with the names the variables expected to be read. The variables are read from the stdin by reading from the stream cin and the values entered by the user are printed to the stdout using the stream cout¹ :

There are a couple of things to explain. Notice that after reading each variable from the standard input stream cin, we call the function getline. By calling getline(cin,buf), a whole line is read from the input stream cin into the string buf² . Then the statement

has the effect of reading three doubles from the stdin and put the rest of the line in the string buf. Since we never use buf, this is a mechanism to discard the rest of the line of input. The reason for doing so will become clear later.

Objects of type string in C++ store character sequences. In order to use them you have to include the header

and, e.g., declare them like

Then you can store data in the obvious way, like buf="Hello World!", manipulate string data using operators like buf=buf1 (assign buf1 to buf), buf=buf1+buf2 (concatenate buf1 and buf2 and store the result in buf), buf1==buf2 (compare strings) etc.

Finally, endl is used to end all the cout statements. This has the effect of adding a newline to the output stream and flush the output³ .

Next, the program initializes the state of the computation. This includes checking the validity of the parameters entered by the user, so that the computation will be possible. For example, the program computes the expression 2.0*PI/omega, where it is assumed that omega has a non zero value. We will also demand that and . An if statement will make those checks and if the parameters have illegal values, the exit statement⁴ will stop the program execution and print an informative message to the standard error stream cerr⁵ . The program opens the file Circle.dat for writing the calculated values of the position and the velocity of the particle.

The line myfile.precision(17) sets the precision of the floating point numbers (like double) printed to myfile to 17 significant digits accuracy. The default is 6 which is a pity, because doubles have up to 17 significant digits accuracy.

If or the corresponding exit statements are executed which end the program execution. The optional error messages are included after the stop statements which are printed to the stderr. The value of the period is also calculated and printed for reference.

The main calculation is performed within the loop

The first statement sets the initial value of the time. The statements between within the scope of the while(condition) are executed as long as condition has a true value. The statement t=t+dt increments the time and this is necessary in order not to enter into an infinite loop. he statements put in place of the dots ......... calculate the position and the velocity and print them to the file Circle.dat:

Notice the use of the functions sin and cos that calculate the sine and cosine of an angle expressed in radians. The header cmath is necessary to be included.

The program is stored in the file Circle.cpp and can be found in the accompanied software. The extension .cpp is used to inform the compiler that the file contains source code written in the C++ language. Compilation and running can be done using the commands:

The switch -o cl forces the compiler g++ to write the binary commands executed by the program to the file⁶ cl. The command ./cl loads the program instructions to the computer memory for execution. When the programs starts execution, it first asks for the parameter data and then performs the calculation. A typical session looks like:

The lines shown above that start with a # character are printed by the program and the lines without # are the values of the parameters entered interactively by the user. The user types in the parameters and then presses the Enter key in order for the program to read them. Here we have used , , , , and .

You can execute the above program many times for different values of the parameters by writing the parameter values in a file using an editor. For example, in the file Circle.in type the following data:

Each line has the parameters that we want to pass to the program with each call to cout. The rest of the line consists of comments that explain to the user what each number is there for. We want to discard these characters during input and this is the reason for using getline to complete reading the rest of the line. The program can read the above values of the parameters with the command:

The command ./cl runs the commands found in the executable file ./cl. The < Circle.in redirects the contents of the file Circle.in to the standard input (stdin) of the command ./cl. This way the program reads in the values of the parameters from the contents of the file Circle.in. The > Circle.out redirects the standard output (stdout) of the command ./cl to the file Circle.out. Its contents can be inspected after the execution of the program with the command cat:

We list the full program in Circle.cpp below:

2.1.1 Plotting Data

We use gnuplot for plotting the data produced by our programs. The file Circle.dat has the time t and the components x, y, vx, vy in five columns. Therefore we can plot the functions and by using the gnuplot commands:

Figure 2.4: The plots (left) and (right) from the data in Circle.dat for , , , , and .

The second line puts the second plot together with the first one. The results can be seen in figure 2.4.

Let’s see now how we can make the plot of the function . We can do that using the raw data from the file Circle.dat within gnuplot, without having to write a new program. Note that . The function atan2 is available in gnuplot⁷ as well as in C++. Use the online help system in gnuplot in order to see its usage:

Therefore, the right way to call the function is atan2(y-y0,x-x0). In our case x0=y0=1 and x, y are in the 2nd and 3rd columns of the file Circle.dat. We can construct an expression after the using command as in page 129, where $2 is the value of the second and $3 the value of the third column:

The second command is broken in two lines by using the character ∖ so that it fits conveniently in the text⁸ . Note how we defined the values of the variables x0, y0 and how we used them in the expression atan2($3-x0,$2-y0). We also plot the lines which graph the constant functions and which mark the limit values of . The gnuplot variable⁹ pi is predefined and can be used in formed expressions. The result can be seen in the left plot of figure 2.4.

The velocity components as function of time as well as the trajectory can be plotted with the commands:

We close this section by showing how to do a simple animation of the particle trajectory using gnuplot. There is a file animate2D.gnu in the accompanied software which you can copy in the directory where you have the data file Circle.dat. We are not going to explain how it works¹⁰ but how to use it in order to make your own animations. The final result is shown in figure 2.5. All that you need to do is to define the data file¹¹ , the initial time t0, the final time tf and the time step dt. These times can be different from the ones we used to create the data in Circle.dat. A full animation session can be launched using the commands:

The first line defines the data file that animate2D.gnu reads data from. The second line sets the range of the plots and the third line defines the time parameters used in the animation. The final line launches the animation. If you want to rerun the animation, you can repeat the last two commands as many times as you want using the same or different parameters. E.g. if you wish to run the animation at “half the speed” you should simply redefine dt=0.05 and set the initial time to t0=0:

2.1.2 More Examples

We are now going to apply the steps described in the previous section to other examples of motion on the plane. The first problem that we are going to discuss is that of the small oscillations of a simple pendulum. Figure 2.6 shows the single oscillating degree of freedom , which is the small angle that the pendulum forms with the vertical direction.

Figure 2.6: The simple pendulum whose motion for is described by the program SimplePendulum.cpp.

The motion is periodic with angular frequency and period . The angular velocity is computed from which gives

We note that the acceleration of gravity is hard coded in the program and that the user can only set the length of the pendulum. The data file SimplePendulum.dat produced by the program, contains two extra columns with the current values of and the angular velocity .

A simple session for the study of the above problem is shown below¹² :

The next example is the study of the trajectory of a particle shot near the earth’s surface¹³ when we consider the effect of air resistance to be negligible. Then, the equations describing the trajectory of the particle and its velocity are given by the parametric equations

Figure 2.7: The trajectory of a particle moving under the influence of a constant gravitational field. The initial conditions are set to , and .

The structure of the program is similar to the previous ones. The user enters the magnitude of the particle’s initial velocity and the shooting angle in degrees. The initial time is taken to be . The program calculates and and prints them to the stdout. The data is written to the file Projectile.dat. The full program is listed below and it can be found in the file Projectile.cpp in the accompanied software:

A typical session for the study of this problem is shown below:

Next, we will study the effect of air resistance of the form . The solutions to the equations of motion

Figure 2.8: The forces that act on the particle of figure 2.7 when we assume air resistance of the form .

Programming the above equations is as easy as before, the only difference being that the user needs to provide the value of the constant . The full program can be found in the file ProjectileAirResistance.cpp and it is listed below:

Figure 2.9: The plots of , (left) and , (right) from the data produced by the program ProjectileAirResistance.cpp for , , , and . We also plot the asymptotes of these functions as .

Figure 2.10: Trajectories of the particles shot with , in the absence of air resistance and when the air resistance is present in the form with .

We also list the commands of a typical session of the study of the problem:

Long commands have been continued to the next line as before. We defined the gnuplot variables v0x, v0y, g and k to have the values that we used when running the program. We can use them in order to construct the asymptotes of the plotted functions of time. The results are shown in figures 2.9 and 2.10.

The last example of this section will be that of the anisotropic harmonic oscillator. The force on the particle is

(2.11)

where the “spring constants” and are different in the directions of the axes and . The solutions of the dynamical equations of motion for , , and are

We have set in the program above. The user must enter the two angular frequencies and and the corresponding times. A typical session for the study of the problem is shown below:

The results for and are shown in figure 2.11.

Figure 2.11: The trajectory of the anisotropic oscillator with and .

2.2 Motion in Space

By slightly generalizing the methods described in the previous section, we will study the motion of a particle in three dimensional space. All we have to do is to add an extra equation for the coordinate and the component of the velocity . The structure of the programs will be exactly the same as before.

The first example is the conical pendulum, which can be seen in figure 2.12. The particle moves on the plane with constant angular velocity . The equations of motion are derived from the relations

(2.13)

where . Their solution¹⁵ is

(2.17)

and when , .

In the program that we will write, the user must enter the parameters , , the final time and the time step . We take . The convention that we follow for the output of the results is that they should be written in a file where the first 7 columns are the values of , , , , , and . The full program is listed below:

In order to compile and run the program we can use the commands shown below:

The results are recorded in the file ConicalPendulum.dat. In order to plot the functions , , , , , we give the following gnuplot commands:

The results are shown in figure 2.13.

Figure 2.13: The plots of the functions of the program ConicalPendulum.cpp for , .

In order to make a three dimensional plot of the trajectory, we should use the gnuplot command splot:

Figure 2.14: The plot of the particle trajectory of the program ConicalPendulum.cpp for , . We can click and drag with the mouse on the window and rotate the curve and see it from a different angle. At the bottom left of the window, we see the viewing direction, given by the angles degrees (angle with the axis) and degrees (angle with the axis).

The result is shown in figure 2.14. We can click on the trajectory and rotate it and view it from a different angle. We can change the plot limits with the command:

We can animate the trajectory of the particle by using the file animate3D.gnu from the accompanying software. The commands are similar to the ones we had to give in the two dimensional case for the planar trajectories when we used the file animate2D.gnu:

The result can be seen in figure 2.15.

Figure 2.15: The particle trajectory computed by the program ConicalPendulum.cpp for , and plotted by the gnuplot script animate3D.gnu. The title of the plot shows the current time and the particles coordinates.

The program animate3D.gnu can be used on the data file of any program that prints t x y z as the first words on each of its lines. All we have to do is to change the value of the file variable in gnuplot.

Next, we will study the trajectory of a charged particle in a homogeneous magnetic field . At time , the particle is at and its velocity is , see figure 2.16.

Figure 2.16: A particle at time is at the position with velocity in a homogeneous magnetic field .

The magnetic force on the particle is and the equations of motion are

We are now ready to write a program that calculates the trajectory given by (2.20) . The user enters the parameters and , shown in figure 2.16, as well as the angular frequency (Larmor frequency). The components of the initial velocity are and . The initial position is calculated from the equation . The program can be found in the file ChargeInB.cpp:

A typical session in which we calculate the trajectories shown in figures 2.17 and 2.18 is shown below:

Figure 2.17: The plots of the functions calculated by the program in ChargeInB.cpp for , , degrees.

Figure 2.18: The trajectory calculated by the program in ChargeInB.cpp for , , degrees as shown by the program animate3D.gnu. The current time and the coordinates of the particle are printed on the title of the plot.

2.3 Trapped in a Box

In this section we will study the motion of a particle that is free, except when bouncing elastically on a wall or on certain obstacles. This motion is calculated by approximate algorithms that introduce systematic errors. These types of errors¹⁶ are also encountered in the study of more complicated dynamics, but the simplicity of the problem will allow us to control them in a systematic and easy to understand way.

2.3.1 The One Dimensional Box

The simplest example of such a motion is that of a particle in a “one dimensional box”. The particle moves freely on the axis for , as can be seen in figure 2.19. When it reaches the boundaries and it bounces and its velocity instantly reversed. Its potential energy is

(2.21)

which has the shape of an infinitely deep well. The force within the box and at the position of the walls.

Figure 2.19: A particle in a one dimensional box with its walls located at and .

Initially we have to know the position of the particle as well as its velocity (the sign of depends on the direction of the particle’s motion) at time . As long as the particle moves within the box, its motion is free and

where the last line gives the testing condition for the wall collision and the subsequent change of the velocity.

The full program that realizes the proposed algorithm is listed below and can be found in the file box1D_1.cpp. The user can set the size of the box L, the initial conditions x0 and v0 at time t0, the final time tf and the time step dt:

In this section we will study the effects of roundoff errors in numerical computations. Computers store numbers in memory, which is finite. Therefore, real numbers are represented in some approximation that depends on the amount of memory that is used for their storage. This approximation corresponds to what is termed as floating point numbers. C++ is supposed to provide at least three basic types of floating point numbers, float, double and long double. In most implementations¹⁷ , float uses 4 bytes of memory and double 8. In this case, float has an accuracy to, approximately, 7 significant digits and double 17. See Chapter 1 of  [8] and  [14] for details. Moreover, float represent numbers with magnitude in the, approximate, range while double in . Note that variables of the integer type (int, long, ...) are exact representations of integers, whereas floating point numbers are approximations to reals.

In the program shown above, we used numbers of the float type instead of double in order to exaggerate roundoff errors. This way we can study the dependence of this type of errors on the accuracy of the floating point numbers used in a program¹⁸ . In order to do that, we declared the floating point variables as float:

We also used numerical constants of type float. This is indicated by the letter f at the end of their names: 2.0 is a constant of type double (the C++ default), whereas 2.0f is a constant of type float. Determining the accuracy of floating point constants is a thorny issue that can be the cause on introducing subtle bugs in a program and the programmer should be very careful about doing it carefully.

Finally we changed the form of the output. Since a float represents a real number with at most 7 significant digits, there is no point of printing more. That is why we used the statements

For purposes of studying the numerical accuracy of our results, we used 9 digits of output, which is, of course, slightly redundant. setw(17) prints the numbers of the next output of the stream myfile using at least 17 character spaces. This improves the legibility of the results when inspecting the output files. The use of setw requires the header iomanip.

The computed data is recorded in the file box1D_1.dat in three columns. Compiling, running and plotting the trajectory using gnuplot can be done as follows:

Figure 2.20: The trajectory of a particle in a box with , , , , . The plot to the right magnifies a detail when which exposes the systematic errors in determining the exact moment of the collision of the particle with the wall at and the corresponding maximum value of , .

The trajectory is shown in figure 2.20. The effects of the systematic errors can be easily seen by noting that the expected collisions occur every units of time. Therefore, on the plot to the right of figure 2.20, the reversal of the particle’s motion should have occurred at , .

The reader should have already realized that the above mentioned error can be made to vanish by taking arbitrarily small . Therefore, we naively expect that as long as we have the necessary computer power to take as small as possible and the corresponding time intervals as many as possible, we can achieve any precision that we want. Well, that is true only up to a point. The problem is that the next position is determined by the addition operation x+v*dt and the next moment in time by t+dt. Floating point numbers of the float type have a maximum accuracy of approximately 7 significant decimal digits. Therefore, if the operands x and v*dt are real numbers differing by more than 7 orders of magnitude (v*dt x), the effect of the addition x+v*dt=x, which is null! The reason is that the floating point unit of the processor has to convert both numbers x and v*dt into a representation having the same exponent and in doing so, the corresponding significant digits of the smaller number v*dt are lost. The result is less catastrophic when v*dt x with , but some degree of accuracy is also lost at each addition operation. And since we have accumulation of such errors over many intervals tt+dt, the error can become significant and destroy our calculation for large enough times. A similar error accumulates in the determination of the next instant of time t+dt, but we will discuss below how to make this contribution to the total error negligible. The above mentioned errors can become less detrimental by using floating point numbers of greater accuracy than the float type. For example double numbers have approximately 17 significant decimal digits. But again, the precision is finite and the same type of errors are there only to be revealed by a more demanding and complicated calculation.

The remedy to such a problem can only be a change in the algorithm. This is not always possible, but in the case at hand this is easy to do. For example, consider the equation that gives the position of a particle in free motion

(2.24)

Let’s use the above relation for the parts of the motion between two collisions. Then, all we have to do is to reverse the direction of the motion and reset the initial position and time to be the position and time of the collision. This can be done by using the loop:

In the above algorithm, the error in the time of the collision is not vanishing but we don’t have the “instability” problem of the dt limit¹⁹ . Therefore we can isolate and study the effect of each type of error. The full program that implements the above algorithm is given below and can be found in the file box1D_2.cpp:

Compiling and running the above program is done as before and the results are stored in the file box1D_2.dat.

2.3.2 Errors

In this section we will study the effect of the systematic errors that we encountered in the previous section in more detail. We considered two types of errors: First, the systematic error of determining the instant of the collision of the particle with the wall. This error is reduced by taking a smaller time step . Then, the systematic error that accumulates with each addition of two numbers with increasing difference in their orders of magnitude. This error is increased with decreasing . The competition of the two effects makes the optimal choice of the result of a careful analysis. Such a situation is found in many interesting problems, therefore it is quite instructive to study it in more detail.

When the exact solution of the problem is not known, the systematic errors are controlled by studying the behavior of the solution as a function of . If the solutions are converging in a region of values of , one gains confidence that the true solution has been determined up to the accuracy of the convergence.

In the previous sections, we studied two different algorithms, programmed in the files box1D_1.cpp and box1D_2.cpp. We will refer to them as “method 1” and “method 2” respectively. We will study the convergence of the results as by fixing all the parameters except and then study the dependence of the results on . We will take , , , , , so that the particle will collide with the wall every 10 units of time. We will measure the position of the particle ²⁰ as a function of and study its convergence to a limit²¹ as .

The analysis requires a lot of repetitive work: Compiling, setting the parameter values, running the program and calculating the value of for many values of . We write the values of the parameters read by the program in a file box1D_anal.in:

Then we compile the program

and run it with the command:

By using the pipe |, we send the contents of box1D_anal.in to the stdin of the command ./box by using the command cat. The result can be found in the last line of the file box1D_1.dat:

The third number in the above line is the value of the velocity. In a file box1D_anal.dat we write and the first two numbers coming out from the command tail. Then we decrease the value in the file box1D_anal.in and run again. We repeat for 12 more times until reaches the value²² . We do the same²³ using method 2 and we place the results for in two new columns in the file box1D_anal.dat. The result is

Convergence is studied in figure 2.21. The 1st method maximizes its accuracy for , whereas for the error becomes % and the method becomes useless. The 2nd method has much better behavior that the 1st one.

We observe that as decreases, the final value of approaches the expected . Why don’t we obtain , especially when is an integer? How many steps does it really take to reach , when the expected number of those is ? Each time you take a measurement, issue the command

which measures the number of lines in the file box1D_1.dat and compare this number with the expected one. The result is interesting:

where the second column has the number of steps computed by the program and the third one has the expected number of steps. We observe that the accuracy decreases with decreasing and in the end the difference is about 5%! Notice that the last line should have given , an error comparable to the period of the particle’s motion.

We conclude that one important source of accumulation of systematic errors is the calculation of time. This type of errors become more significant with decreasing . We can improve the accuracy of the calculation significantly if we use the multiplication t=t0+i*dt instead of the addition t=t+dt, where i is a step counter:

The main loop in the program box1D_1.cpp becomes:

The full program can be found in the file box1D_4.cpp of the accompanying software. We call this “method 3”. We perform the same change in the file box1D_2.cpp, which we store in the file box1D_5.cpp. We call this “method 4”. We repeat the same analysis using methods 3 and 4 and we find that the problem of calculating time accurately practically vanishes. The result of the analysis can be found on the right plot of figure 2.21. Methods 2 and 4 have no significant difference in their results, whereas methods 1 and 3 do have a dramatic difference, with method 3 decreasing the error more than tenfold. The problem of the increase of systematic errors with decreasing does not vanish completely due to the operation x=x+v*dt. This type of error is harder to deal with and one has to invent more elaborate algorithms in order to reduce it significantly. This will be discussed further in chapter 4.

Figure 2.21: The error where is the value calculated by method and the exact value according to the text.

2.3.3 The Two Dimensional Box

A particle is confined to move on the plane in the area and . When it reaches the boundaries of this two dimensional box, it bounces elastically off its walls. The particle is found in an infinite depth orthogonal potential well. The particle starts moving at time from and our program will calculate its trajectory until time with time step . Such a trajectory can be seen in figure 2.23.

If the particle’s position and velocity are known at time , then at time they will be given by the relations

The full program can be found in the file box2D_1.cpp. Notice that we introduced two counters nx and ny of the particle’s collisions with the walls:

A typical session for the study of a particle’s trajectory could be:

Notice the last line of output from the program: The particle bounces off the vertical walls 6 times (nx=6) and from the horizontal ones 13 (ny=13). The gnuplot commands construct the diagrams displayed in figures 2.22 and 2.23.

Figure 2.22: The results for the trajectory of a particle in a two dimensional box given by the program box2D_1.cpp. The parameters are , , , , , , , , .

Figure 2.23: The trajectory of the particle of figure 2.22 until . The origin of the arrow is at the initial position of the particle and its end is at its current position. The bold lines mark the boundaries of the box.

In order to animate the particle’s trajectory, we can copy the file box2D_animate.gnu of the accompanying software to the current directory and give the gnuplot commands:

The last line repeats the same animation at half speed. You can also use the file animate2D.gnu discussed in section 2.1.1. We add new commands in the file box2D_animate.gnu so that the plot limits are calculated automatically and the box is drawn on the plot. The arrow drawn is not the position vector with respect to the origin of the coordinate axes, but the one connecting the initial with the current position of the particle.

The next step should be to test the accuracy of your results. This can be done by generalizing the discussion of the previous section and it is left as an exercise for the reader.

2.4 Applications

In this section we will study simple examples of motion in a box with different types of obstacles. We will start with a game of ... mini golf. The player shoots a (point) “ball” which moves in an orthogonal box of linear dimensions and and which is open on the side. In the box there is a circular “hole” with center at and radius . If the “ball” falls in the “hole”, the player wins. If the ball leaves out of the box through its open side, the player loses. In order to check if the ball is in the hole when it is at position , all we have to do is to check whether .

Figure 2.24: The trajectory of the particle calculated by the program MiniGolf.cpp using the parameters chosen in the text. The moment of ... success is shown. At time the particle enters the hole’s region which has its center at and its radius is .

Initially we place the ball at the position at time . The player hits the ball which leaves with initial velocity of magnitude at an angle degrees with the axis. The program is found in the file MiniGolf.cpp and is listed below:

In order to run it, we can use the commands:

You should construct the plots of the position and the velocity of the particle. You can also use the animation program found in the file MiniGolf_animate.gnu for fun. Copy it from the accompanying software to the current directory and give the gnuplot commands:

The results are shown in figure 2.24.

The next example with be three dimensional. We will study the motion of a particle confined within a cylinder of radius and height . The collisions of the particle with the cylinder are elastic. We take the axis of the cylinder to be the axis and the two bases of the cylinder to be located at and . This is shown in figure 2.26.

The collisions of the particle with the bases of the cylinder are easy to program: we follow the same steps as in the case of the simple box. For the collision with the cylinder’s side, we consider the projection of the motion on the plane. The projection of the particle moves within a circle of radius and center at the intersection of the axis with the plane. This is shown in figure 2.25. At the collision, the component of the velocity is reflected , whereas remains the same. The velocity of the particle before the collision is

Figure 2.25: The elastic collision of the particle moving within the circle of radius and center at the point . We have that . The initial velocity is where . After reflecting the new velocity of the particle is .

The transformation , will be performed in the function reflectVonCircle(vx,vy,x,y,xc,yc,R). Upon entry to the function, we provide the initial velocity (vx,vy), the collision point (x,y), the center of the circle (xc,yc) and the radius of the circle²⁴ R. Upon exit from the function, (vx,vy) have been replaced with the new values²⁵ .

The program can be found in the file Cylinder3D.cpp and is listed below:

Note that the function atan2 is used for computing the angle theta. This function, when called with two arguments atan2(y,x), returns the angle in radians. The correct quadrant of the circle where lies is chosen. The angle that we want to compute is given by atan2(y-yc,x-xc). Then we apply equations (2.29) and (2.31) and in the last two lines we enforce the particle to be at the point , exactly on the circle.

Figure 2.26: The trajectory of a particle moving inside a cylinder with , , computed by the program Cylinder3D.cpp. We have chosen , , , , .

A typical session is shown below:

In order to plot the position and the velocity as a function of time, we use the following gnuplot commands:

We can also compute the distance of the particle from the cylinder’s axis as a function of time using the command:

In order to plot the trajectory, together with the cylinder, we give the commands:

The command set parametric is necessary if one wants to make a parametric plot of a surface . The cylinder (without the bases) is given by the parametric equations with , .

We can also animate the trajectory with the help of the gnuplot script file Cylinder3D_animate.gnu. Copy the file from the accompanying software to the current directory and give the gnuplot commands:

The result is shown in figure 2.26.

The last example will be that of a simple model of a spacetime wormhole. This is a simple spacetime geometry which, in the framework of the theory of general relativity, describes the connection of two distant areas in space which are asymptotically flat. This means, that far enough from the wormhole’s mouths, space is almost flat - free of gravity. Such a geometry is depicted in figure 2.27. The distance traveled by someone through the mouths could be much smaller than the distance traveled outside the wormhole and, at least theoretically, traversable wormholes could be used for interstellar/intergalactic traveling and/or communications between otherwise distant areas in the universe. Of course we should note that such macroscopic and stable wormholes are not known to be possible to exist in the framework of general relativity. One needs an exotic type of matter with negative energy density which has never been observed. Such exotic geometries may realize microscopically as quantum fluctuations of spacetime and make the small scale structure of the geometry²⁶ a “spacetime foam”.

We will study a very simple model of the above geometry on the plane with a particle moving freely in it²⁷ .

Figure 2.28: A simple model of the spacetime geometry of figure 2.27. The particle moves on the whole plane except withing the two disks that have been removed. The neck of the wormhole is modeled by the two circles , , and has zero length since their points have been identified. There is a given direction in this identification, so that points with the same are the same (you can imagine how this happens by folding the plane across the axis and then glue the two circles together). The entrance of the particle through one mouth and exit through the other is done as shown for the velocity vector .

We take the two dimensional plane and cut two equal disks of radius with centers at distance like in figure 2.28. We identify the points on the two circles such that the point 1 of the left circle is the same as the point 1 on the right circle, the point 2 on the left with the point 2 on the right etc. The two circles are given by the parametric equations , , for the right circle and , , for the left. Points on the two circles with the same are identified. A particle entering the wormhole from the left circle with velocity is immediately exiting from the right with velocity as shown in figure 2.28.

Then we will do the following:

1. Write a program that computes the trajectory of a particle moving in the geometry of figure 2.28. We set the limits of motion to be and . We will use periodic boundary conditions in order to define what happens when the particle attempts to move outside these limits. This means that we identify the line with the line as well as the line with the line. The user enters the parameters , and as well as the initial conditions , where . The user will also provide the time parameters and for motion in the time interval with step .
2. Plot the particle’s trajectory with , , in the geometry with .
3. Find a closed trajectory which does not cross the boundaries , and determine whether it is stable under small perturbations of the initial conditions.
4. Find other closed trajectories that go through the mouths of the wormhole and study their stability under small perturbations of the initial conditions.
5. Add to the program the option to calculate the distance traveled by the particle. If the particle starts from and moves in the direction to the , position, draw the trajectory and calculate the distance traveled on paper. Then confirm your calculation from the numerical result coming from your program.
6. Change the boundary conditions, so that the particle bounces off elastically at , and replot all the trajectories mentioned above.

Define the right circle by the parametric equations

(2.32)

and the left circle by the parametric equations

(2.33)

The particle’s position changes at time by

Figure 2.29: The particle crossing the wormhole through the right circle with velocity . It emerges from with velocity . The unit vectors , are computed from the parametric equations of the two circles and .

Crossing the circle is determined by the relation

(2.35)

The angle is calculated from the equation

(2.36)

and the point is mapped to the point where

(2.37)

as can be seen in figure 2.29. For mapping , we first calculate the vectors

(2.38)

so that the velocity

(2.39)

where the radial components are and . Therefore, the relations that give the “emerging” velocity are:

(2.40)

Similarly we calculate the case of entering from and emerging from . The condition now is:

(2.41)

The angle is given by

(2.42)

and the point is mapped to the point where

(2.43)

For mapping , we calculate the vectors

(2.44)

so that the velocity

(2.45)

The emerging velocity is:

(2.46)

Systematic errors are now coming from crossing the two mouths of the wormhole. There are no systematic errors from crossing the boundaries , (why?). Try to think of ways to control those errors and study them.

The closed trajectories that we are looking for come from the initial conditions

(2.47)

and they connect points 1 of figure 2.28. They are unstable, as can be seen by taking .

The closed trajectories that cross the wormhole and “wind” through space can come from the initial conditions

It is easy to compile and run the program. See also the files Wormhole.csh and Wormhole_animate.gnu of the accompanying software and run the gnuplot commands:

You are now ready to answer the rest of the questions that we asked in our list.

2.5 Problems

1. Change the program Circle.cpp so that it prints the number of full circles traversed by the particle.
2. Add all the necessary tests on the parameters entered by the user in the program Circle.cpp, so that the program is certain to run without problems. Do the same for the rest of the programs given in the same section.
3. A particle moves with constant angular velocity on a circle that has the origin of the coordinate system at its center. At time , the particle is at . Write the program CircularMotion.cpp that will calculate the particle’s trajectory. The user should enter the parameters . The program should print the results like the program Circle.cpp does.
4. Change the program SimplePendulum.cpp so that the user could enter a non zero initial velocity.
5. Study the limit in the projectile motion given by equations (2.10) . Expand and keep the non vanishing terms as . Then keep the next order leading terms which have a smaller power of . Program these relations in a file
   ProjectileSmallAirResistance.cpp. Consider the initial conditions and calculate the range of the trajectory numerically by using the two programs
   ProjectileSmallAirResistance.cpp, ProjectileAirResistance.cpp. Determine the range of values of for which the two results agree within 5% accuracy.
6. Write a program for a projectile which moves through a fluid with fluid resistance proportional to the square of the velocity. Compare the range of the trajectory with the one calculated by the program ProjectileAirResistance.cpp for the parameters shown in figure 2.10.
7. Change the program Lissajous.cpp so that the user can enter a different amplitude and initial phase in each direction. Study the case where the amplitudes are the same and the phase difference in the two directions are . Repeat by taking the amplitude in the direction to be twice as much the amplitude in the direction.
8. Change the program ProjectileAirResistance.cpp, so that it can calculate also the case.
9. Change the program ProjectileAirResistance.cpp so that it can calculate the trajectory of the particle in three dimensional space. Plot the position coordinates and the velocity components as a function of time. Plot the three dimensional trajectory using splot in gnuplot and animate the trajectory using the gnuplot script animate3D.gnu.
10. Change the program ChargeInB.cpp so that it can calculate the number of full revolutions that the projected particle’s position on the plane makes during its motion.
11. Change the program box1D_1.cpp so that it prints the number of the particle’s collisions on the left wall, on the right wall and the total number of collisions to the stdout.
12. Do the same for the program box1D_2.cpp. Fill the table on page 290 the number of calculated collisions and comment on the results.
13. Run the program box1D_1.cpp and choose L= 10, v0=1. Decrease the step dt up to the point that the particle stops to move. For which value of dt this happens? Increase v0=10,100. Until which value of dt the particle moves now? Why?
14. Change the float declarations to double in the program box1D_1.cpp. Make sure that all the constants that you use become double precision (e.g. 1.0f changes to 1.0). Compare your results to those obtained in section 2.3.2. Repeat problem 2.13. What do you observe?
15. Change the program box1D_1.cpp so that you can study non elastic collisions , with the walls.
16. Change the program box2D_1.cpp so that you can study inelastic collisions with the walls, such that , , .
17. Use the method of calculating time in the programs box1D_4.cpp and box1D_5.cpp in order to produce the results in figure 2.21.
18. Particle falls freely moving in the vertical direction. It starts with zero velocity at height . Upon reaching the ground, it bounces inelastically such that with a parameter. Write the necessary program in order to study numerically the particle’s motion and study the cases .
19. Generalize the program of the previous problem so that you can study the case . Animate the calculated trajectories.
20. Study the motion of a particle moving inside the box of figure 2.30. Count the number of collisions of the particle with the walls before it leaves the box.

    ---

    
    

    Figure 2.30: Problem 2.20.

    ---

21. Study the motion of the point particle on the “billiard table” of figure 2.31. Count the number of collisions with the walls before the particle enters into a hole. The program should print from which hole the particle left the table.

    ---

    
    

    Figure 2.31: Problem 2.21.

    ---

22. Write a program in order to study the motion of a particle in the box of figure 2.32. At the center of the box there is a disk on which the particle bounces off elastically (Hint: use the routine reflectVonCircle of the program Cylinder3D.cpp).

    ---

    
    

    Figure 2.32: Problem 2.22.

    ---

23. In the box of the previous problem, put four disks on which the particle bounces of elastically like in figure 2.33.

    ---

    
    

    Figure 2.33: Problem 2.23.

    ---

24. Consider the arrangement of figure 2.34. Each time the particle bounces elastically off a circle, the circle disappears. The game is over successfully if all the circles vanish. Each time the particle bounces off on the wall to the left, you lose a point. Try to find trajectories that minimize the number of lost points.

    ---

    
    

    Figure 2.34: Problem 2.24.

    ---

Chapter 3
Logistic Map

Nonlinear differential equations model interesting dynamical systems in physics, biology and other branches of science. In this chapter we perform a numerical study of the discrete logistic map as a “simple mathematical model with complex dynamical properties”  [23] similar to the ones encountered in more complicated and interesting dynamical systems. For certain values of the parameter of the map, one finds chaotic behavior giving us an opportunity to touch on this very interesting topic with important consequences in physical phenomena. Chaotic evolution restricts out ability for useful predictions in an otherwise fully deterministic dynamical system: measurements using slightly different initial conditions result in a distribution which is indistinguishable from the distribution coming from sampling a random process. This scientific field is huge and active and we refer the reader to the bibliography for a more complete introduction  [23, 24, 25, 26, 27, 28, 29, 40].

3.1 Introduction

The most celebrated application of the logistic map comes from the study of population growth in biology. One considers populations which reproduce at fixed time intervals and whose generations do not overlap.

The simplest (and most naive) model is the one that makes the reasonable assumption that the rate of population growth of a population is proportional to the current population: