MacKinnon - Computational Physics

Computational Physics — 3rd/4th Year OptionAngus MacKinnonSeptember 26, 2002Computational PhysicsCourse DescriptionThe use of computers in physics, as well as most other branches of science and engineering, has increasedmany times along with the rapid development of faster and cheaper hardware.

This course aims to give thestudent a thorough grounding in the main computational techniques used in modern physics. It is particularly important in this course that the students should learn by doing. The course is therefore designed suchthat a significant fraction of the students’ time is spent actually programming specific physical problemsrather than learning abstract techniques.The course will cover problems in 4(5) broad sections:Ordinary differential equations, such as those of classical mechanics.Partial differential equations, such as Maxwell’s equations and the Diffusion and Schrödinger equations.Matrix methods, such as systems of equations and eigenvalue problems applied to Poisson’s equationand electronic structure calculations.Monte Carlo and other simulation methods, such as the Metropolis algorithm and molecular dynamics.(If time permits:) Computer Algebra; an introduction using Maple to the uses and abuses of algebraiccomputing in physics.This is not a short course in computing science, nor in programming.

It focuses specifically on methodsfor solving physics problems. Students will be expected to be familiar with basic programming: successfulcompletion of the 1st year computing Lab. is sufficient. There is no requirement that the practical work bedone using Microsoft C++ on the departmental computers, but anyone wishing to use some other programming language or computer should consult the lecturer beforehand. This is to make sure there is both helpavailable from demonstrators and that it will be possible to assess the work satisfactorily.A detailed timetable for the course is given below.Each of the 4(5) sections of the course will be accompanied by a short problem sheet including practicalexercises and longer projects.Students will be expected to do the first project and 2 others. There will be a choice of projects availablein the later sections.This year you will be required to submit your program electronically in addition to including it as anappendix in your project reports.

You should also include an example input file and the correspondingoutput file. This will enable the assessors to check and run your program if necessary. Your programshould be sent as a plain text e-mail attachment(s)1 to Computational-Physics, A. MacKinnonor ph.compphys@ic.ac.uk by the stated deadline. Do not send it to my personal e-mail address. Thefirst project is to be handed in to Martin Morris in the Computing Lab by 5.00pm on the 28th October.These will be assessed and returned by the demonstrators who should provide the students with feedbackwhich will be of use in the longer projects later on. This project will count forof the final result.1 UsingOutlook click Insert then File1During the Lab. sessions staff and demonstrators will be available, but the students are expected todo the practical work and the problems outside these times as well.

I am aware that there is an GeneralRelativity lecture scheduled for 3.00pm on Tuesdays. Nevertheless the Lab. will be reserved for this coursefor the whole afternoon and demonstrators will be available during this period.The report on the 2nd project, which should be chosen from the 2 offered on Partial Differential Equations, should be submitted by 5.00pm on November 18th.The 3rd project may be chosen from those offered on Matrix and Monte Carlo methods (or possibly onComputer Algebra) and the report should be submitted by 5.00pm on Friday December 13th, the last dayof term. This timetable has been designed to make sure that students are working at a steady rate duringthe term.From 10.00am – 12.00 noon on Monday 6th January 2003, there will be a short written test under examconditions.

This will count forof the final result. The 2nd and 3rd projects will also count for.Thus the distribution of marks for the various parts of the course will be approximatelyProject 1Project 2Project 3TestWeekSep. 30Oct. 7Oct. 14Oct. 21Oct. 28Nov. 4Nov. 11Nov. 18Nov. 25Dec. 2Dec. 9Jan. 6LecturesThu. 9.00amFri. 2.00pmFri. 4.00pmFri. 4.00pmFri. 4.00pmFri. 4.00pmFri.

4.00pm10303030PracticalsTue. 2.00–5.00pmTue. 2.00–5.00pmTue. 2.00–5.00pmTue. 2.00–5.00pmTue. 2.00–5.00pmThu. 9.00am Tue. 2.00–5.00pmThu. 9.00am Tue. 2.00–5.00pmThu. 9.00am Tue. 2.00–5.00pmThu. 9.00am Tue. 2.00–5.00pmThu. 9.00am Tue. 2.00–5.00pmTest: Mon. 10.00am–12.00noonDeadlinesMon. 5.00pm 1st ProjectMon. 5.00pm 2nd ProjectFri. 5.00pm 3rd ProjectFor those with surnames beginning with A–K the test will be in lecture theatre 2 and for those withsurnames beginning with L–Z in lecture theatre 3.There will be no examination in the summer.A copy of the timetable as well as contact details for the demonstrators is posted on the Exchange server.This can be accessed from Outlook under //Public Folders/All Public Folders/Physics/Computational Physics.Further information is available on the World Wide Web underhttp://www.cmth.ph.imperial.ac.uk/angus/Lectures/compphys/.Chapter 1Ordinary Differential Equations1.1 Types of Differential EquationIn this chapter we will consider the methods of solution of the sorts of ordinary differential equations(ODEs) which occur very commonly in physics.

By ODEs we mean equations involving derivatives withrespect to a single variable, usually time. Although we will formulate the discussion in terms of linearODEs for which we know the analytical solution, this is simply to enable us to make comparisons betweenthe numerical and analytical solutions and does not imply any restriction on the sorts of problems to whichthe methods can be applied.

In the practical work you will encounter examples which do not fit neatly intothese categories.The work in this section is also considered in chapter 16 of Press et al. (1992) or chapter II of Potter(1973)We consider 3 basic differential equations: $ , ()+* !#"%$ %& !#" ' %& !#".- )+* %&The Decay equation(1.1)The Growth equation(1.2)The Oscillation equation(1.3)which are representative of most more complex cases.Higher order differential equations can be reduced to 1st order by appropriate choice of additional variables. The simplest such choice is to define new variables to represent all but the highest order derivative.For example, the damped harmonic oscillator equation, usually written as/ 1 0 0 32 54 (1.4)17 in the form of a pair of 1st order ODEscan be rewritten in terms of and velocity 6 6 2 4 1 (1.5) / $6 / 16(1.6)Similarly any 8 th order differential equation can be reduced to 8 1st order equations.Such systems of ODEs can be written in a very concise notation by defining a vector, y say, whoseelements are the unknowns, such as and 6 in (1.1).

Any ODE in 8 unknowns can then be written in thegeneral form y " %&; 9 f y :(1.7)3Ordinary Differential Equations48where y and f are –component vectors.Remember that there is no significance in the use of the letter in the above equations. The variable isnot necessarily time but could just as easily be space, as in (1.9), or some other physical quantity.Formally we can write the solution of (1.7) as"<%&; y <" &=$3>3? f " y "<CBD& : CBE&CB(1.8)@%GFH .

?AAlthoughby integrating both sides over the interval(1.8) is formally correct, in practice it is"<%&usually impossible to evaluate the integral on the right–hand–side as it presupposes the solution y . Weywill have to employ an approximation.All differential equations require boundary conditions. Here we will consider cases in which all theboundary conditions are defined at a particular value of (e.g.).

For higher order equations theboundary conditions may be defined at different values of . The modes of a violin string at frequencyobey the equationN with boundary conditions such that 0J L $ * 0 1K 0 M 0I *(1.9)at both ends of the string. We shall consider such problems inchapter 3.1.2 Euler MethodO PRQS=TU$VCQ> ?AWYXZ f " y "< B & : B & B O f " y "< Q & : Q &\[?W"]CQS#TJ&U y "]CQ1&=$ O f " y "<CQ& : CQ1&^[yConsider an approximate solution of (1.8) over a small intervalto obtainor, in a more concise notation,yQS#T yQ $ O f " yQ : Q &; yQ $ O fQ [by writing the integral as(1.10)O(1.11)We can integrate over any larger interval by subdividing the range into sections of width and repeating(1.11) for each part.Equivalently we can consider that we have approximated the derivative with a forward difference y yQS=T $ yQ [=__ QO__We will also come across centred and backward differences,$ y ac `b yQS=T f O yQeT centred __ Q yQ $ yQegT__ bd O backward(1.12)(1.13)respectively.Here we have used a notation which is very common in computational physics, in which we calculateyyat discrete values of given by, and ffy.In what follows we will drop the vector notation y except when it is important for the discussion.Q "<CQ&%Q 8gO Q " Q : CQ&Ordinary Differential Equations51.2.1 Order of AccuracyHow accurate is the Euler method? To quantify this we consider a Taylor expansion of QS=T h Q O _ Q___and substitute this into (1.11) Q O i_ Q O f C0 0 0 _ Q______0 0 O f 0 __ Q__ [[[ Q $lQ [[["<%& around Q(1.14)O %jk"A Q : Q &(1.15) (1.16)O __ Q :__where we have used (1.7) to obtain the final form.