Computational physics

345 Pages·2003·1.852 MB·English

by  Hjorth-Jensen M.

COMPUTATIONAL PHYSICS
M. Hjorth-Jensen
Department of Physics, University of Oslo, 2003

iii
Preface
In 1999, when we started teaching this course at the Department of Physics in Oslo, Compu-
tational Physics and Computational Science in general were still perceived by the majority of
physicistsand scientistsastopicsdealingwithjustmeretoolsand numbercrunching,andnotas
subjects of their own. The computational background of most students enlisting for the course
oncomputationalphysicscouldspanfromdedicatedhackersandcomputerfreakstopeoplewho
basicallyhadneverusedaPC. Themajorityofgraduatestudentshadaveryrudimentaryknowl-
edge of computationaltechniques and methods. Fouryears later most studentshavehad a fairly
uniform introduction to computers, basic programming skills and use of numerical exercises in
undergraduatecourses. Practically everyundergraduatestudentin physicshasnowmadeaMat-
lab or Maple simulation of e.g., the pendulum, with or without chaotic motion. These exercises
underscore the importance of simulations as a means to gain novel insights into physical sys-
tems, especially for those cases where no analytical solutions can be found or an experiment
is to complicated or expensive to carry out. Thus, computer simulations are nowadays an inte-
gral part of contemporary basic and applied research in the physical sciences. Computation is
becoming as important as theory and experiment. We could even strengthen this statement by
sayingthatcomputationalphysics,theoreticalphysicsandexperimentalareallequallyimportant
in our daily research and studies of physical systems. Physics is nowadays the unity of theory,
experiment and computation. The ability "to compute" is now part of the essential repertoire of
research scientists. Several new fields have emerged and strengthened their positions in the last
years, such as computational materials science, bioinformatics, computational mathematics and
mechanics, computationalchemistry and physicsand so forth, just to mention a few. To be able
to e.g., simulate quantal systems will be of great importance for future directions in fields like
materialsscienceand nanotechonology.
This ability combines knowledge from many different subjects, in our case essentially from
thephysicalsciences,numericalanalysis,computinglanguagesandsomeknowledgeofcomput-
ers. These topics are, almost as a rule of thumb, taught in different, and we would like to add,
disconnectedcourses. Onlyatthelevelofthesisworkisthestudentconfrontedwiththesynthesis
of all these subjects, and then in a bewildering and disparate manner, trying to e.g., understand
old Fortran 77 codes inherited from his/her supervisor back in the good old ages, or even more
archaic, programs. Hours may have elapsed in front of a screen which just says ’Underflow’, or
’Bus error’, etc etc, without fully understanding what goes on. Porting the program to another
machinecould evenresultintotallydifferentresults!
The first aim of this course is therefore to bridge the gap between undergraduate courses
in the physical sciences and the applications of the aquired knowledge to a given project, be it
eitherathesisworkoranindustrialproject. Weexpectyoutohavesomebasicknowledgeinthe
physical sciences, especially within mathematics and physics through e.g., sophomore courses
in basic calculus, linear algebraand general physics. Furthermore, having taken an introductory
courseonprogrammingissomethingwerecommend. Assuch,anoptimaltimingfortakingthis
course, wouldbe when you are close to embark on a thesiswork, orif you’vejust started with a
thesis. But obviously,you shouldfeel freeto chooseyourowntiming.
We haveseveralotheraimsas wellin additionto prepareyoufor athesiswork,namely
iv
We would like to give you an opportunity to gain a deeper understanding of the physics
(cid:15)
youhavelearnedinothercourses. Inmostcoursesoneisnormallyconfrontedwithsimple
systems which provide exact solutions and mimic to a certain extent the realistic cases.
Many are however the comments like ’why can’t we do something else than the box po-
tential?’. In several oftheprojects wehopeto present somemore’realistic’cases to solve
by various numerical methods. This also means that we wish to give examples of how
physics can be applied in a much broader context than it is discussed in the traditional
physicsundergraduatecurriculum.
To encourage you to "discover" physics in a way similar to how researchers learn in the
(cid:15)
contextofresearch.
Hopefullyalso tointroducenumericalmethodsand new areas ofphysicsthat can bestud-
(cid:15)
ied withthemethodsdiscussed.
Toteach structuredprogrammingin thecontextofdoingscience.
(cid:15)
The projects we propose are meant to mimic to a certain extent the situation encountered
(cid:15)
duringathesisorprojectwork. Youwilltipicallyhaveatyourdisposal1-2weekstosolve
numerically a given project. In so doing you may need to do a literature study as well.
Finally,wewouldlikeyouto writeareport forevery project.
Theexamreflectsthisproject-likephilosophy. Theexamitselfisaprojectwhichlastsone
(cid:15)
month. Youhavetohandinareportonaspecificproblem,andyourreportformsthebasis
foran oral examinationwithafinal grading.
Our overall goal is to encourage you to learn about science through experience and by asking
questions. Our objective is always understanding, not the generation of numbers. The purpose
of computing is further insight, not mere numbers! Moreover, and this is our personal bias, to
deviceanalgorithmandthereafterwriteacodeforsolvingphysicsproblemsisamarvelousway
of gaining insightinto complicated physical systems. The algorithmyou end up writing reflects
inessentiallyallcases yourownunderstandingofthephysicsoftheproblem.
Mostofyouarebynowfamiliar,throughvariousundergraduatecoursesinphysicsandmath-
ematics, with interpreted languages such as Maple, Mathlab and Mathematica. In addition, the
interestinscriptinglanguagessuch asPythonorPerl hasincreased considerablyinrecent years.
The modern programmer would typically combine several tools, computing environments and
programming languages. A typical example is the following. Suppose you are working on a
project which demands extensive visualizations of the results. To obtain these results you need
however a programme which is fairly fast when computational speed matters. In this case you
would most likely write a high-performance computingprogramme in languages which are tay-
lored for that. These are represented by programming languages like Fortran 90/95 and C/C++.
However,tovisualizetheresultsyouwouldfindinterpretedlanguageslikee.g.,Matlaborscript-
ing languages like Python extremely suitable for your tasks. You will therefore end up writing
e.g., a script in Matlab which calls a Fortran 90/95 ot C/C++ programme where the number
crunchingisdoneandthen visualizetheresultsofsay awaveequationsolverviaMatlab’slarge
v
library of visualizationtools. Alternatively,you could organize everythinginto a Python or Perl
script which does everything for you, calls the Fortran 90/95 or C/C++ programs and performs
thevisualizationinMatlabas well.
Being multilingualisthusafeaturewhichnotonlyappliesto modernsocietybuttocomput-
ingenvironmentsaswell.
However, there is more to the picture than meets the eye. This course emphasizes the use of
programminglanguageslikeFortran90/95andC/C++insteadofinterpretedoneslikeMatlabor
Maple. Computational speed is not the only reason for this choice of programming languages.
The main reason is that we feel at a certain stage one needs to have some insights into the algo-
rithmused, itsstabilityconditions,possiblepitfallslikelossofprecision, ranges ofapplicability
etc. Although we will at various stages recommend the use of library routines for say linear
algebra1, our belief is that one should understand what the given function does, at least to have
a mere idea. From such a starting point we do further believe that it can be easier to develope
morecomplicatedprograms,onyourown. Wedothereforedevotesomespacetothealgorithms
behind variousfunctions presented in the text. Especially,insightintohow errors propagateand
how to avoid them is a topic we’d like you to pay special attention to. Only then can you avoid
problemslikeunderflow,overflowand lossofprecision. Suchacontrolisnotalwaysachievable
withinterpreted languagesand canned functionswhere theunderlyingalgorithm
Needlesstosay,theselecturenotesareupgradedcontinuously,fromtypostonewinput. And
we do always benifit from your comments, suggestions and ideas for making these notes better.
It’sthroughthescientificdiscourseandcritics weadvance.
1Such library functionsare often tayloredto a givenmachine’sarchitecture and should accordinglyrun faster
thanuserprovidedones.

Contents
I Introduction to Computational Physics 1
1 Introduction 3
1.1 Choiceofprogramminglanguage . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Designingprograms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2 Introduction to C/C++andFortran90/95 9
2.1 Gettingstarted . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.1.1 Representation ofintegernumbers . . . . . . . . . . . . . . . . . . . . . 15
2.2 Real numbersand numericalprecision . . . . . . . . . . . . . . . . . . . . . . . 18
2.2.1 Representation ofreal numbers . . . . . . . . . . . . . . . . . . . . . . 19
2.2.2 Further examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.3 Lossofprecision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3.1 Machinenumbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3.2 Floating-pointerror analysis . . . . . . . . . . . . . . . . . . . . . . . . 32
2.4 Additionalfeatures ofC/C++ andFortran 90/95 . . . . . . . . . . . . . . . . . . 33
2.4.1 Operators in C/C++ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.4.2 Pointers and arrays inC/C++. . . . . . . . . . . . . . . . . . . . . . . . 35
2.4.3 Macros inC/C++ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.4.4 Structures in C/C++ andTYPE in Fortran90/95 . . . . . . . . . . . . . 39
3 Numerical differentiation 41
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2 Numericaldifferentiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2.1 Thesecond derivativeof . . . . . . . . . . . . . . . . . . . . . . . . . 45
x
e
3.2.2 Error analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.2.3 Howto makefigures withGnuplot . . . . . . . . . . . . . . . . . . . . . 54
3.3 Richardson’sdeferred extrapolationmethod . . . . . . . . . . . . . . . . . . . . 57
4 Classes,templates andmodules 61
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.2 Afirst encounter, thevectorclass . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.3 Classesand templatesinC++ . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.4 UsingBlitz++ withvectorsandmatrices . . . . . . . . . . . . . . . . . . . . . . 68
vii
viii CONTENTS
4.5 Buildingnew classes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.6 MODULEandTYPE declarations inFortran 90/95 . . . . . . . . . . . . . . . . 68
4.7 ObjectorientinginFortran 90/95 . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.8 Anexampleofuseofclassesin C++and ModulesinFortran 90/95 . . . . . . . . 68
5 Linear algebra 69
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.2 Programmingdetails . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.2.1 Declaration offixed-sized vectorsandmatrices . . . . . . . . . . . . . . 70
5.2.2 Runtimedeclarations ofvectorsand matrices . . . . . . . . . . . . . . . 72
5.2.3 Fortran features ofmatrixhandling . . . . . . . . . . . . . . . . . . . . 75
5.3 LUdecompositionofamatrix . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.4 Solutionoflinearsystemsofequations . . . . . . . . . . . . . . . . . . . . . . . 80
5.5 Inverseofamatrixand thedeterminant . . . . . . . . . . . . . . . . . . . . . . 81
5.6 Project: Matrixoperations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
6 Non-linear equations and roots ofpolynomials 87
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.2 Iterationmethods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
6.3 Bisectionmethod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
6.4 Newton-Raphson’smethod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.5 Thesecant methodand othermethods . . . . . . . . . . . . . . . . . . . . . . . 94
6.5.1 Calling thevariousfunctions . . . . . . . . . . . . . . . . . . . . . . . . 97
6.6 Rootsofpolynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.6.1 Polynomialsdivision . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.6.2 Root findingby Newton-Raphson’smethod . . . . . . . . . . . . . . . . 97
6.6.3 Root findingby deflation . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.6.4 Bairstow’smethod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
7 Numerical interpolation, extrapolationandfitting ofdata 99
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
7.2 Interpolationand extrapolation . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
7.2.1 Polynomialinterpolationand extrapolation . . . . . . . . . . . . . . . . 99
7.3 Qubicsplineinterpolation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
8 Numerical integration 105
8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
8.2 Equalstepmethods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
8.3 Gaussianquadrature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
8.3.1 Orthogonalpolynomials,Legendre . . . . . . . . . . . . . . . . . . . . 112
8.3.2 Mesh pointsand weightswithorthogonalpolynomials . . . . . . . . . . 115
8.3.3 Applicationtothecase . . . . . . . . . . . . . . . . . . . . . . . 116
N = 2
8.3.4 General integrationintervalsforGauss-Legendre . . . . . . . . . . . . . 117
CONTENTS ix
8.3.5 Otherorthogonalpolynomials . . . . . . . . . . . . . . . . . . . . . . . 118
8.3.6 Applicationsto selected integrals . . . . . . . . . . . . . . . . . . . . . 120
8.4 TreatmentofsingularIntegrals . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
9 Outlineofthe Monte-Carlo strategy 127
9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
9.1.1 First illustrationoftheuseofMonte-Carlomethods,crudeintegration . . 129
9.1.2 Second illustration,particles inabox . . . . . . . . . . . . . . . . . . . 134
9.1.3 Radioactivedecay . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
9.1.4 Program exampleforradioactivedecay ofonetypeofnucleus . . . . . . 137
9.1.5 Brief summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
9.2 PhysicsProject: Decay of Bi and Po . . . . . . . . . . . . . . . . . . . . . 140
210 210
9.3 Randomnumbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
9.3.1 Properties ofselected randomnumbergenerators . . . . . . . . . . . . . 144
9.4 Probabilitydistributionfunctions . . . . . . . . . . . . . . . . . . . . . . . . . . 146
9.4.1 Thecentral limittheorem . . . . . . . . . . . . . . . . . . . . . . . . . . 148
9.5 ImprovedMonteCarlo integration . . . . . . . . . . . . . . . . . . . . . . . . . 149
9.5.1 Change ofvariables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
9.5.2 Importancesampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
9.5.3 Acceptance-Rejection method . . . . . . . . . . . . . . . . . . . . . . . 157
9.6 MonteCarlo integrationofmultidimensionalintegrals . . . . . . . . . . . . . . . 157
9.6.1 Brute forceintegration . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
9.6.2 Importancesampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
10 Random walksand theMetropolis algorithm 163
10.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
10.2 Diffusionequationand randomwalks . . . . . . . . . . . . . . . . . . . . . . . 164
10.2.1 Diffusionequation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164
10.2.2 Random walks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167
10.3 Microscopicderivationofthediffusionequation . . . . . . . . . . . . . . . . . . 172
10.3.1 Discretized diffusionequationand Markovchains . . . . . . . . . . . . . 172
10.3.2 Continuousequations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
10.3.3 Numerical simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177
10.4 TheMetropolisalgorithmand detailedbalance . . . . . . . . . . . . . . . . . . 180
10.5 Physicsproject: simulationoftheBoltzmanndistribution . . . . . . . . . . . . . 184
11 Monte Carlomethods instatisticalphysics 187
11.1 Phasetransitionsinmagneticsystems . . . . . . . . . . . . . . . . . . . . . . . 187
11.1.1 Theoretical background . . . . . . . . . . . . . . . . . . . . . . . . . . 187
11.1.2 TheMetropolisalgorithm . . . . . . . . . . . . . . . . . . . . . . . . . 193
11.2 Program example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
11.2.1 Program forthetwo-dimensionalIsingModel . . . . . . . . . . . . . . . 195
11.3 Selected resultsfortheIsingmodel . . . . . . . . . . . . . . . . . . . . . . . . . 199
x CONTENTS
11.3.1 Phase transitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199
11.3.2 Heat capacity and susceptibilityas functionsofnumberofspins . . . . . 200
11.3.3 Thermalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201
11.4 Otherspinmodels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201
11.4.1 Potts model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201
11.4.2 XY-model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201
11.5 Physicsproject: simulationoftheIsingmodel . . . . . . . . . . . . . . . . . . . 201
12 Quantum Monte Carlomethods 203
12.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
12.2 VariationalMonteCarlo forquantummechanicalsystems . . . . . . . . . . . . . 204
12.2.1 FirstillustrationofVMCmethods,theone-dimensionalharmonicoscillator206
12.2.2 Thehydrogenatom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 209
12.2.3 Metropolissamplingforthehydrogenatomand theharmonicoscillator . 211
12.2.4 A nucleon inagaussianpotential . . . . . . . . . . . . . . . . . . . . . 215
12.2.5 Theheliumatom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 216
12.2.6 Program exampleforatomicsystems . . . . . . . . . . . . . . . . . . . 221
12.3 Simulationofmolecularsystems . . . . . . . . . . . . . . . . . . . . . . . . . . 228
12.3.1 TheH molecule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 228
+
12.3.2 Physic2sproject: theH molecule . . . . . . . . . . . . . . . . . . . . . . 230
2
12.4 Many-bodysystems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 230
12.4.1 Liquid He . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 230
4
12.4.2 Bose-Einsteincondensation . . . . . . . . . . . . . . . . . . . . . . . . 232
12.4.3 Quantumdots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233
12.4.4 Multi-electronatoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233
13 Eigensystems 235
13.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 235
13.2 Eigenvalueproblems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 235
13.2.1 Similaritytransformations . . . . . . . . . . . . . . . . . . . . . . . . . 236
13.2.2 Jacobi’smethod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 237
13.2.3 DiagonalizationthroughtheHouseholder’smethodfortri-diagonalization 238
13.3 Schrödinger’sequation(SE)throughdiagonalization . . . . . . . . . . . . . . . 241
13.4 Physicsprojects: Bound statesin momentumspace . . . . . . . . . . . . . . . . 248
13.5 Physicsprojects: Quantummechanical scattering . . . . . . . . . . . . . . . . . 251
14 Differential equations 255
14.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 255
14.2 Ordinarydifferentialequations(ODE) . . . . . . . . . . . . . . . . . . . . . . . 255
14.3 Finitedifferencemethods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 257
14.3.1 ImprovementstoEuler’salgorithm,higher-ordermethods . . . . . . . . 259
14.4 Moreon finitedifferencemethods,Runge-Kuttamethods . . . . . . . . . . . . . 260
14.5 AdaptiveRunge-Kuttaand multistepmethods . . . . . . . . . . . . . . . . . . . 261