computational_physics/cp.tex

   1 \documentclass[pdf]{prosper}
   2
   3 \usepackage{verbatim}
   4 \usepackage[german]{babel}
   5 \usepackage[latin1]{inputenc}
   6 \usepackage[T1]{fontenc}
   7 \usepackage{amsmath}
   8 \usepackage{ae}
   9 \usepackage{aecompl}
  10 \usepackage{color}
  11 \usepackage{graphicx}
  12 \graphicspath{{./img/}}
  13 \usepackage{hyperref}
  14
  15 \title{introduction to computational physics}
  16 \subtitle{basic concepts and approaches}
  17 \author{frank zirkelbach}
  18 \email{frank.zirkelbach@physik.uni-augsburg.de}
  19 \institution{experimantal physics {\footnotesize IV} - university of augsburg}
  20
  21 \begin{document}
  22
  23 \maketitle
  24
  25 \overlays{5}{
  26 \begin{slide}{outline}
  27 \begin{itemstep}
  28   \item motivation
  29   \item history of computing hardware/software
  30   \item warning
  31   \item computational techniques
  32   \item summary
  33 \end{itemstep}
  34 \end{slide}}
  35
  36 \overlays{3}{
  37 \begin{slide}{motivation}
  38 \FromSlide{1}{
  39   \begin{center}
  40   \begin{figure}[h]
  41   \includegraphics[width=8cm]{cp_appl_field.eps}
  42   \end{figure}
  43   \end{center}
  44 }
  45 \FromSlide{2}{
  46 challenge:
  47 \begin{itemize}
  48   \item precise mathematical theory
  49   \item often: solving theory's equations ab-initio is not realistic
  50   \item only a few models can be solved exactly
  51 \end{itemize}}
  52 \FromSlide{3}{
  53 $\Rightarrow$ study and implementation of numerical algorithms
  54 }
  55 \end{slide}}
  56
  57 \overlays{5}{
  58 \begin{slide}{history of computing hardware}
  59   \begin{minipage}[t]{10cm}
  60      \onlySlide*{1}{\begin{center} \includegraphics[height=3cm]{abacus.eps} \end{center}}
  61      \onlySlide*{2}{\begin{center} \includegraphics[height=3cm]{eniac.eps} \hspace{1cm} \includegraphics[height=3cm]{tube.eps} \end{center}}
  62      \onlySlide*{3}{\begin{center} \includegraphics[height=3cm]{z1.eps} \end{center}}
  63      \onlySlide*{4}{\begin{center} \includegraphics[height=3cm]{pdp1.eps} \hspace{1pt} \includegraphics[height=3cm]{transistor.eps} \end{center}}
  64      \onlySlide*{5}{\begin{center} \includegraphics[height=3cm]{pdp8.eps} \hspace{1cm} \includegraphics[height=3cm]{ic.eps} \end{center}}
  65      %\FromSlide{6}{\begin{center} \includegraphics[height=3cm]{} \end{center}}
  66   \end{minipage}
  67   \begin{minipage}[b]{10cm}
  68     \begin{itemstep}
  69       \item $3000 \, bc$: abacus - first calculating device
  70       \item $1945$: eniac - electrical digital computer
  71       \item $1938/41$: z1/3 - featuring memory and programmability
  72       \item $1960$: pdp-1 - transistor based computers
  73       \item $1964$: pdp-8 - integrated circuit computers
  74     \end{itemstep}
  75   \end{minipage}
  76 \end{slide}}
  77
  78 \overlays{6}{
  79 \begin{slide}{history of computing hardware}
  80   \begin{minipage}[t]{10cm}
  81      \onlySlide*{1}{\begin{center} \includegraphics[height=3cm]{4004.eps} \end{center}}
  82      \onlySlide*{2}{\begin{center} \includegraphics[height=3cm]{cray2.eps} \hspace{1pt} \includegraphics[height=3cm]{cray2_i.eps} \end{center}}
  83      \onlySlide*{3}{\begin{center} \includegraphics[height=3cm]{apple2.eps} \includegraphics[height=3cm]{c64.eps} \end{center}}
  84      \onlySlide*{4}{\begin{center} \includegraphics[height=3cm]{intel1.eps} \includegraphics[height=3cm]{intel2.eps} \end{center}}
  85      \onlySlide*{5}{\begin{center} \includegraphics[height=3cm]{mips.eps} \hspace{1pt} \includegraphics[height=3cm]{ppc.eps} \end{center}}
  86      \onlySlide*{6}{\begin{center} \includegraphics[height=3cm]{cluster1.eps} \hspace{1cm} \includegraphics[height=3cm]{cluster2.eps} \end{center}}
  87   \end{minipage}
  88   \begin{minipage}[b]{10cm}
  89     \begin{itemstep}
  90       \item $1970$: intel 4004 - first single chip $\mu$-processor
  91       \item $1977/85$: cray1/2 - vector supercomputer
  92       \item $1977/82/85$: 6502/6510/m68k - first pc
  93       \item $1978/82/85 $: 8086/80286/80386
  94       \item $1985$: mips - first risc design
  95       \item $1990/2000$: massive parallel computing
  96     \end{itemstep}
  97   \end{minipage}
  98 \end{slide}}
  99
 100 \overlays{11}{
 101 \begin{slide}{history of computing software}
 102   \begin{itemstep}
 103     \item $1946$: plankalk"ul - high-level programming language
 104     \item $1950$: assembler - translating instruction mnemonics
 105     \item $1954$: fortran - {\scriptsize formula translation}
 106     \item $1963$: basic - {\scriptsize beginner's all purpose symbolic instruction code}
 107     \item $1964$: os/360 - batch processing operating system
 108     \item $1969$: unix - multics port to pdp-8, pdp-11/20
 109     \item $1972$: c programming language - thompson, ritchie
 110     \item $1978/84/85$: apple os/atari, amiga os/mac os
 111     \item $1981/85/92/95$: ms-dos/windows 1.0/3.x/95
 112     \item $1983$: gnu project - unix-like free software development
 113     \item $1991$: linux - open-source kernel
 114   \end{itemstep}
 115 \end{slide}}
 116
 117 \overlays{7}{
 118 \begin{slide}{warning - numerical errors}
 119   \begin{itemstep}
 120     \item machine accuracy $\epsilon_m$
 121           \begin{itemize}
 122             \item ieee 64-bit floating point format: $v = -1^s 2^{-e} m$ \\
 123                   \begin{tabular}{lll}
 124                   $s$: & signe & 1 bit \\
 125                   $m$: & mantissa & 52 bit \\
 126                   $e$: & exponent & 11 bit \\
 127                   \end{tabular}
 128             \item $\epsilon_m$: smallest floating point with $1 + \epsilon_m \neq 1$ \\
 129                   $\epsilon_m \approx 2.22 \times 10^{-16}$ \hspace{2pt} (roundoff error)
 130           \end{itemize}
 131     \item truncation error $\epsilon_t$
 132           \begin{itemize}
 133             \item discrete approximation of continuous quantity
 134             \item persists even on hypothetical perfect computer ($\epsilon_m = 0$)
 135             \item machine independent, characteristic of used algorithm
 136           \end{itemize}
 137   \end{itemstep}
 138 \end{slide}}
 139
 140 \overlays{4}{
 141 \begin{slide}{warning - recursive functions}
 142   \begin{itemstep}
 143     \item avoid recursive functions!
 144           \verbatiminput{fak1.c}
 145     \item better:
 146           \verbatiminput{fak2.c}
 147   \end{itemstep}
 148 \end{slide}}
 149
 150 \begin{slide}{computational techniques}
 151   \begin{minipage}{5.5cm}
 152     \begin{itemize}
 153       \item rough discretization
 154       \item solution of linear algebraic equations
 155       \item interpolation and extrapolation
 156       \item integration of functions
 157       \item evaluation of (special) functions
 158       \item monte carlo methods
 159     \end{itemize}
 160   \end{minipage}
 161   \begin{minipage}{5.5cm}
 162     \begin{itemize}
 163       \item eigensystems
 164       \item spectral applications
 165       \item modeling of data
 166       \item ordinary differential equations
 167       \item two point boundary value problems
 168       \item partial differential \\
 169             equations
 170     \end{itemize}
 171   \end{minipage}
 172 \footnote{http://www.nr.com/}
 173 \end{slide}
 174
 175 \overlays{3}{
 176 \begin{slide}{first steps: rough discretization}
 177   \begin{itemstep}
 178     \item example: homogenous field of force $\vec{F} = (0,-mg)$ \\
 179           \begin{tabular}{ll}
 180           equation of motion: & $\vec{F} = m \vec{a} = m \frac{d^2 \vec{r}}{dt^2}$ \\
 181           initial condition: & $\vec{r}(t=0) = \vec{r_0} = (x_0,y_0)$ \\
 182                              & $\frac{d \vec{r}}{dt}|_{t=0} = (v_{x_0},v_{y_0})$ \\
 183           \end{tabular}
 184     \item algorithm using discretized time ($T = N \tau$):
 185           \begin{tabular}{lll}
 186           $x^1 = x_0;$ & $y^1 = y_0;$ & \\
 187           $v^1_x = v_{x_0};$ & $v^1_y = v_{y_0};$ & \\
 188           loop: & $x^2 = x^1 + \tau v^1_x;$ & $y^2 = y^1 + \tau v^1_y;$ \\
 189                 & $v^2_x = v^1_x;$ & $v^2_y = v^1_y - g \tau;$ \\
 190                 & $x^1 = x^2;$ & $y^1 = y^2$ \\
 191                 & $v^1_x = v^2_x;$ & $v^1_y = v^2_y;$ \\
 192           \end{tabular}
 193     \item euler's method for solving o.d.e.
 194   \end{itemstep}
 195 \end{slide}}
 196
 197 \overlays{10}{
 198 \begin{slide}{euler's method: error estimation}
 199   \begin{itemstep}
 200     \item truncation error $\epsilon_t$:
 201           \begin{itemize}
 202             \item $x_{t+\tau} = x(t) + v(t,x) \tau + O(\tau^2)$
 203             \item period $T$: $O(\tau^{-1})$ steps $\Rightarrow \epsilon_t \sim O(\tau)$
 204           \end{itemize}
 205     \item machine accuracy:
 206           \begin{itemize}
 207              \item every floating point step: error of $O(\epsilon_m)$
 208              \item $O(\tau^{-1})$ steps $\Rightarrow$ error of $O(\frac{\epsilon_m}{\tau})$
 209           \end{itemize}
 210     \item optimum:
 211           \begin{itemize}
 212             \item $\epsilon \sim \frac{\epsilon_m}{\tau} + \tau$
 213             \item 64-bit: $\epsilon_m \sim 10^{-16} \Rightarrow \tau \sim 10^{-8}$
 214             \item 32-bit: $\epsilon_m = 1.19 \times 10^{-7} \Rightarrow \tau \sim 3 \times 10^{-4}$
 215           \end{itemize}
 216   \end{itemstep}
 217 \end{slide}}
 218
 219 \begin{slide}{euler's method: accuracy}
 220   \includegraphics[width=10cm]{euler.eps}
 221 \end{slide}
 222
 223 \overlays{9}{
 224 \begin{slide}{monte carlo methods}
 225   \begin{itemstep}
 226     \item algorithms for solving computational problems using random numbers
 227     \item deterministic pseudo-random sequences
 228     \item applications:
 229           \begin{itemize}
 230             \item monte carlo integration
 231             \item metropolis algorithm
 232             \item simulated annealing
 233           \end{itemize}
 234     \item advantages:
 235           \begin{itemize}
 236             \item more efficient than other methods
 237             \item no need fo simplifying assumptions
 238           \end{itemize}
 239   \end{itemstep}
 240 \end{slide}}
 241
 242 \overlays{5}{
 243 \begin{slide}{random number generator}
 244 linear congruential generator:
 245   \begin{itemstep}
 246     \item $I_{j+1} = ( a I_{j} + c ) \, mod \, m$ \\
 247           $a$: multiplier, $c$: increment \\
 248           $m$: modulus, $I_0$: seed
 249     \item minimal standard by park and miller: \\
 250           $a = 7^5 = 16807, \quad m = 2^{31} - 1 = 2147483647, \quad c = 0$
 251     \item always seed the rng
 252   \end{itemstep}
 253 \FromSlide{4}{
 254 $\Rightarrow$ sequence of integers $\in [0,m[$ \\
 255 }
 256 \vspace{2pt}
 257 \FromSlide{5}{
 258 division by modulus $\Rightarrow$ uniform deviates : \\
 259 \[
 260   p(x)dx = \left\{
 261   \begin{array}{ll}
 262     dx & 0 \leq x < 1 \\
 263     0  & \textrm{sonst}
 264   \end{array} \right.
 265 \]
 266 }
 267 \end{slide}}
 268
 269 \overlays{8}{
 270 \begin{slide}{special deviates}
 271   \begin{itemstep}
 272     \item transformation method:
 273           \begin{itemize}
 274             \item arbitrary probability distribution $\rho(y)$
 275             \item trafo: $p(x) dx = \rho(y) dy \Rightarrow x = \int_{- \infty}^y \rho(y) dy$
 276             \item get inverse of $x(y) \Rightarrow y(x)$
 277           \end{itemize}
 278     \item rejection method: \\
 279           \begin{minipage}{5cm}
 280           \begin{itemize}
 281             \item $p(x) \in [a,b]$ mit $p(x) \geq 0 \quad \forall x \in [a,b]$
 282             \item uniformly distributed $x \in [a,b]$ und $y \in [0,p_m]$
 283             \item if $y \leq p(x)$ use $x$, else reject $x$
 284           \end{itemize}
 285           \end{minipage}
 286           \begin{minipage}{5cm}
 287           \includegraphics[width=5cm]{rej_meth.eps} \\
 288           \end{minipage}
 289   \end{itemstep}
 290 \end{slide}}
 291
 292 \overlays{5}{
 293 \begin{slide}{monte carlo integration}
 294   basics:
 295   \begin{itemstep}
 296     \item $I = \int_{\Omega} f d \Omega$
 297     \item instead of regular $x_i$, choose them at random
 298     \item $I \approx \Omega <f> \pm \Omega \sqrt{\frac{<f^2> - <f>^2}{N}}$ \\
 299           $<f> = \frac{1}{N} \sum_{i=1}^{N} f(\vec{x_i})$ \\[6pt]
 300           $<f^2> = \frac{1}{N} \sum_{i=1}^{N} f^2(\vec{x_i})$
 301   \end{itemstep}
 302 \FromSlide{4}{
 303   example: gambling for $\pi$ \\
 304 }
 305 \FromSlide{5}{
 306   \[
 307   \begin{array}{l}
 308   \pi = \int_{-1}^1 \int_{-1}^1 p(x,y) dx dy \approx \frac{4}{N} \sum_{i=1}^N p(x_i,y_i) \\[6pt]
 309   \textrm{with } p(x,y) = \left\{
 310     \begin{array}{ll}
 311     1 & x^2 + y^2 \leq 1 \\
 312     0 & \textrm{sonst}
 313     \end{array} \right.
 314   \end{array}
 315   \]
 316 }
 317 \end{slide}}
 318
 319 \overlays{5}{
 320 \begin{slide}{metropolis algorithm}
 321 ising model:
 322   \begin{itemstep}
 323     \item $d$-dimensional periodic lattice
 324     \item two possible states for magnetic moment at site $i$: \\
 325           $\mu_i = \mu S_i \qquad S_i = \pm 1 \quad \forall i$
 326     \item nearest neighbors moments interact, \\
 327           interaction strength $\frac{J_{ij}}{\mu^2}$
 328   \end{itemstep}
 329 \FromSlide{4}{
 330 $\Rightarrow$ hamiltonian: $H = - \sum_{(i,j)} J_{ij} S_i S_j$ \\
 331 }
 332 \FromSlide{5}{
 333 partition function: \\
 334 \[
 335 Z = \sum_{i=1}^N e^{\frac{-E_i}{k_B T}} = Tr(e^{-\beta H})
 336 \]
 337 }
 338 \end{slide}}
 339
 340 \overlays{2}{
 341 \begin{slide}{metropolis algorithm}
 342   \begin{itemstep}
 343     \item importance sampling: \\
 344           $<A> = \sum_i p_i A_i \approx \frac{1}{N} \sum_{i=1}^N A_i$ , with \\[6pt]
 345           $\qquad p_i = \frac{e^{- \beta E_i}}{Z}$
 346     \item detailed balance \\[6pt]
 347           sufficient condition for equilibrium: \\
 348           \[
 349           W(A \rightarrow B) p(A) = W(B \rightarrow A) p(B)
 350           \]
 351           $\Rightarrow \frac{W(A \rightarrow B)}{W(B \rightarrow A)} = \frac{p(B)}{p(A)} = e^{\frac{- \Delta E}{k_B T}}$ \\[6pt]
 352           with $\Delta E = E(B) - E(A)$
 353   \end{itemstep}
 354 \end{slide}}
 355
 356 \overlays{5}{
 357 \begin{slide}{metropolis algorithm}
 358   \begin{itemstep}
 359     \item choose $W$: \\
 360           \[
 361           W(A \rightarrow B) = \left\{
 362           \begin{array}{ll}
 363             e^{- \beta \Delta E} & : \Delta E > 0 \\
 364             1 & : \Delta E < 0
 365           \end{array} \right.
 366           \]
 367     \item algorithm:
 368           \begin{itemize}
 369              \item visit every lattice site
 370              \item calculate $\Delta E$ for spin flip
 371              \item flip spin if $r \leq e^{\frac{-\Delta E}{k_B T}}$
 372           \end{itemize}
 373   \end{itemstep}
 374 \end{slide}}
 375
 376 \begin{slide}{summary}
 377   \begin{itemize}
 378     \item importance of computational physics
 379     \item things to keep in mind when doing computational physics
 380     \item euler's method for solving o.d.e.
 381     \item introduction to monte carlo methods
 382   \end{itemize}
 383 \end{slide}
 384
 385 \end{document}