Table of Contents

1. Cover
2. Related Titles
3. Title
4. Authors
5. Copyright
6. Dedication
7. Preface
8. Acknowledgments
9. 1: Introduction
   1. 1.1 Computational Physics and Computational Science
   2. 1.2 This Book’s Subjects
   3. 1.3 This Book’s Problems
   4. 1.4 This Book’s Language: The Python Ecosystem
   5. 1.5 Python’s Visualization Tools
   6. 1.6 Plotting Exercises
   7. 1.7 Python’s Algebraic Tools
10. 2: Computing Software Basics
    1. 2.1 Making Computers Obey
    2. 2.2 Programming Warmup
    3. 2.3 Python I/O
    4. 2.4 Computer Number Representations (Theory)
    5. 2.5 Problem: Summing Series
11. 3: Errors and Uncertainties in Computations
    1. 3.1 Types of Errors (Theory)
    2. 3.2 Error in Bessel Functions (Problem)
    3. 3.3 Experimental Error Investigation
12. 4: Monte Carlo: Randomness, Walks, and Decays
    1. 4.1 Deterministic Randomness
    2. 4.2 Random Sequences (Theory)
    3. 4.3 Random Walks (Problem)
    4. 4.4 Extension: Protein Folding and Self-Avoiding Random Walks
    5. 4.5 Spontaneous Decay (Problem)
    6. 4.6 Decay Implementation and Visualization
13. 5: Differentiation and Integration
    1. 5.1 Differentiation
    2. 5.2 Forward Difference (Algorithm)
    3. 5.3 Central Difference (Algorithm)
    4. 5.4 Extrapolated Difference (Algorithm)
    5. 5.5 Error Assessment
    6. 5.6 Second Derivatives (Problem)
    7. 5.7 Integration
    8. 5.8 Quadrature as Box Counting (Math)
    9. 5.9 Algorithm: Trapezoid Rule
    10. 5.10 Algorithm: Simpson’s Rule
    11. 5.11 Integration Error (Assessment)
    12. 5.12 Algorithm: Gaussian Quadrature
    13. 5.13 Higher Order Rules (Algorithm)
    14. 5.14 Monte Carlo Integration by Stone Throwing (Problem)
    15. 5.15 Mean Value Integration (Theory and Math)
    16. 5.16 Integration Exercises
    17. 5.17 Multidimensional Monte Carlo Integration (Problem)
    18. 5.18 Integrating Rapidly Varying Functions (Problem)
    19. 5.19 Variance Reduction (Method)
    20. 5.20 Importance Sampling (Method)
    21. 5.21 von Neumann Rejection (Method)
    22. 5.22 Nonuniform Assessment
14. 6: Matrix Computing
    1. 6.1 Problem 3: N–D Newton–Raphson; Two Masses on a String
    2. 6.2 Why Matrix Computing?
    3. 6.3 Classes of Matrix Problems (Math)
    4. 6.4 Python Lists as Arrays
    5. 6.5 Numerical Python (NumPy) Arrays
    6. 6.6 Exercise: Testing Matrix Programs
15. 7: Trial-and-Error Searching and Data Fitting
    1. 7.1 Problem 1: A Search for Quantum States in a Box
    2. 7.2 Algorithm: Trial-and-Error Roots via Bisection
    3. 7.3 Improved Algorithm: Newton–Raphson Searching
    4. 7.4 Problem 2: Temperature Dependence of Magnetization
    5. 7.5 Problem 3: Fitting An Experimental Spectrum
    6. 7.6 Problem 4: Fitting Exponential Decay
    7. 7.7 Least-Squares Fitting (Theory)
    8. 7.8 Exercises: Fitting Exponential Decay, Heat Flow and Hubble’s Law
16. 8: Solving Differential Equations: Nonlinear Oscillations
    1. 8.1 Free Nonlinear Oscillations
    2. 8.2 Nonlinear Oscillators (Models)
    3. 8.3 Types of Differential Equations (Math)
    4. 8.4 Dynamic Form for ODEs (Theory)
    5. 8.5 ODE Algorithms
    6. 8.6 Runge–Kutta Rule
    7. 8.7 Adams–Bashforth–Moulton Predictor–Corrector Rule
    8. 8.8 Solution for Nonlinear Oscillations (Assessment)
    9. 8.9 Extensions: Nonlinear Resonances, Beats, Friction
    10. 8.10 Extension: Time-Dependent Forces
17. 9: ODE Applications: Eigenvalues, Scattering, and Projectiles
    1. 9.1 Problem: Quantum Eigenvalues in Arbitrary Potential
    2. 9.2 Algorithms: Eigenvalues via ODE Solver + Search
    3. 9.3 Explorations
    4. 9.4 Problem: Classical Chaotic Scattering
    5. 9.5 Problem: Balls Falling Out of the Sky
    6. 9.6 Theory: Projectile Motion with Drag
    7. 9.7 Exercises: 2- and 3-Body Planet Orbits and Chaotic Weather
18. 10: High-Performance Hardware and Parallel Computers
    1. 10.1 High-Performance Computers
    2. 10.2 Memory Hierarchy
    3. 10.3 The Central Processing Unit
    4. 10.4 CPU Design: Reduced Instruction Set Processors
    5. 10.5 CPU Design: Multiple-Core Processors
    6. 10.6 CPU Design: Vector Processors
    7. 10.7 Introduction to Parallel Computing
    8. 10.8 Parallel Semantics (Theory)
    9. 10.9 Distributed Memory Programming
    10. 10.10 Parallel Performance
    11. 10.11 Parallelization Strategies
    12. 10.12 Practical Aspects of MIMD Message Passing
    13. 10.13 Scalability
    14. 10.14 Data Parallelism and Domain Decomposition
    15. 10.15 Example: The IBM Blue Gene Supercomputers
    16. 10.16 Exascale Computing via Multinode-Multicore GPUs
19. 11: Applied HPC: Optimization, Tuning, and GPU Programming
    1. 11.1 General Program Optimization
    2. 11.2 Optimized Matrix Programming with NumPy
    3. 11.3 Empirical Performance of Hardware
    4. 11.4 Programming for the Data Cache (Method)
    5. 11.5 Graphical Processing Units for High Performance Computing
    6. 11.6 Practical Tips for Multicore and GPU Programming
20. 12: Fourier Analysis: Signals and Filters
    1. 12.1 Fourier Analysis of Nonlinear Oscillations
    2. 12.2 Fourier Series (Math)
    3. 12.3 Exercise: Summation of Fourier Series
    4. 12.4 Fourier Transforms (Theory)
    5. 12.5 The Discrete Fourier Transform
    6. 12.6 Filtering Noisy Signals
    7. 12.7 Noise Reduction via Autocorrelation (Theory)
    8. 12.8 Filtering with Transforms (Theory)
    9. 12.9 The Fast Fourier Transform Algorithm
    10. 12.10 FFT Implementation
    11. 12.11 FFT Assessment
21. 13: Wavelet and Principal Components Analyses: Nonstationary Signals and Data Compression
    1. 13.1 Problem: Spectral Analysis of Nonstationary Signals
    2. 13.2 Wavelet Basics
    3. 13.3 Wave Packets and Uncertainty Principle (Theory)
    4. 13.4 Short-Time Fourier Transforms (Math)
    5. 13.5 The Wavelet Transform
    6. 13.6 Discrete Wavelet Transforms, Multiresolution Analysis
    7. 13.7 Principal Components Analysis
22. 14: Nonlinear Population Dynamics
    1. 14.1 Bug Population Dynamics
    2. 14.2 The Logistic Map (Model)
    3. 14.3 Properties of Nonlinear Maps (Theory and Exercise)
    4. 14.4 Mapping Implementation
    5. 14.5 Bifurcation Diagram (Assessment)
    6. 14.6 Logistic Map Random Numbers (Exploration)
    7. 14.7 Other Maps (Exploration)
    8. 14.8 Signals of Chaos: Lyapunov Coefficient and Shannon Entropy
    9. 14.9 Coupled Predator–Prey Models
    10. 14.10 Lotka–Volterra Model
    11. 14.11 Predator–Prey Chaos
23. 15: Continuous Nonlinear Dynamics
    1. 15.1 Chaotic Pendulum
    2. 15.2 Visualization: Phase-Space Orbits
    3. 15.3 Exploration: Bifurcations of Chaotic Pendulums
    4. 15.4 Alternate Problem: The Double Pendulum
    5. 15.5 Assessment: Fourier/Wavelet Analysis of Chaos
    6. 15.6 Exploration: Alternate Phase-Space Plots
    7. 15.7 Further Explorations
24. 16: Fractals and Statistical Growth Models
    1. 16.1 Fractional Dimension (Math)
    2. 16.2 The Sierpiński Gasket (Problem 1)
    3. 16.3 Growing Plants (Problem 2)
    4. 16.4 Ballistic Deposition (Problem 3)
    5. 16.5 Length of British Coastline (Problem 4)
    6. 16.6 Correlated Growth, Forests, Films (Problem 5)
    7. 16.7 Globular Cluster (Problem 6)
    8. 16.8 Fractals in Bifurcation Plot (Problem 7)
    9. 16.9 Fractals from Cellular Automata
    10. 16.10 Perlin Noise Adds Realism
    11. 16.11 Exercises
25. 17: Thermodynamic Simulations and Feynman Path Integrals
    1. 17.1 Magnets via Metropolis Algorithm
    2. 17.2 An Ising Chain (Model)
    3. 17.3 Statistical Mechanics (Theory)
    4. 17.4 Metropolis Algorithm
    5. 17.5 Magnets via Wang–Landau Sampling
    6. 17.6 Wang–Landau Algorithm
    7. 17.7 Feynman Path Integral Quantum Mechanics
    8. 17.8 Feynman’s Space–Time Propagation (Theory)
    9. 17.9 Exploration: Quantum Bouncer’s Paths
26. 18: Molecular Dynamics Simulations
    1. 18.1 Molecular Dynamics (Theory)
    2. 18.2 Verlet and Velocity–Verlet Algorithms
    3. 18.3 1D Implementation and Exercise
    4. 18.4 Analysis
27. 19: PDE Review and Electrostatics via Finite Differences and Electrostatics via Finite Differences
    1. 19.1 PDE Generalities
    2. 19.2 Electrostatic Potentials
    3. 19.3 Fourier Series Solution of a PDE
    4. 19.4 Finite-Difference Algorithm
    5. 19.5 Assessment via Surface Plot
    6. 19.6 Alternate Capacitor Problems
    7. 19.7 Implementation and Assessment
    8. 19.8 Electric Field Visualization (Exploration)
    9. 19.9 Review Exercise
28. 20: Heat Flow via Time Stepping
    1. 20.1 Heat Flow via Time-Stepping (Leapfrog)
    2. 20.2 The Parabolic Heat Equation (Theory)
    3. 20.3 Assessment and Visualization
    4. 20.4 Improved Heat Flow: Crank–Nicolson Method
29. 21: Wave Equations I: Strings and Membranes
    1. 21.1 A Vibrating String
    2. 21.2 The Hyperbolic Wave Equation (Theory)
    3. 21.3 Strings with Friction (Extension)
    4. 21.4 Strings with Variable Tension and Density
    5. 21.5 Vibrating Membrane (2D Waves)
    6. 21.6 Analytical Solution
    7. 21.7 Numerical Solution for 2D Waves
30. 22: Wave Equations II: Quantum Packets and Electromagnetic
    1. 22.1 Quantum Wave Packets
    2. 22.2 Time-Dependent Schrödinger Equation (Theory)
    3. 22.3 Algorithm for the 2D Schrödinger Equation
    4. 22.4 Wave Packet–Wave Packet Scattering
    5. 22.5 E&M Waves via Finite-Difference Time Domain
    6. 22.6 Maxwell’s Equations
    7. 22.7 FDTD Algorithm
    8. 22.8 Application: Wave Plates
    9. 22.9 Algorithm
    10. 22.10 FDTD Exercise and Assessment
31. 23: Electrostatics via Finite Elements
    1. 23.1 Finite-Element Method
    2. 23.2 Electric Field from Charge Density (Problem)
    3. 23.3 Analytic Solution
    4. 23.4 Finite-Element (Not Difference) Methods, 1D
    5. 23.5 1D FEM Implementation and Exercises
    6. 23.6 Extension to 2D Finite Elements
32. 24: Shocks Waves and Solitons
    1. 24.1 Shocks and Solitons in Shallow Water
    2. 24.2 Theory: Continuity and Advection Equations
    3. 24.3 Theory: Shock Waves via Burgers’Equation
    4. 24.4 Including Dispersion
    5. 24.5 Shallow-Water Solitons: The KdeV Equation
    6. 24.6 Solitons on Pendulum Chain
33. 25: Fluid Dynamics
    1. 25.1 River Hydrodynamics
    2. 25.2 Navier–Stokes Equation (Theory)
    3. 25.3 2D Flow over a Beam
    4. 25.4 Theory: Vorticity Form of Navier–Stokes Equation
34. 26: Integral Equations of Quantum Mechanics
    1. 26.1 Bound States of Nonlocal Potentials
    2. 26.2 Momentum–Space Schrödinger Equation (Theory)
    3. 26.3 Scattering States of Nonlocal Potentials
    4. 26.4 Lippmann–Schwinger Equation (Theory)
35. Appendix A: Codes, Applets, and Animations
36. Bibliography
37. Index
38. End User License Agreement

Guide

1. Cover
2. Table of Contents
3. Begin Reading

List of Illustrations

1. 1 Introduction
   1. Figure 1.1 (a) A representation of the multi-disciplinary nature of computational physics as an overlap of physics, applied mathematics and computer science, and as a bridge among them. (b) Simulation has been added to experiment and theory as a basic approach in the search for scientific truth. Although this book focuses on simulation, we present it as part of the scientific process.
   2. Figure 1.2 Screen dumps of two x–y plots produced by EasyVisual.py using the Visual package. Plot (a) uses default parameters, while plot (b) uses user-supplied options.
   3. Figure 1.3 (a) Output from the program 3GraphVisual.py that places three different types of 2D plots on one graph using Visual. (b) Three frames from a Visual animation of a quantum mechanical wave packet produced with HarmosAnimate.py.
   4. Figure 1.4 Some 3D shapes created with single commands in VPython.
   5. Figure 1.5 Matplotlib plots. (a) Output of EasyMatPlot.py showing a simple, x–y plot. (b) Output from GradesMatPlot.py that places two sets of data points, two curves, and unequal upper and lower error bars, all on one plot.
   6. Figure 1.6 (a,b) Columns show two separate outputs, each of two figures, produced by MatPlot2figs.py. (We used the slider button to add some space between the plots.)
   7. Figure 1.7 (a) A 3D wire frame. (b) A surface plot with wire frame. Both are produced by the program Simple3dplot.py using Matplotlib.
   8. Figure 1.8 A 3D scatter plot produced by the program Scatter3dPlot.py using Matplotlib.
   9. Figure 1.9 A comparison of a saw tooth function to the sum of its Fourier components, with the number of included waves varied interactively by a Matplotlib slider. FourierMatplot.py, which produced this output, was written by Oscar Restrepo.
   10. Figure 1.10 (a) A Mayavi surface plot of the function as seen in the screen viewer. (b) A rotatable visualization of a spherical harmonic in which the radial distance represents the value of the function.
   11. Figure 1.11 (a) A Mayavi visualization in which tubes are used to connect a set of data points. (b) A Mayavi visualization using arrows (glyphs) to represent a vector field.
   12. Figure 1.12 (a) A Mayavi visualization of the contours of the scalar field . (b) A volume rendering of the same scalar field.
   13. Figure 1.13 (a) A Mayavi visualization of the same scalar field using cut planes. (b) A visualization of the same scalar field combining cut planes and contours.
2. 2 Computing Software Basics
   1. Figure 2.1 A schematic view of a computer’s kernel and shells. The hardware is in the center surrounded by increasing higher level software.
   2. Figure 2.2 A flowchart illustrating a program to compute projectile motion. Plot (a) shows the basic components of the program, and plot (b) shows are some of its details. When writing a program, first map out the basic components, then decide upon the structures, and finally fill in the details. This is called top-down programming.
   3. Figure 2.3 The limits of single-precision floating-point numbers and the consequences of exceeding these limits (not to scale). The hash marks represent the values of numbers that can be stored; storing a number in between these values leads to truncation error. The shaded areas correspond to over- and underflow.
3. 3 Errors and Uncertainties in Computations
   1. Figure 3.1 A schematic of the N steps in a random walk simulation that end up a distance from the origin. Notice how the Δx's for each step add vectorially.
   2. Figure 3.2 The first four spherical Bessel functions j_l(x) as functions of x. Notice that for small x, the values for increasing l become progressively smaller.
   3. Figure 3.3 A log–log plot of relative error vs. the number of points used for a numerical integration. The ordinate value of ∼ 10⁻¹⁴ at the minimum indicates that ∼ 14 decimal places of precision are obtained before round-off error begins to build up. Notice that while the round-off error does fluctuate indicating a statistical aspect of error accumulation, on average it is increasing but more slowly than did the algorithm's error decrease.
   4. Figure 3.4 The error in the summation of the series for e^−x vs. N for various x values. The values of x increase vertically for each curve. Note that a negative initial slope corresponds to decreasing error with N, and that the dip corresponds to a rapid convergence followed by a rapid increase in error. (Courtesy of J. Wiren.)
4. 4 Monte Carlo: Randomness, Walks, and Decays
   1. Figure 4.1 (a) A plot of successive random numbers (x, y) = (r_i, r_i+1) generated with a deliberately “bad” generator, (b) A plot generated with the built in random number generator. While the plot (b) is not the proof that the distribution is random, the plot (a) is a proof enough that the distribution is not random.
   2. Figure 4.2 A plot of a uniform pseudorandom sequence r_i, vs. i. The points are connected to make it easier to follow the order. While this does not prove that a distribution is random, it at least shows the range of values and that there is fluctuation.
   3. Figure 4.3 (a) A schematic of the N steps in a random walk simulation that end up a distance R from the origin. Note how the Δx’s for each step add vectorially. (b) A simulated walk in 3D from Walk3D.py.
   4. Figure 4.4 (a) The steps taken in seven 2D random walk simulations. (b) The distance covered in two walks of N steps using different schemes for including randomness. The theoretical prediction (4.20) is the straight line.
   5. Figure 4.5 Two self-avoiding random walks that simulate protein chains with hydrophobic (H) monomers in light gray, and polar (P) monomers in black. The dark dots indicate regions where two H monomers are not directly connected.
   6. Figure 4.6 Circle: A sample containing N nuclei, each of which has the same probability of decaying per unit time, Graphs: Semilog plots of the number of nuclei vs. time for five simulations with differing initial numbers of nuclei. Exponential decay would be a straight line with bumps, similar to the initial behavior for N = 100 000.
5. 5 Differentiation and Integration
   1. Figure 5.1 A trajectory of a projectile experiencing air resistance. Forward-difference approximation (slanted dashed line) and central-difference approximation (horizontal line) for the numerical first derivative at time t. (A tangent to the curve at t yields the correct derivative.) The central difference is seen to be more accurate.
   2. Figure 5.2 The integral dx is the area under the graph of f(x) from a to b. Here we break up the area into four regions of equal widths h and five integration points.
   3. Figure 5.3 Different shapes used to approximate the areas under the curve. (a) Straight-line sections used for the trapezoid rule. (b) Two parabolas used in Simpson’s rule.
   4. Figure 5.4 Log–log plots of the error in the integration of exponential decay using the trapezoid rule, Simpson’s rule, and Gaussian quadrature vs. the number of integration points N. Approximately 15 decimal places of precision are attainable with double precision (a), and seven places with single precision (b). The algorithms are seen to stop converging when round-off error (the fluctuating and increasing part near the bottom) starts to dominate.
   5. Figure 5.5 (a) Stones into a pond as a technique for measuring its area. The ratio of “hits” to total number of stones thrown equals the ratio of the area of the pond to that of the box. (b) The evaluation of an integral via a Monte Carlo (stone throwing) technique of the ratio of areas.
   6. Figure 5.6 The area under the curve f(x) is the same as that under the horizontal line whose height .
   7. Figure 5.7 A ring of radius a carries a current l. Find vector potential at point P.
   8. Figure 5.8 The von Neumann rejection technique for generating random points with weight W(x). A random point is accepted if it lies below the curve of W(x) and rejected if it lies above. This generates a random distribution weighted by whatever W(x) function is plotted.
6. 6 Matrix Computing
   1. Figure 6.1 Two weights connected by three pieces of string and suspended from a horizontal bar of length L. The lengths are all known, but the angles and the tensions in the strings are to be determined.
   2. Figure 6.2 A free body diagram for one weight in equilibrium. Balancing the forces in the x and y directions for all weights leads to the equations of static equilibrium.
   3. Figure 6.3 (a) Row-major order used for matrix storage in Python, C and Java. (b) Column-major order used for matrix storage in Fortran. The table at the bottom shows how successive matrix elements are actually stored in a linear fashion in memory.
7. 7 Trial-and-Error Searching and Data Fitting
   1. Figure 7.1 A graphical representation of the steps involved in solving for a zero of f(x) using the bisection algorithm. The bisection algorithm takes the midpoint of the interval as the new guess for x, and so each step reduces the interval size by one-half. Four steps are shown for the algorithm.
   2. Figure 7.2 A graphical representation of the steps involved in solving for a zero of f(x) using the Newton–Raphson method. The Newton–Raphson method takes the new guess as the zero of the line tangent to f(x) at the old guess. Two guesses are shown.
   3. Figure 7.3 Two examples of how the Newton–Raphson algorithm may fail if the initial guess is not in the region where f(x) can be approximated by a straight line. (a) A guess lands at a local minimum/maximum, that is, a place where the derivative vanishes, and so the next guess ends up at x = ∞. (b)The search has fallen into an infinite loop. The technique know as “backtracking” could eliminate this problem.
   4. Figure 7.4 A function of the reduced magnetism m at three reduced temperatures t. A zero of this function determines the value of the magnetism at a particular value of t.
   5. Figure 7.5 Three fits to data. Dashed: Lagrange interpolation using an eight-degree polynomial; Short dashes: cubic splines fit ; Long dashed: Least-squares parabola fit.
   6. Figure 7.6 A reproduction of the experimental measurement by Stetz et al. (1973) giving the number of decays of a π meson as a function of time since its creation. Measurements were made during time intervals (box sizes) of 10-ns width. The dashed curve is the result of a linear least-squares fit to the log N(t).
   7. Figure 7.7 A linear least-squares best fit of a straight line to data. The deviation of theory from experiment is greater than would be expected from statistics, which means that a straight line is not a good theory to describe these data.
   8. Figure 7.8 A linear least-squares best fit of a parabola to data. Here we see that the fit misses approximately one-third of the points, as expected from the statistics for a good fit.
8. 8 Solving Differential Equations: Nonlinear Oscillations
   1. Figure 8.1 A mass m (the block) attached to a spring with restoring force F_k(x) driven by an external time-dependent driving force (the hand).
   2. Figure 8.2 (a) The potentials of an harmonic oscillator (solid curve) and of an anharmonic oscillator (dashed curve). If the amplitude becomes too large for the anharmonic oscillator, the motion becomes unbound. (b) The shapes of the potential energy function V(x) ∝ |x|^p for p = 2 and p = 6. The “linear” and “nonlinear” labels refer to the restoring force derived from these potentials.
   3. Figure 8.3 A sequence of uniform steps of length h taken in solving a differential equation. The solution starts at time t = 0 and is integrated in steps of h until t = T.
   4. Figure 8.4 Euler’s algorithm for integration of a differential equation one step forward in time. This linear extrapolation with the slope evaluated at the initial point is seen to lead to an error Δ.
   5. Figure 8.5 The rk2 algorithm for integration of a differential equation uses a slope (bold line segment) evaluated at the interval’s midpoint, and is seen to lead to a smaller error than Euler’s algorithm in Figure 8.4.
   6. Figure 8.6 The logarithm of the relative error in the solution of an ODE obtained with rk4 using a differing number N of time steps over a fixed time interval. The logarithm approximately equals the negative of the number of places of precision. Increasing the number of steps used for a fixed interval is seen to lead to smaller error.
   7. Figure 8.7 The position vs. time for oscillations within the potential for four different initial amplitudes. Each is seen to have a different period.
9. 9 ODE Applications: Eigenvalues, Scattering, and Projectiles
   1. Figure 9.1 Computed wave function and the square-well potential (bold lines). The wave function computed by integration in from the left is matched to the one computed by integration in from the right (dashed curve) at a point near the right edge of the well. Note how the wave function decays rapidly outside the well.
   2. Figure 9.2 Two guesses for the energy that are either too low or too high to be an eigenvalue. We see that the low-E guess does not oscillate enough to match onto a dying exponential, and that the high-E guess oscillates too much.
   3. Figure 9.3 A photograph of a pinball machine in which multiple scatterings occur from the bumpers.
   4. Figure 9.4 Scattering from the potential which in some ways models a pinball machine. The incident velocity v is in the y direction, and the impact parameter (x value) is b. The velocity of the scattered particle is v’ and its scattering angle is θ.
   5. Figure 9.5 The trajectories of a projectile fired with initial velocity V₀ in the θ direction. The lower curve includes air resistance.
   6. Figure 9.6 (a) The gravitational force on a planet at distance r from the sun. The x and y components of the force are indicated. (b) Output from the applet Planet showing the precession of a planet’s orbit when the gravitational force ∝ 1/r⁴ (successive orbits do not lie on top of each other).
   7. Figure 9.7 A snapshot from the animated output of the code UranusNeptune.py (on Instructor’s disk) showing: (a) the orbits of Uranus (inner circle) and of Neptune (outer circle) with the Sun in the center. The arrows indicate the Uranus–Neptune force that causes a perturbation in the orbits. (b) The perturbation in the angular position of Uranus as a result of the presence of Neptune.
10. 10 High-Performance Hardware and Parallel Computers
    1. Figure 10.1 The logical arrangement of the CPU and memory showing a Fortran array A(N) and matrix M(N, N) loaded into memory.
    2. Figure 10.2 Typical memory hierarchy for a single-processor, high-performance computer (B = bytes, k, M, G, T = kilo, mega, giga, tera).
    3. Figure 10.3 The elements of a computer’s memory architecture in the process of handling matrix storage.
    4. Figure 10.4 (a) Multitasking of four programs in memory at one time. On a SISD computer the programs are executed in round robin order. (b) Four programs in the four separate memories of a MIMD computer.
    5. Figure 10.5 (a) A generic view of the Intel core-2 dual-core processor, with CPU-local level-1 caches and a shared, on-die level-2 cache (courtesy of D. Schmitz), (b) The AMD Athlon 64 X2 3600 dual-core CPU (Wikimedia Commons).
    6. Figure 10.6 Two views of parallel computing (courtesy of Yuefan Deng).
    7. Figure 10.7 The theoretical maximum speedup of a program as a function of the fraction of the program that potentially may be run in parallel. The different curves correspond to different numbers of processors.
    8. Figure 10.8 A master process and two slave processes passing messages. Notice how in this not-well-designed program there are more sends than receives, and consequently the results may depend upon the order of execution, or the program may even lock up.
    9. Figure 10.9 A graphical representation of weak vs. strong scaling. Weak scaling keeps each node doing the same amount of work as the problem is made bigger. Strong scaling has each node doing less work (running for less time) as the number of nodes is made bigger.
    10. Figure 10.10 A representation of the lattice in a 2D rectangular space upon which Laplace’s equation is solved using a finite difference approach. The lattice sites with white centers correspond to the boundary of the physical system, upon which boundary conditions must be imposed for a unique solution. The large circles in the middle of the lattice represent the algorithm used to solve Laplace’s equation in which the potential at the point (x, y) = (i, j)Δ is set to the average of the potential values at the four nearest-neighbor points.
    11. Figure 10.11 The building blocks of Blue Gene (from Gara et al. (2005)).
    12. Figure 10.12 The single-node memory system (as presented by Gara et al. (2005)).
    13. Figure 10.13 (a) A 3D torus connecting 2 × 2 × 2 compute nodes. (b) The global collective memory system. (c) The control and Gb-ethernet memory system (from Gara et al. (2005)).
    14. Figure 10.14 A schematic of an exascale computer in which, in addition to each chip having multiple cores, a graphical processing unit is attached to each chip (adapted from Dongarra (2011)).
11. 11 Applied HPC: Optimization, Tuning, and GPU Programming
    1. Figure 11.1 Running time vs. dimension for an eigenvalue search using Tune.java and tune.f90.
    2. Figure 11.2 The cache manager’s view of RAM. Each 128-B cache line is read into one of four lines in cache.
    3. Figure 11.3 A schematic of a GPU (b) and its interactions with a CPU (a). The GPU is seen to have streaming multiprocessors, special function units and several levels of memory. Also shown is the communication between memory levels and with the host CPU.
    4. Figure 11.4 A schematic of how the contiguous elements of a matrix must be tiled for parallel processing (from Dongarra, 2011).
    5. Figure 11.5 (a) A 1 × 1 grid with a 1D block containing 32 threads, and the corresponding threadIdx.x that goes from 0 to 31. (b) A 1D 4 × 1 grid, that has 4 blocks with 8 threads each.
12. 12 Fourier Analysis: Signals and Filters
    1. Figure 12.1 (a) A sawtooth function that repeats infinitely in time. (b) The Fourier spectrum of frequencies contained in this function (natural units). Also see Figure 1.9 in which we show the result of summing a finite number of terms of this series.
    2. Figure 12.2 A plot of the functions sin(πt/2) and sin(2πt). If the sampling rate is not high enough, these signalsmay appear indistinguishable in a Fourier decomposition. If the sample rate is too low and if both signals are present in a sample, the deduced low-frequency component may be contaminated by the higher frequency component signal.
    3. Figure 12.3 (a) A function that is a signal plus noise s(t) + n(t); (b) the autocorrelation function vs. time deduced by processing this signal; (c) the power spectrum obtained from autocorrelation function; (d) the signal plus noise after passage through a low-pass filter.
    4. Figure 12.4 A schematic of an input signal f(t) passing through a filter h that outputs the function g(t).
    5. Figure 12.5 (a) An RC circuit arranged as a low-pass filter. (b) An RC circuit arranged as a high-pass filter.
    6. Figure 12.6 A delay-line filter in which the signal at different times is scaled by different amounts c_i.
    7. Figure 12.7 The rectangle function rect(ω) that is constant for a finite frequency interval. The Fourier transform of this function is sinc(t).
    8. Figure 12.8 The basic butterfly operation in which elements y_p and y_q on the left are transformed into y_p + Zy_q and y_p – Zy_q on the right.
    9. Figure 12.9 The butterfly operations performing a FFT on the eight data on the left leading to eight transforms on the right. The transforms are different linear combinations of the input data.
    10. Figure 12.10 A modified FFT in which the eight input data on the left are transformed into eight transforms on the right. The results are the same as in the previous figure, but now the output transforms are in numerical order whereas in the previous figure the input signals were in numerical order.
13. 13 Wavelet and Principal Components Analyses: Nonstationary Signals and Data Compression
    1. Figure 13.1 The input time signal (13.1) we wish to analyze. The signal is seen to contain additional frequencies as time increases. The boxes are possible placements of windows for short-time Fourier transforms.
    2. Figure 13.2 Four possible mother wavelets that can be used to generate entire sets of daughter wavelets. (a–d) Morlet (real part), Mexican hat, Daub4 e6 (explained later), and Haar. The daughter wavelets are generated by scaling and translating these mother wavelets.
    3. Figure 13.3 (a) A wave packet in time corresponding to the functional form (13.4) with ω₀ = 5 and N = 6. (b) The Fourier transform in frequency of this same wave packet.
    4. Figure 13.4 Four wavelet basis functions (daughters) generated by scaling (s) and translating (τ) an oscillating Gaussian mother wavelet. (a–d) (s = 1, τ = 0), (s = 1/2, τ = 0), (s = 1, τ = 6), and (s = 2, τ = 60). Note how s < 1 is a wavelet with higher frequency, while s > 1 has a lower frequency than the s = 1 mother. Likewise, the τ = 6 wavelet is just a translated version of the τ = 0 one directly above it.
    5. Figure 13.5 A schematic representation of the steps followed in performing a wavelet transformation over all time displacements and scales. The upper signal is first analyzed by evaluating its overlap with a narrow wavelet at the signal’s beginning. This produces a coefficient that measures the similarity of the signal to the wavelet. The wavelet is successively shifted over the length of the signal and the overlaps are successively evaluated. After the entire signal is covered, the wavelet is expanded and the entire analysis is repeated.
    6. Figure 13.6 Comparison of an input and reconstituted signal (13.23) using Morlet wavelets (the curves overlap nearly perfectly). As expected for Fourier transforms, the reconstruction is least accurate near the endpoints.
    7. Figure 13.7 The continuous wavelet spectrum obtained by analyzing the input signal with Morelet wavelets. Observe how at small values of time τ there is predominantly one frequency present, how a second, higher frequency (smaller scale) component enters at intermediate times, and how at larger times a still higher frequency components enter (figure courtesy of Z. Dimcovic).
    8. Figure 13.8 A graphical representation of the relation between time and frequency resolutions (the uncertainty relation). Each box represents an equal portion of the time-frequency plane but with different proportions of time and frequency.
    9. Figure 13.9 A multifrequency dyadic (power-of-2) filter tree used for discrete wavelet transformations. The L boxes represent low-pass filters and the H boxes represent high-pass filters. Each filter performs a convolution (transform). The circles containing “↓ 2” filter out half of the signal that enters them, which is called subsampling or factor-of-2 decimation. The signal on the left yields a transform with a single low and two high components (less information is needed about the low components for a faithful reproduction).
    10. Figure 13.10 An input signal at the top is processed by a tree of high- and low-band filters. The outputs from each filtering are down-sampled with half the data kept. The process continues until there are only two data of high-band filtering and two data of low-band filtering. When complete, the total number of output data equals the total number of signal data. The process of wavelet analysis is thus equivalent to a series of filterings.
    11. Figure 13.11 In successive passes, the filtering of the original signal at the top goes through the pyramid algorithm and produces the outputs shown. The sampling is reduced by a factor of 2 in each step. Note that in the upper graphs, we have connected the points to emphasize their continuous nature while in the lower graphs we plot the individual output points as histograms.
    12. Figure 13.12 (a) The Daub4 e6 wavelet constructed by inverse transformation of the wavelet coefficients. This wavelet has been found to be particularly effective in wavelet analyses. (b) The sum of Daub4 e10 and Daub4 1e58 wavelets of different scale and time displacements.
    13. Figure 13.13 (a) The normalized data and eigenvectors of covariance matrix. (b) The normalized data using the PCA eigenvectors as basis.
14. 14 Nonlinear Population Dynamics
    1. Figure 14.1 The insect population x_n vs. the generation number n for the two growth rates: (a) μ = 2.8, a single attractor; (b) μ = 3.3, a double attractor.
    2. Figure 14.2 The insect population x_n vs. the generation number n for two different growth rates: (a) μ = 3.5, a quadruple attractor; (b) μ = 3.8, a chaotic regime.
    3. Figure 14.3 Bifurcation plot of attractor population x^∗ vs. growth rate μ for the logistic map. The inset shows some details of a three-cycle window. (The colors or gray scales indicate the regimes over which we distributed the work on different CPUs when run in parallel.)
    4. Figure 14.4 Bifurcation plot for the Gaussian map f(x) = exp(−bx²) + μ (courtesy of W. Hager).
    5. Figure 14.5 Lyapunov coefficient (bottom) and bifurcation values (top) for the logistic map as functions of the growth rate μ. Note how the Lyapunov coefficient, whose value is a measure of chaos, changes abruptly at the bifurcations.
    6. Figure 14.6 Shannon entropy (top) and Lyapunov coefficient (bottom) for the logistic map. Note the close relation between the thermodynamic measure of disorder (entropy) and the nonlinear dynamics measure of chaos (Lyapunov).
    7. Figure 14.7 (a) The time dependencies of the populations of prey p(t) (solid curve) and of predator P(t) (dashed curve) fromthe Lotka–Volterra model. (b) A phase–space plot of Prey population p as a function of predator population P. The different orbits correspond to different initial conditions.
    8. Figure 14.8 A chaotic attractor for the 4D Lotka–Volterra model projected onto three axes.
    9. Figure 14.9 (a) The Lotka–Volterra model of prey population p(t) (solid curve) and predator population P(t) (dashed curve) vs. time when a prey population limit is included, (b) Prey population p as a function of predator population P.
    10. Figure 14.10 Lotka–Volterra model with predation efficiency and prey limitations. From left to right: overdamping, b = 0.01; damped oscillations, b = 0.1, and limit cycle, b = 0.3.
15. 15 Continuous Nonlinear Dynamics
    1. Figure 15.1 (a) A pendulum of length / driven through resistive air (dotted arcs) by an external sinusoidal torque (semicircle). The strength of the external torque is given by f and that of air resistance by a. (b) A double pendulum with neither air resistance nor a driving force. In both cases there is a gravitational torque.
    2. Figure 15.2 The data screen (a) and the output screen (b) of the applet HearData that converts data into sounds. Columns of [t_i, x(t_i)] data are pasted into the data window, processed into the graph in the output window, and then converted to sound data that are played by Java.
    3. Figure 15.3 Three potentials and their behaviors in phase space. The different orbits below the potentials correspond to different energies, as indicated by the limits of maximum displacements within the potentials (dashed lines). (a) A repulsive potential leads to open orbits characteristic of nonperiodic motion. The trajectories cross at the hyperbolic point in the middle, an unstable equilibrium point. (b) The harmonic oscillator leads to symmetric ellipses; the closed orbits indicate periodic behavior, and the symmetric trajectories indicate a symmetric potential. (c) A nonharmonic oscillator. Notice that the ellipse-like trajectories at the bottom are neither ellipses nor symmetric.
    4. Figure 15.4 (a) Phase-space trajectories for a pendulum including “over the top” motions. (b) The corresponding θ dependence of the potential.
    5. Figure 15.5 (a) Position vs. time for two initial conditions of a chaotic pendulum that end up with the same limit cycle. (b) A phase space plot of position versus velocity for the limit cycle shown in (a) (courtesy of W. Hager).
    6. Figure 15.6 Position vs. time, phase-space plot, and Fourier spectrum for a chaotic pendulum with ω₀ = 1, a = 0.2, f = 0.52, and ω = 0.666 and three different initial conditions. Column (a) displays three dominant cycles, the column (b) only one, while column (c) has multiple cycles. (Examples of chaotic behavior can be seen in Figure 15.7).
    7. Figure 15.7 Some examples of complicated behaviors of a realistic pendulum. From top to bottom in each column of the figure, the first three rows show θ(t) behavior, a broadband Fourier spectrum, a phase-space diagram containing regular patterns with dark bands, and a broad Fourier spectrum. These features are characteristic of chaos. In the bottom two rows we see how the behavior changes abruptly after a slight change in the magnitude of the force and that for f = 0.54 there occur the characteristic broad bands of chaos.
    8. Figure 15.8 (a,b) Output from the applet Pend1 that produces an animation of a chaotic pendulum, along with the corresponding position vs. time and phase-space plots. (b)The resulting Fourier spectrum produced by Pend1. (c,d) outputs from two runs of the applet Pend2 producing an animation of two chaotic pendula along with the corresponding phase-space plots, and the final output with limit cycles (dark bands).
    9. Figure 15.9 A bifurcation diagram for the damped pendulum with a vibrating pivot (see also the similar diagram for a double pendulum, Figure 15.11). The ordinate is | dθ/dt |, the absolute value of the instantaneous angular velocity at the beginning of the period of the driver, and the abscissa is the magnitude of the driving force f. Note that the heavy line results from the overlapping of points, not from connecting the points (see enlargement in the inset).
    10. Figure 15.10 Photographs of a double pendulum built by a student in the OSU Physics Department. The upper pendulum consists of two separated shafts so that the lower one can rotate completely around. Both pendula can go over their tops. The first two frames show the pendulum released from rest and then moving quickly. The photography with faster shutter speeds stops the motion in various stages (photograph by R. Landau).
    11. Figure 15.11 (a) Phase-space trajectories for a double pendulum with m₁ = 10m₂ and with two dominant attractors. (b) A bifurcation diagram for the double pendulum displaying the instantaneous velocity of the lower pendulum as a function of the mass of the upper pendulum (both plots courtesy of J. Danielson).
    12. Figure 15.12 A 3D plot of a Lorenz attractor output from the program LorentzAtract.py (on the Instructor’s disk).
16. 16 Fractals and Statistical Growth Models
    1. Figure 16.1 (a) A statistical fractal Sierpiński gasket containing 10 000 points. Note the self-similarity at different scales. (b) A geometric Sierpiński gasket constructed by successively connecting the midpoints of the sides of each equilateral triangle. The first three steps in the process are labeled as A, B, C.
    2. Figure 16.2 (a) A fractal fern generated by 30 000 iterations of the algorithm 16.14). Enlarging this fern shows that each frond with a similar structure. (b) A fractal tree created with the simple algorithm (16.17).
    3. Figure 16.3 A simulation of the ballistic deposition of 20 000 particles onto a substrate of length 200. The vertical height increases in proportion to the length of deposition time, with the top being the final surface.
    4. Figure 16.4 Examples of the use of “box” counting to determine fractal dimension. (a) In the top “boxes” are circles and the perimeter is being covered. In the bottom an entire figure is being covered. In (b) a “coastline” is being covered by boxes of two different sizes (scales). The fractal dimension can be deduced by recording the number of box of different scale needed to cover the figures.
    5. Figure 16.5 Fractal dimensions of a line, box, and coastline determined by box counting. The slope at vanishingly small scale determines the dimension.
    6. Figure 16.6 (a) A view that might be seen in the undergrowth of a forest or after a correlated ballistic deposition. (b) The probability of particle i + 1 sticking in one column depends upon the distance d from the previously deposited particle i.
    7. Figure 16.7 (a) A globular cluster of particles of the type that might occur in a colloid, (b) The applet Dla2en.html lets you watch these clusters grow. Here the cluster is at the center of the circle, and random walkers are started at random points around the circle.
    8. Figure 16.8 Number 8 by the American painter Jackson Pollock. (Used with permission, Neuberger Museum, State University of New York.) Some researchers claim that Pollock’s paintings exhibit a characteristic fractal structure, while some other researchers question this (Kennedy, 2006). See if you can determine the fractal dimensions within this painting.
    9. Figure 16.9 The rules for two versions of the Game of Life. The rules, given graphically on the top row, create the gaskets below (output obtained from the applet JCellAut in the auxiliary files).
    10. Figure 16.10 The coordinates used in adding Perlin noise. (a) The rectangular grid used to locate a square in space and a corresponding point within the square. As shown with the arrows, unit vectors g_i with random orientation are assigned at each grid point. (b) A point within each square is located by drawing the four p_i. The g_i vectors are the same as on the left.
    11. Figure 16.11 The mapping used in adding Perlin noise. (a) The numbers s, t, u, and v are represented by perpendiculars to the four vertices, with lengths proportional to their values. (b) The function 3p² − 2p³ is used as a map of the noise at a point like (x, y) to others close by.
    12. Figure 16.12 Perlin noise mapping, (a) The point (x,y) is mapped to point (s_x,x_y). (b) Using (16.33). Then three linear interpolations are performed to find c, the noise at (x, y).
    13. Figure 16.13 After the addition of Perlin noise, the random scatterplot in (a) becomes the clusters on (b).
    14. Figure 16.14 (a) The output from the Pov-Ray ray-tracing program that took as input the 2D coherent random noise plot in Figure 16.13 and added height and fog. (b) An image of a surface of revolution produced by Pov-Ray in which the marble-like texture is created by Perlin noise.
17. 17 Thermodynamic Simulations and Feynman Path Integrals
    1. Figure 17.1 The 1D lattice of N spins used in the Ising model of magnetism. The interaction energy between nearest-neighbor pairs E = ±J is shown for aligned and opposing spins.
    2. Figure 17.2 An Ising model simulation on a 1D lattice of 100 initially aligned spins (on the left). Up spins are indicated by circles, and down spins by blank spaces. Although the system starts with all up spins (a “cold” start), the system is seen to form domains of up and down spins as time progresses.
    3. Figure 17.3 Simulation results from a 1D Ising model of 100 spins. (a) Energy and specific heat as functions of temperature; (b) magnetization as a function of temperature.
    4. Figure 17.4 The energy, specific heat, and magnetization as a function of temperature from a 2D Ising model simulation with 40 000 spins. Evidence of a phase transition at the Curie temperature kT = 2.5 is seen in all three functions. The values of C and E have been scaled to fit on the same plot as M (courtesy of J. Wetzel).
    5. Figure 17.5 Wang–Landau sampling used in the Ising model 2D Ising model on an 8 × 8 lattice. (a) Logarithm of the density of states log g(E) ∝ S vs. the energy per particle. (b) The histogram H(E) showing the number of states visited as a function of the energy per particle. The aim of WLS is to make this function flat.
    6. Figure 17.6 The numbering scheme used in our WLS implementation of the 2D Ising model with an 8 × 8 lattice of spins.
    7. Figure 17.7 In the Feynman path-integral formulation of quantum mechanics a collection of paths connect the initial space–time point A to the final point B. The solid line is the trajectory followed by a classical particle, while the dashed lines are additional paths sampled by a quantum particle. A classical particle somehow “knows” ahead of time that travel along the classical trajectory minimizes the action S.
    8. Figure 17.8 (a) The probability distribution for the harmonic oscillator ground state as determined with a path-integral calculation (the classical result has maxima at the two turning points). (b) A space–time quantum path resulting from applying the Metropolis algorithm.
    9. Figure 17.9 (a) A path through a space–time lattice that starts and ends at x = x_a = x_b. The action is an integral over this path, while the path integral is a sum of integrals over all paths. The dotted path BD is a transposed replica of path AC. (b) The dashed path joins the initial and final times in two equal time steps; the solid curve uses N steps each of size ε. The position of the curve at time t_j defines the position x_j.
    10. Figure 17.10 The analytic and quantum Monte Carlo solution for the quantum bouncer. The continuous line is the Airy function squared and the dashed line |ψ₀(z)|² after a million trajectories.
18. 18 Molecular Dynamics Simulations
    1. Figure 18.1 The molecule–molecule effective interaction arises from the many-body interaction of the electrons and nucleus in one molecule (circle) with the electrons and nucleus in another molecule (other circle). Note, the size of the nucleus at the center of each molecule is highly exaggerated, and real electrons have no size.
    2. Figure 18.2 The Lennard–Jones effective potential used in many MD simulations. Note the sign change at r = 1 and the minimum at r 1.1225 (natural units). Also note that because the r-axis does not extend to r = 0, the infinitely high central repulsion is not shown.
    3. Figure 18.3 (a) Two frames from an animation showing the results of a 1D MD simulation that starts with uniformly spaced atoms. Note the unequal spacing resulting from an image atom moving in from the left after an atom left from the right. (b) Two frames from the animation of a 2D MD simulation showing the initial and an equilibrated state. Note how the atoms start off in a simple cubic arrangement but then equilibrate to a face-centered-cubic lattice. In both the cases, the atoms remain confined as a result of the interatomic forces.
    4. Figure 18.4 The infinite space generated by imposing periodic boundary conditions on the particles within the simulation volume (shaded box). The two-headed arrows indicate how a particle interacts with the nearest version of another particle, be that within the simulation volume or an image. The vertical arrows indicate how the image of particle 4 enters when the actual particle 4 exits.
    5. Figure 18.5 The kinetic, potential, and total energy as a function of time or number of steps for a 2D MD simulation with 36 particles (a), and 300 particles (b), both with an initial temperature of 150 K. The potential energy is negative, the kinetic energy is positive, and the total energy is seen to be conserved (flat).
    6. Figure 18.6 (a) The temperature after equilibration as a function of initial kinetic energy for a 2D MD simulation with 36 particles. The two are nearly proportional. (b) The pressure vs. temperature for a simulation with several hundred particles. An ideal gas (noninteracting particles) would yield a straight line (courtesy of J. Wetzel).
    7. Figure 18.7 A simulation of a projectile shot into a group of particles. The energy introduced by the projectile is seen to lead to evaporation of the particles (courtesy of J. Wetzel).
19. 19 PDE Review and Electrostatics via Finite Differences and Electrostatics via Finite Differences
    1. Figure 19.1 (a) The shaded region of space within a square in which we determine the electric potential by solving Laplace’s equation. There is a wire at the top kept at a constant 100 V and a grounded wire (dashed) at the sides and bottom. (b) The computed electric potential as a function of x and y. The projections onto the shaded xy plane are equipotential (contour) lines.
    2. Figure 19.2 The analytic (Fourier series) solution of Laplace’s equation summing 21 terms. Gibbs overshoot leads to the oscillations near x = 0, and persist even if a large number of terms is summed.
    3. Figure 19.3 The algorithm for Laplace’s equation in which the potential at the point (x, y) = (i, j)Δ equals the average of the potential values at the four nearest-neighbor points. The nodes with white centers correspond to fixed values of the potential along the boundaries.
    4. Figure 19.4 (a) A simple model of a parallel-plate capacitor within a box. A realistic model would have the plates close together, in order to condense the field, and the enclosing grounded box so large that it has no effect on the field near the capacitor. (b) A numerical solution for the electric potential for this geometry. The projection on the xy plane gives the equipotential lines.
    5. Figure 19.5 A guess as to how charge may rearrange itself on finite conducting plates.
    6. Figure 19.6 (a) A visualization of the computed electric potential for a capacitor with finite width plates. (b) A visualization of the charge distribution along one plate determined by evaluating ∇² U(x, y) (courtesy of J. Wetzel). Note the “lightening rod” effect of charge accumulating at corners and points.
    7. Figure 19.7 (a) The geometry of a capacitor formed by placing two long, square cylinders within each other. (b) The geometry of a capacitor formed by placing two long, circular cylinders within each other. The cylinders are cracked on the side so that wires can enter the region.
    8. Figure 19.8 (a) Computed equipotential surfaces and electric field lines for a realistic capacitor. (b) Equipotential surfaces and electric field lines mapped onto the surface for a 3D capacitor constructed from two tori.
20. 20 Heat Flow via Time Stepping
    1. Figure 20.1 A metallic bar insulated along its length with its ends in contact with ice. The bar is displayed in dark gray and the insulation is of lighter gray.
    2. Figure 20.2 The algorithm for the heat equation in which the temperature at the location x = iΔx and time t = (j + 1)Δt is computed from the temperature values at three points of an earlier time. The nodes with white centers correspond to known initial and boundary conditions (the boundaries are placed artificially close for illustrative purposes).
    3. Figure 20.3 A numerical calculation of the temperature vs. position and vs. time, with isotherm contours projected onto the horizontal plane.
    4. Figure 20.4 Temperature vs. position and time when two bars at differing temperatures are placed in contact at t = 0. The projected contours show the isotherms.
21. 21 Wave Equations I: Strings and Membranes
    1. Figure 21.1 (a) A stretched string of length L tied down at both ends and under high enough tension so that we can ignore gravity. The vertical disturbance of the string from its equilibrium position is y(x, t). (b) A differential element of the string showing how the string’s displacement leads to the restoring force.
    2. Figure 21.2 The solutions of the wave equation for four earlier space–time points are used to obtain the solution at the present time. The boundary and initial conditions are indicated by the white-centered dots.
    3. Figure 21.3 The vertical displacement as a function of position x and time t for a string initially placed in a standing mode when friction is included. Notice how the wave moves up and down with time (courtesy of J. Wiren).
    4. Figure 21.4 The vertical displacement as a function of position and time of a string initially plucked simultaneously at two points, as shown by arrows. Note that each initial peak breaks up into waves traveling to the right and to the left. The traveling waves invert on reflection from the fixed ends. As a consequence of these inversions, the t 12 wave is an inverted t = 0 wave.
    5. Figure 21.5 (a) A uniform string suspended from its ends in a gravitational field assumes a catenary shape. (b) A force diagram of a section of the catenary at its lowest point. Because the ends of the string must support the entire weight of the string, the tension now varies along the string.
    6. Figure 21.6 A small part of an oscillating membrane and the forces that act on it.
    7. Figure 21.7 The standing wave pattern on a square box top at three different times.
22. 22 Wave Equations II: Quantum Packets and Electromagnetic
    1. Figure 22.1 The position as a function of time of a localized electron confined to a square well (computed with the code SqWell.py available in the instructor’s site). The electron is initially on the left with a Gaussian wave packet. In time, the wave packet spreads out and collides with the walls.
    2. Figure 22.2 The probability density as a function of time for an electron confined to a 1D harmonic oscillator potential well. (a) A conventional surface plot, (b) a color visualization.
    3. Figure 22.3 The probability density as a function of position and time for an electron incident upon and passing through a slit. Significant reflection is seen to occur.
    4. Figure 22.4 The probability density as a function of x and y at three different times for an electron confined to a 2D parabolic tube (infinite in the y direction). The electron is initially placed in a Gaussian wave packet in both the x and y directions, and it is to be noted how there is spreading of the wave packet in the y direction, but not in the x direction.
    5. Figure 22.5 (a) Six frames from an animation of the two-particle density ρ(x₁, x₂, t)as a function of the position of particle 1 with mass m and of the position of particle 2 with mass 10m. (b) This same collision as seen with the single-particle densities ρ(x₁, t) and ρ(x₂, t). The numbers in the left-hand corners are the times in units of 100Δt. Note that each plot ends at the walls of the containing box, and that particle 1 “bounces off” a wall between times 36 and 86.
    6. Figure 22.6 A time sequence of two single-particle wave packets scattering from each other. The particles have equal mass, a mean kinetic energy equal to a quarter of the well’s depth, and the wavefunction has been symmetrized.
    7. Figure 22.7 A single EM pulse traveling along the z-axis. The coupled E and H pulses are indicated by solid and dashed curves, respectively, and the pulses at different z values correspond to different times.
    8. Figure 22.8 The algorithm for using the known values of E_x and H_y at three earlier times and three different space positions to obtain the solution at the present time. Note that the values of E_x are determined on the lattice of filled circles, corresponding to integer space indices and half-integer time indices. In contrast, the values of H_y are determined on the lattice of open circles, corresponding to half-integer space indices and integer time indices.
    9. Figure 22.9 The E field (light) and the H field (dark) at the initial time (a) and at a later time (b). Periodic boundary conditions are used at the ends of the spatial region, which means that the large z wave continues into the z = 0 wave.
    10. Figure 22.10 E and H fields at t = 100 for a circularly polarized wave in free space.
    11. Figure 22.11 One frame from the program quarterwave.py (on Instructor’s site) showing a linearly polarized EM wave entering a quarter-wave plate from the left and leaving as a circularly polarized wave on the right (the arrow on the left oscillates back and forth at 45° while the one on the right rotates).
23. 23 Electrostatics via Finite Elements
    1. Figure 23.1 A finite element solution to Laplace’s equation for two metal plates with a charge density between them. The dots are the nodes x_i, and the lines connecting the nodes are the finite elements.
    2. Figure 23.2 Basis functions used in finite-elements solution of the 1D Laplace equation. (a) A set of overlapping basis functions φ_i. Each function is a triangle from x_i−1 to x_i+1. (b) A Piecewise-linear function. (c) A piecewise-quadratic function.
    3. Figure 23.3 Exact (line) vs. FEM solution (points) for the two-plate problem for N = 3 and N = 11 finite elements (N = 3 displaced upwards for clarity). On this scale, the N = 11 solution IS identical to the exact one.
    4. Figure 23.4 (a) Decomposition of a domain into triangular elements. Smaller triangles are used in regions of rapid variation and larger triangles are used in regions of slow variation. Discretization errors occur at boundaries. (b) A decomposition of a rectangular domain into 32 right triangles on a mesh with 25 nodes.
    5. Figure 23.5 (a) Linear triangular elements in the x–y plane. (b) Linear triangular element (master element) in the η plane.
24. 24 Shocks Waves and Solitons
    1. Figure 24.1 (a) A re-creation on the Union Canal near Edinburgh of Russel’s soliton (www.ma.hw.ac.uk/solitons/press.html, Nature v. 376, 3 August 1995, p. 373). (b) Two shallow-water solitary waves crossing each other computed with the code Soliton.py. The taller soliton on the left catches up with and overtakes the shorter one at t 5. The waves resume their original shapes after the collision.
    2. Figure 24.2 The formation of a tsunami. A single two-level waveform at time zero progressively breaks up into eight solitons (labeled) as time increases. The tallest soliton (1) is narrower and faster in its motion to the right. You can generate an animation of this with the program SolitonAnimate.py.
    3. Figure 24.3 A visualization showing the wave height vs. position for increasing times showing the formation of a shock wave (sharp edge) from an initial sine wave.
    4. Figure 24.4 A 1D chain of pendulums coupled with a torsion bar on top. The pendulums swing in planes perpendicular to the length of the bar.
    5. Figure 24.5 The dispersion relation for a linearized chain of pendulums.
    6. Figure 24.6 A circular ring soliton at times 8, 20, 40, 60, 80, and 120. This has been proposed as a model for an elementary particle.
25. 25 Fluid Dynamics
    1. Figure 25.1 Side view of the flow of a stream around a submerged beam (a) and around two parallel plates (b). Both beam and plates have length L along the direction of flow. The flow is seen to be symmetric about the center-line and to be unaffected at the bottom and at the surface by the submerged object.
    2. Figure 25.2 The boundary conditions for two thin submerged plates. The surrounding box is the integration volume within which we solve the PDEs and upon whose surface we impose the boundary conditions. In practice, the box is much larger than L and H
    3. Figure 25.3 Boundary conditions for flow around the beam in Figure 25.1. The flow is symmetric about the centerline, and the beam has length L in the x direction (along flow).
    4. Figure 25.4 Two visualizations of the stream function u for Reynold’s number R = 5. The one on the left uses contours, while the one on the right uses color.
    5. Figure 25.5 (a) The velocity field around the beam as represented by vectors. (b) The vorticity as a function of x and y for the same velocity field. Rotation is seen to be largest behind the beam.
26. 26 Integral Equations of Quantum Mechanics
    1. Figure 26.1 A dark particle moving in a dense medium in which it interacts with all particles present. The nonlocality of the potential felt by the dark particle at r arises from the interactions at all r′.
    2. Figure 26.2 The grid of momentum values on which the integral equation is solved.
    3. Figure 26.3 (a) A projectile (dark particle at r) scattering from a dense medium. (b) The same process viewed in the CM system where the projectile and target always have equal and opposite momenta.
    4. Figure 26.4 Three different paths in the complex k′ plane used to evaluate line integrals when there are singularities. Here the singularities are at k and −k, and the integration variable is k′.
    5. Figure 26.5 The energy dependence of the cross section for l = 0 scattering from an attractive delta-shell potential with λb = 15. The dashed curve is the analytic solution (26.39), and the solid curve results from numerically solving the integral Schrödinger equation, either via direct matrix inversion or via LU decomposition.

List of Tables

1. 1 Introduction
   1. Table 1.1 Topics for one-quarter (10 Weeks) scientific computing course.
   2. Table 1.2 Topics for two-quarters (20 Weeks) computational physics course.
   3. Table 1.3 Some common Matplotlib commands.
2. 2 Computing Software Basics
   1. Table 2.1 The IEEE 754 Standard for Primitive Data Types.
   2. Table 2.2 Representation scheme for normal and abnormal IEEE singles.
   3. Table 2.3 Representation scheme for IEEE doubles.
3. 3 Errors and Uncertainties in Computations
   1. Table 3.1 Approximate values for spherical Bessel functions (from maple).
4. 4 Monte Carlo: Randomness, Walks, and Decays
   1. Table 4.1 A table of a uniform, pseudo-random sequence r_i generated by Python’s random method.
5. 5 Differentiation and Integration
   1. Table 5.1 Elementary weights for uniform-step integration rules.
   2. Table 5.2 Types of Gaussian integration rules.
   3. Table 5.3 Points and weights for 4-point Gaussian quadrature (for checking computation).
6. 6 Matrix Computing
   1. Table 6.1 The operators of NumPy and their effects.
7. 7 Trial-and-Error Searching and Data Fitting
   1. Table 7.1 Experimental values for a scattering cross section (f(E) in the theory), each with absolute error ±σ_i, as a function of energy (x_i in the theory).
   2. Table 7.2 Temperature vs. distance as measured along a metal rod.
   3. Table 7.3 Distance vs. radial velocity for 24 extragalactic nebulae.
8. 8 Solving Differential Equations: Nonlinear Oscillations
   1. Table 8.1 Comparison of ODE solvers for different equations.
9. 10 High-Performance Hardware and Parallel Computers
   1. Table 10.1 Computation of c = (a + b)/(d × f).
   2. Table 10.2 Computation of matrix [C] = [A] + [B]-
   3. Table 10.3 Vector processing of matrix [A] + [B] = [C].
   4. Table 10.4 A typical organization of a program containing both serial and parallel tasks.
10. 12 Fourier Analysis: Signals and Filters
    1. Table 12.1 Reordering for 16 data complex points.
11. 13 Wavelet and Principal Components Analyses: Nonstationary Signals and Data Compression
    1. Table 13.1 PCA data.
12. 18 Molecular Dynamics Simulations
    1. Table 18.1 Parameter and scales for the Lennard-Jones potential.
13. 19 PDE Review and Electrostatics via Finite Differences and Electrostatics via Finite Differences
    1. Table 19.1 Three categories of PDE based on the value of their discriminant d.
    2. Table 19.2 The relation between boundary conditions and uniqueness for PDEs.
14. Appendix A: Codes, Applets, and Animations
    1. Table A.1 Python codes, Chapters 1–26.
    2. Table A.2 Animations (A player such as VLC or QuickTime required formpeg and avi, and a Web browser for gifs).