### Create a StudySoup account

#### Be part of our community, it's free to join!

Already have a StudySoup account? Login here

# COMPUTATIONAL PHYSICS PHYS 5322

TTU

GPA 3.77

### View Full Document

## 16

## 0

## Popular in Course

## Popular in Physics 2

This 75 page Class Notes was uploaded by Harry Jerde on Thursday October 22, 2015. The Class Notes belongs to PHYS 5322 at Texas Tech University taught by Volobouev in Fall. Since its upload, it has received 16 views. For similar materials see /class/226440/phys-5322-texas-tech-university in Physics 2 at Texas Tech University.

## Similar to PHYS 5322 at TTU

## Popular in Physics 2

## Reviews for COMPUTATIONAL PHYSICS

### What is Karma?

#### Karma is the currency of StudySoup.

#### You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!

Date Created: 10/22/15

Monte Carlo Simulations Igor Volobouev Computational Physics PHYS 5322 Texas Tech University Monte Carlo Methods 9 A very large class of computational methods which use randomness to generate typical situations in physics engineering mathematics chemistry geology economics linguistics etc Used to solve problems in which exhaustive consideration of problem configurations is impractical or impossible the complete phase space is too large Often work very well Usually the most complicated question is understanding of convergence towards the desired answer The Metropolis Algorithm Let say we want to compute a thermodynamic quantity A in a canonical ensemble negw 56 This can be done approximately by sampling the states E Flat in S using a weight factor exp3E5IZ simple sampling Does not work too well because the density of states typically increases with E I Mth probability exp3E5IZ importance sampling Much more ef cient However it is not obvious how to do this without owing Z e p oblem was solved in 1953 perhaps earlier by Nicholas Metropolis and his colleagues The paper titled Equation of State Calculation by Fast Computing Machinesquot J Chemical Physics vol 21 no 6 1953 explains the method Markov Chains Markov Chain is a sequence of random variables in which the probability of the current observation depends on the past through the most recent observation only PX0X1v X 2v PX0X1 Possible X values are called states A Markov chain is timehomogeneous if the transition probabilities do not change over time We will only consider timehomogeneous chains in what follows Random walk is an example of a Markov chain PX0X11X1 PX0X11X1 05 Markov Chain Properties A Markov chain is called periodic if it can only return to its present state at times tk t2k etc Example random walk has period k2 If there is no such k the chain is aperiodic The chain is irreducible if we can get from any state to any other possibly in several steps E The chain is called recurrent if we are sure to come back to any state and transient otherwise A recurrent state is positive recurrent if the expected time till return is finite otherwise it is null recurrent Ergodic Theorem Stationary distribution for a Markov chain probability distribution over the states so that an ensemble of chains started from this distribution does not change the distribution Markov chain is said to possess a limiting distribution if an ensemble of chains arrives to a stationary distribution no matter what the initial distribution is Ergodic theorem a Markov chain which is aperiodic irreducible and positive recurrent has a limiting distribution For such a chain the ensemble average of some quantity when N a w is the same as the time average when t a w E Similar in spirit to the uniqueness theorem for the Poisson s equation Metropolis Idea Set up a Markov chain whose limiting distribution is exp 3ESZ In particular impose the detailed balance restriction the average number of transitions from state A to state B is the same as in the reverse direction for any A and B PX0AIX4BPB PX0BIX4APA If we want to get the correct limiting distribution PAPB should be exp 3EAEB This means that the Markov chain should be set up so that PXoAIX4BPXOBIX4A eXp3EAEB Do it like this PX0AX4B 1 if EAsEB and PX0AIX4B eXP3EAEB if EAgtEB MetropolisHastings Method The Metropolis algorithm can be generalized to any probability distribution notjust exp3ESZ Suppose we want to simulate a random variable X having density function Suppose we can simulate Ywith conditional density qyx Then the algorithm for sampling from fX is 1 At k1 stage simulate yw from qyx 2 Accept this new point and set xw yw with probability P mfwmm may gm MUM l m 3 Return to step 1 The lsing Model A simple model of a ferromagnetic E JZSISJ LLHZSI w x In the absence of spin interactions 3 HIMWm Mean field theory H O s tanh zJlts Critical temperature To is determined from 6 z 1 2 Critical exponent s as Tc 7 T Monte Carlo for the lsing Model Initialize all spins randomly at high T or to s 1 or all iat low T For each spin in the lattice ICalculate Em the energy required to flip the spin a If Emp s 0 flip the spin E If Ehp gt 0 generate a random number r between 0 anJ 1 Flip ifrs exp5Emp EAfter cycling many times overall spins record various thermodynamic quantities such as ltsgt Several example implementations are available on the web See for example httpwww nha ihu 39 39 39 quot quot 39 html Convergence of the Metropolis Algorithm Small phase space Markov chains can be studied using the transition matrix This matrix contains the transition probabilities from any state A to any other state B The eigenvalues of the transition matrix determine the convergence speed or mixing time The eigenvalue of l the largest corresponds to the limiting distribution The second largest eigenvalue determines the longest relaxation time The convergence acceleration methods do exist similar in spirit to the multigrid methods for solving PDEs For introduction see for example the Cluster Acceleration section of the handout article Build Your Own Ising Model Please implement your own version ofthe 2 dimensional Ising model Visualization of the 2d plots with Python has been demonstrated by the WaveSolver2dexamplepyquot code in Lab 6 Make your code flexible so that you can easily calculate and plot average spin heat capacitance etc as a function of temperature We will continue working on this project on Thursday Root Finding and Basic Function Minimization Igor Volobouev Computational Physics PHYS 5322 Texas Tech University Root Finding Problem 1d fx 0 Usually easy unless you need to solve gazillion equations in a limited amount oftime and fX is not cheap to evaluate Spefial case polynomials d ffc 0 Usually hard except for the important special case of linear systems Success in the nonlinear case depends crucially on having a good first guess for the solution Bracketing the Root The root is bracketed on a b if fa and fb have opposite signs Ability to bracket the root is the main difference between 1d and multivariate cases quoti ii 2 quotik a a 3 b a A quotA iquot Bisection quot irst bracket the root hen divide the interval in two and keep the root bracketed Repeat until the desired tolerance is reached You must think carefully what tolerance is appropriate for your problem Converges to a root a discontinuity point or a singularity The rate of convergence is usually expressed as em constant X an m Other Methods Without Derivatives A variety of methods speed up the convergence by making assumptions about the function behavior Superlinear convergence m gt 1 is achieved when the assumptions turn out to be correct ee Numerical Recipes for details NewtonRhapson Method Based on Taylor series near the root X fx fxf39x a Z fx f 39x Converges quadratically near the root very fast Unfortunately it can shoot off to outer space or enter a nonconvergent cycle away from the root It is usually combined with bisections to ensure convergence Roots of Polynomials The number of roots is known in advance Direct formulabased methods are appropriate for polynomial up to 3rd degree Assume simpler form near the root and iterate Muller s Laguerre s methods Then deflate the polynomial and repeat for another root JenkinsTraub and eigenvalue methods are known to be very robust but usually somewhat slower Linear Systems E Equations of the type G11 G21 0131 041 Axb B Common way to solve LU decomposition ALU 0 0 0quot 311 512 513 lml G11 G12 G22 0 0 0 B22 523 324 G21 G22 G32 G33 OJ L0 0 533 534 G31 G32 G42 G43 G44 0 0 0 344 G41 G42 ThusAxLUxLUXb Let y UX Solve for y in L y b by forward substitution Solve for X in U X y by backward substitution CL13 C523 a3 G43 G14 l 124 G34 Q44 E This is the technology for example behind the Iinsolve module in SciPy LAPACK is often used for CFortan calculations Nonlinear Systems Very tough no general solution methods known A multivariate extension ofthe NewtonRhapson method works if a good initial guess can be provided e 5 g J 5 J 3 for The last equation can be solved eg by LU decomposition Degenerate Jacobian spells trouble Function Minimization A particular optimization problem which happens very often in physics Dichotomies local global 1d multivariate using not using derivatives Bracketing the Minimum Three points are needed to bracket the minimum a lt b lt c so that fb lt fa and f b lt f c f x Golden Section Search Similar to nding roots by bisection just keep reducing the bracket Next point to be tried is the one which is fraction 3 V52 038197 into the larger of the two intervals measuring from the central point of the triplet The interval is thus reduced by factor 061803 on each step Brent s Method Assume the function is parabolic near the minimum Pick the next point at the minimum of the parabola Converges quickly if the starting bracket is not too far away Derivative information can be easily included Simplex Method Multivariate Also called NelderMead method Simplex operations are illustrated below for a tetrahedron mp1 a m Iugh hm contraction multiple columctmn Faster Methods Parabolic Assumptions Conjugate gradient methods QuasiNewton methods DavidonFletcher Powell DFP and BroydenFletcher GoldfarbShanno BFGS Steepest descent does not work particularly well due to its very slow progress in narrow valleys Finding the Global Minimum Very dif cult Thing to try first just apply local minimum codes with many starting points Use quasiMC to choose initial points Simulated annealing might or might not help analogs of phase space and energy required and they are not always easy to construct Basic ODE Solving Igor Volobouev Computational Physics PHYS 5322 Texas Tech University Newton s Laws 11 mii Particle model physics is simple and works for a lot of problems Usually force is a function of time coordinate andjor Lelocity but not acceleration F Ft c17 With the computer we can handle much more complex Fthan analytically Simple Model no x in F If the force is independent from particle coordinates the equation becomes the first order ODE for velocity d6 m F t7 dt Terminology order is the degree ofthe derivative F can be a highly complex function of v yet this remains the first order equation Example bicycle motion on a flat terrain F 7 lesz v 2 Numerical Approach Let say we want to know vt Taylor series give or 1 dz 17 7At 17 2 At dt t 2 alt2 tr If At is small then 17TAt177 N TFTEAt dt tr m Euler s Algorithm Euler s algorithm moves the system forward in small steps 2 q Mm m 121 1 N Selfstarting from an initial condition M hoice of a good At value is empirical and it depends on the target precision and the typical time scales in the system Precision of the Euler s Algorithm 9 Euler s method ignores terms ofthe order At2 in the Taylor s expansion If we run the simulation for time T then there will be TAt such terms Therefore the algorithmic uncertainty of the Euler s method is proportional to At This kind of proportionality is usually expressed as 0At big Oh notation Other computational schemes use more terms in the Taylor s expansion which results in 0A V convergence with N gt 1 These schemes perform substantially better in case F is a sufficiently smooth function Second Order ODEs Newton s equation for coordinates d2 alt2 Second orderquot due to second derivative m 160116 d x alt2 Needs more initial conditions 2 00 9800 Standard Form for ODEs Write all ODEs in the following form dIt a t dt f y Trick generalize coordinates 0 7 0 17 0 y 7x y 7 dt 1 dt 7y if F0Iy0y1 dz m First Principle of Efficient Programming Don t Repeat Yourself Euler s algorithm worksjust ne for the standard ODE form It is to our bene t to separate the ODE solution algorithm and the physical model There will be many types of physical systems and a few different ODE solvers We want to recombine them easily without duplicating code DRY All the physics is contained in F Callables and Functors Fcan be a complicated function Best way to represent such a function in Python is to create a callable object Any object which can be called like a function is a callable To make such an object implement a class which includes the special function callself Callables can perform a lot of work during their initialization stage which can later be used to save computing time no need to repeat the initialization at every step In C similar objects are called functors They must provide the function operator Let s Look at Some Code Examples Now Go to the course web page httphiqhenerqvphvsttueduiqvComputationalthsicsLecturesLe cture2Examgles Model for the Air Drag Stoke s drag viscous resistance Quadratic drag a 1 A FD ECpszv Drag dependence on altitude via air density p poeXPy yo y0 kBTmg m10gtlt104m Realistic Physics of Baseball Escalation of the Air Drag Model Velocitydependent drag coefficient Variable surface roughness knuckleball Headtail wind Magnus force Magnus Effect Important for taking into account the effects of ball spin Combination of the Bernoulli effect and formation of boundary layers FM Sow X vbau Air Drag from First Principles Can we figure out all these forces from first principles ie Newton s laws and the laws of thermodynamics Possibly Newton s laws in application to fluids are expressed in terms of NaVIerStokes equations The Computational Fluid Dynamics CFD discipline attempts to solve these equations numerically A lot of applied CFD research is carried out at Boein Lockheed SA etc Wind tunnels are sti essential The Cla Mathematics Institute offers 1000 00 to anybody who can prove or disprove general existence and smoothness of NavierStokes solutions Let s Get to Work Please code up the air drag models and run the Euler solver on them Do your results look like something you d expect Discuss with your neighbor By next time please read Appendix A from Giordano amp Nakanishi Errors and Uncertainties in Computations Igor Volobouev Computational Physics PHYS 5322 Texas Tech University Taxonomy of Errors Wrong or Insufficient Models Hardware Errors 39 Uncertainties in the Data pproximation Errors RoundOff Errors Bugs Unlike in math there is no way to formally prove that a program is correct at least not yet 1996 explosion of the Ariane 5 rocket cost 500 M an unintended overflow of a short integer variable in the rocket guidance software T88ical commercial software has 17 bugs per 1 0 lines of code You will have to deal with bugs in your code all the time Fortunately things will gradually Improve quot 39 better 39 39 39 Testing unit integration functional breakdown Better programming languages Wrong or Insufficient Models Learning how to get the models right is the one of the main purposes of this course Bodies move in elliptical orbits in a central gravitational field Such a nice abstraction Now lets try to use satellite orbits so that we can determine our own coordinates with 1 m precision Suddenly we have to take into account that E Earth is neither round nor uniform it rotates and has tides E There are many other bodies in the solar syste Sun exerts radiation pressure I There is some residual atmosphere even at 20000 km altitude E Clocks run at different rates on the satellites and on the ground I Signals move in curved trajectories just as Einstein predicted i Atmosphere ionosphere in particular is a refractive m dium If you ve done all you can and the model still does not work perhaps you have discovered new physics Hardware Errors Usually easy to deal with lThe computer actually tells you when your hard drive starts dying Just pay attention However certain failure mechanisms are difficult to detect and control Design mistakes pentium FDIV bug Electromigration Power glitches Single event upsets Redundancy helps when calculations are missioncritical Example triple or quadruple computer redundancy is standard on space systems and ICBMs You might want to repeat a longrunning calculation Uncertainties in the Data Scientific data is almost always imprecise EWe will usually assume that no gross mistakes were during data collection The area of knowledge which studies methods of extracting meaningful results from gar age contaminated data is called robust statistics One of the early decisions you have to face when you build your computational model what to do With measurement errors These are the common chOIces H Ignore at your peril Assume uncertainties are Gaussian use linear error propagation techniques Model the measurement process in detail Approximation Errors Bread and butter of applied mathematics studies We will see many examples in this course Most often there is a tradeoff between the time needed to run the algorithm and the accuracy achieved 9quot i xn z i Cn equot gxN quot0 quot0 RoundOff Errors Numbers are represented in a computer with limited range and accuracy Hypothetical example let say computer keeps four decimal places Then 13 03333 and 23 06667 Now 213 23 00001 0 Representing numbers in double precision floating oint eveg such number in a computer is enco ed using 6 bits is a reasonable trade off between accuracy range and computation speed This may not be optimal for all types of applications Arbitraryprecision arithmetic software ackages do exist and can be used if necessary They incur a huge hit in computation speed FPU becomes useless Floating Point Numbers Most computers conform to the IEEE 754 standard This is great it allows you to get the same result using the same algorithm on different computers and coded in different programming languages 9 haulquot Sign in hit lwrm xFP 15 x 1 f x 2840 mantissa1f1b51 gtlt2391 b50 gtlt2392 b0 gtlt23952 IEEE 754 Representation for Doubles Number Type s e and f Value Normal 0 lt e lt 2047 HY xlfxzmm Subnormal e O f O 1s xofxzlm Signed zero e O f O 1s x00 s0 e2047 f0 nf s1 e2047 f0 nf Not a number e 2047 f O NaN Properties of IEEE 754 Doubles 9 Precision 1 part in 252 approximately 16 decimal places Range from 49 x 10324 to 18 x 10308 Underflow happens when result magnitude is smaller than the smallest positive number Usually converted to O Overflow happens when result magnitude is larger than the largest positive number Usually converted to lnf 00 oxen eoeo eoeo etc produce NaNs Model for Computer Accuracy Machine precision am This is the smallest positive number so that 1 am a6 1 for the computer An arbitrary number X is represented on the computer as XFP X1 5 Isl lt 5m where the actual value for the error a is not known RoundOff Disasters If you subtract two large numbers and end up with a small result precision of the result will be worse than the precision of the terms This effect is called subtractive cancellation E Beware rapidly changing functions which map large numbers into small This effect is known as error ampli cation Subtractive Cancellation a b c in the computer turns into aFP bFPCFP b15b C15c Then aFPa 1 ab ba ac ca But aFPa 1 ea by definition If a is small then b z c and aa z ba ab 80 z ba maxeb C where ba is a large number Error Amplification Let y fx In the computer this turns into yFP fx1ex Using Taylor series at x pr W mm x y lt1 f39ltxgtltex xgtygt Then 8y ex f39x xy so that the error amplification factor is large in case f39x xy is large RoundOff Error Accumulation a b x c in the computer turns into aFPbFPxCFPb1 512 x0 5c Then aFPa 1 ab 1 0 1 ab so That is relative errors are added Assume that all errors are uncorrelated and ofthe order am Then the overall error from N multiplications grows as amxN Error Minimization Approximation errors in algorithms often behave as GNk for large IV where N is the number of steps and d and k are some empirical constants If we assume that roundoff errors from different steps are uncorrelated then the total error behaves as stotal GNk BN where 3 is the error from one step usually on the order of several am The first term dominates at small N and the second at large N As a consequence there is an optimal value of N at which atotal is minimized trapezoid lerrorl a 1039 39 10 100 N Fig 54 Log log plot of the error in integration of exponential decay using the trapezoid rule Simpson s rule and Gaussian quadrature versus the number of integration points N Approxi mately seven decimal places of precision are attainable with single precision shown here and 15 places with double precision From the book Computational Physics by RH Landau et al Empirical Approach Best way to study convergence compare with a known answer for some parameter values in your problem If there is no reference case study the behavior of your result R as a function of N Plot RN 7R2N R N as a function of logmN RZN approximates the exact answer If you obtain a straight line dropoff you can educe the value of k from the slope If your code is fast run until you see the onset of the roundoff error dominat39on Be sensible the error model just described might not work for your particular algorithm Also converged correct logic Let s Play I Hopefully this lecture will help you scratch your head in a more systematic manner when results of your calculations do not look right I Let s do some exercises now h l Iquot I I nhvs HI quot II 1quot 1 alPhysicsLabslab1html I Reading assignment for the next class chapters 1 and 2 from Giordano amp Nakanishi Don t forget the programming homework you already have Patterns and Fractals Igor Volobouev Computational Physics PHYS 5322 Texas Tech University Important Reminder H You need to define and let me know the topic of yourfinal project by March 13 see the course Syllabus E Suggestions for selecting a topic Apply some of the techniques learned in this course to a complex problem in your primary field of study Take any of the computational methods studied in this course ODEPDE solving optimization numerical integration Monte Carlo data analysis etc to its limits E Search the web for quotcomputational physics projects Several journals publish articles on various modern physics topics including computational at a very accessible level I Physics Today I Computing in Science amp Engineering I American Journal of Physics European Journal of Physics I Reviews of Modern Physics Reports on Progress in Physics E There are many more specialized journals International Journal of Modern Physics C Computer Physics Communications ACM Transactions on Modeling and Computer Simulation SIAM Journal on Scientific Computing Mathematics and Computers in Simulation etc httparxivorg B Talk to me if you feel stuckoverwhelmed Mandelbrot Set Defined as the set of values I in the complex plane for which the iteration Xn1 Xn2 0 remains bounded One of the first fractals studied boundary is the fractal not the set itself Exhibits selfsimilarity around certain points in the structure Fractal Dimension 50 51 A There is no unique way to AA define the dimensionality of A a fractal AA 5 Amount of stuff M Ld A A A AA Density is then 0 L A 569 feta Mathematicians call this d the HausdorffBesicovitch dimension 2 Sierpinski triangle d 1 585 The Koch Curve Example fractal curve discussed by GampN Fractal dimension is defined from L6r NSLs L1d d n4n3 1262 Gen 1 Fractal Dimension on a Computer Box counting method 2d EUse finer and finer grids N2 squares Of course the smallest square should still be larger than the fractal resolution pixel size ECount the number of squares SN needed to cover the fractal This number should scale as SN Nd Average over many realizations of the fractal Plot SN vs N using loglog plot Statistical SelfSimilarity Exact selfsimilarity a set is composed of n scaled copies of itself translated and possibly rot Certain types of fractals exhibit statistical self similarity This means that certain average quantities such as distribution momen ts are invariant to a change of scales Fractals in Nature It seems occurrences of fractals in nature are very common E Largescale galaxy struc ures Earthquake magnitudes es Networks of rivers or blood vessels 39 39 lant sha es pidemics forest fires suf x 3 v E Mechanisms for Producing Fractal Patterns Phase transitions E Density fluctuations occur at all scales near the t critical pom Umbrella term selforganized criticality Refers to the tendency of large systems to drive themselves to a critical state with a wide range of length and time sca es Avalanches in a pile of sand is a canonical example Universality Behavior of systems near the critical point is unusually independent from the microscopic system details Both magnetization of a erromagnetic near its Curie point and the density ofa liquid near its critical point behave as T TC3 with the same exponent 3 0325 n two different systems have the same behavior near the critical point they belong to the same universality class System universality properties are closely related to its selfsimilarity under scale transformations coarse graining The way in which the system maps into itself under these transformations is called renormalization group

### BOOM! Enjoy Your Free Notes!

We've added these Notes to your profile, click here to view them now.

### You're already Subscribed!

Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'

## Why people love StudySoup

#### "I was shooting for a perfect 4.0 GPA this semester. Having StudySoup as a study aid was critical to helping me achieve my goal...and I nailed it!"

#### "Selling my MCAT study guides and notes has been a great source of side revenue while I'm in school. Some months I'm making over $500! Plus, it makes me happy knowing that I'm helping future med students with their MCAT."

#### "I was shooting for a perfect 4.0 GPA this semester. Having StudySoup as a study aid was critical to helping me achieve my goal...and I nailed it!"

#### "Their 'Elite Notetakers' are making over $1,200/month in sales by creating high quality content that helps their classmates in a time of need."

### Refund Policy

#### STUDYSOUP CANCELLATION POLICY

All subscriptions to StudySoup are paid in full at the time of subscribing. To change your credit card information or to cancel your subscription, go to "Edit Settings". All credit card information will be available there. If you should decide to cancel your subscription, it will continue to be valid until the next payment period, as all payments for the current period were made in advance. For special circumstances, please email support@studysoup.com

#### STUDYSOUP REFUND POLICY

StudySoup has more than 1 million course-specific study resources to help students study smarter. If you’re having trouble finding what you’re looking for, our customer support team can help you find what you need! Feel free to contact them here: support@studysoup.com

Recurring Subscriptions: If you have canceled your recurring subscription on the day of renewal and have not downloaded any documents, you may request a refund by submitting an email to support@studysoup.com

Satisfaction Guarantee: If you’re not satisfied with your subscription, you can contact us for further help. Contact must be made within 3 business days of your subscription purchase and your refund request will be subject for review.

Please Note: Refunds can never be provided more than 30 days after the initial purchase date regardless of your activity on the site.