Computational Physics PHYS 4007
Popular in Course
Popular in Physics 2
This 282 page Class Notes was uploaded by Iva Cormier on Sunday October 11, 2015. The Class Notes belongs to PHYS 4007 at East Tennessee State University taught by Donald Luttermoser in Fall. Since its upload, it has received 123 views. For similar materials see /class/221411/phys-4007-east-tennessee-state-university in Physics 2 at East Tennessee State University.
Reviews for Computational Physics
Report this Material
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/11/15
PHYS40075007 Computational Physics Course Lecture Notes Section VII Dr Donald G Luttermoser East Tennessee State University Version 41 Abstract These Class notes are designed for use of the instructor and students of the course PHYS 40075007 Computational Physics taught by Dr Donald Lutternioser at East Tennessee State University VII Methods of Differentiation and Integration A Derivatives 1 The derivative of a function is de ned by a b C 01 df06 E hm f06 0 x dx hgt0 h In computational work functions and their independent variables are given by tabulated data and or derived data Since there is a subtraction in Eq VII l subtraction cancellation can lead to rather large errors in the deter mination of a derivative via numerical techniques The computer s nite word length can cause the numera tor to uctuate between 0 and the machine precision em as the denominator approaches zero In this section we will discuss the various techniques used to calculate derivatives numerically and estimating the error in the derivative 2 Method 1 Forward Difference 3 0 Write out the function as a Taylor series at a position of one step forward from the current position h2 h M h M was 51 g gwc 39 39 VII2 where h is the step size as shown in Figure VII l We obtain the forward dz erence derivative algorithm by solving Eq VII l for fx Joex 2 W7 VH3 VII1 VII 2 PH Y5 4007 5007 Computational Physics h 2 HM E f 06 7 VII4 where the subscript c denotes the computed expression i The approximation of Eq VII 3 has an error pro portional to h as shown in Eq VII 4 ii We can make the approximation error smaller and smaller by making h smaller and smaller iii We can t make it too small however since all precision will be lost if the subtraction cancellation error becomes larger than the step size c Consider for example my a bx VII5 the exact derivative is f 2bx VII 6 and the computed derivative is W m m h me This would only be a good approximation if h ltlt 206 2m bh VII7 3 Method 2 Central Difference a An improved approximation to the derivative starts with the basic de nition Eq Vll l For this technique in stead of making a step of h forward we form a central di erence by stepping forward by h 2 and stepping back ward by h 2 06 h2 x h2 h EDCfxh Vii 8 Donad G Luttermoser ETSU VII 3 y AA xhl2 x xhl2 xh K Figure VII 1 Forward difference solid line and central difference dashed line methods for numerical integration i The symbol DC represents center difference ii Carrying out the Taylor series for x I h gives 1 f 2 fx i f m VII9 iii The important difference from Eq VII 3 is that when x h2 is subtracted from x h2 all terms containing an odd power of h in the Taylor series cancel iv Therefore the central difference algorithm be comes accurate to one order higher than h gt hg v If f3h224 ltlt f2h2 then the error with the central difference method should be smaller than the forward difference b For the polynomial given in Eq VII 5 the central dif VII 4 PH Y5 4007 5007 Computational Physics ference gives the exact answer regardless of the size of h m m h2 m h2 h 2m VII10 x 4 Errors in Taking Derivatives a One should always try to keep calculation errors em at a minimum This typically occurs when 60 m eaer in Eq V 33 b Because differentiation subtracts 2 numbers that are close in value the limit of roundoff error is essentially machine precision fx m fx h mg m in VII11 h h gt em m 6 VII12 c The approximation error with the forward difference al gorithm Eq VII 4 is an 9h term while that with the central difference algorithm Eq VII 9 is an 9h2 term M N 102 a1er N 2 VII13 C f3 h eaim m 24 VII14 d The value of h for which roundoff and approximation er rors are equal is therefore in N fd ifem h N apply 7 2 7 em C 1432 I g gaipm 24 VII 16 2 m 24 m gt hi 6 3 6 VII17 W7 Cd f3 Donald G Luttermoser ETSU e As an example for the ex and cosx functions f m fa m f 3 in single precision with em 10 7 one would get hfd m 00005 and had 001 This makes the central difference algorithm better for calculating derivatives since a larger h would mean a smaller error gt here the error in the central difference method is 20 times smaller than the er ror in the forward difference method 5 The Method in Calculating Second Derivatives a Taking second derivatives will involve an additional sub traction step as compared to the rst derivative calcula tion leading to a larger subtraction error b We can remedy this with a little algebra in the central difference method Taking the second derivative of Eq VII 8 and then rewriting the rst derivatives with a for ward and backward difference equation we get ism W mm W m h2gtl VII19gt Note however that one must keep the second I Z derivative formula in this form to minimize can cellation error E Integration 1 The method of numerical integration is sometimes referred to as numerical quadrature gt summing boxes a The de nition of an integral given by Riemann is himh from m VII20gt VII 5 VII 6 PH Y5 4007 5007 Computational Physics b If we ignore the limit the integral just becomes a sum mation of boxes or quadrilaterals lying below the curve b N a my day m mi w VII21 where N is the number of points in the interval 1 b and f is evaluated at those interval points 2quot f E The wi s are integration weights which are proportional to h the distance between points 239 and 239 l c The different integration schemes presented here will all make use of Eq VII 21 d In the simplest integration schemes the integrand is ap proximated by a few terms in the Taylor series expansion of f and the terms are integrated gt typically adding more terms in the series yield higher precision e This is the basis of the NewtonCotes methods gt the total interval is divided into equal subintervals with the integrand evaluated at equally spaced points 06 Two such Newton Cotes methods include the i Trapezoid rule a rst order method ii Simpson doh l rule a second order method f More accurate integration methods involve the use of non equally spaced intervals gt these are the methods of Gaus sian quadrature i Gaussian quadrature methods are typically more precise than Newton Cotes methods as long as there are no singularities e denominators going to in nity non continuous functions in the integrand or its derivative Donald G Luttermoser ETSU ii You are better off to remove such singularities algebraically before attempting Gaussian quadra ture For example 1 x fa die L01f xdx 01fxdx VII22 iii Regions where a function is slowly varying re quire fewer integration points and regions with rapid variations in the function will have many in tegration points in order not to miss any oscillation as can be seen for such functions evenly spaced integration points will not represent the true inte gral accurately 2 Trapezoid Rule a b C As an introduction to the numerically integration of a function consider the generic integral I b mm VII23 The most straight forward way to solving such an integral function is to evaluate f at a few points and t a simple curve eg piecewise linear through these points One way to do this is to t trapezoids of equal width h under the curve represented by f and add up the total areas of these trapezoids Let h b a VII24 N l7 x ai 1h i12N VII 25 where a and b correspond to the initial and nal end points and N the total number of points in the interval VII 7 VII 8 PH Y5 4007 5007 Computational Physics 1 b square bracket means we include the endpoints 239 6 there are N 2 points in between the endpoints Note that the trapezoid rule requires an odd number of points N i Straight lines connect the points and this piecewise linear function serves as our tting curve ii The integral of this tting function is easy to com pute since it is the sum of the areas of trapezoids The area of a single trapezoid is 1 Ti 596241 06ifi1 fi VII26 iii The true integral is estimated as the sum of the areas of the trapezoids so Nil IRITT1T239HTN1 Z VII27 il Notice that the last term in the sum is N 1 since there is one fewer panel than grid points iv The general formula simpli es if we take equally spaced grid points as given by Eq VII 25 v Then the area for an individual trapezoid e the area of the one trapezoid within that interval becomes 1 Ti hfil fi 7 VII28 or in our original integral notation 41 fltxgtdx 2 W Ehfi hfi VII 29 As such our weights in Eq VII 21 for the indi vidual points are U h Donald G Luttermoser ETSU d e f g We now need to add up all of the trapezoids in the subin tervals across the entire interval 1 b gt the trapezoid rule 1 1 TWO Ehfl hf2 hf th71 Eth 1 N71 hf1 fN h fi VII 30 OIquot 1695degfihf2hfgth71ng VII31 Note that since each internal point gets counted twice it has a weight of h whereas the endpoints get counted just once and thus have weights of h 2 LUZ39 hh VII 32 2 2 lt gt Our approximation error also called the truncation error or the algorithmic error here for the trapezoid rule can be written either as can I ITxh T12b luffc VII33 for some value x C that lies in 1 b or as Gaper 112 hglfb fa 9h4 VII 34 This error is proportional to h2 and Eq VII 34 warns you that the trapezoidal rule will have dif culties if the derivative diverges at the end points VII 9 VII 10 PH Y5 4007 5007 Computational Physics 3 Bomb erg Integration 3 m But how many panels ta trapezoids should one use ta how small should h be One way to decide is to re peat the calculation with a smaller interval If the answer doesn t vary signi cantly then we accept it as correct i This doesn t work with pathological functions or in unusual scenarios but this won t happen very often ii In the trapezoidal rule if the number of panels is a power of two we can halve the interval size without having to recompute all the points iii De ne the sequence of interval sizes as l hn VII35 hl b a hg b a iv For n l we have only one panel so hvw1 aMKwfwt 2 VII36 v There is a simple recursive formula for calculating Thn1 from Thn 1 2n1 Thn1 hn1 fla 2i lhn1l 7 VII37 By using this recursive method we can keep adding pan els until the answer appears to converge can greatly improve this process by using a method called Romb erg integration However we Donald G Luttermoser ETSU VII 11 i The method builds a triangular table of the form R11 R21 M2 R31 R 2 ii The formula for the rst column is just the recur sive trapezoidal rule Rn1 1mm VII38 iii The successive columns to the right are com puted using the Richardson extrapolation for mula Rnlml 7 nlm VII 39 1 4m llRm 1m Rnm iv The most accurate estimate for the integral is Rnm the value at the bottom right corner of the table v To understand why the Romberg scheme works consider the approximation 16 truncation error capprxhn I Thn then using Eq VII 34 swam gulfb fal 0 VII40 Since hnH hnZ 1 Ehilfb W W VII41 Eappnc hn1 VII 12 PH Y5 4007 5007 Computational Physics vi Consider the second column of the Romberg ta ble The approximation error for Rn 12 is I Rm 1 2 I IThn1 lIThn1 Thnl ewxmn ieapmlthn1gt swim g hi W fal 00 90 VII42 vii Notice that the hi term serendipitously cancels out leaving us with a truncation error that is of order hi The next third column of the Romberg table removes this term and so forth c The following is an example of an integration subroutine using the Romberg method If you are interested in using this routine let me know and send you a copy of the routine shown below written in Fortran 77 Note that func is an externally supplied subroutine that contains the function to be integrated subroutine rombfa b N func param R NR l Routine to compute integrals by Romberg algorithm real RNRNR param external func h b a This is the coarsest panel R11 h2 funcaparam funcbparam up 1 up is the number of panels do i 2 N h h 2 Donald G Luttermoser ETSU do k 1 npl 2 This loop goes k135 sum sum funcakh param enddo npl Compute the first column entry Ri1 Ri1 O5Ri 11 hsum m 1 do j 2 i m 4 m Rij Rij 1 Rij 1 Ri 1j 1m 1gt enddo enddo return end 4 Simpson s Rule Simpson eh 3 b C 01 Instead of tting two adjacent points with trapezoids we will now t three adjacent points with parabolas x m ago x 7 VII 43 for each interval still keeping the intervals evenly spaced The area of each section is then the integral of this parabola Iih 0063 062 Th A ax2 x7dx 7x Ii VII44 This is equivalent to integrating the Taylor series up to the quadratic term Hence the Simpson rule is a second order polynomial method In order to relate the parameters oz and 7 to the func tion we consider an interval from 1 to 1 20 2 37 lmy x 7 dx VII 45 Compute the other columns Ri2Rii VII 13 VII 14 PH Y5 4007 5007 Computational Physics i Note however for the function itself f1 a a flt1gtlt 1gt f0 ND 7 W N a 71 0 VII46 ii Using the results of Eq VII 46 in Eq VII 45 yields 1ax2 x7dxm 3 3 VII47 iii Because 3 values of the function are needed we evaluate the integral over two adjacent intervals gt evaluate the functions at the two endpoints and the middle zih xih Ii AH mm fxdxIiihfxdx h 4h h gfiil Efi gfz H VII 48 2 e Simpson s rule requires the elementary integration to be over pairs of intervals gt this requires the number of intervals to be even and hence the number of points N to be odd f To integrate across the entire range 1 b we add up con tributions from each pair of subintervals counting all but the rst and last endpoints twice ffltxgtdx w gf1f2fgf4m fN71fN VII49 Donald G Luttermoser ETSU i From this integral we see that the total set of weights is h 4h 2h 4h 2h 4h h 10 i u ii The sum of your weights provides a useful check on your integration w N 1 h VII51 il Remember N must be odd 5 Errors in the NewtonCotes Methods a As we have said above and as was the case for differ entiation our Newton Cotes methods are essentially just Taylor series expansions of a function We will highlight here what is derived in more detail in the textbook b The approximation ie truncation or algorithmic error 6mm can be estimated by the next higher term in the Taylor series that is not used in the evaluation of the integral 3 EaPpDGt O 7 b a 5 EaPPDFS O a where the subscripts t and s in the subscripts refer to trapezoid and Simpson s rule respectively VII53 c The relative error e is just these approximation errors divided by the value of the function as VII54 VII 15 VII 16 PH Y5 4007 5007 Computational Physics 1 Assume that after N steps the relative roundoff error is random so em VN em VII 55 where em is the machine precision N 10 7 for single pre cision 10 15 for double precision e v Assume that the total error is given by Eq V 33 Etot Em apprx VII56 i We want to determine the value of N that min imizes the total error This will occur when the two relative errors are of equal magnitude which we approximate even further by assuming that the two errors are equal Eappmet 6T0 6appnc f ii Now let s make the following assumptions n f m 1 VII58 f l b 1 hi VII 59 a gt N f For the trapezoid rule Eq VII 57 becomes Mb agt3 1 VN em m fN2 N2 VII 60 l i Then the optimum numbers of N steps for the trapezoid rule are N i l i 110 725 631 for single precision h 110 1 5 5 106 for double precision VII62 Donald G Luttermoser ETSU VII 17 ii The corresponding errors are 60 Wem 3 X 10 6 for single precision 10 12 for double precision VII 63 g For Simpson s rule Eq VII 57 becomes N fl4b a5 1 N 6m N W m l gtN m VII 65 i Then the optimum numbers of N steps for Simp son s rule are N i l i 110 729 36 for single precision T h T 110 15 9 2154 for double precision VII 66 ii The corresponding errors are i 6 X 10 7 for single precision 6T0 N N 6m 7 l 5 X 10 14 for double precision VII 67 h These results are illuminating because they show that i Simpson s why you little rule is an improve ment over the trapezoid rule ii It is possible to obtain an error close to the ma chine precision with Simpson s rule and with other higher order integration algorithms iii Obtaining the best numerical approximation to an integral is not obtained by letting N gt 00 but with a relatively small N g 1000 VII 18 PH Y5 4007 5007 Computational Physics 6 There are higher order Newton Cotes methods the 3rd degree 38th method and the 4th degree Milne method See Table 41 in your textbook on page 48 for the elementary weights that one would use for these two methods 7 Gaussian Quadrature a Often one cannot choose an evenly spaced interval of grid points to do a numerical integration From Eq VII 21 we have my day m w1f061 waxN VII68 One then can ask the question is there an optimal choice for the grid points or nodes x and the weights 10 to solve this integral b This question leads us to formulate a new class of integra tion formulas known collectively as Gaussian quadra ture i In this class we will use only the most common formula namely GaussLegendre quadrature c There are many other kinds of Gaussian quadra ture that treat speci c types of integrands such as the Gauss Laguerre formularization which is optimal for integrals of the form 5 e xfm dx 0 The derivation of the other Gaussian formulas see Table 42 in your textbook is similar to our analysis of Gauss Legendre quadrature ii The theory of Gaussian integration is based on the following theorem 0 Let qx be a polynomial of degree N such that b qxpxxkdx 0 VII 69 a Donald G Luttermoser ETSU where k 12 N 1 and px is a speci ed weight function 0 Call 061 062 xN the roots of the polynomial Using these roots as grid points plus a set of weights 101102 wN we construct an integra tion formula of the form Ab fltxgtpltxgt dx z w1fx1 mom VII70 0 There is a set of w s for which the integral for mula will be exact if f is a polynomial of de gree lt 2N iii The weights can be determined from the Three Point GaussianLegendre Rule For example consider the interval 1 l with px 1 This gives us a Gaussian Legendre formula For inte grals in the interval 1 b it is easy to transform them as ib a meldx 2 1fzdz VII 71 using the change of variable 06 b a b 12 o The rst step is to nd polynomial We want a three point rule so that qx is a cubic 106 co 0106 02062 03x3 VII 72 0 From the theorem of Eq VIII 69 we know that lqux 1Qdx1x2qxdx0 VII73 VII 19 VII 20 PH Y5 4007 5007 Computational Physics 0 Plugging in and doing each integral we get the equations 2 2 2 2 2 2 7 7 7 7 7 0 VII 74 00362 361363 3605C2 l o The general solution is co 001 a02 003 513 where a is some constant This arbitrary constant cancels out in the second step 0 Using this in Eq VIII 72 gives a polynomial solution of 5 3 106 i063 Ex VII 75 Notice that this is just the Legendre polynomial 0 Next we need to nd the roots of qx 13306 This cubic is easy to factor The roots are 961 35x2 0 and x3 Using the grid points in Eq VIII 70 gives 1 M dx m w1fMw2f0w3fW VII76 0 Finally to nd the weights The above formula must be exact for 106 x5 We can use this to work out values of wl wg and 103 It turns out to be su icient to consider just 106 and x2 which produces 3 equations 2 LU1 LU2 LU3 0 35w1 35w3 VII77 2 3 3 3w13w3 o This linear system of equations is easy to solve giving wl 59w2 89w3 59 Donald G Luttermoser ETSU iv c You will typically never have to write your own Gaussian quadrature subroutines since they are included in math libraries and built into some programming languages 6 9 the IDL functions NT2D and NT3D the QGAUS An alternative way of nding the weights is to use the identity 2 1 xddxPNxz2 where N 3 in our example This formula may be derived from the recurrence relation for Leg endre polynomials w VII78 There are various advantages and disadvantages in using Gaussian integration Advantage A very high order accuracy is ob tained for just a few points often this method yields excellent results using fewer than 10 points Disadvantages l The node points and weights must be computed or obtained from tables This step is nontrivial if you want to use many node points Using more than N 20 points is rarely worth it since badly behaved functions will spoil the results in any case 2 Unlike Newton Cotes integration the method does not lead itself to iteration nor is it easy to estimate the error subroutine in Numerical Recipes VII 21 PHYS40075007 Computational Physics Course Lecture Notes Appendix B Dr Donald G Luttermoser East Tennessee State University Version 41 Abstract These Class notes are designed for use of the instructor and students of the course PHYS 40075007 Computational Physics I taught by Dr Donald Lutternioser at East Tennessee State University B Scienti c Computing Using C A Tutorial Introduction to C 1 Here is everyone s rst C program It print the words hello 77 world on the screen 4 helloc jgc 31193 from KampR include ltstdiohgt int mainvoid printfquothello worldnquot return 0 a A C program functions executable code variables data b Every program must have a function called main execu tion starts there main may call other functions These functions can be taken from one of the following sources i As written by the programmer in the same le as mam ii As written by the programmer in separate les iii Called from prede ned libraries c C has no concept of a Program as in Fortran just func tions d The C view is that main is called by the operating system the operating system may pass arguments to main as well as receive return values 69 return 0 however above we choose to make main take no arguments void Appendix B 2 PHYS 40075007 Computational Physics e is a comment and is ignored by the compiler the C standard says that comments cannot be nested Also comments must be terminated explicitly ie new line does not terminate them this is a common source of compiler errors and for the unwary can be very dif cult to trace For your own sanity the header comment should include i Name of the program to correspond to the le name ii Authors name even if copiedl iii Date iv Brief indication of function and purpose of the program 69 assignment number or project and what it does f Program layout is very important I suggest two spaces indentation for each block I suggest that you avoid tabs certainly avoid having the program pressed up against the right margin with all the white space on the left hand side g The include ltstdiohgt is a directive to the C preproces sor to include the contents of a le called stdioh stdioh contrains declarations of STanDard Input Output func tions include means that the full contents of the le stdioh are inserted where the include directive appears only then is the program passed to the C compiler proper include is C s way of importing h In the de nition of a function encloses the argument list in this case main expects no arguments void explic itly denotes has no arguments Donald G Luttermoser ETSU Appendix B 3 i enclose the statements in a function printf is a library function for printing output on the users screen in this case the string of characters between quotes n is C notation for newline printf will not insert a newline by itself k n represents a single character neccessary since typing the ltreturngt key in the middle of a string is not really practical Other control characters are t tab a alert bell double quote r carriage return backslash itself l 1 return 0 in C programs can return values to the operating system in UNIX 0 signi es that the program terminated normally other values indicate error conditions This is useful if you intend putting your program in a UNIX shell script likewise DOS bat les B Variables and Arithmetic 1 Convert degrees Fahrenheit to degrees Celsius 5 TC TF 32 Then print a list as follows 0 17 20 6 40 4 60 15 300 148 u ftocc Fahr to Celsius table 39gc 31089 Copied from KampR p9 include ltstdiohgt int mainvoid nquot fahr celsius fahr fahr step return 0 2 Dissection a All variables must be declared Usually at the beginning of their function b Built in types in C int integer can be 16 bits or 32 bits oat oating point char single text character really it is a small integer taking on values in 0255 short short integer possibly 8 bits maybe 16 long long integer 32 bits double double precision oating point more signi cant digits larger exponent c Assignment statement upper 300 1 While loop whiefahr lt upper executable statements body in here Operation of while Donald G Luttermoser E TSU Appendix B 5 i Condition is tested ii If condition is true body is executed iii Go back to l execution con tinues at the statement following the loop iv However if condition is false 3 Note the textual layout style used as mentioned earlier my sug gestion is to indent each block by 2 characters the Kernighan amp Ritchie KampR book on C The C Programming Language considered by many to be the C bible use the next tab but with that you very quickly get to the right hand side of the page on the other hand we have to be able see what is part of a block and what isn t 4 For me it is essential that the closing brace lines up with while 5 lnteger arithmetic causes truncation hence 59 gt 0 6 printf a b C d e Has multiple arguments the rst is always a string con stant then this may contain codes 6g d to say how subsequent arguments are to be printed d print the second argument as a decimal integer printf is not part of the C language but it s in the standard library which is available with all C compilers Field widths can be speci ed eg 6d Uses right justi cation unless speci ed otherwise Appendix B 6 7 8 a b PHYS 40075007 Computational Physics ln ftocc we have a problem with integer arithmetic not really satisfactory truncation 69 59 0 Therefore we need a oating point version of ftocc u ftoclc Fahr to Celsius table floating point version jgc 31089 Copied from KampR p12 include ltstdiohgt int mainvoid float fahr celsius int lower upper step lower 0 lower limit of table fahr upper 300 upper limit step 20 step size fahr whilefahr lt upper celsius 5090 fahr320 printfquot30f 61fnquot fahr celsius fahr fahr step return 0 lower Dissection Mixing of operands is allowed in C 69 fahr lower whiefahr lt upper but this does not mean that there is no type checking in the examples given the mts are con verted to floats before the operation is done In general you should be very careful with this feature Actually it might be better to make type conversion explicit this can be done with a type cast thus fahr floatower f conversion speci cation in printf i 30f 3 characters wide no decimal point no fraction digits Donald G Luttermoser ETSU Appendix B 7 ii 61f 6 characters wide decimal point plus 1 fractional digit c Other conversions i d print as decimal integer no width speci ed ii 6d decimal integer at least 6 chars wide iii f oating point no width spec iv 6f oat 6 chars wide v 2f 2 chars after decimal point width uncon strained vi X hexadecimal vii o octal viii s char string ix 00002 for itself C The for Statement 1 Another version of ftocc ftoc2c u ft062C Fahr to Celsius table demo of for statement jgc 31089 Copied from KampR p13 include ltstdiohgt int mainvoid int fahr forfahr0fahrlt300fahrfahr20 Appendix B 8 PHYS 40075007 Computational Physics printfquotd Z61fnquotfahr5090fahr32 return 0 2 Major differences a In printf celsius is replaced by a complex expression that evaluates to celsius this should be nothing new to those who have used Fortran the general rule is A variable of some type can be replaced by an arbitrarily complicated expression that evaluates to a value of that type b for statement and loop There are three statements con tained Within the in a for statement i fahr 0 initialization ii fahr lt 300 loop continuation control evaluate the condition if true execute the body other wise jump out iii fahr fahr 20 increment do this AFTER the rst and subsequent loops BEFORE evaluation of control condition iv The for body may be a single statement or a block D Symbolic Constants and the Preprocessor 1 The general form of a symbolic de nition is define ltsymboic namegt ltrepacement textgt Donald G Luttermoser ETSU Appendix B 9 2 The C preprocessor replaces all occurrences of ltsymboiC namegt with ltrepacement textgt just like an editor global replace com mand 3 Final version ftoc3c u ftocSc Fahr to Celsius table demo of Symbolic constants jgc 31089 Copied from KampR p15 include ltstdiohgt define LOWER 0 lower limit of table define UPPER 300 define STEP 20 int mainvoid int fahr float celsius for fahr LOWER fahr lt UPPER fahr fahr STEP celsius 5090 fahr32 printfquot3d 61fnquot fahr celsius return 0 4 Dissection a The preprocessor is an extremely important part of C When you compile a C program two things happen First the preprocessor runs through the le second actual com pilation Therefore in the example given above all oc currences of the string LOWER get replaced with 0 the compiler never sees the LOWER and LOWER is not a variable b Adopt a convention to use upper case for symbolic con stants be very careful what else you use upper case for Appendix B 10 PHYS 40075007 Computational Physics c There is no after the define statement 1 Nor does a define statement contain an E Character Input and Output 1 There is a standard input output library which deals primarily with streams of characters or text streams cf UNIX les a A text stream a sequence of characters divided into lines b A line zero or more characters followed by a newline character 2 The standard input output library MUST make each input output stream conform to this model no matter what is in the physical le 3 C getcharo read next input char from the standard input stream ie usually the keyboard 4 putcharc putprint char c on the screen standard output stream 5 There are equivalent functions for les 6 Buffered and Echoed lnput a On most computers input from the keyboard is echoed and bu ered b Echoed input When you type a character on the key board the computer input output system immediately Donald G Luttermoser ETSU Appendix B 11 echoes it to the screen terminal immediately means BE FORE it is presented to the reading program see buffered below Incidentally this means that the terminal it self does not display the typed character just what is echoed from the host computer c Buffered input While the user is typing the com puter stores all the typed characters in a buffer array and presents the array line of characters to the reading program ONLY AFTER A CARRIAGE RETURN En ter has been typed 7 Terminal input output and 1 0 Redirection a O and UNIX has a fairly uni ed view of text input output Reading from the keyboard is like reading from a le device called stdin However there are special stdin routines that hide that fact 6 9 getchar below is exactly equivalent to getcstdin Likewise putcharc is equiva lent to putcc stdout where stdout is the standard output device ie the screen b Incidentally putcc stdprn is a handy way of writing to a PO printer as is putcc stdaux to the auxiliary com munications port where stdprn is the standard printer device and stdaux is the standard auxiliary device BUT THESE ARE NOT STANDARD C UNIX c So when we talk about le l O below we include ter minal lO All the programs below can be tested using the terminal Note ltCtrgt z is EOE end of le for a keyboard PHYS 40075007 Computational Physics Appendix B 12 If you want to test the programs using les you can use 01 input output redirection i First you must get the compiler to compile to an exe le say Cioexe ii Second create a text le say testdat Then en ter at the Unix prompt Cio lt testdat which tells do to read from testdat gt the lt redi rects the program to read from the le instead of the keyboard iii If you want to send the output to a le say testoutdat use Cio lt testdat gt testoutdat File Copying Make a program with the following algorithm read a char while char is not an eof char output the char just read read a char 4 ciolc copy input to output version 1 gc J 31089 copied from KampR p16 include ltstdiohgt int mainvoid int c c getchar whilec EOF Donald G Luttermoser ETSU Appendix B 13 putcharc c getchar return 0 Dissection a denotes not equal to b int C instead of Char Cl gt getchar must be able to indicate errors 69 when there is no more input end of le EOF EOF is de ned in stdioh using define normally it requires more bits than a plain character 9 Another version using a common C cliche J ci02c copy input to output version 2 demo of C file reading clicheidiom c 31089 copied from KampR p17 include ltstdiohgt int mainvoid int c whilecgetchar EOF putcharc return 0 Dissection a whiecgetchar EOF This is quite a common id iomcliche in C it is made possible because assignment statements in C have a value the value is the last value assigned Experienced C programmers will all understand what it does However do NOT take this as general per mission to write dense cryptic code b Note the importance of brackets Appendix B 14 10 a b PHYS 40075007 Computational Physics Character Counting Objective count the characters in the input stream stop at EOF u ctch1c counts input chars version 1 jgc 31089 copied from KampR p18 include ltstdiohgt int mainvoid long nc nc0 whilegetchar EOF nc printfquotldnquot nc return 0 Dissection operator increment by one decrement by one Note As single statements nc and nc are equiva lent but give different values in expressions 69 int i n n 6 i n POST increment gives i 7 and n 7 the increment is done BEFORE the assignment wheras n 6 i n POST increment gives i 6 and n 7 long integer is at least 32 bits long d tells printf to expect a long Donald G Luttermoser E TSU 11 12 Appendix B 15 Line Counting 4 include ltstdiohgt count lines in input int mainvoid int 6 n1 n1 0 while cgetchar EOF ifc n nl printfquotdnquot n1 return 0 I u Dissection a if statement tests condition in if true executes statement or group following otherwise false skip them b means is equal to CAUTION in place of can be syntactically correct and so not trapped by compiler a nasty hard to nd error results c Character constants 69 n represents an integer value equal to the value of the character in the machine s charac ter set In ASCII n 10 decimal A 65 decimal Word Counting Word any sequence of chars that does not contain a white space char 216 blank tab newline u include ltstdiohgt define IN 1 inside a word define OUT 0 outside a word count words lines and chars in input Appendix B 16 PHYS 40075007 Computational Physics int mainvoid int 6 n1 nw me state state OUT n1 nw me O whilecgetchar EOFH nc ifc n nl ifc II c n II c t state OUT else ifstate OUT state IN nw printfquotd Zd dnquot n1 nW n6 return 0 Dissection a n nw nc 0 recall an assignment has a value te nc 0 has the value 0 b rneans OR ampamp rneans AND c Expressions with and ampamp are evaluated left to right evaluation is stopped as soon as the result is known e 9 int i1 j2 k3 ifi1 ll J3 k3 the evaluation of the condition can and does stop as soon as j3 evaluates to FALSE why 1 if expression statement else statement As usual statement can be a group in brackets Donald G Luttermoser E TSU 13 Appendix B 17 Arrays Below is a program to count the occurrences of each nu meric digit of white spaces and of all other characters together without using 12 named counters include ltstdiohgt counts digits whitespaces others int mainvoid int c i nwhite nother int ndigit10 nwhite mother 0 fori0 ilt10 i ndigiti O whilecgetchar EOF ifcgt 0 ampamp Clt 9 ndigit c O else ifc II c n II C t nwhite else nother printfquotdigits quot fori0 ilt10 i printfquot Zdquot ndigiti printfquot white space Zd other dnquot nwhite nother return 0 Dissection 3 b c d 8 int ndigit10 an array of 10 ints The for loops must go 0 l 9 since in C array subscripts must start at 0 The test if Cgt 0 ampamp ampamp denotes logical AND C 390 assumes that 0 l 9 have consecutive increasing values C 390 is an integer expression gt it s OK for a subscript Appendix B 18 PHYS 40075007 Computational Physics f ifCondition1 statement else if condition2 statement statement This is the model for a multiway decision you can have any number of else if cond stmt groups between rst if and nal else the nal else catches anything that didn t satisfy any of the previous condi tions g lndenting style Again please note that we don t want to run off the right hand edge of the paper F Functions 1 We have already encountered functions that have been written for us eg printf getchar putchar 2 Here is a program that reads ints and oats and does some arith metic with them We ll use it as a case study of many of the aspects of functions 4 math1c jgc 6196 demo of functions in same file include ltstdiohgt int getIntvoid int i Donald G Luttermoser ETSU Appendix B 19 scanfquotdquot ampi return i float getFloatvoid float f scanfquotfquot ampf return f float addffloat a float b float c c a b return c int main float x y z w int i j k printfquotenter first intquot i getInt printfquotenter second intquot j getInt kij printfquotd d dnquot i j k printfquotenter first floatquot x getFloat printfquotenter second floatquot y getFloat z x y printfquotf f fnquot x y z w addfx y printfquotf f fnquot x y w printfquotf f fnquot x y addfx y return 0 Dissection a A function de nition consists of Appendix B 20 PHYS 40075007 Computational Physics return type function nameparameter list if any de nitions of internal variables locals statements b The functions may be all in the same le along with main or may be spread across many les c The variables a b C are local to addf they are invisible elsewhere d a b c are called parameters or formal parameters or colloquially formats e In the call addfx y X y are called arguments or actual parameters or colloquially actuals f return ltexpresstongt returns a value to main and passes control back to main g Unless the function returns a value some don t return is not essential but if the function declaration indicates that the function returns a value then there must be an appropriate return statement h The calling function can ignore the returned value 3 Here is another version with different layout here we ve kept main at the beginning and functions at the end consequently we ve had to declare functions getlnt etc before they are called u math2c jgc 6196 demo of functions in same file now with functions after main and with declarations include ltstdiohgt Donald G Luttermoser ETSU Appendix B 21 must declare functions before calls int getIntvoid float getFloatvoid float addffloat a float b int main float x y z w int i j k printfquotenter first intquot i getInt printfquotenter second intzquot 39 getInt ki39 printfquotd d dnquot i j k c printfquotenter first floatquot x getFloat printfquotenter second floatquot y getFloat z x y printfquotf f fnquot x y z w addfx y printfquotf f fnquot x y w function call can replace variable addf returns a float value printfquotf f fnquot x y addfxy you can mix types if the function has been declared and if types are compatible here it is possible to coerce i to float essentially the compiler inserts a cast floati usually it is best for the programmer to do this explicitly printfquotd f fnquot i y addfiy return 0 int getIntvoid int i scanfquotdquot ampi Appendix B 22 PHYS 40075007 Computational Physics return i float getFloatvoid float f scanfquotfquot ampf return f float addffloat a float b return a b Dissection a int getlntvoid declares the type of getlnt in C par lance it is called the prototype for getlnt b getlnt has type void gt int c addf has type float float gt float 1 Parameter names in prototypes are neither signi cant nor neccessary But can provide good documentation e Prototypes and the way functions are declared are the biggest change between ANSI C and the earlier versions Easier for compiler to check errors 4 Normally it is a good idea to split a program into modules 69 the arithmetic functions we ve written may be tested an stable whilst the main program is subject to change it s nonsensical to have to re compile the functions each time so we ll put them in a separate le funsc 5 Actually there are bigger and better reasons for such modularity but we ll come to them later Thus math3c Donald G Luttermoser ETSU Appendix B 23 J math3c jgc 7196 demo of functions in separate file module funsc and with declarations in funsh include ltstdiohgt include quotfunshquot int main float x y z w int i j k printfquotenter first intquot i getInt printfquotenter second intquot j getInt kij printfquotd d dnquot i j k printfquotenter first floatquot x getFloat printfquotenter second floatquot y getFloat z x y printfquotf f fnquot x y z w addfxy printfquotf f fnquot x y w printfquotf f fnquot x y addfxy return 0 and funsc J funsc jgc 7196 demo of functions in separate file called by math3c include quotfunshquot int getIntvoid Appendix B 24 PHYS 40075007 Computational Physics int i scanfquotdquot ampi return i float getFloatvoid float f scanfquotfquot ampf return f float addffloat a float b return a b and we need funsh u funsh jgc 7196 demo of functions in separate file called by math3c include ltstdiohgt int getIntvoid float getFloatvoid float addffloat a float b Dissection a When using a library or external module you must get into the habit of producing a header le that con tains declarations of the functions and declarations of other things required by the functions eg symbolic con stants b If you do not include declarations the compiler reverts to pre ANSl habits it assumes implicitly that all functions return an int and that actual and formal arguments agree Donald G Luttermoser E TSU c d e f g 6 4 Appendix B 25 in this pre ANSl mode the line in math2c above printf d f f n i y addfiy would fail miserably Fortunately gcc GNU C compiler seems to issue warnings in such cases other compilers may not be so helpful gcc Wall is better Note that on all Unixes except Linux the standard C compiler executable is ran with the CC command On Linux the standard is the GNU C compiler which is ran with the gcc command The prototypes of standard library functions are contained in header les 69 stdioh is just an ordinary text le containing prototypes of printf etc in addition stdioh contains defines h les are Lot compiled they get compiled as part of any le in which they are included recall OFF the C PreProcessor executes all commands before the com piler proper When you have separate modules or compilation units obViously the compilation amp linking procedure must be modi ed The separate units can be compiled to object as before gcc C math3c gcc C funsc yields math3o gt funso Now link gcc 0 math math3o funso yields executable math Actually you can do all this in one line gcc 0 math math3c funsc Raise an integer m to a positive integer power n powerm n tstpow tester for power Appendix B 26 PHYS 40075007 Computational Physics jgc 31089 copied from KampR p 2425 include ltstdiohgt int powerint m int n int mainvoid int i fori0 ilt10 i printfquotd Zd dnquot i power2i power3i return 0 4 power raise base to nth power n not negative jgc 31089 copied from KampR p25 int powerint base int n int i p P1 fori1 iltn i p p base return p G Parameters i Pass by Value 1 Call by value means that the arguments are copied into tempo rary storage for the function essentially the formal parameters are local variables that get created then initialised to the value of the arguments thus above in the call power2 i rst a local int called base is created and initialised to value 2 likewise a local n gets created and initialised to the value that i has Contrast Subroutine function parameters in Fortran where the subroutine function has access to the original argument variable pass by reference Donald G Luttermoser ETSU Appendix B 27 3 We have no explicit call by reference in C must be done using pointers But conversely in C arrays are passed by reference and you cannot pass them by value 4 Another version of power to demonstrate call by value L power raise base to nth power n nonneg ver2 to show use of call by value int powerint base int n int p for p1 ngt0 n p p base return p Dissection Note that parameter n can be modi ed without affecting the argument n is local and like all local variables temporary H Character Arrays 1 To read a set of lines nd and print the longest Pseudo code while there is another line if it is longer than he current longest save it save its length print the longest line 2 Solution Break into functions following the pseudo code spec i cation 1 getlme fetch the next line work out its length 2 copy use for saving a line cannot assign strings in C L lilc find and print longest line jgc 41089 Copied from KampR p29 PHYS 40075007 Computational Physics Appendix B 28 include ltstdiohgt define MAXLINE 1000 size of buffer for line int getlinechar lineint maxline void copychar tochar from int mainvoid int len current line length int max max length so far char lineMAXLINE current input line char longestMAXLINE buffer for longest line max0 whilelengetlinelineMAXLINEgt0 iflengtmax maxlen copylongestline ifmaxgt0 there was at least one line printfquotslongest return 0 u getline read a line into s RETURN length int getlinechar s int lim int ci fori0iltlim1ampampcgetcharEOFampampC n i sic ifc n sic i Si 0 return i u copy copy quotfromquot into quottoquot void copychar tochar from int i Donald G Luttermoser E TSU i0 whileto i from i 0 i Dissection 3 b C d It is not neccessary to declare the length of s in getline storage has already been allocated in main void copy states explicitly that copy does not return a value allows compiler to check inconsistent usage Null terminated strings convention in C heon stored as heon 0 7chars NB the issue of buffer over ow is ignored like many of these programs we have not provided full bullet proo ng I External Variables and Scope 1 Local variables 69 Iine Iongest in main are private or LOCAL to main 2 The i in copy is unrelated to i in getline ie it has LOCAL SCOPE 3 Such variables are also called AUTOMATIC into existence when the function is called they disappear when the function is exited 4 Contrast STATIC to call Side issue STATIC when applied to external variables or functions also has another effect see KampR Chap 46 namely Appendix B 29 they only come still local but values are retained from call Appendix B 30 PHYS 40075007 Computational Physics that their names are invisible outside the le in which they are declared 5 Automatic assumed unless static int C C declared static 6 Declaration auto int c equivalent to int c gt hardly ever used since defaults to auto 7 External Variables a Similar to variables in a Fortran COMMON block b They are globally accessible GLOBAL SCOPE c Remain in existence permanently Their LIFETIME or DURATION or SPAN is for the duration of the execution of the program gt Scope and lifetime are important issues in programming languages 1 Thus externalsglobals can be used to communicate data between functions but there are all sorts of reasons why you should never use them mostly to do with reduction of coupling between software units e Rules of use of externals i An external variable must be DEFINED exactly once OUTSIDE any function ii It must be DECLARED in each function that ac cesses it and the declaration must announce that it is external Donald G Luttermoser ETSU Appendix B 31 extern int max Appendix Example B71 ic rewritten to use exter nal variables u lillc find and print longest line version 2 to demo extern variables jgc 41089 Copied from KampR p32 include ltstdiohgt define MAXLINE 1000 size of buffer for line external definitions follow NB outside any function int max max length so far char lineMAXLINE current input line char longestMAXLINE buffer for longest line int getlinevoid void copyvoid int mainvoid int len current line length extern int max extern char longest max0 whilelengetlinegt0 iflengtmax maxlen copyltgt ifmaxgt0 there was at least one line printfquotslongest return 0 u getline read a line into external line RETURN length int getlinevoid int ci extern char line Appendix B 32 PHYS 40075007 Computational Physics fori0iltMAXLINE1 ampampcgetcharEOFampampC n i lineic ifc n lineic i linei 0 return i l u copy copy extern line into extern longest void copyvoid int i extern char linelongest i0 whilelongestilinei 0 i Dissection i Externals max ine longest are DEFINED out side any function ii Externals are DECLARED before they are used in a function iii Exception to 2 if de nition appears in the source le before use then extern need not be stated ex plicitly but it is much better for documentation to state it explicitly iv As they are typed getline and Copy could be placed in separate les from main the LINKER would connect the external functions V The external data could likewise be de ned in a Donald G Luttermoser ETSU Appendix B 33 separate le which could be compiled and linked even though it only contains data vi Empty argument list must use explicit void If you de ned the function Copy as copy the compiler would assume pre ANSl C the major difference with pre ANSl was much more relaxed argument parameter consistency checking In fact there was no consistency checking gt the called function just assumed that the caller had passed the correct argument types so you had all sorts of nasties like integers being interpreted as oating point vii External variables are against most of the prin ciples of software engineering modularity informa tion hiding lack of coupling f I recommend that you adopt the following standard all uses of globals extern must be accompanied with a jus ti cation contained in a comment placed nearby it is surprising how with a little thought they can be completely avoided and you really appreciate it when it comes to testing J De nition Declaration 1 We will adopt the same convention as KampR for the use of these terms Appendix B 34 PHYS 40075007 Computational Physics 2 DEFINITION the variable is created or function body given including allocation of storage 3 DECLARATION The properties of the variable or function are announced to the program but no storage allocated K Creation of Executables 1 Basics Compiling and Linking a Source programs like helloc math3c etc are Lot directly executable There are a number of stages in creation of the executable and executing it b The rst stage of creating an executable is to compile helloc into object code gcc C helloc Doing this and puts the hello object code into a le helloo This object code is essentially machine code with code to initialize main a call to printf and code to return 0 But the code for printf is not present c The second stage is to link the machine code in helloo with appropriate object code for printf gcc 0 hello helloo This produces the executable le hello no extension 1 Even though we use gcc or CC gcc actually invokes a command called Id to do the linking loading loading refers to loading of external code e Here the extraction of the code for printf from a library is kept implicit however a library is nothing more than an object le with appropriate indexes to each function Donald G Luttermoser E TSU Appendix B 35 f Finally to execute hello you type g hello at the Unix prompt This reads the contents of hello into memory and starts execution at an appropriate start ad dress The situation may be made more clear cut if we look at the two module program math3c funsc introduced above Recall gcc C math3c gcc C funsc yields math3o gt funso gcc o math math3o funso yields executable math However you can do all this in one line gcc o math math3c funsc This process is diagrammed in Figure B 1 2 Other Libaries a Not all system functions are in the default libraries that are searched by gCCld For example math functions like sqrt take math4c if you try to link using gcc ansi pedantic Wall Wstrict prototypes o math math4c funsc The linker will complain that the call to sq rt is unresolved To correct this you must explicitly mention the maths library thus gcc ansi pedantic Wall Wstrict prototypes o math math4c funsc lm u math4c jgc 7196 extension of math3c but showing use of math library Appendix B 36 PHYS 40075007 Computational Physics library compiler gcc c math3c gt matho gcc c funs C gt funs ol V V V I linker V exe cut able math Figure B 1 Compiling and Linking Donald G Luttermoser ETSU Appendix B 37 function sqrt ie need 1 include ltmathhgt 2 put lm in linker command include ltstdiohgt include ltmathhgt include quotfunshquot int main float x y z w v int i j k printfquotenter first floatquot x getFloat printfquotenter second floatquot y getFloat w addfx y printfquotf f fnquot x y w make sure we have a nonnegative number w fabsw v sqrtw printfquotsquareroot of f is fnquot w v return 0 b Also libraries exist to carry out a variety of X Window functions 3 Static vs Shared Libraries The code for addf getlnt and getFloat is linked statically ie the object code for each of the functions is copied into the le math For common functions like printf this can become wasteful hence shared libraries in which the linker inserts just a pointer to shared code that held in the operating system or maybe somewhere on disk 4 Make a When you get to more complicated systems of multi module programs make can become useful Appendix B 38 b PHYS 40075007 Computational Physics Here is a make le for math4c makefile for math4 jgc 7196 here we define a macro CC as gcc the GNU Ccompiler CCgcc ansi pedantic Wall Wstrictprototypes how to create math the executable 1 math4o and funso must be up to date math math4o funso 2 create prog by linking math4o funso note must use TAB here spaces no good CC o math math4o funso lm specify what funso depends on don t forget the TAB for the cc command funso funsc funsh CC c funsc dependencies for math4o math4o math4c CC c math4c cleanup directory of object files invoke as make clean ie a separate shell command typically after you re finished developing amp executable math is OK note that the rm command is on the line after the clean target and as usual TABbed clean rm f funso math4o Donald G Luttermoser E TSU C d e f g Appendix B 39 The standard practice is to name this le Makefile then invoking make causes math to be made all programs compiled and linked Make takes are of dependencies 69 if math4o and funso already exist make will not recompile them On the other hand if they exist but if say math4c has changed since the last compilation to math4o then make will recompile math4c but not funcc In the example above make Clean will delete all the o les If you want to call the make le something else 239 6 math4m maybe you want to keep many make les in one directory though advanced make users have other solutions to that problem then you can invoke makem as make f makem Make is very much part of UNIX and is used for very much more than just C programs Another point because of its ddly syntax nobody ever creates a make le from scratch nd a template make le one that works that is close to what you require and change bits to suit Be careful with the places where the TAB syntax matters 5 Warnings for gcc and cc a b As I have said earlier gcc and CC without warnings switched on is positively dangerous I have suggested the following for either gcc or CCI Appendix B 40 c 01 PHYS 40075007 Computational Physics gcc ansi pedantic Wa Wstrict prototypes The following list explains the various switches on GNU CC 69 gcc pedantic pedanticerrors Wall enable a host of warning messages Wpointerarith issue warnings for nonANSI standard C issue errors for nonANSI standard C warn of dependency on quotsizeofquot function or void Wcastqual warn of pointer cast that removes type qualifier Wcastalign warn of pointer cast that increases alignment eritestrings make string constants type quotconst charquot and warn about nonconst char mix Wconversion warn when prototype causes different type conversion Waggregatereturn warn of functions returning structures Wstrictprototypes warn if argument types missing Wmissingprototypes warn if global function has no prototype Wredundantdecls warn of multiple declarations in same scope Wnestedexterns warn of externs within functions Winline warn if function cannot be inlined Werror make all warnings errors ansi ANSI standard C It is recommended that the following be used with gcc ansi pedantic Wa Wpointer arith Wcast qua erite stri ngs Wconversion Wno strict prototypes Wmissing prototypes Wnested externs Wi n i ne Wno error Jamie Blustein ltjamie csduwocagt has notes about us ing gcc at httpwwwcsduwocajamieRefsFootnotesgcc html Donald G Luttermoser ETSU Appendix B 41 L Interactions with the Operating System 1 Since Unix and C are closely related indeed Unix is written in C C was built to allow easy access to the operating system 2 The function int system const Char s which resides in ltstdibhgt allows a C program to pass a string to the operating system to be processed For instance system quotdatequot causes the Unix program date to be run it prints the date and time of day on the standard output system returns a system dependent integer status from the command executed In the Unix system the status returned is the value returned by exit 3 The function void exit int status also in ltstdibhgt causes nor mal program termination Open les are flushed open streams are closed and control is returned to the environment How sta tus is returned to the environment is implementation dependent but zero is taken as successful termination The values EXITSUC CESS and EXITFALU RE may also be used Note that if one used exit to terminate a program one would not use return in these cases 4 The function getenv const Char name in ltstdibhgt returns the environment string associated with name or NU LL if no string exists 5 There are various time functions de ned in lttimehgt to return time information from the operating system see Appendix B10 in Kernigham amp Ritchie Appendix B 42 6 PHYS 40075007 Computational Physics The function int remove const char filename ltstdiohgt re moves the named le It returns non zero if the attempt fails This is the function called by the Unix command rm The function int rename const char oldname const char new name ltstdiohgt changes the name of a le It returns non zero if the attempt fails This is the function called by the Unix command mv The following commands can give useful information about the processes that are running on the machine They are all in ltstdiohgt a int getpid and int getppid returns the process s ID and the parent s process ID b int chdir char pathname sets the process s current work The process must have execute permission from the directory mg directory to the directory stored in pathname to succeed chdir returns 0 if successful 1 if not c int nice int delta adds delta to a process s current priority value Only a super user may specify a delta that leads to a negative priority value Legal priority values lie between 20 and 19 the lower the value the larger amount of time the CPU spends on the process If a delta is speci ed that takes a priority value beyond a user s limit the priority value is truncated to the limit If nice succeeds it returns a new nice value otherwise it returns 1 Note that this can cause problems since a nice value of 1 is legal 1 int getuid and int getgid returns the user s and the user s group 1D PHYS40075007 Computational Physics Course Lecture Notes Section XI Dr Donald G Luttermoser East Tennessee State University Version 41 Abstract These Class notes are designed for use of the instructor and students of the course PHYS 40075007 Computational Physics taught by Dr Donald Lutternioser at East Tennessee State University XI Fourier Analysis and NonLinear Oscillations A Introduction 1 This section of the notes covers Fourier Transform Methods or spectral methods a For some problems the Fourier transform is simply an ef cient computational tool for accomplishing certain com mon manipulation of data b In other cases we have problems for which the Fourier transform or the related power spectrum is itself of intrinsic interest 2 A physical process can be described either in a time domain by values of some quantity f as a function of time t e 9 f t or else in a frequency domain where the process is speci ed by giving its amplitude F generally a complex number indicating phase also as a function of frequency V or w 27w where w is the angular frequency eg Fy or a In quantum mechanics the time domain is described by the wave function lx t as a function of both displace ment x and time t and the frequency domain is described by the amplitude function Ak where the independent variable k is typically a function of energy E which of course is related to the frequency of the particle b In optics the time domain is given by the flow of photons as wave trains and the frequency domain is given by the spectrum of these photons amplitude as a function of frequency Xl 2 PHYS 40075007 Computational Physics B Fourier Analysis 1 We will use the example of the Schrodinger equation from quan tum mechanics for our description of Fourier analysis In quan tum one often encounters wave functions that take the form 1 00 lx t foo Ak WW dk XI1 a At t 0 this equation takes on a form that may be familiar out 127T LAkeik1dk XI2 b Eq XI 2 reveals that the amplitude function Ak is the Fourier transform of the wave function lxt at t 0 gt the amplitude function is related to the wave function at t 0 by a Fourier integral 2 Fourier analysis the generation and deconstruction of Fourier series and integrals are the mathematical methods that underlies the construction of wave packets by superposition a Mathematicians commonly use Fourier analysis to rip func tions apart representing them as sums or integrals of sim ple component functions each which is characterized by a single frequency b This method can be applied to any function that is piecewise continuous ie that has at most a nite number of nite discontinuities In quantum mechanics wave functions must be continuous so as such satis es this condition gt prime candidates for Fourier analysis c Whether we represent f Via a Fourier series or Fourier integral depends on whether or not this function is pe Donald G Luttermoser E TSU 01 e f 3 Fourier Series Fourier series are not mere mathematical de vices they can be generated in the laboratory or telescope gt a spectrometer decomposes an electromagnetic wave into spectral lines each with a different frequency and amplitude intensity Thus a spectrometer decomposes a periodic function in a fashion riodic gt any function that repeats itself is said to be periodic More precisely if there exists a nite number L such that at L at then is periodic with period L We can write any function that is periodic or that is de ned on a nite interval as a Fourier series However if f is non periodic or is de ned on the in nite interval from oo to 00 we must use a Fourier integral analogous to the Fourier series a b c Suppose we want to write a periodic piecewise continuous function f as a series of simple functions Let L denote the period of f x and choose as the origin of coordinates the midpoint of the interval de ned by this period L 2 g at g L 2 If we let an and bn denote real expansion coe icients we can write the Fourier series of this function as fxao 1 x m an cos lt27TnZgt bn s1n lt27TnZgt XI 3 We calculate the coe icients in Eq XI 3 from the func tion f as 1 a0 Z LL52 at dx XI 4 XI 4 PHYS 40075007 Computational Physics an 2L2 x cos Zuni day 71 12 L L2 L XI 5 2 L2 53 bn LLLZ x s1n Zuni day 71 12 XI 6 1 Notice that the summation in Eq XI 3 contains an in nite number of terms In practice we retain only a nite number of terms gt this approximation is called trun cation i Truncation is viable only if the sum converges to whatever accuracy we want before we chop it off ii Truncation is not as extreme an act as it may seem If is normalizable then the expan sion coe icients in Eq XI 3 decrease in magni tude with increasing n ie lanl gt0 and lbnl gt0 asn gt 00 iii Under these conditions which are satis ed by physically admissible wave functions the sum in Eq XI 3 can be truncated at some nite maxi mum value mum of the index n Trial and error is typically needed to determine the value of mm that is required for the desired accuracy iv If f is particularly simple all but a small nite number of coe icients may be zero One should always check for zero coe icients rst before evalu ating the integrals in Eqs XI 4 5 6 Donad G Luttermoser ETSU 4 01 The Power of Parity One should pay attention as to whether one is integrating an odd or an even function Trigonometric functions have the well known parity properties odd even sin x sin 06 cos x cos 06 As such if f is even or odd then half of the expansion coef cients in its Fourier series are zero a If is odd fx then an0 71 12 L x 7 bn s1n lt27mfgt b If is even fx then If bn0 71 12 00 x XI10 ll x 721a cos lt27mfgt c If is either an even or and odd function it is then said to have de nite parity The Complex Fourier Series lf does not have a de nite parity we can expand it in a complex Fourier series a To derive this variant on the Fourier series in Eq XI 3 we just combine the coef cients an and bn so as to introduce the complex exponential function eigmxL 1122 Z OneiZNnxL TLOO XI11 b Note carefully that in the complex Fourier series in Eq XI ll the summation runs from oo to 00 The expan PHYS 40075007 Computational Physics sion coef cients on for the complex Fourier series are 1 L2 7i rrnz on ZL2fxe 2 Make XI12 Exercise Derive Eqs XI 11 and XI 12 and thereby determine the relationship of the coefficients on of the complex Fourier series of a function to the coefficients an and bn of the corresponding real series Fourier Integrals Any normalizable function can be expanded in an in nite number of sine and cosine functions that have in nitesirna y differing arguments Such an expansion is called a Fourier integral a A function f can be represented by a Fourier inte gral provided the integral 3 00 f dx exists gt all wave functions satisfy this condition for they are normalizable b The Fourier integral has the form mo Lgke dk which is the inverse Fourier transform XI13 c The function gk plays the role analogous to that of the expansion coef cients on in the complex series Eq XI ll The relationship of gk to f is more clearly exposed by the inverse of Eq XI 13 gov L mar dx which is the famed Fourier transform equation In rnath XI 14 ematical parlance f and gk are said to be Fourier transforms of one another 1 More precisely gk is the Fourier transform of at and is the inverse Fourier transform of Donald G Luttermoser E TSU e When convenient we will use the shorthand notation Wk f lx 0 and in 0 71 Ak XI 15 to represent Eqs Xl 14 and XI 13 respectively in the realm of quantum mechanics f Many useful relationships follow from the intimate rela tionship between and most important is the BesselParseval relationship L 39lede L 9k2dk XI16 C The Time Domain versus the Frequency Domain For our purposes the 1 In the time domain ht a physical process is described by some quantity h as a function of time t 2 In the frequency domain Hf the process is speci ed by giving its amplitude H generally a complex number indicating phase also as a function of frequency f where oo lt f lt oo 00 Essentially ht and H f are two different representations of the same function related by the transform equations Hf ht em dt ht L 1 H f Wt df a If t is measured in seconds then f is measured in HZ ls b If h is a function of position x in meters then H is a XI17 XI18 function of wavenumber k in cycles per meter 4 Sometimes angular frequency in radians second is used instead of frequency where w 2n f Then Eqs XI 17 and XI 18 are rewritten as L Mt e2 t dt L i 27F PHYS 40075007 Computational Physics XI19 Hw 67m dw XI 20 a To introduce symmetry between these two equations of ten the 1 27F coef cient is split between the two integrals introducing a 1 27F coef cient to each equation as done in Eqs Xl l21314 b For this section we will follow the f notation 5 can be made If As we have seen the following statements about these functions Mt is real Mt is imaginary Mt 39 Mt is odd M M M t is real and even t is real and odd t is imaginary and even H Hf te even Hf 239e odd is real and even is imaginary and odd is imaginary and even Mt is imaginary and odd H f is real and odd These symmetries will be useful in order to develop computa tional ef ciency in coding 6 Mt ltgt Mat 1 Whatb Mt to gt gt gt 727rifot Hm 1 H lal HUM HO HO Useful scalings and shifting equations XI21 XI22 no scaling fa time scaling frequency scaling XI 23 27rifto fa time shifting XI 24 freq shifting XI 25 Donald G Luttermoser E TSU D Fourier Transform of Discretely Sampled Data 1 In the most common situations function ht is sampled ie measurements taken at evenly spaced intervals in time a Let T denote the time interval between consecutive sam ples such that hn hnT n 3 2 10123 XI 26 Often T is called the sampling interval b The reciprocal of the time interval is called the sampling rate If T is measured in seconds then the sampling rate is measured in HZ cycles per second 2 Sampling Theorem and Aliasing a For any sampling interval there is a special frequency f5 called the Nyquist critical frequency given by fa E i XI27 i If a sine wave of the Nyquist critical frequency is sampled at its positive peak value then the next sample will be at the negative trough value the sample after that at the positive peak again and so on ii Expressed otherwise Critical sampling of a sine wave is two sample points per cycle iii One frequently chooses to measure time in units of the sampling interval T In this case the Nyquist critical frequency is just the constant 12 XI 10 PH Y5 4007 5007 Computational Physics b The Nyquist critical frequency is important for two dis tinct reasons the rst good the second bad i The sampling theorem If a continuous function ht sampled at interval 739 happens to be band width limited to frequencies smaller than f5 then the function ht is completely determined by its sample hn Explicitly hm T i0 hnsin27rfct 717 W700 7T0 n7 XI 28 ii Sampling a continuous function that is Lot band width limited to less than the Nyquist critical fre quency will miss information outside then range of fc lt f lt fc gt this is called aliasing Any fre quency component outside of the frequency range f5 f5 is aliased te falsely translated into that range by the very act of discrete sampling The ef fects of this are shown in Figure XI l 3 Discrete Fourier Transform a Suppose we have N consecutive sampled values thMtk tkEkT k012N 1XI29 For description purposes let s assume that ht is an even function b Now determine estimates of the frequency from fc to fc at the distinct points de ned by n N N n27 77 Xl 3 f 7 n 2 7 7 2 The extreme values of n correspond to the lower and upper limits of the Nyquist critical frequency range Donald G Luttermoser ETSU XI 11 Figure XI 1 The continuous function ht shown in a is sampled every 739 seconds The actual Fourier transform of this function7 Hf7 is shown If we ignore the points outside the Nyquist critical frequency7 an chased and incorrect Fourier transform is computed as shown in X1 12 PH Y5 4007 5007 Computational Physics c The discrete values for H f are now determined with 00 N71 N71 ht62mf tdt Z hkeZTrzfntkT 7 Z hkeZTrzknN 00 k0 k0 XI31 d The summation shown in Eq XI 31 is called the dis crete Fourier transform of the N points hk If we de ne Hn with then have N71 Hn E Z hkem mN XI32 k0 Eq XI 31 can then be written as where fn is given by Eq XI 30 i Remember that N is a measure of the period of the function H Therefore Hn H NW with n 1 2 ii To simplify the calculation one generally lets the n in Hn vary from 0 to N 1 one complete period Then n and k in hk vary exactly over the same range iii When this convention is followed you must re member that zero frequency corresponds to n 0 positive frequencies 0 lt f lt fc correspond to val ues 1 g n g N 2 1 while negative frequencies fclt flt 0correspondto N21 3 ng N 1 The value N 2 corresponds to both f fc and f fc Donald G Luttermoser E TSU e The discrete inverse Fourier transform then takes the form 1 N71 hk 7 Z HneimnN XI34 N n0 f Finally in analogy to Eq XI l6 we can write the dis crete form of Parseval s theorem as N71 2 1 N71 2 Z Ihkl N Z lHnl XI35gt k0 n0 E Fast Fourier Transform 1 We can ask the question how much time is required to compute a function of the form N71 Hn Z Wnkhk k0 XI36 where the vector of Ms is multiplied by a matrix whose 71 Mm is the constant W to the power n X k with W given by W E chiN XI37 hence a Fourier series expression a This matrix multiplication requires N 2 complex multipli cations plus a smaller number of operations to generate the required powers of W b As such discrete Fourier transforms appear then to be a 9N2 process 2 We can speed the calculations up to order 9Nlog2 N opera tions with an algorithm called the Fast Fourier Transform or FFT for short a The difference between N log2 N and N 2 calculations is immense XI 13 X1 14 PHYS 40075007 Computational Physics b With N 106 this corresponds to 003 seconds and 20 minutes on a 1 GHZ processor respectively 3 FFTs were rst developed for computational coding in the mid 1960s by J W Cooley and J W Tukey The earliest discoveries of the FFT was made by Danielson and Lanczos in 1942 The DanielsonLanczos Lemma is as follows a A discrete Fourier transform of length N can be rewrit ten as the sum of two discrete Fourier transforms each of length N 2 b One is formed from the even numbered points of the orig inal N the other from the odd numbered points The mathematical proof is Fk eQTrijkij j0 N271 I N271 I 2wzk2Nf2j Z 2wzk21Nf2j1 j0 j0 N27 1 I N271 Z 27rzkjN2f2j Wk Z emiwN2 f2 1 139 F0 F WkF XI38 O c Note that k in the equations above varies from 0 to 71 not just to N2 Never the less the transforms F the even sum and F k0 the odd sum are periodic in k with length N 2 As such each is repeated through 2 cycles to obtain Fk d For this to be the most effective N should be an integer multiple of 2 If it is not one should pad the vectors with zeros until the next power of 2 is reached Donald G Luttermoser E TSU 4 X1 15 Virtually all mathematics software packages have FFTs built into the command structure or available in a math library FFTs are most often used in the convolution and deconvolution of spectral data PHYS40075007 Computational Physics Course Lecture Notes Section V Dr Donald G Luttermoser East Tennessee State University Version 41 Abstract These Class notes are designed for use of the instructor and students of the course PHYS 40075007 Computational Physics taught by Dr Donald Lutternioser at East Tennessee State University V Error Analysis and Uncertainties A Errors and Uncertainties 1 Types of Errors a Blunders Can control these with sleep gt look for b C inconsistencies in data points i Typographical errors in program or data ii Running wrong program iii Using wrong data le iv Using wrong equations for analysis Random Errors Can t control these gt harder to de tect i Electronic uctuations due to power surges ii Cosmic ray damage on detectors and or chips iii Somebody pulled the plug Systematic Errors Very hard to detect i Faulty calibration of equipment ii Bias from observer or experimenter iii These must be estimated from analysis of exper imental conditions and techniques iv In some cases corrections can be made to data to V 2 PHYS 40075007 Computational Physics compensate for these errors where type and extent of error is known V In other cases the uncertainties resulting from these errors must be estimated and combined with uncertainties from statistical uctuations d Approximation Errors i These arise from simplifying the mathematics so that the problem can be solved on a computer gt in nite series being replaced by a nite sum 00 xn I 7 V l 6 20 quotl N mm 22 EUR e1 any N v2 ii In Eq V 2 506N is the total absolute er TOT iii This type of error also is called algorithmic error the remainder or truncation error iv To have a small approximation error N gt x in Eq V Z e Roundoff Errors i lnaccuracies in stated numbers on a computer due to a nite typically small number of bits being used on a computer to represent numbers ii The total number of possible machine numbers Donald G Luttermoser ETSU V 3 is much less than the total number of real num bers which of course is in nite iii The more calculations a computer does the larger the roundo error iv This error can cause some algorithms to become unstable v If the roundo error value gt value of the number being represented gt the result is referred to as garbage vi For instance a computer may compute the fol lowing as 2 06666666 06666667 00000001 y 0 vii We need to worry about signi cant gures on a computer 0 On a 32 bit CPU single precision is stored in w word gt 4 bytes REALM or REAL o On a 32 bit CPU double precision is stored in M words gt 8 bytes REAL8 Trailing digits beyond these bytes are lost if you ask a code to print more digits out the printed digits should be considered as garbage f What kind of error do we make in representing a real number in a oating point number system V 4 PHYS 40075007 Computational Physics i Absolute Error true value approximate value ii Relative Error true value approximate value true value Relative error is not de ned if the true value is zero 2 Accuracy versus Precision a Accuracy is how close an experiment comes to the true value i It is a measure of the correctness of the result ii For an experimenter it is a measure of how skilled the experimenter is iii For a programmer it is a measure on how good they are at programming and the assumptions used in the algorithm b Precision of an experiment is a measure of how exactly the result is determined without reference to what the results means i It is a measure of the precision of the instruments being used in the experiment ii The precision of an experiment is dependent on how well we can overcome or analyze random er TOYS iii In programming it is a measure of how many bits are used to store numbers and perform calculations Donald G Luttermoser ETSU 3 Uncertainties 3 b C The term error signi es a deviation of the result from some true value i However since we often cannot know the true value of a measurement prior to the experiment we can only determine estimates of the errors inherent to the experiment ii The difference between two measurements is called the discrepancy between the results iii The discrepancy arises due to the fact that we can only determine the results to a certain uncer tainty There are two classes of uncertainties i The most prominent type Those which result from uctuations in repeated measurements of data from which the results are calculated ii The secondary type Those which result from the fact that we may not always know the appropriate theoretical formula for expressing the result Probable Error The magnitude of the error which we estimate we have made in our determination of the results i This does not mean that we expect our results to be wrong by this amount ii Instead it means that if we are wrong in our re sults it probably won t be wrong by more than the probable error V6 PH Y5 4007 5007 Computational Physics iii As such probable error will be synonymous with uncertainty in a measurement or calculation 4 Implications for Numerical Computing 3 0 Most numbers that we use in oating point computations must be presumed to be somewhat in error They may be rounded from the values we have in mind by i lnput conversion errors ii lnexact arithmetic from the limited signi cance of the input numbers due to limited RAM sizes of the computer How is the value of function f affected by errors in the argument x i Suppose x has a relative error of 6 so that we ac tually use the value 061 6 ii Then the value of the function is f06l iii If f is differentiable the absolute error in the value of f caused by the error in x can be approx imated by x 606 m emfx V B iv The relative error is then f06 06 x fx m 606 V 4 M M l l v As we had above suppose 61 then the ab solute error is approximately exex and the relative error is 661 Donald G Luttermoser ETSU vi These errors can be serious if x is large Thus even if the routine EXP were perfect there could be serious errors in using EXPY for ex when x is large vii Another signi cant example is cosine 06 when x is near 7T 2 The absolute error is approximately 7r exs1nx m e 1 viii This is not troublesome but the relative error is very much so sinx 7T 1 606 6 7 7 2 0 cosx Very small changes in 06 near 7T 2 cause very large relative changes in cos x gt we say the evaluation is unstable there The accurate values cos 157079 063267949 x 10 5 cos 157078 163267949 x 10 5 demonstrate how a small change in the argument can have a profound effect on the function value c An important part of computational physics is deriving solutions to problems in terms of series When one is working to high accuracies or if the series converges slowly it is necessary to add up a great many terms to sum the series V7 V 8 PHYS 40075007 Computational Physics B Useful Theorems in Computational Physics 1 Intermediate Value Theorem Let ba a continuous func tion on the closed interval 1 b If for some number oz and for some 061062 6 1b we have fx1 3 oz 3 fx2 then there is some point c E 1 b such that oz f C V5 Here the notation 1 b means the interval consisting of the real numbers 06 such that 1 S 06 S b 2 Rolle s Theorem Let be continuous on the closed nite interval 1b and differentiable on 1b lf f1 fb 0 there is a point c E 1 b such that fC 0 V5 Here the notation 1 b means the interval consisting of the real numbers at such that 1 lt 06 lt b 3 MeanValue Theorem for Integrals Let 906 be a non negative function integrable on the interval 1 b If f is con tinuous on 1 b then there is a point c E 1 b such that rwwmwmfwx mmdx Va more to come in Vll of the notes 4 MeanValue Theorem for Derivatives Let be continu ous on the nite closed interval 1 b and differentiable on 1 b Then there is a point c E 1 b such that fo wamp more to come in Vll of the notes Donald G Luttermoser ETSU 5 Taylor s Theorem with Remainder Let have the continuous derivative of order n l on some interval 1 b con taining the points x and 060 Set x fxo wm xo 39ltngt fxo WW 2 06 me n x me Ram V9 Then there is a number 0 between 06 and 060 such that i fn1c n1 Rn1 i x0 Let f be a continuous function on the nite closed interval 1 b Then f assumes its maximum and minimum values on 1 b 16 there are points 061 062 E 1 b such that f061 S 1 S 1W2 V41 for all 06 E 1 b Integration by Parts Let and 996 be real values func tions with derivatives continuous on 1 b Then Abf tgt dz ftgt 2 fogt dt V42 Fundamental Theorem of Integral Calculus Let be continuous on the interval 1 b and let Fm f mm for all x E1b v13 Then F is differentiable on 1 b and V 14 V9 V 10 PHYS 40075007 Computational Physics C The Mathematics of Errors and Uncertainties 1 Subtractive Cancellation a For the following let the actual numbers be represented by unmarked variables and those on the computer be des ignated with a c subscript b The representation of a simple subtraction is then ab Cgt acbcCCJ ac b1 61 01 EC V 16 gt 1 6129 EEC a a a where the error of the number variable is given by computer number actual number 6 V l8 actual number gt c From Eq V l7 the average error in a is a weighted average of the errors in b and c d We can have some cases however when the error in 1 increases when b m 0 because we subtract off and thereby lose the most signi cant parts of both numbers gt this leaves the least signi cant parts e If you subtract two large numbers and end up with a small one there Will be less signi cance in the small one i For example if a is small it must mean that b m c and so 1ea V49 1 6a m gm EC V ZO Donald G Luttermoser ETSU V 11 ii This shows that even for small relative errors in b and c the uncertainty in ac will be large since it is multiplied by ba remember b m 6 gtgt 1 iii This subtraction cancellation pokes its ugly head in the series for 6 2 Multiplicative Errors a Error in computer multiplication arises in the following way abgtltcgtacbcgtltcc V 21 i 16bl c a i 1ea m l 61 65 V 22 Since 61 and EC can have opposite signs the error in ac is sometimes larger and sometimes smaller than the individ ual errors in be and cc b Often we can estimate an average roundoff error for a series of multiplications by assuming that the computer s representation of a number differs randomly from the ac tual number c In these cases we can use the random walk technique i Let R be the average total distance covered in N steps each of length 7 then R m W7 v 23 ii Each step of a multiplication has a roundoff error of length em the machine s precision d 8 PHYS 40075007 Computational Physics iii In the random walk analogy the average relative error em arising after a large number N of multi plicative steps is then em xNem V 24 We will be making use of Eq V 24 many times in this course When roundoff errors do not occur randomly careful anal ysis is needed to predict the dependence of the error on the number of steps N i If there are no cancellations of error the relative error may increase like N em ii Some recursive algorithms where the production of errors is coherent 69 upward recursive Bessel functions the error increases like N l em This is something to keep in mind when you hear about computer calculations requiring hours to complete i A fast computer may complete 1010 oating point operations per second Hence a program running for 3 CPU hours performs about 1014 operations ii Then after 3 hours we can expect roundoff errors to have accumulated to a relative importance of 107 em iii For the error to be smaller than the answer this demands that em lt 104 Donald G Luttermoser ETSU iv As such this long of a calculation with 32 bit arithmetic hence inherently possesses only 6 to 7 places of precision probably contains much noise 3 De nitions from Statistics and Probability Theory a b C The mean u of the parent population e measure ments is de ned as the limit of the sum N determinations x of the quantity 06 divided by the number N determina tions 1 N a E N V 25 The mean is therefore equivalent to the centroid or aver age value of the quantity x The median 112 of the parent population is de ned as that value for which in the limit of an in nite number of determinations x of the quantity x half of the observa tions will be less than the median and half will be greater than the median i In terms of the parent distribution this means that the probability P is 50 that any measure ment x will be large or smaller than the median Px 3 MW 13062 2 MW 50 V 26 ii Much computer time is wasted by guring out median values As such fast sorting routines have been developed in many programming languages 69 SORT in IDL Then the median is just the element in an array that is at the midway point The most probable value pm of the parent popula tion is that value for which the parent distribution has V 14 PHYS 40075007 Computational Physics greatest value Pom 2 m um V27 d The deviation 12 of any measurement mi from the mean u of the parent distribution is de ned as i Deviations are generally de ned with respect to the mean for computational purposes ii If u is the true value of the quantity then 12 is the true error in 062 iii The average deviation 8 for an in nite number of observations must vanish by virtue of the de nition of the mean see Eq V 25 lim d lim N 400 N 400 1 N 1310ltNn21xigt MM M0 iv The average deviation oz therefore is de ned as the average of the magnitudes of the deviations which are given by the absolute values of the devi ations 1 N o E NT ugt V 29 v The average deviation is a measure of the disper sion of the expected observations about the mean Donald G Luttermoser ETSU e The variance 02 is de ned as the limit of the average of the squares of the deviations of the mean 1 N 2 1 N 2 2 NT06i H gt lejgo H V BO i The standard deviation 0 is the square root of the variance ii The standard deviation is thus the root mean square of the deviations where we de ne the root mean square to be the square root of the mean or average of the square of an argument iii In computational work the standard deviation 0 is considered an appropriate measure of the w taiinty of a measurement or a calculation 4 Errors in Algorithms An algorithm is often characterized by its step size h or by the number of steps N it takes to reach its goal If the algorithm is good it should give an exact answer in the limit h gt 0 or N gt 00 Here we present methods for determining the error in your code a 0 Let s assume that an algorithm takes a large number of N steps to get a good answer and that the approximation error approaches zero like oz 6313er z W In Eq V Bl oz and are empirical constants that would change for different algorithms and may be constant only for N gt oo V 16 PH Y5 4007 5007 Computational Physics c Meanwhile the roundoff error keeps accumulating as you take more steps in a random fashion following em 2 W em v 32 where em is the machine precision see Eq V 24 d The total error is just the sum of these two errors 6tot 6apprx 67 0 7 OZ 2 W V N m e Say we have a test problem where an analytic solution ex ists By comparing a computed solution with the analytic solution we can deduce the total error em i If we then plot loget0t against logN we can use the slope of this graph e the power of N to de duce which term in Eq V 34 dominates the total error gt if the slope 12 roundoff error domi nates if not the slope and the approximation error dominates ii Alternatively we can start at very large values of N where there should be no approximation error and continuously lower N to deduce the slope in the approximation error 5 Optimizing with Known Error Behavior For this section we will assume that the approximation error has oz 1 2 SO 1 Eappm Z W v35 a The extremum occurs when demtdN 0 which gives 1 2N73 ENilgem 0 7 Donald G Luttermoser ETSU V 17 N lZem 4N 3 N52 7 V36 b For single precision on a 32 bit chip processor em 2 10 7 so the minimum total error is reached when 4 N52 N 10 7 N z 1099 and 1 Etot 2 W m m 8 x 104 33 x 10 2 4 x 106 c This shows for a typical algorithm most of the error is due to roundoff Also even though that this is the minimum error the best we can do is to get some 40 times machine precision in single precision 1 As such to reduce roundoff error we should come up with a code that requires less steps N to achive convergence 6 Sections 312 and 313 of your textbook go into further details on optimizing error behavior and setting up an empirical error analysis You should look these sections over in detail D Propagation of Uncertainties 1 We will demonstrate the propagation of errors with a simple ex ample Suppose we wish to nd the volume V of a box of length L width W and height H a We measure Lo We and H0 and determine lo Lo We Ho b c PH Y5 4007 5007 Computational Physics How do the uncertainties in L0 W0 and H0 affect the uncertainty in V0 Let L W H and V be the actual to true value then AL LO L AW WO W and AH Ho H The error in V is approximately the sum of the products of the errors in each dimension times the effect that di mension has on the nal value of V 8V 8V 8V of for our example AV WOHOALLO Ho AWLOWO AH V 39 Dividing this equation by Eq V 37 gives AV AL AW AH V0 LO W0 H0 V40 2 In general we do not know the actual errors in the determina tion of any of the parameters The following describes how we estimate the uncertainties 3 b c 01 Let s de ne an fuv V41 Assume the most probable value for x is given by T op V42 lndividual results for x are found by individual mea surements of other parameters that is u 1 giving 06 fuZv V43 Using Eq V BO the variance in x is then 02 lim i296 m v44 Naoo Donald G Luttermoser ETSU e Following Eq V 38 the deviation for measurement in x is 806 806 i i i 7 7 i 7 7 06 06 u ultaugtv vltavgt f Combining these two equations gives 1 8x 806 2 2 i 7 i 7 7 7 I 19310 N Z llul Bu l U2 U av l l 806 2 806 2 7 2 7 7 2 7 uZ u Bu v2 1 av 2uz vi U 203 3 3 v 46 where a 2 gig 2 uz mm mi V47 g If measurements in u and v are uncorrelated one should get as many negative values as positive values for the terms in the series of Eq V 47 As such this summation in Eq V 47 gt 0 So 806 2 806 2 27 2 7 2 7 ax7aultaugt 0v BU V48 3 Speci c Formulas a Addition and Subtraction Let x be the weighted sum of u and v 06 W i be V 49 i Then the partial derivatives are simply the weight ing coef cients which are constant 37 a 279 ib v50 b C PH Y5 4007 5007 Computational Physics ii Eq V46 then becomes 2 2 Oz a 03 bgag 2ab03w V 51 Multiplication and Division Now let x be the weighted product of u and v V 52 06 iauv i The partial derivatives of each variable contain the values of the other variable 8x 806 am erm ii Eq V46 yields a 1212203 61211203 Zaguvaiv a V54 which can be expressed more symmetrically as 2 2 2 2 J2 u U UU x2 V 55 iii Similarly if x is obtained through division x i V 56 v the variance for x is given by 2 2 2 2 a a a a JJJ 2 u V 57 062 u2 v2 uv Powers Let x be obtained by raising the variable u to a power 06 auib V 58 i The partial derivative of x with respect to u is iabuib 1 ibg Bu V 59 Donald G Luttermoser ETSU ii Eq V46 becomes amp b7 06 u V60 d Exponentials Let x be obtained by raising the natural base to a power proportional to u 06 aeib V 61 i The partial derivative of x with respect to u is 806 iabeibu ibx V62 Bu ii Eq V46 becomes a i u V 63 06 b0 e Logarithms Let x be obtained by taking a log of u 06 alnibu V 64 i The partial derivative of x with respect to u is e 3 ii Eq V46 becomes a I 4 V66 a au which is essentially the inverse of Eq V 63 PHYS40075007 Computational Physics Course Lecture Notes Section VIII Dr Donald G Luttermoser East Tennessee State University Version 41 Abstract These Class notes are designed for use of the instructor and students of the course PHYS 40075007 Computational Physics taught by Dr Donald Lutternioser at East Tennessee State University VIII Matrices and Solutions to Linear Equations A Introduction Setting Up the Problem 1 There may be times when you have a system of N linear equations with N unknowns 111061 112062 ame b1 VIII 1 121061 122062 angN b2 VIII 2 aN11 1N2532 39 39 39 aNNxN bN VIII 3 a In many cases the a and b values are known so your problem is to solve for all of the x values b To solve this problem we must set the problem up as a matrix equation 111 112 MN 061 b1 121 122 a2N 062 i 52 VIII 4 am am INN 06N 5N AX B VIII5 c The solution for the X vector is then found by inverting the A matrix mle A lB VIII6 X A lB VIII7 2 Before discussing the techniques for carrying out such an inver sion we need to go over some basic linear algebra and discuss various types of matrices that one might encounter in physics VIII 2 PH Y5 4007 5007 Computational Physics 3 Also since matrices play a big role in quantum mechanics we will use the formalism that is used in QM to describe vectors and matrices B Linear Algebra 1 In classical mechanics vectors are typically de ned in Cartesian coordinates as A A a aim oxyy azz VIII8 Note that one also can use the 239 j k notation for the unit vectors a b c 2 A vector space consists of a set of vectors agt W W together with a set of real or complex scalars a b c Vectors are added via the component method such that a i r a i my a i 5sz 012 i 52 VIII9 However in quantum mechanics often we will have more than 3 coordinates to worry about indeed sometimes there may be an in nite amount of coordinates As such we will introduce a new notation the so called bra and ket notation to describe vectors la 0 Note that the gtlt in the bra de nition means take the complex conjugate multiply all 239 1 terms by l in vector oz ket bra VIII10 a 1 which are subject to 2 operations a Vector addition The sum of any 2 vectors is another vector la W 7gt VIII11 V Donald G Luttermoser ETSU VIII 3 i Vector addition is commutative la W W lagt VIII12 ii Vector addition is associative la l M la W 7gt VIII13 iii There exists a zero or null vector 0 with the property that la I0 la VIII 14 for every vector a iv For every vector a there is an associated in verse vector a such that la I oz VIII 15 b Scalar multiplication The product of any scalar with any vector is another vector ala VIII 16 i Scalar multiplication is distributive with respect to vector addition altagt If ala aw VIII 17 and with respect to scalar addition a ba ala bla VIII 18 ii It is also associative aba aba VIII 19 VIII 4 PH Y5 4007 5007 Computational Physics iii Multiplications by the null and unit vector are magt 0t Hagta VHL2m Note that a 1agt c A linear combination of the vectors a W W is an expression of the form ala bl cly VIII 21 i A vector is said to be linearly independent of the set a W W if it cannot be written as a linear combination of them eg unit vectors 1 Q and ii A collection of vectors is said to span the space if every vector can be written as a linear combination of the members of this set iii A set of linearly independent vectors that spans the space is called a basis gt 9339 and 2 or i jk de ne the Cartesian basis which is a 3 D orthogonal basis iv The number of vectors in any basis is called the dimension of the space Here we will introduce the nite bases analogous to unit vectors lelgt7 62gt7 7 lengt7 of any given vector IMMMmenmm wman which is uniquely represented by the ordered n tuple of its components la lt gt a1 a2 an VIII 24 Donald G Luttermoser ETSU VIII 5 v It is often easier to work with components than with the abstract vectors themselves Use whatever method to which you are most comfortable 3 ln 3 dimensions we encounter 2 kinds of vector products the dot product and the cross product The latter does not generalize in any natural way to n dimensional vector spaces but the former does and is called the inner product a The inner product of 2 vectors agt and is a com pleX number which we write as mm with the following properties lt lagt lt0z gt VIII25 ltalagt 2 0 VIII26 ltaagt0 ltgt agt0gt VIII 27 04bl gtclvgt blta gtcltoz7gt VIII28 ltal gt glam VIII 29 b A vector space with an inner product is called an inner product space c Because the inner product of any vector with itself is a non negative number Eq VIII 26 its square root is real we call this the norm think of this as the length of the vector a E 0401 VIII 30 d A unit vector whose norm is l is said to be normalized e Two vectors whose inner product is zero are called or thogonal gt a collection of mutually orthogonal nor malized vectors dzlog 62 VIII 31 VIII 6 PH Y5 4007 5007 Computational Physics is called an orthonormal set where 61 is the Kro necker delta f Components of vectors can be written as 12 clla VIII 32 g For vectors that are co linear and proportional to each other the Schwarz inequality can be applied to these vectors ltCv gt2 S ltaagtlt l gt VIII33 and we can de ne the complex angle between a and W by the formula al Wla ltalagtlt l gt39 cos 6 VIII 34 4 A linear transformation T the hat on an operator from this point forward will imply that the operator is a linear transfor mation don t confuse it with the hat of a unit vector takes each vector in a vector space and transforms it into some other vector a gt o Toz with the proviso that the operator is linear Taoz m ama bT VIII35 a We can write the linear transformation of basis vectors as fled ETz jlez39 j 172aquot7n39 VIII35 i1 This is also the de nition of a tensor as such the oper ator T is also a tensor b If a is an arbitrary vector la a1le1 anlen Z aj jgt VIII 37 Donald G Luttermoser ETSU o o 8 then no W TL 1 zwmw F 121 J VIII 38 T takes a vector with components 11 a2 an into a vec tor with components a Z Ejaj j1 If the basis is orthonormal it follows from Eq VIII 36 that Ti eilTleJ VIII40 or in matrix notation T11 T12 quot39 Tln T T T quot39 T VIII41 T791 ng Tim The sum of 2 linear transformations is 3 Tagt 5p Tloz VIII42 or again in matrix notation U S T ltgt Uij 317 TH VIII 43 The product of 2 linear transformations is the net effect of performing them in succession rst T the S In matrix notation U ST ltgt Uik Z SlyTM VIII 44 j1 this is the standard rule for matrix multiplication to nd the ikth element of the product you look at the ith row of S and the 16th column of T multiply corresponding entries and add 137 lei Vlll 7 VIII 8 PHYS 40075007 Computational Physics f The transpose of a matrix is the same set of elements in T but with the rows and columns interchanged T11 T21 39 39 39 Tnl T12 T22 39 39 39 Tn T VIII45 Tln T2n 39 39 39 Tnn Note that the transpose of a column matrix is a row ma trixl g A square matrix is symmetric if it is equal to its trans pose re ection in the main diagonal upper left to lower right leaves it unchanged it is antisymmetric if this operation reverses the sign SYMMETRIC T T ANTISYMMETRIC T T VIII46 h The complex conjugate T is obtained by taking the complex conjugate of every element Th T1962 Tin 1i x x Ti T21 T22 T2n a2 aquot 7 VIII47 ng Tquot a Tril TLTL i A matrix is real if all its elements are real and imaginary if they are all imaginary REAL Tquot T lMAGlNARY Tquot T VIII 48 A square matrix is Hermitian or selfadj oint as de ned by Tl E if it is equal to its Hermitian conjugate if Hermitian conjugation introduces a minus sign the ma trix is skew Hermitian or antiHermitian HERMITIAN Tl T SKEW HERMITIAN Tl T VIII49 Donald G Luttermoser ETSU k III II 0 With this notation the inner product of 2 vectors with respect to an orthonormal basis can be written in matrix form aw al b VIII50 Matrix multiplication is not in general commutative ST y Ts commutator the difference between 2 orderings is called the s T E ST TS VIII51 It can also be shown that one can write the following commutator relation A3 0 A 1 O 21013 VIII52 The transpose of a product is the product of the transpose m reverse order SNT TS VIII53 and the same goes for Hermitian conjugates STl TlSl VIII54 The unit matrix is de ned by 1 0 0 1 quot39 VIII55 j j i In other words 127 52739 VIII 56 The inverse of a matrix written T4 is de ned by T lT TT l 1 VIII57 VIII 9 Vlll 10 PH Y5 4007 5007 Computational Physics i A matrix has an inverse if and only if its deter minant is nonzero in fact 1 1 det T 7 where C is the matrix of cofactors VIII58 ii The oofaotor of element Ti is lij times the determinant of the submatrix obtained from T by erasing the ith row by the jth column iii As an example for taking the inverse of a matrix let s assume that T is a 3x3 matrix of form T11 T12 T13 T T21 T22 T23 VIII 59 T31 T32 T33 lts determinant is then T11 T12 T13 detT ITI T21 T22 T23 T31 T32 T33 11 T22 T23 12 T21 T23 T32 T33 T31 T33 T21 T22 T 13 T31 T32 T11 T22T33 T23T32 T12 T21ng T23T31 T13 T21T32 T22T31 VIII 60 Donald G Luttermoser ETSU iv For this 3x3 matrix the matrix of cofactors is given by T22 T23 T21 T23 T21 T22 T32 T33 T31 T33 T31 T32 C T12 T13 T11 T13 T11 T12 T32 T33 T31 T33 T31 T32 T12 T13 T11 T13 T11 T12 T22 T23 T21 T23 T21 T22 Vlll Gl v The transpose of this cofactor matrix is then see Eq VIII45 T22 T32 T12 T32 T12 T22 T23 T33 T13 T33 T13 T23 T21 T31 T11 T31 T11 T21 T23 T33 T13 T33 T13 T23 T21 T31 T11 T31 T11 T21 T22 T32 T12 T32 T12 T22 Vlll GZ vi A matrix without an inverse is said to be singu lar vii The inverse of a product assuming it exists is the product of the inverses m reverse order ST 1 T IS l VIII63 p A matrix is unitary if its inverse is equal to its Hermitian conjugate UNITARY Ul U l VIII64 q The trace of a matrix is the sum of the diagonal elements TrT E g VIII65 il Vlll 11 Vlll 12 5 A vector under a linear transformation that obeys the following PH Y5 4007 5007 Computational Physics and has the property T1quotT1T2 T1quotT2T1 equation TIM Ala plex number A is called the eigenvalue a b c 01 Notice that any nonzero multiple of an eigenvector is still an eigenvector with the same eigenvalue ln matrix form the eigenvector equation takes the form Ta Aa VIII 68 for nonzero a or T A1a 0 VIII 69 here 0 is the zero matrix whose elements are all zero If the matrix T A1 had an inverse we could multiply both sides of Eq VIII 69 by T A1 1 and conclude that a 0 But by assumption a is not zero so the matrix T A1 must in fact be singular which means that its determinant vanishes T11 X T12 quot39 Tln T T A T n detT 1 21 l 22 l 0 Tnl Tn 39 VIII70 Expansion of the determinant yields an algebraic equation for A om Cn1 1 0M co 0 VIII71 VIII67 is called an eigenvector of the transformation and the com Donald G Luttermoser ETSU Vlll 13 where the coef cients C depend on the elements of T This is called the characteristic equation for the matrix its solutions determine the eigenvalues Note that it is an nth order equation so it has n complex roots i Some of these root may be duplicates so all we can say for certain is that an n x n matrix has at least one and at most 71 distinct eigenvalues ii In the cases where duplicates exist such states are said to be degenerate iii To construct the corresponding eigenvectors it is generally easiest simply to plug each A back into Eq VIII 68 and solve by hand for the components of a see Examples VIII 2 and VIII 3 6 In many physical problems involving matrices in both classical mechanics and quantum mechanics it is desirable to carry out a real orthogonal similarity transformation or a unitary transfor mation to reduce the matrix to its diagonal form 239e all non diagonal elements equal to zero a If eigenvectors span the space we are free to use them as a basis Tlf1gt A1f1gt Tlf2gt A2f2gt b The matrix representing T takes on a very simple form in this basis with the eigenvalues strung out along the main Vlll 14 PH Y5 4007 5007 Computational Physics diagonal and all other elements zero A1 0 0 T 2 VIII72 0 0 An c The normalized eigenvectors are equally simple 1 0 0 l 0 am 0 a2 0 a 0 0 0 l VIII73 d A matrix that can be brought to diagonal form Eq VIII 72 by change of basis is said to be diagonalizable e In a geometrical sense diagonaliZing a matrix is equiva lent to rotating the bases of a matrix about some point in the space until all of the off diagonal elements go to zero If D is the diagonalized matrix of matrix M the operation that diagonalizes M is D SMS l VIII74 where matrix S is called a similarity transformation Note that the inverse of the similarity matrix can be con structed by using the eigenvectors in the old basis as the columns of 84 84 am VIII75 f There is great advantage in bringing a matrix to diagonal form it is much easier to work with Unfortunately Donald G Luttermoser ETSU Vlll 15 not every matrix can be diagonalized the eigenvec tors have to span the space for a matrix to be diagonalizable 7 The Hermitian conjugate of a linear transformation called a Hermitian transformation is that transformation Tl which when applied to the rst member of an inner product gives the same result as if T itself had been applied to the second vector Tlozl mfg VIII76 for all vectors a and a Note that the notation used in Eq VIII 76 is commonly used but incorrect actually means and Tlod means the inner product of the vector Tla b Note that we can also write mm alTb Tlalb aw VIII77 c In quantum mechanics a fundamental role is played by Hermitian transformations Tl The eigenvectors and eigenvalues of a Hermitian transformation have 3 cru cial properties i The eigenvalues of a Hermitian transforma tion are real ii The eigenvectors of a Hermitian transfor mation belonging to distinct eigenvalues are orthogonal iii The eigenvectors of a Hermitian transfor mation span the space Vlll 16 PH Y5 4007 5007 Computational Physics Example VIIIil Given the following two matrices 112 20 2 A203B010 22 222 232 compute a A B 13 A132 0 A2 BL C1 A2 9 N2 f N2 3 TTB2 h detB and B l Check that BB 1 1 Does A have an inverse Solution a Sum the respective elements of the matrix Solution b Multiply rows of A by columns of B 20 10132 2022 AB 4032 009 2206 42022 0 226 204 3 132 32 432 9 6 22 62 6 22 6 Solution c A B AB BA we already have AB 202 20 2 210 2i BA 020 000 030 i64i i0 4i l94 0 0 0 2 0 3 632 32 12 7 3 1 32 32 0 0 0 A B 4 32 9 6 22 2 0 3 62 6 22 6 6 32 32 12 Donald G Luttermoser ETSU Vlll 17 3 1 32 32 2 32 9 3 22 632 62 6 Solution d Transpose of A ip A about the diagonal A l 2 22 l 0 2239 2 3 2 Solution e Complex conjugate of A A multiply each 2 term by l in Ai l 2 2239 l 239 0 3 22 2 Solution f Hermitian of A AlEAquot Solution g Trace of B TrB Z l 2 2239 l 0 2239 239 3 2 3 ZBZi212 391 Solution h Determinant of B detB 22 0 00 Solution i lnverse of B 0 20 24 0 1 V111 18 PH Y5 4007 5007 Computational Physics where 10 00 01 32 i2 i3 2 0 z O i 2 i 20 C 32 2392 23 32336 0 1 24 20 Z 1 0 0 0 01 then 2 ampi 34 a 0 3 0 i 62 40 1 6i062 2i0 2i BB 1 a 000 030 000 2i0 2i 39 12 104 1 3 0 0 1 0 0 g 0 3 0 0 1 0 0 0 3 0 0 1 1fdetA y 0 then A has an inverse detA 106i 14 6ii 4i 0 6i 46i40 Assuch Example VIII72 Find the eigenvalues and normalized eigenvectors of the following matrix 1 1 M e lt 0 1 gt Can this matrix be diagonalized Solution 0mM AD Donald G Luttermoser ETSU Vlll 19 only one eigenvalue From Eq VIII 68 we get 1 1 11 i 1 l1 0 l 12 i 12 i We get two equations from this eigenvector equation a1a2 a1 12 12 The second equation tells us nothing but the rst equation shows us that 12 0 We still need to gure out the value for 11 We can do this by normalizing our eigenvector a agtz 2 2 1 alal lail Z lall2 la2l2 lall2 or 11 1 Hence our normalized eigenvector Since these eigenvectors do not span the space as described on page VIII 4 A2cii this matrix cannot be diagonalized Example VIII73 Find the eigenvalues and eigenvectors of the following matrix 2 0 2 M 2239 239 2239 1 0 1 Vlll 20 PH Y5 4007 5007 Computational Physics Solution The characteristic equation is Z A 0 2 M 1 2239 i A 2239 1 0 l A i A 2239 2239 i A 2A 0 l A02 1 0 2 AW W l A 0 2l0 i Ml 2 A AAA22 2A 22 2M2A2A2iAiA2 A2 A32i 2A 31i2 z 0 To nd the roots to this characteristic equation factor out a A and use the quadratic formula solution equation 0 31i z39 A21z z A1 0 1z l li2 4z39 A2 2 i 1z l l2z 1 4z 2 i 1z l 2239 2 However note that l 239 2239 As such the equation above becomes 1z l 1 12 273 i 1z ll z39 2 1z l z 2 71 A 2 2 1 39 1 39 239 A3 z z ziz7 2 72 Donald G Luttermoser ETSU Vlll 21 so the roots of A e the eigenvalues are 0 l and 239 Now let s call the components of the rst eigenvector agt a1 a2 13 which corresponds to eigenvalue A1 0 The eigenvector equation becomes 2 0 2 a1 a1 0 2i i 2i a2 0 a2 0 1 0 1 a3 a3 0 which yield 3 equations 211 213 0 2z39a1 2112 22113 0 l1 13 0 The rst equation gives 13 11 the second gives 12 0 and the third is redundant with the rst equation We can nd the values for 11 and 13 by norrnaliZing 3 2 1 Ola lail lall2 lt122 la32 lall2 lall2 2 a12 or 11 a3 Hence our eigenvector for A1 is l agtai 0 forA10 For the second eigenvector let s call it b we have 2 0 2 b1 b1 b1 2239 239 2239 b2 1 b2 b2 l 0 1 b3 b3 b3 which yield the equations 2b1 2b3 b1 2ib1ib2 Zibg b2 b1 b3 b 7 Vlll 22 PHYS 40075007 Computational Physics with the solutions b3 12b1 and b2 1 i2b1 Normalizing gives 3 1 lt l gti21bi2 bl2b22b32 Ib1vlt1ygtlt1gigt W W i W 4 2 2 2 1 2 7 7 7b 4lbil 4b1 41 7 2 77 b jllla or b1 So b2 1 and b3 giving our nal eigenvector for A2 as 2 l gtbg 1 2 forA21 l Finally the third eigenvector call it W c is 2 0 2 C1 C1 ZC1 2i i 2i C2 i C2 ZCQ l 0 1 C3 C3 ZC3 which gives the equations 2C1 2C3 ZC1 223901 202 2203 202 C1 C3 ZC3 with the solutions 03 01 0 with 02 undetermined Once again we can normalize our eigenvector to determine this undetermined 02 coef cient 1 3 2 WM lCil Z W2 C22 lcsl2 C22 a Donald G Luttermoser ETSU or Cg l which gives our third eigenvector 7gtc l forA3z39 0 C Computational Linear Algebra 1 Here we will be solving the set of equations as described in Eqs VIII123 a b However for the present let N number of unknowns and M number of equations If N M there s a good chance we can obtain a unique solution i However if one or more of the M equations is a linear combination of the others gt row degen eracy occurs gt we will not be able to obtain a unique solution ii Or if all equations contain variables in exactly the same linear combination gt column degeneracy occurs gt no unique solution can be found iii For square N M matrices row degeneracy implies column degeneracy and vise versa iv If a set of equations are degenerate the matrix is said to be singular Vlll 24 c d 8 PH Y5 4007 5007 Computational Physics Numerically other things can go wrong i If some equations are close to being a linear com bination of the other equations in the set roundoff errors may render them linear dependent at some stage of the calculations ii Accumulated roundoff errors in the solution can swamp the true solution gt this can occur if N is too large gt direct substitution of the so lution back into the original equations can verify this Linear sets with 2 lt N lt 50 can be routinely solved in single precision without resorting to sophisticated meth ods gt in double precision N gt 200 without worrying about roundoff error As we have discussed solution to a set of linear equations involve inverting matrices To write the most e icient ma trix inverter one needs to know how numbers are stored Assuming we have matrix 111 112 am 121 122 a2N A 0M1 0M2 aMN i Column storage which IDL calls row major stor age a11121 aM1a12a22 aM2 a1Na2N aMN gt Fortran and IDL use this method Donald G Luttermoser ETSU 2 The basic process of solving linear systems of equations is to eliminate variables until you have a single equation with a single ii Row storage which IDL calls column major storage a11a12a1Na21a22agNaM1aM2aMN gt C and C use this method iii The techniques we will be discussing here are de signed with column storage in mind unknown 3 For nonlinear problems an iterative scheme is developed that solves a linearized version of the equations 4 Equations that do Lot depend upon time are called autonomous systems that is w M Q d f lf initial conditions are of the form for all 239 l g 239 g N and t the solution points in the N dimensional space of the variables are called steady state If we start at steady state we stay there forever Locating steady states for linear equations or ODE s is important since they are used in stability analysis prob lems 9 It is easy to see that Xquot 061 le is a steady state if and only if VIII 79 vmwa Vlll 25 Vlll 26 e 0 PH Y5 4007 5007 Computational Physics 01quot for all 239 VIII 80 since this implies that dXdt 0 Hence locating steady states reduces to the problem of solving N equations in the N unknowns This problem is also called nding roots of f 96 5 We shall now discuss the various numerical techniques used in solVing sets of linear equations D Gaussian Elimination 1 The problem of solVing 0 is diVided into two important classes a b C Techniques used for linear equations such as Gaussian elimination and matrix inversion Techniques used for nonlinear equations such as New ton s Method With Gaussian elimination we set up the N linear equa tions with N unknowns in the form of Eqs VIII 123 111061 112062 aiNxN b1 0 121x1 12ng angN b2 0 Imam aNgxg aNNxN bN 0 VIII 81 or it matrix form A a b 0 VIII 82 Donald G Luttermoser ETSU where 111 112 061 b1 Aa2ia22quotxxgbb2 VIII 83 1 One can then go through a process of eliminating variables 8 by adding or subtracting one equation to or from the other equations Example VIIIi4 Take the equations 2x1x2 4 4x1 x2 2 We want to eliminate 961 from the second equation we multiply the rst equation by 2 and subtract the rst equation from the second equation giving 3x2 6 or 062 2 This step is known as forward elimination For larger sets of equations the forward elimination procedure eliminates 961 from the second equation then eliminates x1 and 062 from the third equation and so on The last equation will only contain the variable 061V which can then be solved We then carry out a backsubstitution am is then plugged back into the N 1 equation to solve for xN1 In the example above substitute x2 2 into the rst equation giving 2x124 or 0611 This method of solving systems of linear equations is called Gaussian elimination A portion of a Fortran 77 code Vlll 27 Vlll 28 PHYS 40075007 Computational Physics that might perform such a Gaussian elimination would be written as gtllt Forward elimination DO K 1 N 1 Go to column k operate DO K1 N on the rows below column k COEFF AK AKK DOJKLN AUJAUJCOEFFAU D ENDDO ALKCOEFF BUBUCOEFFEXK ENDDO ENDDO Then the backsubstitution is performed Via gtllt Backsubstitution XN BN ANN Start from bottom and work DO N 1 1 1 work upward Note This loop SUM B goes from n 1 to 1 in steps of 1 DO J 1 N Skip lower triangular part SUM SUM AJXJ ENDDO X SUMA ENDDO f Gaussian elimination is a simple procedure yet it has its pitfalls Consider the set of equations 8x1 2 x3 5 x1 x2 3 x1 x3 4 In the limit 5 gt 0 the solution is 961 1062 2063 3 For these equations the forward elimination step would start by multiplying the rst equation by 15 and sub Donald G Luttermoser ETSU tracting it from the second and third equations giving 8581 x2 x3 5 l lax2 15x3 3 55 l x2 l 15x3 4 55 i Of course if 5 0 we have big problems since the 15 factors blow up ii Even if 5 y 0 but is small we are going to have serious roundoff problems In this case 15 gtgt 1 so the equations above become 8581 x2 x3 5 l x2 15x3 55 l x2 15x3 55 At this point it is clear that we may not proceed since the second and third equations are now iden tical gt 3 unknowns with only 2 equations g Fortunately there is a simple x we can just interchange the order of the equations before doing the forward elim ination x1 x2 3 8x1 2 x3 5 x1 x3 4 i The next step of forward elimination gives x1 x2 3 l 5x2 x3 5 35 x2 x3 4 3 ii Roundoff eliminates the 5 terms giVing x1 x2 3 x2 x3 5 x2 x3 1 VIII 29 Vlll 30 h PH Y5 4007 5007 Computational Physics iii The second step of forward elimination removes 062 from the third equation using the second equa tion x1 x2 3 x2 x3 5 2x3 6 iv You can easily substitute back giving 961 1062 2x3 13 Algorithms that rearrange the equations when they spot small diagonal elements are said to pivot The price of pivoting is just a little extra bookkeeping in the program but it is essential to use pivoting in all but the smallest matrices Even with pivoting one cannot guarantee being safe from roundoff problems when dealing with very large matrices The program below performs Gaussian elimination with pivoting subroutine geaabbnnpx Perform Gaussian elimination to solve aax bb Matrix aa is physically np by np but only n by n is used n lt np parameter nmax 100 real aanpnpbbnpxnp real anmaxnmax bnmax integer indexnmax real scalenmax if np gt nmax then print ERRUR Matrix is too large for ge routine stop end if do i1n bi bbi Copy vector do j1n aij aaij Copy matrix end do end do Vlll 31 Donald G Luttermoser ETSU Forward elimination do i1n indexi i scalemax 0 do j1N scalemax amaxlscalemaxabsaij end do sca1ei end do scalemax do k1N 1 ratiomax 0 do ikn ratio absaindexiksca1eindexi if ratio gt ratiomax then ji ratiomax ratio end if end do indexk indexj indexj indexk indexk indexk do ik1n coeff aindexikaindexkk do jk1n aindexij aindexij coeffaindexkj end do aindexik coeff bindexi bindexi aindexikbindexk end do end do Back substitution xn bindexnaindexnn do in11 sum bindexi do ji1n sum sum aindexijxj end do xi end do sumaindexii return end 2 Working with Matrices 3 It is easy to obtain determinants of a matrix using Gaus sian elimination After completing forward elimination Vlll 32 0 PH Y5 4007 5007 Computational Physics one simply computes the product of the coef cients of the diagonal elements Take the equations in Example VIII 4 Aim With forward elimination these equations become the matrix is 4 6 2061 062 3x2 The products of the coef cients of the diagonal elements of this matrix is 2 3 6 which is the determinant of A above However it should be noted that this method is slightly more complicated when pivoting is used If the number of points is odd the determinant is the negative of the product of the coef cients of the diagonal elements Matrix Inverse and Gaussian Elimination i Recall the linear equation in Eq VIII 82 A x b 0 VIII84 where we solved for the vector X by Gaussian elim ination Note however that we could also have solved it with a little matrix algebra x A lb VIII85 where A 1 is the matrix inverse of A ii It shouldn t surprise you that the inverse of a ma trix is computed by repeated applications of Gaus sian elimination or a variant called LU decompo sition Donald G Luttermoser ETSU iii As we have already discussed the inverse of a matrix is de ned by A A4 I VIII86 where I is the identity matrix 1 0 0 0 1 0 I 0 0 1 VIII 87 iv De ning the column vectors 1 0 0 1 0 e1 0 e2 0 i a eN 0 a 39 39 1 VIII 88 we may write the identity matrix as a row vector of column vectors I e1 e2 eN VIII 89 v If we solve the linear set of equations AX1 81 the solution vector x1 is the rst column of the inverse matrix A l vi If we proceed this way with the other e s we will compute all of the columns of A4 In other words our matrix inverse equation Eq VIII 86 is solved by writing it as A x1 x2 xN e1 e2 eN VIII 91 Vlll 33 Vlll 34 PH Y5 4007 5007 Computational Physics vii After computing the X s we build A 1 as A 1X1X2 xN VIII92 viii It is usually not necessary to write your own routines to do matrix inverse since virtually all pro gramming languages has routines that will do this for you For instance Matlab has the built in func tion invA IDL has NVERTA Fortran uses the Numerical Recipes LU DCMP subroutine which can be freely downloaded LINPACK also has inverting matrices routines and C has similar routines to Fortran ix A handy formula to remember involves the in verse of a 2 x 2 matrix A4 1 111022 112021 For larger matrices the formulas quickly become 122 112 VIII93 a21 111 very messy c Singular and IllConditioned Matrices i A matrix that has no inverse is said to be singu lar e 9 1 1 2 2 And remember a singular matrix has a determi A nant of zero ii Sometime a matrix is not singular but is so close to being singular that roundoff errors may push it over the edge A trivial example would be 1 5 1 l 2 2 l a Donald G Luttermoser ETSU where 5 ltlt 1 iii The condition of a matrix indicates how close it is from being singular a matrix is said to be ill conditioned if it is almost singular iv Formally the condition criterion is de ned as the normalized distance between a matrix and the nearest singular matrix All of the programming languages mentioned above also have the ability of returning this normalized distance with either the inverse function or a separate function E LU Decomposition 1 Suppose we are able to write a matrix as the product of 2 matri ces L U A VIII 94 where L is lower triangular has elements only on the diagonal and below and U is upper triangular has elements only on the diagonal and above 2 In the case of a 3 x 3 matrix A we would have 0111 0 0 511 512 513 an 112 113 021 022 0 0 522 523 121 122 123 031 O32 O33 0 0 533 131 132 133 VIII95 a We can use a decomposition such as Eq VIII 94 to solve the linear set Ax LU x L Ux b VIII96 by rst solving for the vector y such that L y b VIII97 Vlll 35 Vlll 36 PH Y5 4007 5007 Computational Physics and then solving U X y VIII 98 b The advantage to this method is that the solution of a triangular set of equations is quite trivial Thus Eq Vlll 97 can be solved by forward substitution as follows i yl VIII 99 0111 1 221 yi bi Zaijyj i2737quot aN aii j1 c Then Eq VIII 98 can then be solved by backsubstttutmg exactly as in M LN VIII100 NN 1 N x 7 yi ijxj iN 1N 21 it ji1 1 Equations VIII 99 and VIII 100 total for each right hand side b N 2 executions of an inner loop containing one multiply and one add If we have N right hand sides which are the unit column vectors which is the case when we are inverting a matrix then taking into account the leading zeros reduces the total execution count of Eq VIII 99 from N3 to N3 while Eq VIII 100 is un changed e Notice that once we have the LU decomposition of A we can solve with as many right hand sides as we then care to one at a time This is a distinct advantage over the Gaussian elimination scheme described earlier Donald G Luttermoser ETSU 3 Performing the LU Decomposition a b C d 8 How do we solve for L and U given A First we write out the ijth component of Eqs VIII 94 and VIII 95 That component is always the sum beginning with ai1 1j 12739 The number of terms in the sum depends however on whether 239 or j is the smaller number We have in fact the 3 cases ilt j iai151j Oli2 2j 39 39 39 aii ij aij VIII101 ij iai151j 01239252139 quot39 Olii jj aij VIII102 1 gt1 3 Otii ij 01239232139 39 39 39 Otij jj aij VIII103 Eqs VIII 101 VIII 103 total N2 equations for the N2 N unknown a s and s the diagonal being represented twice Since the number of unknowns is greater than the number of equations we are invited to specify N of the unknowns arbitrarily and then try to solve for the others In fact it is always possible to take 01 E 1 239 1N VIII 104 Often the Crout algorithm is used to solve the set of N2 N equations 69 Eqs VIII 101103 for all the a s and s This is done by just arranging the equations in a certain order i Set 04 12 1 N Eq VIII104 Vlll 37 Vlll 38 PH Y5 4007 5007 Computational Physics ii For each j 12 3 N do these 2 procedures First for 2 l 2 j use Eqs VIII 101 Vlll 102 and VIII 103 to solve for ij namely 23971 52739 au Z aik kj VIII 105 121 When 2 1 in Eq VIII 28 the summation term is taken to mean zero iii Second for 2 j lj 2 N use Eq Vlll 103 to solve for 02 namely 1 7394 Olij 7 a27 Z Oak8M 52239 21 Be sure to do both procedures before going on to the next 339 iv In brief Crout s method lls in the combined ma trix of 02 s and s 511 512 513 021 522 523 0131 0132 533 by columns from left to right and Within each col umn from top to bottom f Pivoting is absolutely essential for the stability of Crout s method P02222021 p2210t239ng interchange of rows can be implemented ef ciently and this is enough to make the method stable The Numem39cal Rec239pe s subroutine LUD CM P is an LU decomposition routine using Crout s method with partial pivoting I recommend its use whenever you need to solve a linear set of equations Donald G Luttermoser ETSU Vlll 39 9quot H N H I39m39u39 39 Hulu39l39b Figure VIII 1 A two mass coupled harmonic oscillator with the origin set at the position of the left wall ie Case A F Coupled Harmonic Oscillators 1 A canonical example of a system of linear equations is the case of a coupled harmonic oscillator as shown in Figure VIII 1 Each spring has an unstretched length of L1 L2 and L3 in this example and a spring constant of 61162 and k3 In between each spring is an object of mass m1 and m2 Finally the distance between the non moving wall is Lw 2 The equation of motion for block i is 106 i U dUi i dt i 27 7 7 VIII1 7 d mi 0 where E is the net force on block i 3 At the steady state the velocities 1 are zero and the net forces E are zero gt static equilibrium 4 When working with coupled oscillators one must de ne a frame of reference from which the measurements are made For in stance one could de ne the reference point to be the left wall of the system Case A as shown in Figure VIII 1 then the net Vlll 40 PH Y5 4007 5007 Computational Physics Figure VIII 2 A two mass coupled harmonic oscillator with the origins for each coordinate set at the massls equilibrium positions 16 Case B force equations become F1 m1ab391 11 L1 k2gg2 x1 L2 VIII 108 F2 mgab39g 22 x1 L2 x2 L3 5 To solve these equations we can use the matrix techniques that have been described earlier in this section 6 9 Gaussian elimina tion or LU decomposition by writing these equations in matrix form F k k k x l 1k 2 4623163 1 VIII110 k1L1k2L2 k2L2 k3L3 Lw or F KAx b VIII 111 6 One could also choose the equilibrium positions of each block as the reference Case B see Figure VIII 2 and write the net force equations as Donald G Luttermoser ETSU VIII 41 F1 m1ab391 k1x1 k2061 x2 VIII112 F2 ngug k3x2 k2ltx2 x1 VIII113 7 In matrix form Case B takes the form of F1 i k1 62 kg 531 F2 7 k2 k2 kg 962 VIII 114 or in shorthand notation F KB X VIII 115 As can be seen the unstretohed lengths of the springs do not enter into the second case since measurements are being made from equilibrium positions PHYS40075007 Computational Physics Course Lecture Notes Section VII Dr Donald G Luttermoser East Tennessee State University Version 41 Abstract These Class notes are designed for use of the instructor and students of the course PHYS 40075007 Computational Physics taught by Dr Donald Lutternioser at East Tennessee State University VII Methods of Differentiation and Integration A Derivatives 1 The derivative of a function is de ned by a b C 01 df06 E hm f06 0 x dx hgt0 h In computational work functions and their independent variables are given by tabulated data and or derived data Since there is a subtraction in Eq VII l subtraction cancellation can lead to rather large errors in the deter mination of a derivative via numerical techniques The computer s nite word length can cause the numera tor to uctuate between 0 and the machine precision em as the denominator approaches zero In this section we will discuss the various techniques used to calculate derivatives numerically and estimating the error in the derivative 2 Method 1 Forward Difference 3 0 Write out the function as a Taylor series at a position of one step forward from the current position h2 h M h M was 51 g gwc 39 39 VII2 where h is the step size as shown in Figure VII l We obtain the forward dz erence derivative algorithm by solving Eq VII l for fx Joex 2 W7 VH3 VII1 VII 2 PH Y5 4007 5007 Computational Physics h 2 HM E f 06 7 VII4 where the subscript c denotes the computed expression i The approximation of Eq VII 3 has an error pro portional to h as shown in Eq VII 4 ii We can make the approximation error smaller and smaller by making h smaller and smaller iii We can t make it too small however since all precision will be lost if the subtraction cancellation error becomes larger than the step size c Consider for example my a bx VII5 the exact derivative is f 2bx VII 6 and the computed derivative is W m m h me This would only be a good approximation if h ltlt 206 2m bh VII7 3 Method 2 Central Difference a An improved approximation to the derivative starts with the basic de nition Eq Vll l For this technique in stead of making a step of h forward we form a central di erence by stepping forward by h 2 and stepping back ward by h 2 06 h2 x h2 h EDCfxh Vii 8 Donad G Luttermoser ETSU VII 3 y AA xhl2 x xhl2 xh K Figure VII 1 Forward difference solid line and central difference dashed line methods for numerical integration i The symbol DC represents center difference ii Carrying out the Taylor series for x I h gives 1 f 2 fx i f m VII9 iii The important difference from Eq VII 3 is that when x h2 is subtracted from x h2 all terms containing an odd power of h in the Taylor series cancel iv Therefore the central difference algorithm be comes accurate to one order higher than h gt hg v If f3h224 ltlt f2h2 then the error with the central difference method should be smaller than the forward difference b For the polynomial given in Eq VII 5 the central dif VII 4 PH Y5 4007 5007 Computational Physics ference gives the exact answer regardless of the size of h m m h2 m h2 h 2m VII10 x 4 Errors in Taking Derivatives a One should always try to keep calculation errors em at a minimum This typically occurs when 60 m eaer in Eq V 33 b Because differentiation subtracts 2 numbers that are close in value the limit of roundoff error is essentially machine precision fx m fx h mg m in VII11 h h gt em m 6 VII12 c The approximation error with the forward difference al gorithm Eq VII 4 is an 9h term while that with the central difference algorithm Eq VII 9 is an 9h2 term M N 102 a1er N 2 VII13 C f3 h eaim m 24 VII14 d The value of h for which roundoff and approximation er rors are equal is therefore in N fd ifem h N apply 7 2 7 em C 1432 I g gaipm 24 VII 16 2 m 24 m gt hi 6 3 6 VII17 W7 Cd f3 Donald G Luttermoser ETSU e As an example for the ex and cosx functions f m fa m f 3 in single precision with em 10 7 one would get hfd m 00005 and had 001 This makes the central difference algorithm better for calculating derivatives since a larger h would mean a smaller error gt here the error in the central difference method is 20 times smaller than the er ror in the forward difference method 5 The Method in Calculating Second Derivatives a Taking second derivatives will involve an additional sub traction step as compared to the rst derivative calcula tion leading to a larger subtraction error b We can remedy this with a little algebra in the central difference method Taking the second derivative of Eq VII 8 and then rewriting the rst derivatives with a for ward and backward difference equation we get ism W mm W m h2gtl VII19gt Note however that one must keep the second I Z derivative formula in this form to minimize can cellation error E Integration 1 The method of numerical integration is sometimes referred to as numerical quadrature gt summing boxes a The de nition of an integral given by Riemann is himh from m VII20gt VII 5 VII 6 PH Y5 4007 5007 Computational Physics b If we ignore the limit the integral just becomes a sum mation of boxes or quadrilaterals lying below the curve b N a my day m mi w VII21 where N is the number of points in the interval 1 b and f is evaluated at those interval points 2quot f E The wi s are integration weights which are proportional to h the distance between points 239 and 239 l c The different integration schemes presented here will all make use of Eq VII 21 d In the simplest integration schemes the integrand is ap proximated by a few terms in the Taylor series expansion of f and the terms are integrated gt typically adding more terms in the series yield higher precision e This is the basis of the NewtonCotes methods gt the total interval is divided into equal subintervals with the integrand evaluated at equally spaced points 06 Two such Newton Cotes methods include the i Trapezoid rule a rst order method ii Simpson doh l rule a second order method f More accurate integration methods involve the use of non equally spaced intervals gt these are the methods of Gaus sian quadrature i Gaussian quadrature methods are typically more precise than Newton Cotes methods as long as there are no singularities e denominators going to in nity non continuous functions in the integrand or its derivative Donald G Luttermoser ETSU ii You are better off to remove such singularities algebraically before attempting Gaussian quadra ture For example 1 x fa die L01f xdx 01fxdx VII22 iii Regions where a function is slowly varying re quire fewer integration points and regions with rapid variations in the function will have many in tegration points in order not to miss any oscillation as can be seen for such functions evenly spaced integration points will not represent the true inte gral accurately 2 Trapezoid Rule a b C As an introduction to the numerically integration of a function consider the generic integral I b mm VII23 The most straight forward way to solving such an integral function is to evaluate f at a few points and t a simple curve eg piecewise linear through these points One way to do this is to t trapezoids of equal width h under the curve represented by f and add up the total areas of these trapezoids Let h b a VII24 N l7 x ai 1h i12N VII 25 where a and b correspond to the initial and nal end points and N the total number of points in the interval VII 7 VII 8 PH Y5 4007 5007 Computational Physics 1 b square bracket means we include the endpoints 239 6 there are N 2 points in between the endpoints Note that the trapezoid rule requires an odd number of points N i Straight lines connect the points and this piecewise linear function serves as our tting curve ii The integral of this tting function is easy to com pute since it is the sum of the areas of trapezoids The area of a single trapezoid is 1 Ti 596241 06ifi1 fi VII26 iii The true integral is estimated as the sum of the areas of the trapezoids so Nil IRITT1T239HTN1 Z VII27 il Notice that the last term in the sum is N 1 since there is one fewer panel than grid points iv The general formula simpli es if we take equally spaced grid points as given by Eq VII 25 v Then the area for an individual trapezoid e the area of the one trapezoid within that interval becomes 1 Ti hfil fi 7 VII28 or in our original integral notation 41 fltxgtdx 2 W Ehfi hfi VII 29 As such our weights in Eq VII 21 for the indi vidual points are U h Donald G Luttermoser ETSU d e f g We now need to add up all of the trapezoids in the subin tervals across the entire interval 1 b gt the trapezoid rule 1 1 TWO Ehfl hf2 hf th71 Eth 1 N71 hf1 fN h fi VII 30 OIquot 1695degfihf2hfgth71ng VII31 Note that since each internal point gets counted twice it has a weight of h whereas the endpoints get counted just once and thus have weights of h 2 LUZ39 hh VII 32 2 2 lt gt Our approximation error also called the truncation error or the algorithmic error here for the trapezoid rule can be written either as can I ITxh T12b luffc VII33 for some value x C that lies in 1 b or as Gaper 112 hglfb fa 9h4 VII 34 This error is proportional to h2 and Eq VII 34 warns you that the trapezoidal rule will have dif culties if the derivative diverges at the end points VII 9 VII 10 PH Y5 4007 5007 Computational Physics 3 Bomb erg Integration 3 m But how many panels ta trapezoids should one use ta how small should h be One way to decide is to re peat the calculation with a smaller interval If the answer doesn t vary signi cantly then we accept it as correct i This doesn t work with pathological functions or in unusual scenarios but this won t happen very often ii In the trapezoidal rule if the number of panels is a power of two we can halve the interval size without having to recompute all the points iii De ne the sequence of interval sizes as l hn VII35 hl b a hg b a iv For n l we have only one panel so hvw1 aMKwfwt 2 VII36 v There is a simple recursive formula for calculating Thn1 from Thn 1 2n1 Thn1 hn1 fla 2i lhn1l 7 VII37 By using this recursive method we can keep adding pan els until the answer appears to converge can greatly improve this process by using a method called Romb erg integration However we Donald G Luttermoser ETSU VII 11 i The method builds a triangular table of the form R11 R21 M2 R31 R 2 ii The formula for the rst column is just the recur sive trapezoidal rule Rn1 1mm VII38 iii The successive columns to the right are com puted using the Richardson extrapolation for mula Rnlml 7 nlm VII 39 1 4m llRm 1m Rnm iv The most accurate estimate for the integral is Rnm the value at the bottom right corner of the table v To understand why the Romberg scheme works consider the approximation 16 truncation error capprxhn I Thn then using Eq VII 34 swam gulfb fal 0 VII40 Since hnH hnZ 1 Ehilfb W W VII41 Eappnc hn1 VII 12 PH Y5 4007 5007 Computational Physics vi Consider the second column of the Romberg ta ble The approximation error for Rn 12 is I Rm 1 2 I IThn1 lIThn1 Thnl ewxmn ieapmlthn1gt swim g hi W fal 00 90 VII42 vii Notice that the hi term serendipitously cancels out leaving us with a truncation error that is of order hi The next third column of the Romberg table removes this term and so forth c The following is an example of an integration subroutine using the Romberg method If you are interested in using this routine let me know and send you a copy of the routine shown below written in Fortran 77 Note that func is an externally supplied subroutine that contains the function to be integrated subroutine rombfa b N func param R NR l Routine to compute integrals by Romberg algorithm real RNRNR param external func h b a This is the coarsest panel R11 h2 funcaparam funcbparam up 1 up is the number of panels do i 2 N h h 2 Donald G Luttermoser ETSU do k 1 npl 2 This loop goes k135 sum sum funcakh param enddo npl Compute the first column entry Ri1 Ri1 O5Ri 11 hsum m 1 do j 2 i m 4 m Rij Rij 1 Rij 1 Ri 1j 1m 1gt enddo enddo return end 4 Simpson s Rule Simpson eh 3 b C 01 Instead of tting two adjacent points with trapezoids we will now t three adjacent points with parabolas x m ago x 7 VII 43 for each interval still keeping the intervals evenly spaced The area of each section is then the integral of this parabola Iih 0063 062 Th A ax2 x7dx 7x Ii VII44 This is equivalent to integrating the Taylor series up to the quadratic term Hence the Simpson rule is a second order polynomial method In order to relate the parameters oz and 7 to the func tion we consider an interval from 1 to 1 20 2 37 lmy x 7 dx VII 45 Compute the other columns Ri2Rii VII 13 VII 14 PH Y5 4007 5007 Computational Physics i Note however for the function itself f1 a a flt1gtlt 1gt f0 ND 7 W N a 71 0 VII46 ii Using the results of Eq VII 46 in Eq VII 45 yields 1ax2 x7dxm 3 3 VII47 iii Because 3 values of the function are needed we evaluate the integral over two adjacent intervals gt evaluate the functions at the two endpoints and the middle zih xih Ii AH mm fxdxIiihfxdx h 4h h gfiil Efi gfz H VII 48 2 e Simpson s rule requires the elementary integration to be over pairs of intervals gt this requires the number of intervals to be even and hence the number of points N to be odd f To integrate across the entire range 1 b we add up con tributions from each pair of subintervals counting all but the rst and last endpoints twice ffltxgtdx w gf1f2fgf4m fN71fN VII49 Donald G Luttermoser ETSU i From this integral we see that the total set of weights is h 4h 2h 4h 2h 4h h 10 i u ii The sum of your weights provides a useful check on your integration w N 1 h VII51 il Remember N must be odd 5 Errors in the NewtonCotes Methods a As we have said above and as was the case for differ entiation our Newton Cotes methods are essentially just Taylor series expansions of a function We will highlight here what is derived in more detail in the textbook b The approximation ie truncation or algorithmic error 6mm can be estimated by the next higher term in the Taylor series that is not used in the evaluation of the integral 3 EaPpDGt O 7 b a 5 EaPPDFS O a where the subscripts t and s in the subscripts refer to trapezoid and Simpson s rule respectively VII53 c The relative error e is just these approximation errors divided by the value of the function as VII54 VII 15 VII 16 PH Y5 4007 5007 Computational Physics 1 Assume that after N steps the relative roundoff error is random so em VN em VII 55 where em is the machine precision N 10 7 for single pre cision 10 15 for double precision e v Assume that the total error is given by Eq V 33 Etot Em apprx VII56 i We want to determine the value of N that min imizes the total error This will occur when the two relative errors are of equal magnitude which we approximate even further by assuming that the two errors are equal Eappmet 6T0 6appnc f ii Now let s make the following assumptions n f m 1 VII58 f l b 1 hi VII 59 a gt N f For the trapezoid rule Eq VII 57 becomes Mb agt3 1 VN em m fN2 N2 VII 60 l i Then the optimum numbers of N steps for the trapezoid rule are N i l i 110 725 631 for single precision h 110 1 5 5 106 for double precision VII62 Donald G Luttermoser ETSU VII 17 ii The corresponding errors are 60 Wem 3 X 10 6 for single precision 10 12 for double precision VII 63 g For Simpson s rule Eq VII 57 becomes N fl4b a5 1 N 6m N W m l gtN m VII 65 i Then the optimum numbers of N steps for Simp son s rule are N i l i 110 729 36 for single precision T h T 110 15 9 2154 for double precision VII 66 ii The corresponding errors are i 6 X 10 7 for single precision 6T0 N N 6m 7 l 5 X 10 14 for double precision VII 67 h These results are illuminating because they show that i Simpson s why you little rule is an improve ment over the trapezoid rule ii It is possible to obtain an error close to the ma chine precision with Simpson s rule and with other higher order integration algorithms iii Obtaining the best numerical approximation to an integral is not obtained by letting N gt 00 but with a relatively small N g 1000 VII 18 PH Y5 4007 5007 Computational Physics 6 There are higher order Newton Cotes methods the 3rd degree 38th method and the 4th degree Milne method See Table 41 in your textbook on page 48 for the elementary weights that one would use for these two methods 7 Gaussian Quadrature a Often one cannot choose an evenly spaced interval of grid points to do a numerical integration From Eq VII 21 we have my day m w1f061 waxN VII68 One then can ask the question is there an optimal choice for the grid points or nodes x and the weights 10 to solve this integral b This question leads us to formulate a new class of integra tion formulas known collectively as Gaussian quadra ture i In this class we will use only the most common formula namely GaussLegendre quadrature c There are many other kinds of Gaussian quadra ture that treat speci c types of integrands such as the Gauss Laguerre formularization which is optimal for integrals of the form 5 e xfm dx 0 The derivation of the other Gaussian formulas see Table 42 in your textbook is similar to our analysis of Gauss Legendre quadrature ii The theory of Gaussian integration is based on the following theorem 0 Let qx be a polynomial of degree N such that b qxpxxkdx 0 VII 69 a Donald G Luttermoser ETSU where k 12 N 1 and px is a speci ed weight function 0 Call 061 062 xN the roots of the polynomial Using these roots as grid points plus a set of weights 101102 wN we construct an integra tion formula of the form Ab fltxgtpltxgt dx z w1fx1 mom VII70 0 There is a set of w s for which the integral for mula will be exact if f is a polynomial of de gree lt 2N iii The weights can be determined from the Three Point GaussianLegendre Rule For example consider the interval 1 l with px 1 This gives us a Gaussian Legendre formula For inte grals in the interval 1 b it is easy to transform them as ib a meldx 2 1fzdz VII 71 using the change of variable 06 b a b 12 o The rst step is to nd polynomial We want a three point rule so that qx is a cubic 106 co 0106 02062 03x3 VII 72 0 From the theorem of Eq VIII 69 we know that lqux 1Qdx1x2qxdx0 VII73 VII 19 VII 20 PH Y5 4007 5007 Computational Physics 0 Plugging in and doing each integral we get the equations 2 2 2 2 2 2 7 7 7 7 7 0 VII 74 00362 361363 3605C2 l o The general solution is co 001 a02 003 513 where a is some constant This arbitrary constant cancels out in the second step 0 Using this in Eq VIII 72 gives a polynomial solution of 5 3 106 i063 Ex VII 75 Notice that this is just the Legendre polynomial 0 Next we need to nd the roots of qx 13306 This cubic is easy to factor The roots are 961 35x2 0 and x3 Using the grid points in Eq VIII 70 gives 1 M dx m w1fMw2f0w3fW VII76 0 Finally to nd the weights The above formula must be exact for 106 x5 We can use this to work out values of wl wg and 103 It turns out to be su icient to consider just 106 and x2 which produces 3 equations 2 LU1 LU2 LU3 0 35w1 35w3 VII77 2 3 3 3w13w3 o This linear system of equations is easy to solve giving wl 59w2 89w3 59 Donald G Luttermoser ETSU iv c You will typically never have to write your own Gaussian quadrature subroutines since they are included in math libraries and built into some programming languages 6 9 the IDL functions NT2D and NT3D the QGAUS An alternative way of nding the weights is to use the identity 2 1 xddxPNxz2 where N 3 in our example This formula may be derived from the recurrence relation for Leg endre polynomials w VII78 There are various advantages and disadvantages in using Gaussian integration Advantage A very high order accuracy is ob tained for just a few points often this method yields excellent results using fewer than 10 points Disadvantages l The node points and weights must be computed or obtained from tables This step is nontrivial if you want to use many node points Using more than N 20 points is rarely worth it since badly behaved functions will spoil the results in any case 2 Unlike Newton Cotes integration the method does not lead itself to iteration nor is it easy to estimate the error subroutine in Numerical Recipes VII 21 PHYS40075007 Computational Physics Course Lecture Notes Section VI Dr Donald G Luttermoser East Tennessee State University Version 41 Abstract These Class notes are designed for use of the instructor and students of the course PHYS 40075007 Computational Physics taught by Dr Donald Lutternioser at East Tennessee State University VI Methods of Data Fitting A Methods of Fitting Local Data 1 There are two main types of data tting a 0 Local data tting means that you are just tting a few points in the complete population of data points pre sented We will look at three types of local data tting in the class Interpolation Lagrangian Interpolation and Cubic Splines The second type of data tting is global data tting Here we t the entire population of data points to an assumed analytic function We will investigate two types of global data tting Fitting data to Gaussian pro les or distributions and to linear equations using the method of linear least squares 2 Linear Interpolation 3 b C There are times when we may wish to know the value of an unmeasured datum that lies between two measured data in a population The easiest way to do this is to draw a straight line between the two known data points and t the unmeasured point on that straight line Let s say that we have two measured data P1061y1 and P2x2y2 and we wish to determine the y Value of the data that lies on x such that 061 lt 06 lt 062 To determine this we just set up the following ratio y y1y2 y1 x x1 xg x1 m VI 1 3 PH Y5 4007 5007 Computational Physics 1 Then solving for y we get yy1y2 y1 x2 06 061 111 m 06 061 VI 2 e There may be times where you need a regularly spaced grid of data points e x x constant where j z39 l to pass to another function such as an integration routine see Vll of the notes You would then need to interpolate your data to determine your data s y values at the regularly spaced x values from your original uneven data Such is the time when you would need a linear interpolation procedure The IDL function INTERPOL will carry out such an operation for you Lagrange Interpolation a Suppose that the data on either side of the unknown point we are interested in is not following a linear distribution but instead is following a curved path Then a straight line is not the best representation of your data b In such cases it is better to t to a polynomial instead of a straight line This leads us to an interpolation scheme developed by Lagrange c One can t an n l order polynomial 2 a0 1106 12062 17406714 06 z 062 VI 3 to n values of the function f evaluated at points 06 d The formula is written as the sum of polynomials 1055 f11x f2206 f 39 39 39 annW a VI 4 Donald G Luttermoser ETSU 8 f g h where x x x x1x x2 06 06 A 7 7 7 x j l1xi 96139 96239 56106239 062 562 an VI 5 and note that 00 Z A 1 VI 6 For 3 points Eq VI 4 provides a second degree poly nomial while for eight points it gives a seventh order polynomial It is easy to see that for 2 points Lagrange interpolation just reduces to the linear interpolation de scribed above The Lagrange interpolation has Q constraints that the data be evenly spaced The difference between the value of the polynomial eval uated at some 06 and that of the actual function is equal to the remainder x x1x x2x xn lt n R z n 1 gtltcgt where C lies somewhere in the interval of interpolation VI 7 but is otherwise undetermined Lagrangian interpolation is typically used to t only local data though it can be use to t the data globally How ever there are better ways to t nonlinear global data gt polynomial least square tting which we are not going to have time to cover in this class 4 Interpolation with Cubic Splines a A sophisticated variation of the n 4 Lagrange interpo lation known as cubic splines sometimes offers superior ts to nonlinear data m C d e f g h PH Y5 4007 5007 Computational Physics In this method cubic polynomials are t to the data points representing some polynomial function in each in terval with the constraint that the rst and second deriva tives of the polynomials be continuous from one interval to the next This continuity of slope and curvature makes the spline t particularly eye pleasing The basic approximation of splines is the representation of the function f in the subinterval mlxZH with a cubic polynomial x 2 1206 for x g x g gel1 VI 8 mac 12 Mac xi m xi f3x gel3 v19 where fa f 2 and f3 are the rst second and third derivatives of the function f evaluated at 06 The computational chore is to determine these derivatives in terms of the N tabulated values of flu The matching of f from one interval to the next at the nodes provides the equations fixZ1 fi1xi1 i l N l VT10 Also the rst and second derivatives must be continuous at the nodes n mmm Emm m WHD Finally an additional equation is needed to determine all of the constants For this the third derivatives are found Donad G Luttermoser ETSU by approximating them from the second derivatives f l f5 06i1 06239 ff 2 VI12 i In practice it is better to use a cubic spline routine sup plied by the programming language 69 SPLINE in IDL or math library B Data Fitting Gaussian Distributions 1 The most important probability distribution for use in statistical analysis of data is the Gaussian or normal error distribution a Mathematically it is an approximation to the binomial distribution for the special limiting case where the num ber of possible different observations 71 become in nitely large and the probability of success p for each is nitely large so that np gtgt 1 b Physically it is useful because it seems to describe the distribution of random observations for most experiments c In spectra it describes the shape of many spectral lines This type of data tting can be considered both a global and a local data tting technique depending on what is being tted In the case of spectra one often nds either absorption lines and or emission lines superimposed on a gradually varying continuum Often it is useful just to isolate the individual lines from the continuum and t these lines to the proper pro le distribution More times than not such lines follow a Gaussian distribution Figure Vl l shows the Gaussian pro le and labels the various components to this distribution The Gaussian distribution func PH Y5 4007 5007 Computational Physics 04 Figure VI 1 The Gaussian distribution function tion is de ned by U PG06MU exp QC VI13 a It is a continuous function describing the probability that from a parent distribution with a mean u and a standard deviation 0 the value of the random observation would be x b The probability function is properly de ned such that the probability dPgx 110 that the value of a random ob servation will fall within the interval day around 06 is given by dPgxia Pgxua dx VI 14 C Gaussian distributions are bell shaped and characterized by the full width at half maximum FWHM which is gen erally called the half width F gt the range of 06 between Donald G Luttermoser ETSU the values which PG is half the maximum value 1 1 PG Hiira aagt EPGCU MUGU 2 v115 T 23540 VI16 d A tangent drawn along the portion of steepest decent of the curve intersects the curve at the e l2 points x Mia and intersects the x axis at points x u l 20 see dashed lines in Figure Vl l P001 i 0 u 0 6 12 P001 u 0 VI17 The mean and standard deviation of the Gaussian distribution of Eq VI 13 are given by the u and 0 parameters respectively The probable error is de ned to be the location on the Gaussian distribution where half of the observations will fall be tween uzl PE Mathematically it can be shown that PE 067450 02865 F VI l8 The integral probability Adm u a is de ned as the proba bility that the value of any random measurement at will have a deviation from it less than 20 where z is the dimensionless range Z lac ua It is found with Mza 1 Z 1 2 A P d 7 51 d Cxnuaa Liza Cxnuaa x m L26 x7 VI 19 and normalized such that Agz oo l VI 20 Note that the smaller 0 is the more precise our observations and or calculations will be The standard deviation 0 is often referred to as the best estimate of x representing the actual value PH Y5 4007 5007 Computational Physics 8 Gaussian Distributions in Nature Absorption and Emission Lines 3 b C d In many cases in nature and in the laboratory 6 9 stellar spectra and emission discharge tubes dark and or bright lines are seen in the spectra objects emit Often the shapes of these lines in either wavelength fre quency or energy space follow a Gaussian shape If a spectral line has a rest also called lab wavelength A0 we can describe a line pro le A such that it is normalized over the entire pro le f MAM 1AA 1 VI21 where AA A A0 VI22 and A represents the wavelength at a speci c point along the pro le i In practice 239 6 computer programming one does not carry the integration from 00 to 00 but in stead just carries the integration across the line pro le 239e lAmax 30 to 50 ii Also note that MAM 1AA u dy VI23 where V cA is the frequency of the line pro le If an object that is emitting a spectrum is moving with respect to the observer the rest wavelength ie line center gets shifted via the Doppler effect AAiA Aoi A0 7 A0 70 z VI 24 Donald G Luttermoser ETSU 8 f where here A is the observed line center A0 is the line center as measured in the lab 1 is the radial velocity of the moving object 0 is the speed of light in the vacuum and Z is a parameter called the redshift i If a line is redshifted A gt A0 and UT gt 0 gt the object is moving away from us ii If a line is blueshifted A lt A0 and UT lt 0 gt the object is moving toward us iii Note that Eq VI 24 is valid only when UT ltlt 6 ie the motion is nonrelativistic We also can write a more general expression for the red shift AA 1vTc 1 z Ac il vTc where here we have included the relativistic ie UT E c v125 correction for the Doppler Effect i Rewriting this relativistic formula we can express velocity as a function of redshift UT 2 12 1 g m VI26 ii You will note that there is no way a velocity cal culated from the relativistic form of the Doppler effect can exceed the of light gt the special theory of relativity is not violated If an absorptionemission line has a Gaussian distribu tion its pro le function follows the form 1 AA VI 10 PH Y5 4007 5007 Computational Physics i The parameter Aka is referred to as the Gaussian or Doppler see below width From a comparison of Eqs VI 27 and VI 13 we see immediately that MS xia VI28 ii With this bit of information we can show that the line pro le half width is AA12 1665 AAG VI 29 iii The Gaussian pro le has a peak value of Aka 9 Sometimes line pro les are a conglomerate of other line pro les just as the total error can be the sum of different individual random errors a In such cases we must convolve the individual pro les together with the convolution integral mm mm x we dx VI30gt b The convolution of two Gaussian pro les with AAGl and Akgg at the same A0 yields another Gaussian pro le with AVG A m AVG 10 Spectrographs are a hardcopy version of the Fourier trans form function gt it breaks white light up into its component colors a From this analogy the spectrograph acts as a convolution integral convolving the instrument pro le with the pro le from the emission source b Since good spectrographs are designed to have Gaus sian instrument pro les the Gaussian from the source is Donald G Luttermoser ETSU 11 c d 8 VI 11 convolved with the Gaussian of the spectrograph Spectrographs with AA12 lt l A l A 10 10 m are said to be high resolution or high dispersion spectrographs Spectrographs with AA12 gt 1 A are said to be low resolution or low dispersion spectrographs AA12 effectively determines whether two spectral lines separated by AA can be split into individual lines i If AA lt AAlg the two spectral lines will appear as one in a spectrum ii If AA gt AAlg the two spectral lines will appear as two individual lines Line pro les emitted from a gas tend to have Gaussian shapes due to the fact that the atoms and molecules that make up the gas are randomly moving in the gas gt they are said to be following a Maxwellian or Maxwellian Boltzmann distribution a b c As individual particles absorb and emit photons the pho tons are shifted via the Doppler effect based on their rel ative thermal motion with respect to the observer As such Gaussian pro les are often referred to as Doppler broadening of a spectral line Since the particles are following a Maxwellian distribution in velocity space they have a Gaussian distribution in wavelength or frequency space where the Gaussian or Doppler width is equal to the thermal Doppler width vZkT mg AAD mm A0 VI31 VI 12 12 PH Y5 4007 5007 Computational Physics where k 138062 X 10 23 JK is Boltzmann s constant T is the temperature of the gas 0 is the speed of light and m5 is the mass of the radiating particle 1 As a result of this one can determine the temperature of a gas by measuring the half width of a line arising from a gas after of course one has subtracted off the instrument pro le e In practice life isn t quite this easy as is shown in the Astrophysics course ASTR 34l5 Most of the math oriented programming languages have func tions that will t data to a Gaussian pro le In IDL there are 2 intrinsic functions dealing with tting Gaussians a The GAUSSFIT function computes a non linear least squares t to a function f with from 3 to 6 unknown param eters f is a linear combination of a Gaussian and a quadratic the number of terms is controlled by the key word parameter NTERMS The syntax for this function is Result GAUSSFTX Y A ESTIMATESarray NTERMSinteger3 to Note that the square brackets here means optional parameters i X An n element vector of independent variables data ii Y A vector of dependent variables same vector length as X Donald G Luttermoser ETSU iii A A named optional variable that will contain the coef cients A of the t see below iv ESTIMATES Set this keyword equal to an ar ray of starting estimates for the parameters of the equation lf NTERMS keyword is speci ed the ES TIMATES array should have NTERMS elements If NTERMS is not speci ed the ESTIMATES array should have 6 elements If the ESTIMATES ar ray is not speci ed estimates are calculated by the GAUSSFIT routine v NTERMS Set this keyword to an integer value between 3 and 6 to specify the function to be used for the t The values correspond to the functions shown below In all cases 06 A1 2 A2 o NTERMS6 an AOe ZQ2 A3 A496 A5962 NTERMS6 is the default setting Here A0 is the height of the Gaussian A1 is the center of the Gaussian A2 is the width of the Gaussian A3 is the constant term A4 is the linear term and A5 is the quadratic term NTERMS5 x A052 A3 A406 NTERMS4 my Aoam A3 VI 13 VI 14 PH Y5 4007 5007 Computational Physics NTERMS3 my A062 Note that this last value just ts the data to a Gaussian pro le b Result GAUSSINTX evaluates the integral of the Gaus sian probability function and returns the result The Gaussian integral is de ned as GAUSSINTX 2 Wm 7139 700 If X is double precision the result is double precision oth erwise the argument is converted to oating poin and the result is oating point The result has the same structure as the input argument X C Data Fitting The Method of Least Squares 1 Introduction Here we will discuss tting data to a straight line See your textbook for a more general discussion of Least Squares for higher order equation ts a If the data measurements are a linear function of some other variable we can use the method of linear least squares t b The data will contain errors in their measurement as such we will need to calculate a goodness of t c Let s say that y in a linear sense then y a1 I2 I VI 32 then our task is the gure out 11 to the y intercept 12 to the slope and their uncertainties from the data Donald G Luttermoser ETSU on Note that some authors use the y A Bx notation and others use y b mac notation These of course are all the same equation with 11 A b e y intercepts in the three equations and a2 B m e the slopes The method of least squares also is called linear regres sion 2 The Goodness of Fit X2 3 b c The chisquared test also is called the Method of Max imum Likelihood It is the assumption that the observed set of measurements is more likely to have come from the parent distribution than some other distribution or ran dom event The X2 is determined with X2 ltgt2 T 02 391 02 N l H l g y a1 1206i2 VI 33 2 H MZ l l H 2 where N represents the total number of data measure ments The best t to the assumed functional dependence will occur when X2 is at a minimum 8X2 7 8am T where M is the maximum number of terms in our poly 0 m 1 M VI34 nomial t VI 15 VI 16 PH Y5 4007 5007 Computational Physics 1 Since we are only interested in straight lines here M 2 so a 2 a N 1 2 EX 871i 712yi a1 a2xi N 2 ZgQJi ai a i 0 VI 35 2 1 i 97 722yi a1 a2xi N Z 02yi a1 a2xi 0Vl 36 il i Note that if the uncertainties of each measurement are the same we could just label 02 a and take it outside the summation sign At this point we will not make that assumption and leave the uncertainty in the summation e These equations can be rearranged to yield a pair of si multaneous equations N N N yi a1 06239 Z i Z 7 Z 12 i1022 i1 039 i1 039 N 1 N Z 01 Z 7 02 Z 7 i1 0392 il 0239 N 06 y N 06 062 2 117 D1242 21 0392 il 2 21 0392 N 362 N x2 11 Z 7 12 Z i VI38 239 0239 i10i f Solving Eqs VI 37 and VI 38 for 11 and 12 will give the coef cients for which X2 the sum of the square of the deviations of the data points from the calculated t is a minimum One way to solve these equations is to use the method of determinants a1 1 mata mama A EWWi Mad0 Donald G Luttermoser ETSU 3 VI 17 v139 A 12 1 ma 2yiUz392 A 296202 290iyi022 1 N 1 inyi A2 2 i1aii1ai I v141 Combining Statistical and Experimental Errors 3 b In the equations above yz corresponds to individual mea surements each with their own statistical errors gt statistical uctuations in the collection of nite numbers over a nitely long interval of time gt 0iyi Note that there also can be uncertainties in our indepen dent variable 062 i For example yz might be displacement measure ments and 062 the times the measurements are made ii Since the time measurement also has an uncer tainty 02062 we have 2 uncertainties to worry about iii 02xi is referred to as the experimental error VI 18 PH Y5 4007 5007 Computational Physics c As such the ms in Eqs VI 39 through VI 41 are the combination of these 2 uncertainties which is found with 022 0306239 fl2a VI42 4 If we known the statistical and experimental errors or have de termined an approximate a from the sample variance from your tted function the theory of least squares then gives an expres sion for the variance or uncertainty in the deduced parameters 2 SM 2 E al A 7 02 A v143 5 A measure of the dependence of the parameters on each other is given by the correlation coef cient pa1a2 COVlta1 a2 VI 44 05100 I v1 45 COVG1702 A l a Here cova1 a2 is the covariance of a1 and a2 and vanishes if a1 and a2 are independent b The correlation coe icient pa1 a2 lies in the range 1 g p s 1 i Positive p indicates that the errors in a1 and a2 are likely to have the same sign ii Negative p indicate opposite signs 6 The preceding analytic solutions for the parameters are of the form found in statistics books but are not optimal for numerical calculations because subtractive cancellation can make the an swers unstable We can improve the stability of the calculations Donald G Luttermoser ETSU by rearranging some terms and introducing the averages of the measurements y to the solution 31 l1 y a 1287y VI46 N 062 112 N xi EZ St El 5 Smltw47gt 21 0i i1 0i 1 N 1 N i 7 Z 7 Z VI48 06 Ngx y N221 7 Math oriented programming languages always have least square functions or subroutines available either as a built in command or in a math library a The books titled Numerical Recipes show how to code your own subroutine for least squares should you need to do it See FIT and LFIT in any of the Numerical Recipes books b IDL has built in functions LNFT and REGRESS See the IDL Help utility to see the details of these two functions 8 Though we won t do it here there are least square methods to t non linear functions In IDL see the internal functions CU RVEFT Non linear least squares t to a function GAUSSFT Fits sum of a Gaussian and quadratic POLYFT Polynomial least squares t SFT Determine a polynomial t to a surface VI 19 PHYS40075007 Computational Physics Course Lecture Notes Section II Dr Donald G Luttermoser East Tennessee State University Version 41 Abstract These Class notes are designed for use of the instructor and students of the course PHYS 40075007 Computational Physics taught by Dr Donald Lutternioser at East Tennessee State University II How To Select a Programming Language A Which is the Best Programming Language for Your Work 1 Selection of a programming language for a given task can be a daunting task One must weight the availability of the language the functions it is capable of and the speed and debugging ca pabilities offered by the compiler a b C d 8 Different elds of study typically have their favorite pro gramming language Many mathematicians use Mathematica Maple and Math lab ln physics and astronomy most of the large codes are written in Fortran 77 with the remaining small percentage written in either the newer Fortran 90 or C very few are written in C Astronomers who use a lot of NASA space data typically use the Interactive Data Language IDL for their work Since we don t have time to cover each and every pro gramming language on the market we will intercompare the three most commonly used in physics and astronomy Fortran 77 C and IDL i The details of Fortran 77 and an overview of For tran 9095 can be found in Appendix A of these notes ii The details of C and an overview of C can be found in Appendix B PHYS 40075007 Computational Physics iii The details of IDL can be found in the next sec tion ie lll of these notes f Any coding introduced in these notes will be written in ei ther IDL or Fortran 77 hence Appendix A and lll should be studied carefully so you fully understand what is being written in these code segments Most programming languages have analogous elements in their programming style which can be reduced to 8 sections as shown in Table 11 1 partially lifted from Appendix A of your textbook No matter which programming language you choose you need to come up with a style and be consistent throughout the code a You do this not only to make it easier for other users to gure out what you are doing it makes it easier on you when you go back to a code after you haven t looked at in a long time b Always have numerous comments throughout the code describing What is being done c For my style I use lower case letters except for the rst letter of words at the beginning of a sentence for com ments and upper case letters for the coding itself Though in Table 11 1 1 use lower case letters for demonstration purposes only 1 Try to use a lot of subroutines or procedures and func tions gt make your code as modular as possible PAN DORA has over 4500 subroutines in it Donald G Luttermoser E TSU Table II 1 Analogous Elements in Fortran C and IDL Fortran C IDL program f function fX subroutine fX Xy X do 2 k1kmaX1 do 2 kkmaX1 1 realer X y integer i j realer y100100 integer ii100100 data mn00 character jan read 9 X1 write 9 7Enter radius7 write 9 7radius 7 radius open7 le Xdat status new write7 f102 Din read7100 Din if den eq 0 then write 9 7Trouble7 endif if den eq 0 then X 0 else X X1den Program Structure main double fdouble X void fX Operations WWWJN X yy for k1kltkmaXk for kkmaXkgt1k7 Data Type Declarations double X y int i j double y 100 100 int ii 100 100 de ne mm 00 char jan Input and Output to Screen scanf lf ampX1 printf Enter radius printf radius f radius Input and Output to Files FILE fout out fopen Xda 77 W fprintffout lf Din fscanf n lf ampDin 7 Control Structure if if den printf Trouble Control Structure if else if den 0 xQ else X X1den no main program structure function f X pro f X XAy X for k0kmaX 11 do for kkmaX 10 1 do X 00d0 amp y 00d0 i0ampj0 y dblarr100100 ii intarr100100 defsysv 7lmn 00 readonly jan 7 NULL string read X1 print 7Enter radius7 print 7radius 7 radius openw 7 7Xdat7 printf 7 Din format f1027 readf 7 Din if den eq then begin print 7Trouble7 endif if den eq then begin X 0 endif else begin X X1 den endelse ll 4 PH Y5 4007 5007 Computational Physics Table ll 1 cont Analogous Elements in Fortran C and IDL Fortran C IDL Control Structure for do 2 k1kmax1 for k1kltkmaxk for k0kmaX 11 do begin term rterm term rterm term rterm 2 continue sum sumterm endfor Switch goto 127 choice switch choice case choice of 1 series rst case 1 1 series rst 2 series second series rst 2 series second break else case 2 endcase series second break B The Fortran Programming Language 1 Fortran was rst developed in 1957 by IBM and was the rst high level programming language a Its name is derived from FORmula TRANslation b Fortran has evolved considerably since 1957 gt Fortran l gt Fortran 11 gt Fortran 111 gt Fortran 66 Fortran 1V gt Fortran 77 Fortran V gt Fortran 90 gt Fortran 95 c Fortran is a sequential programming language there is a logical ow to the coding the rst line of the code is operated upon rst the 2nd line second through the last line last 1 Most large number crunching codes in science are written in Fortran 77 which came on the market in 1978 Donald G Luttermoser E TSU 2 Fortran has a number of features that make it useful to the sci enti c community a Portability Portability means that it is easy to move from one machine to another While no large program can move from one machine to another totally without dif culty problems are usually minor provided the For tran 77 standard has been strictly adhered to b Availability Fortran is probably more widely used than any other programming language in the academic com munity It is certainly the most widely used by scientists and engineers Compilers are available for all machines from micro to mainframe likely to be encountered in a university c Ef ciency Partly because of its simple static use of store and the un complicated nature of its constructions and partly because of a massive investment in compiler development numerically intensive programs written in Fortran are generally faster than those in any other high level language The technique of optimization by which a compiler alters the machine code it generates without changing the end result to speed the calculation is more developed for Fortran than any other language 1 Fortran Libraries Very considerable effort has gone over many years into the production of libraries of rou tines that may be called from Fortran programs Chief amongst these are libraries of numerical analysis routines such as the very large NAG library and routines to do graph plotting such as the GINO Dl BOOO and UNIRAS libraries 3 PHYS 40075007 Computational Physics Fortran is designed to be modular a One has a main program that drives the code and calls subroutines and functions b A good Fortran code typically has a short main program which call many subroutines c Subroutines and functions can call other subroutines and functions C Running Fortran on the PCs 1 Some of you may have copied over the free ware GNU Fortran 77 compiler If you plan on using this software read the readmetxt for information on installing and running this compiler note that this g77 compiler must be ran from DOS mode This package contains the SLATEC mathematical library which is similar to the LAPACK library see Appendix A For those of you running a Fortran 7790 code from the PCs in Brown Hall Digital Visual Fortran 5 has been installed on these machines Simply double click on the appropriate menu item in the Start menu or the Developer Studio icon on the main screen One then creates a Fortran code in the editor window and used the various buttons and pull down menus to compile and run the code You will have to retrieve the LAPACK or BLAS libraries from the netlib web site if you wish to use them Then use the Help utility search on math library for instructions on using external libraries To compile and run on a PC using Microsoft39s Development Studio version of Digital s Visual Fortran do the following tasks Donald G Luttermoser E TSU 3 Initially there are one of two things you need to do to get a Fortran code running on a PC with MS Windows If this is the rst time trying to write compile and run a code select the File menu item on the toolbar then follow the next recipe i Select New which will bring up a new GUI ii In this new GUI select QuickWin Application from the Projects folder iii Enter the project name choose whatever name you feel is appropriate do not put spaces in the name make sure the Create a new Workspace but ton has been selected and nally pick a location for the workspace directory that is unique to you make one with the Create New Folder button iv Select New from the File menu in the upper left of the main GUI if your Fortran code has not been written yet and select File from the menu items If you are importing a le which was created somewhere else select Open instead of New v This brings up an editor in which you can type in your code Once done select Save As from the File menu item and save your le to an appro priate lename usually the name of the Workspace name that you selected above If your le had previously been saved to the le name that you selected just select the Save menu item under File b PHYS 40075007 Computational Physics In the main program just after the rst PROGRAM startement include the the line USE DFLIB Then after the declaration and common block statements vi include the line IOLDCT SETTEXTCOLORRGB00FFFFFF Take these lines out if you port your code to a Unix machine vii These two lines are needed in order for Microsoft Windows to bring up a text window and to change the text color to white so that you can read the lines output to this window with WRITE and PRINT statements If you have already set up your workspace and need to do further editing of your program do the following i Highlight the File menu item on the toolbar select the Open Workspace menu item This will bring up a new GUI box ii Find and select workspace name which will have a dsw as the le name suf x iii The Fortran le or les you will need should already be present in this workspace If not open the needed le or les in the File menu item iv Once you have made any editing changes save the le then select the Build item on the top menu bar and compile your les rst any supple mentary Fortran les you have in the workspace then the le that contains the main program sec Donald G Luttermoser ETSU ll 9 tion Then go back to the Build item and se lect Build runnameexe where rmmame is the name of the executable le that the compiler made c Finally to run your code click the exclamation mark 1 above the editing window to run the code A textbox window should appear that contains any output that you made to the terminal usually a PRINT or WRITE6 100 where 100 is an arbitrary FORMAT statement lDl i Note that Microsoft Windows will immediate close the textbox window when the code nishes execut ing ii To prevent this from happening issue the follow ing commands just prior to the last executed line in your code typically a STOP statement PRINT Hit ltreturngt to Clear window READ 110 ASK 110 FORMATA where ASK needs to be de ned as a CHARACTER variable and the 110 number in the FORMAT state ment would be some number that you have not yet used as an ID in that portion of your code either the main program or a subroutine wherever STOP is located This way Microsoft won t close the window until you hit ltreturngt on your keyboard hence allowing you to look at the output D Running Fortran on the Linux Workstation 1 To write a Fortran code under Linux use the following recipe a Change directories to whatever subdirectory you want as your working directory ll 10 2 b c The Linux operating system uses the GNU Fortran compiler called g77 Read the man pages on fester for a listing of compiler PHYS 40075007 Computational Physics Let s say we wish to create a code called atlasf or we already had one by that name lssue the command at the Unix prompt emacs atlasf This will bring up the emacs editor window and start writ ing your code When done click the Save button fol lowed by the Exit button if you don t need to make any more modi cations Let s say that we edit a second le which contains some operating system speci c command like the GETNODEO function that is available on some Fortran compilers called atlasunixf options Compiling and linking the code are performed with one command again at the Unix prompt a b C g7 o atexe W atlasf atlasunix1c Here the 0 ag tells the linking software to save the executable code in a le called atexe if 0 was not in cluded the executable would be saved in a le named aout The W capital W tells the compiler and linker to issue extra warnings such as oating point over ows and underflows during operation see Appendix A Note that GNU Fortran does allow some Fortran 90 con structs when using the ff90 option switch and VAX For tran constructs when using the fvxt option switch See the man pages for further information Donald G Luttermoser E TSU 4 ll 11 Unlike the PC environment whatever output is sent to the ter minal screen will remain after the code nishes execution so you don t have to include the READ command that we described in the PC section of this section E The C Programming Language 1 C like Fortran 77 is referred to as a sequential programming language gt the ow of the operations carried out by the com puter line by line from the beginning of the main program to the end of it Like Fortran C is considered a higher level language but it has some lower level language functionality gt it can access the operating system much easier and more e iciently than Fortran 7790 a A C program functions executable code variables b c 01 data Every program must have a function called main execu tion starts there main may call other functions These functions can be taken from one of the following sources i As written by the programmer in the same le as main ii As written by the programmer in separate les iii Called from prede ned libraries C has no concept of a Program as in Fortran just func tions The C view is that main is called by the operating system the operating system may pass arguments to main as well ll 12 PH Y5 4007 5007 Computational Physics as receive return values 69 return 0 however above we choose to make main take no arguments void e is a comment and is ignored by the compiler the C standard says that comments cannot be nested Also comments must be terminated explicitly 239e new line does not terminate them this is a common source of compiler errors and for the unwary can be very dif cult to trace For your own sanity the header comment should include i Name of the program to correspond to the le name ii Authors name even if copiedl iii Date iv Brief indication of function and purpose of the program 69 assignment number or project and what it does 3 Appendix B of these course notes contains some details and a tutorial in programming with the C language One thing to note is that C is not designed to be a number crunching program ming language It has very limited mathematical operators and one must load additional math libraries at compile time to do any math beyond simple arithmetic see below If you are going to be crunching numbers then Fortran is typically the code of choice 4 Many of you are already familiar with the objectoriented pro gramming language known as C The University has a C Donald G Luttermoser E TSU compiler which you are free to use They do not have a C com piler As such if you wish to program in C you will have to carry out your calculations on the Linux workststion see below F The Origin and History of the C Programming Language 1 10 C is a member of the ALGOL family of programming languages Ancestry ALGOL 60 CPL BCPL B and ALGOL 68 CPL was developed at Cambridge University around 1963 for compiler writing BCPL was derived from CPL as a simple systems language by J Richards at Cambridge 1967 B derived from BCPL Bell Labs ATamp T late 1960s for developing the Unix operating system incidentally for a PDP 7 computer which had 4K 18 bit wordsl developed by Ken Thompson 1970 C 1972 by Thompson amp Dennis Ritchie Successor C Stroustrup ATamp T 1986 superset object oriented stronger type checking it should be noted that the development of C in uenced ANSI C ANSI C 1988 1991 ISO C 1994 ISO C Amendment 1 ll 13 ll 14 PHYS 40075007 Computational Physics G Running C on the Linux Workstation 1 When compiling C code on fester issue the command This will take the C code stored in filenamec compile it and link gcc o filename filenamec it storing the executable in le filename 3 b c d e 2 If your C code uses functions from the standard C math library The gcc command invokes the GNU C compiler Note that if the output option 0 filename was not given the executable would have been stored in le aout However your compilation may not go well An Error indicates a serious flaw in the source code which the com piler cannot circumvent In the event of any error the compiler can produce no usable object code Warning as the name suggests is the compilers way of saying are you sure about this examples are variables de ned but never used fairly benign functions called without any prototype declaration present de nitely not benign unless the types of the actual parameters are identical to the types of the formal parameters you could be in for a nasty shock when the software executes I maintain that C compiler warnings must never be ig nored and insist that students heed this principle issue the command on fester where the lm tells the compiler to load the math library to the gcc o filename m filenamec linker as well PHYS40075007 Computational Physics Course Lecture Notes Section I Dr Donald G Luttermoser East Tennessee State University Version 41 Abstract These Class notes are designed for use of the instructor and students of the course PHYS 40075007 Computational Physics taught by Dr Donald Lutternioser at East Tennessee State University Choosing Hardware and an Operating System A Introduction 1 With the widespread availability of fast computers at affordable prices it has become important to understand the capabilities and details of numerical modeling 2 The way I see it there are four different classes of computer hardware From most to least expensive they are a Super Computers A mainframe with many 64 bit CPU chips vector processors Cray used to be the biggest maker of super computers but they have now gone out of business due to the speed and low cost of the PCs There are still a few Crays in operation around the country b Mainframes Depending on the age of the machine ei ther a 32 or 64 bit CPU running some multi user oper ating system The old VAXes various models 750 780 8600 etc and the IBM 360 are two examples of main frames It is hard to nd machines like this any more since they have become virtually extinct due to the ad vent of workstations and fast PCs c Workstations A machine that runs a multi user oper ating system like Unix or OpenVMS They are essentially scaled down versions of the older mainframe computers They typically have a 64 bit CPU processor or series of processors though some older machines 69 VAXsta tions DECstations and Sparc Stations have only 32 bit chips in them On these machines numerous users can operate the machine remotely and at the console at the ll PHYS 40075007 Computational Physics same time Examples of such machines are the Sun work stations eg the Ultra Stations AlphaStations originally made by Digital Equipment Company which is no longer in business and Silicon Graphics workstations Note that the PC world has been calling their PCs running Win dows NT2000XP workstations but in reality they are not since multi users cannot log in simultaneously not at least without a lot of system administration work 1 Personal Computers PCS Typically a 32 bit though some are now 64 bit CPU processor or series of proces sors running an operating system OS of either some version of Microsoft Windows Linux with X Windows or Mac OS X on the Macintoshes Besides the Macs ex amples of such machines include Dell and HP PCs These machines are designed for single user mode ie one at a time 3 The size of the CPU can be as important as the clock speed of the CPU a Early PCs only had only an 8 bit CPU chip They evolved to 16 bit and now to the current 32 bit Intel and Macin tosh have recently developed a 64 bit chip for the PC class machine and they are currently available Microsoft has recently developed a 64 bit version of Windows for these chips 64 bit versions of Linux have been available for the past few years The Power Mac G5 PC has dual 64 bit CPUs in them running Mac s version of Unix called Mac OS X b Their are two reasons why the size of the CPU is impor tant Donald G Luttermoser E TSU i The maximum and minimum size of a number is limited by the chip size ii The amount of precision possible for a number is also limited by the size though some compilers are sophisticated enough to increase precision for a given chip beyond what the chip can normally provide 4 Computers work with numbers and any character for that mat ter in binary format Two bits Binary ding 0 E 013 and 1 E on a There are 2N integers that can be represented with N bits b The sign of the integer is handled with the rst bit 0 E positive number which leaves N 1 bits to represent the value of the integer c Therefore N bit integers will be in the range absolute value wise of 0 2N4 i Hence an 8 bit machine can handle integers in the range 128 127 ii A 16 bit machine can handle integers in the range 32768 32767 iii A 32 bit machine can handle integers in the range 2147483648 2147483647 iv Finally a 64 bit machine can handle integers in the range 9223372 x 1018 9223372 x 1018 Note I am using real notation here due to the d 8 PHYS 40075007 Computational Physics large size of the number on a computer such a number would have no decimal point Also note that in computer programming languages numbers in scienti c notation use the E D for double precision see below notation 9223372E18 9223372E18 v Calculating or expressing integers that are smaller if negative or larger if positive than the given chip size would result in an underflow negative or overflow of the computer register and typically crash one s program gt the bigger the chip the better it is for numerical work Since binary strings of numbers are not easy for people to work with a compiler is used to translate the binary numbers of the machine to either octal base 8 instead of base 2 decimal base 10 what we normally are used to or hexadecimal base 16 numbers English letters and punctuation marks are processed on many computers as an 8 bit number called ASCII Amer ican Standard Code for Information Interchange format note that the rst sign bit is not used in ASCII proto col Each character on the keyboard including their cap ital representations has an ASCII value associated with it ranging from 0 to 127 eg 0 zero has an ASCII value of 48 A capital a 65 and a 97 Many of the ASCII code numbers between 0 and 127 are nonprint able control characters I will often refer to text les in the Microsoft world that is those les that can be printed or viewed with an editor as ASCII les as they are called in the VMSUnix world Donald G Luttermoser E TSU 5 Numbers and character strings are stored on computers in words where the word length is often expressed in bytes 1 byte 1 B E 8 bits 8 b 1 1 a Conventionally storage size is measured in bytes or kilo bytes KB megabytes MB and gigabytes GB b However here kilo does not mean 1000 instead it is equal to 1 KB 1 K 210 bytes 1024 bytes 1 2 c Confusingly many machines measure their memory size in units of 1 2 kilobytes called blocks or more precisely 512 B 29 bytes 512 bytes 1 3 note that your textbook has a typo in it for Eq 23 512K should read 512B 1 As mentioned above for ASCII format 1 byte is the amount of memory required to store a single character This adds up to a typical typed page requiring N 3 KB Data types on computers There are two basic types of data that are operated on and stored by computers gt integers and real numbers Depending on the programming language one is using there are various types of integers and reals a Integers are stored as N71 I sign gtlt 0112 1 4 n0 where an takes on either a 1 value the bit is set or 0 value the bit is not set For example on an 8 bit b PHYS 40075007 Computational Physics machine we can write 57 as sign bit set to 1 X 1 X 20 0 X 21 0 X 22 1x231gtlt241x250gtlt260gtlt27l 1008163200 57 or in binary format where the bits are listed in the opposite order of the summation above that is the least signif icant digit 20 is recorded on the far right side of the binary number just as it is in decimal notation 57 1 00111001 where the rst bit is the sign bit 1 negative Integers can come in a variety of avors in various pro gramming languages i Logical 1 bit word values FALSE 0 TRUE 1 maximum value 1 ii Short Integer 8 bit or one byte word length range l 128 127 maximum value 27 1 number of digits 3 iii Character Like a short integer typically con taining the ASCII code of the character however the programmer in a higher level language like For tran or IDL uses the string notation characters surrounded by single and or double quotation marks and the compiler changes these characters to the short integer ASCII code then to binary to which the machine then operates iv Integer 16 bit or two byte word length range 32768 32767 maximum value 215 1 number of digits 5 Donald G Luttermoser E TSU l 7 v Long Integer sometimes just called Long or Double Precision Integer 32 bit or four byte word length range 2147483648 2147483647 maximum value 231 1 number of digits 10 vi DoubleLong Integer sometimes just called Double Long or Quad Precision lnteger 64 bit or eight byte word length range 922E18 922E18 maximum value 263 1 number of digits 20 Note that these types of integers are only found on 64 bit architectures and even then only found in some programming languages 69 DEC Fortran 90 c Unlike integers real numbers are indicated with a deci mal points located in the number 1 The computer can handle real numbers in two different ways Fixed point not to be confused with xed or in teger data types and floating point notation i In xed point real notation the number x is rep resented as x x sign x 01712 an112n 1 020 014124 15 ii The rst bit is used to store the sign of the number and the remaining N 1 bits are used to store the oz values note that nm N 2 The particular values for N m and n are machine dependent For instance the number 9875 could be represented in l 8 PHYS 40075007 Computational Physics real xed point binary as 0 1001111 gt sign bit x 1 x 23 0 x 22 0 x 21 1 x 20 1 x 2 1 1 x 2 2 1 x 24 8001050250125 9875 on an 8 bit machine iii All xed point real numbers have the same abso lute error of Z m l the term left off the right hand side of Eq 1 5 iv The correspondingly disadvantage is that small numbers those which the rst string of oz values are zeros have large relative errors v Relative errors tend to be more important than absolute errors we will be covering errors in more detail later in the course xed point real numbers are used mainly in special applications like busi ness e In scienti c work the programming language compilers like Fortran use oating point numbers for reals as such reals are often referred to as oats in some languages i In oating point notation the number x is stored as a sign a mantissa and an exponential eld exp d ii The number is reconstituted as 06mm 18 X mantissa gtlt 2exp d bm 1 6 Here the mantissa contains the signi cant gures of the number s is the sign bit still the rst bit Donald G Luttermoser E TSU of the binary number and the actual exponent of the number has a bias added to it 0 Since we have a sign bit the mantissa will always be positive 0 The bias guarantees that the number stored as the exponent eld is always positive of course the actual exponent can be negative 0 The use of the bias is rather indirect For exam ple a single precision 32 bit word may use 8 bits for the exponent in Eq 1 6 and represent it as an integer This 8 bit integer exponent has a range of 0255 Numbers with actual negative ex ponents are represented by a bias equal to 127 a xed number for a given machine Consequently the exponent has the range 127128 even though the value stored in the exponent in Eq 1 6 is a positive number 0 Of the remaining bits one is use for the sign and the other 23 for the mantissa where mantissa m1 X 2 1 m2 X 2 2 quot39 W23 X 2 23 7 17 with the m stored liked the oz in Eq l 5 c As an example the number 05 is stored as 0 0111 1111 1000 0000 0000 0000 0000 000 on a 32 bit machine where the bias is 0111 11112 12710 f PHYS 40075007 Computational Physics iii Typically the largest possible oating point num ber for a 32 bit machine is 0 1111 1111 1111 1111 1111 1111 1111 111 which has the value 1 for all its bits except the sign and adds up to 2128 34 X 1038 Whereas the smallest possible oating point number for a 32 bit machine is 0 0000 0000 1000 0000 0000 0000 0000 000 which has the value 0 for almost all the bits and adds up to 2128 29 x 10 As built in by the use of the bias the smallest number possible to store is the inverse of the largest Just as was the case for integers reals come in a variety of avors in various programming languages For example i SinglePrecision Real sometimes called Floats or Realgtk4 numbers 32 bit 4 byte word size numbers 6 7 decimal places of precision 1 part in 223 and a range 117549435 X 10 38 340282347x 1038 ii DoublePrecision Real sometimes called just Double Precision or Real8 64 bit 2 32 bit words gt 11 bits used for the exponent 1 for the sign and 52 bits for the mantissa 8 byte word size numbers about 16 decimal places of precision 1 part in 252 and a range 22250738585072014 X 10308 17976931348623157 x 10308 iii QuadPrecision Real Realgtk6 128 bit 4 32 bit words gt 15 bits used for the exponent 1 Donald G Luttermoser E TSU for the sign and 112 bits for the mantissa 16 byte word size numbers about 36 decimal places of precision 1 part in 2112 and a range 3362 X 10 4932 1189 gtlt 104932 Currently this precision is only available on 64 bit Sun workstations iv Complex Numbers Some programming lan guages in particular Fortran and IDL have com plex data types A complex number has the form 2 06 iy 1 8 V l and 06 corresponds to the real part of the number and y to the imaginary part where 239 E These numbers are stored as two element arrays real imaginary of real numbers B The Inner Workings of a Computer 1 A computer does mathematics which is the same thing as run ning programs whether your program is designed to math prob lems or not Via machine language a Machine language is built into the CPU of the machines each CPU family of chips has its own speci c machine language or architecture where machine language manip ulates numbers on the bit level i The 32 bit Pentium chips have one type of machine language primarily based on the 0150 Complex Instruction Set Computer ii Whereas the 64 bit chips in workstation class ma chines have their own machine language based on the R150 Reduced Instruction Set Computer ar chitecture b PHYS 40075007 Computational Physics iii Subtle variations exist in the various vendor RlSC chips as such each type of workstation or PC has to have a speci c operating system on it that talks to the CPU in its speci c machine language As such a discussion of computer hardware needs to in clude a discussion of operating systems However rst we need to de ne a few terms about the inner workings of a computer Refer to Chapter 18 in your textbook i CPU The central processing unit is the fastest part of the computer The CPU consists of a num ber a very high speed memory units called registers containing the instructions sent to the hardware to do things like fetch store and operate on data ii FPU The oating point or arithmetic unit is a piece of hardware designed for oating point arithmetic It may be a separate unit called a math coprocessor or it might be built directly into the CPU iii Cache A small very fast bit of memory some times called the high speed bu er that holds in structions addresses and data in their passage be tween the very fast CPU registers and the slower main RAM memory The main memory is also called dynamic RAM DRAM while the cache is static RAM SRAM iv Cache and data lines The data transferred to and from the cache or CPU are grouped into cache lines or data lines The time it takes to bring data Donald G Luttermoser E TSU from memory into cache is called latency v RAM Random access memory or central mem ory is in the middle memory hierarchy This is where your program resides while it is being pro cessed The contents of RAM is lost upon comple tion of the jobs or the turning off of the machine The RAM of your computer is analogous to the memory centers of your brain vi ROM Read only memory contains data that is hard coded on the chip typically non erasable Information on this chip tells the machine its iden tity senses the deVices hooked to the motherboard called peripherals and tells the machine where to nd the boot software The ROM is analogous to your DNA vii Pages Central memory is organized into pages which are blocks of memory of xed length The operating system labels and organizes its memory pages much like we do with the pages of a book gt they are numbered and kept track of in a table of contents viii Hard disk Finally at the bottom of the mem ory hierarchy is the permanent storage on magnetic disks or optical deVices They are slow but can store large amounts of data The coding of the op erating system compilers and any other documen tation whether in ASCII format or binary format is stored here PHYS 40075007 Computational Physics ix Backup devices Since hard disks work very hard they are the rst thing that is likely to go bad on a computer Once the head crashes on a disk drive the data stored there is typically lost forever As such it is good practice to store the contents of your hard drive e hard disk on a more per manent and safe medium magnetic tapes and CDROMs X Virtual memory Virtual memory permits your program to use more pages of memory than will physically t into RAM at one time Pages not currently in use are stored in slower memory typ ically hard disks in a region called swap space and brought into fast memory only when needed The amount of memory space needed for a program to run is limited by the machine s RAM plus the swap space size In Unix the amount of swap space available and total can be found with the swap on some machines called swapon command The speed of the CPU is determined by the clock chip that is attached to it The faster the clock measured in MHZ or GHZ the faster your CPU CPU clock speed doesn t tell the whole story in computer speed the speed of the cache and the speed of the bus that talks to the hard disk and other peripherals is also important in turn around time in running a program Donald G Luttermoser E TSU C A Brief History of Operating Systems 1 A simple de nition of a operating system is the suite of programs that make the hardware usable The operating system manages the CPU disks and IO devices Manipulation of the operating system was typically one of the hardest aspects of learning to use computers which is why the GUI Graphic User Interface mentality has taken hold of the modern computers the user no long needs to talk to the op erating system the GUI does it for you Which is actually a pain in the butt if you need to actually talk to the operating system We now follow the history of operating systems and programming languages that are important in the scienti c community starting in the 1940s a Early computers had no operating system 69 the IBM 604 b The computer was told what to do with a low level set of commands gt assembly language c The IBM 604 could undertake 60 program steps before using punch cards as a backing store 1 Often the user had to manipulate toggle switches to input the code for the mainframe The 1950s saw the advent and rapid changes in the capability of operating systems a Jobs were batched so that the time between jobs was min imized 01 Q 0 C PHYS 40075007 Computational Physics The user was now distanced from the mainframe and en tered a program via a card reader This era saw the rapid development in higher level pro gramming languages such as Algol 60 and the rst version of Fortran The 1960s saw the advent of multiprogramming and the concept of time sharing 3 b Mum user operating systems were developed which al lowed users to actually log into the computer Since scientists were the primary users of mainframes an of cial programming language for the sciences was de clared Fortran N which was also called Fortran 66 The 1970s saw numerous multi user operating systems come on line Two of the most successful were Unix and VMS Virtual Memory System a b C 01 Unix developed by Bell labs in 1970 and further modi ed by the University of California at Berkeley Unix is written in the C programming language is capable on running on a variety of different architectures VMS was released in 1978 by the Digital Equipment Corpo ration DEC to run on their VAX mainframe computers Whereas Unix commands are often cryptic VMS com mands are based on English words and hence was con sidered much more user friendly than Unix As such it was typically the preferred OS by scientists whereas Unix was favored by the computer scientists Donald G Luttermoser E TSU T 00 e The 1970s also saw a new Fortran standard come on line gt Fortran 77 This version of Fortran was much more versatile and powerful than the earlier version of this pro gramming language f Finally the 1970s saw the beginnings of the PC market open up through the release of the Apple He and IBM PC running DOS computers The 1980s saw an explosion of computers in science through the release of the microcomputer which lasted only a few years and workstation Both of these types of machines were are scaled down versions of mainframes with the workstations being de signed to t on a desk a The most popular workstations at that time was the VAX station running VMS and the Sun workstation running Sun OS b During this decade Unix s popularity grew leaps and bounds whereas VMS declined c Later this decade DEC came out with their version of Unix to run on the VAX chip called Ultrix Machines that had Ultrix installed on it were called DECstations 1 Apple came out with a new 32 bit chip that they placed in a new machine called the Macintosh In the 1990s there was a PO explosion which left the workstations in the dust a Microsoft began to dominate this market with their vari ous avors of the Windows OS note that Windows origi nally was a GUI front end that sat on top of DOS 16 bit b C d e PHYS 40075007 Computational Physics PCS were replaced by 32 bit PCS during this decade PCS started becoming very fast when the Intel Pentium chip was developed DEC came out with a new 64 bit chip the VAX chip is 32 bit called the Alpha which gave rise to the AlphaStations This chip required new versions of DEC S OSS VMS gt OpenVMS and Ultrix gt OSFl which later was renamed Digital Unixin 1998 then Tru64 Unix in 2000 when Compaq bought out DEC Sun then came out with the Spare station 32 bit CPU and a new OS called Solaris to run on it This decade saw the last of the mainframes as the work station and PC explosion made this type of computing obsolete Because of this and the fact that their Alpha Stations were not as popular as their VAXstations DEC went belly up at the end of this decade Finally this decade saw the emergence of a free version of Unix called Linux which was designed to run on PC class machines In the 2000s the bulk of computing is done by PCS In the phys ical science community most of the number crunching is still performed on workstations but many scientists are now switch ing to PCS due their fast speeds and low costs a Sun came out with a new 64 bit chip and workstation called the Ultra stations Their latest and quickest work station is the multi processor Blade computers b For workstations the biggest sellers are Suns Silicon Graph Donald G Luttermoser E TSU ics IBM Workstations and HP Workstations c In 2002 Hewlett Packard bought Compaq and has an nounced that they are dropping the AlphaStation from their line This means the death of both Tru64 Unix and OpenVMS d In 2003 Apple came out with its 64 bit dual processor Power Mac G5 running Mac OS X which is a Unix based operating system e In 2004 The Athlon 64 chip by Advanced Micro Devices AM D and the Intel Xeon 64 bit chip have appeared in PCs These machines run either the 64 bit version of Mi crosoft Windows XP or the 64 bit version of Linux D The Unix Operating System 1 Since most of you are well aware of the workings on the various avors of Microsoft s OS ie Windows at this point I will highlight the use of the Unix operating system since you will have to do a little work on one of the Department s Linux computers in Brown Hall 260 a Unix is the operating system of choice for number crunching 64 bit RISC based and CISC based workstations which are popular in the scienti c community astronomy and physics in particular and in technical corporations b There are various flavors of Unix that exist on the market each design for use by speci c computer platforms i Soaris is the avor of Unix that exist on Sun Mi crosystems platforms 69 Sparc Stations and Utra Station You may also run across an earlier Sun operating system called the SunOS This operating PHYS 40075007 Computational Physics system is no longer supported by Sun ii lrix is the version of Unix used on Silicon Graphics workstations iii HPUX is the Unix avor on Hewlett Packard s RISC based workstations iv Digital Equipment Corporation s DEC version of Unix was called Tru64 Unix previously known as Digital Unix and previous to that OSFl before they went out of business v There are at least two platform independent ver sions of Unix SCO Santa Cruz Operating system Unix and BS D OS Berkeley Systems Development Operating System Unix vi Mac OS X of Apple is based on the BSD version of Unix vii Finally the Linux version of Unix is designed to run on a variety of different platforms from Intel microprocessors to SunUltra and Alpha chip archi tectures Linux is one of the few versions of Unix that can be downloaded for free 2 The History of Unix a Brian Kernighan and Dennis Ritchie both work at Bell Labs and were involved in the development of the Unix operating system and the C language Donald G Luttermoser E TSU b c d e f g h 3 Note that the Unix operating system is set up in a hierarchal In 1968 Ken Thompson of Bell Labs wrote the original version of Unix in PDP 7 a DEC mainframe which pre ceded the VAX assembly language In 1969 a language called TMG was ported to Unix the compiler was written in assembly Using this Thompson created B a pared down version of BCPL B was essentially C without types Dennis Ritchie added character integer and oating point types to B anticipating oating point operations in new versions of the PDP This was the beginning of the C programming language In 1973 the Unix microkernel was rewritten in C The KampR White book was rst published in 1978 Brian Kernighan wrote the body of book and Ritchie wrote the technical appendices The ANSI C standard was released in 1983 system as shown in Figure 21 of your textbook a The user and programmer talks to the shell 69 the Bourne shell sh the Korne shell ksh the C shell csh Note that these shells are not designed well for user inter face at the keyboard updown side arrows won t recall previous commands or edit them As such new shells have been developed that allow keyboard editing csh gt tcsh sh gt bash b C 01 PHYS 40075007 Computational Physics The shell talks to the Unix utilities 16 commands like cp and ls The utilities talk to the kernel 16 the actual operating system The kernel talks to the hardware You will be receiving a special handout that will instruct you in logging into one of the Linux workstations in Brown Hall 260 The node name which is important in the Unix world of the Linux machine that you will be using for this class is fester and its complete lnternet name is festeretsuedu Note that you will have to log into this machine at the console since most of the Internet software has been turned off due to outside hackers E The Structure of Unix 1 Unix consist of two main parts a b The kernel controls i Computer memory ii Resource allocation iii Peripheral systems iv Filestore organization v File security and access The shell provides i User interface or command processor Donald G Luttermoser E TSU 00 ii Utility programs c Examples of shell are the Bourne shell sh C shell csh Bourne Again shell bash Korne shell ksh and the tcsh tcsh Why use Unix a Unix is portable b Unix is a Multi user Operating System c Unix is a Multi tasking Operating System 1 Long running programs can be sent to batch which will run a process in background mode freeing up the terminal for coincidence interactive use e Users have the ability to write scripts to control the run ning of complicated processes f Unix supports the X Window protocol X Windows was ini tially created by computer scientists from MIT to allow a windowing environment for multi user operating systems such as Unix and VMS Getting Started in Unix a Unix is designed to be user interactive that is the user can talk to the operating system directly through ter minal windows Much of the work you will do on the Linux workstation will involve use of such a terminal Window b Users and passwords The following commands are useful for checking on current users and changing your password PHYS 40075007 Computational Physics r00t dev etc usr users sbin man lib bin lutter mwc jmo manl man3 atlas bin idl iraf tex man2 man4 atlas hst mcmath Figure l 1 Directory Tree Structure of a Unix File System passwd who whoami sets users password allows you to see who else is logged on showing the users on the system their terminal ID numbers and when they logged in to the system tells you your own username and when you logged in c Getting Help man Displays the online manual page for a command 1 Other useful utilities are available through pull down menus on the X Window toolbar 69 Email web browser etc 4 Unix File System a The Unix le structure is a hierarchy with the root direc tory at the top see Figure 1 1 b When you log in to a Unix system you are assigned to a login home directory For example user5lutter Donald G Luttermoser E TSU C d The directory where you are working is your current direc tory The pathname of the le gives you the exact location of the le The full path is the absolute pathname relative to the root directory For example usersluttertdlmcmath is the directory where I might keep my observed spectra taken with the McMath Telescope The relative path is the pathname relative to your cur rent directory For example From the login directory userslutter the relative path is tdlmcmath 5 File Basics and Filenames 3 Conventions i Filenames can be any length and almost any char acter is valid Avoid the following N amp l quot ltgt l ii Unix is case sensitive iii There is no division into lename and le ex tension although there are certain naming conven tions For example C f f90 pro tex C language source code Fortran 77 source code Fortran 90 source code lDL source code TEX or ETEX source code Wildcards can be used 7 represents any one character represents any number of characters including none iv l 26 PHYS 40075007 Computational Physics b Directory Navigation pwd prints the working directory to the screen so that you know where you are Cd changes directory takes you back to your home directory Cd changes to the root directory Cd N changes to your login directory c Listing Files Is lists the unhidden les in the current directory 1 Command switches or options are placed after the com mand and before any arguments Is a lists les including hidden les note that a hidden le begins with a Is I lists the contents of a directory in long form 6 File Management Utilities a Copying Files Cp makes a copy of a le The syntax of the command is Cp source le target le mv moves or renames a le The syntax of the command is mv source le target le i Cp makes a new copy of a le leaVing the existing source le unchanged ii The new copy will overwrite a le of the same name in the target directory iii If there is more than one source le then the target le must be a directory iv mv creates a new copy of a le without retaining the source le Donald G Luttermoser E TSU v mv is used to move les from one directory to another vi mv is used to rename les b Removing Files rm deletes unwanted les The syntax of the command is rm lename rm od deletes a group of les ending in old rm i od Unix asks to con rm a deletion before deleting les ending in old c Creating and Removing Directories mkdir creates new subdirectories The syntax of the command is mkdir newdirectory rmdir deletes directories The syntax of the command is rmdir directory i mkdir can create a subdirectory of the current work ing directory ii mkdir can create directories within directories other than the current For example mkdir usersdata newdirectory iii To remove a directory 1 you must be the owner 2 the directory must be empty 3 it must not be the current directory 7 PHYS 40075007 Computational Physics Other Unix Commands and Tools a Displaying Contents of Files cat concatenates and displays the contents of a le The syntax of the command is cat lename If the le is more than 24 lines the display will scroll off the screen displays output from le one screenful at a time The syntax of command is more lename more b Pipes and Redirection The shell normally expects to re ceive input commands from stdin the keyboard Out put is normally sent to stdout the screen and any error messages to stderr the screen pipe allows standard output from one command to be used as standard input for another command For example cat lename more The output of the cat command is piped through the more command gt redirects stdout to a real le For example Is I gt lelist lt redirects stdin to a real le For example new le lt lelist Redirection always sends output to a le whereas pipe sends output to another command There is no limit to the number of pipelines that can be set up Donald G Luttermoser E TSU c Searching for Strings on searches for text in either a single le or group grep of les The syntax of the command is grep search string lename grep i ignores the case of the search string grep lists only the names of les containing the search string WC counts words and lines in les The syntax of the command is WC options lenames WC C displays only the number of characters WC displays only the number of lines WC W displays only the number of words i grep searches one line at a time ii grep looks for strings of text and does not limit itself to whole words iii grep can be used with wildcards File Permissions i The default setting for a new les is r W which means that only the user can read and write a le ii The permissions column consists of ten charac ters The rst character denotes the type of le For example ordinary le d directory link iii The next three show access mode for the user owner of the le e f PHYS 40075007 Computational Physics iv The next three show the access for the group v The last three show access for all others vi The protection symbols mean r read allows user to read the contents of a le allows user to modify a le allows user to execute or run a program le w write X execute User Access The different classes are de ned as u user owner of le g group le belongs to 0 other users a all users Changing Access Privileges Chmod changes the access mode of les you own to restrict or allow access The syntax of the command is Chmod classes 0pcmti0ns lenames operation is one of to add take away or set permissions For example Chmod ugw examplel F Summary of Unix Commands There are a large number of commands in the Unix operating system Table 1 2 lists some of the more common ones Note that with some of these commands l have created new commands using the alias command that make a bit more user friendly output of the original Unix commands These new com mands are listed in Table 1 2 in column labeled DGL Com mand Note that the alias utility is only found in the csh tcsh bash and Korne shells Initially I have set all of you up on Donald G Luttermoser ETSU l 31 faster with either the tcsh shell or bash shell 3 Note that commands that are listed with Unix ags re letters after the initial command word require a minus sign Note 4 77 that some avors of Unix do not require the quali er a Most processes are run in foreground this means that one does not get the Unix prompt back until the process is complete i For typical Unix commands and processes this is the way to go since these commands and processes take a fraction of a second to run ii However for long running processes editing a le for instance it is often wise to put a process in background mode To do this just include an am persand sign at the end of a command 69 emacs myfiletxt amp iii The job will then run in background mode and the Unix prompt will immediately come back to the screen awaiting further input b You can check on background jobs by entering the jobs command i To bring a background job back to foreground mode either enter fg to bring the last entered back ground job to foreground or enter fg PID where PID is the process ID number which can be ob tained with the ps command described above ii If you ever need to stop a job while it is in back ground enter kill 9 PID PHYS 40075007 Computational Physics Table l 1 Unix Commands and My Aliases chmod z39jk le cp le le df emacs le grep str le head le logout lpr le ls ls fl ls 7a man command mkdir dlmame more le mv le le passwd printenv var ps ps iaf pwd rm le rmdir dlmame source le tail le who Change the read write execute protection of a le Copy le le to le le Give information about the mounted disks Edit le le with the emacs editor the shortcut command puts the emacs editor into a pretty screen with big fonts Search files for pattern Displays beginning of le default is ten lines Log off of the terminal session Print the le le on the laser printer Give listing of a directory Give a long listing of a directory Give listing of a directory including hidden les Display information man pages of the Unix command command make directory dlmame Type the contents of le le to the screen one page at a time Move or rename le le to le Change your password to a new password Display the value of the environment variable var Display current process status just gives information on the user s processes Print working directory Delete le le remove directory dz39mame7 if empty Run a Unix shell script called le Displays end of le default is ten lines Displays names and other information about users on system Unix DGL Command Description Command cat le concatenate and display le contents 7 cd dlmame Change to directory dlmame in dlmame cd Go back to the parent directory of the current up directory cd N Go back to home 239e7 login directory out see chmhelp sd t le rename le le where del le run le Donald G Luttermoser E TSU l 33 Tablel 1 continued Unix DGL Command Description Command 7 Print a help le on the chmod command chmhelp Convert a dvi le le to a postscript le one dvips le that can be printed on the laser printer 7 Preview a 111239 le on an X Windows terminal xdvi le 7 Type a beginner s reference guide to Unix help 7 Start lDL idl 7 Run the LATEX le le through the LATEX word latex le processor 7 Delete old duplicate versions of les pur C v Your can also suspend a foreground job by entering ltCtrgt Z le pressing down on the control Ctrl key when you hit the 2 key note that this doesn t work from inside an emacs process ltCtrgt Z will save and exit the le you are editing 4 Finally we have only scratched the surface of Unix in this section of the notes This however should be enough to enable you to work on the Linux workstation G The Emacs Editor 1 Emacs is a user friendly editor that exists on most Unix worksta tions a Although emacs is not shipped with Unix itself it is freely available to the Internet community b It was written and still evolving by the GNU Project a group of computer scientists who don t like the large amounts of money that vendors are charging for software they write software and give it out for free l 34 PHYS 40075007 Computational Physics 2 Unix s standard editor is vi and it is very user unfriendly in stead use emacs 3 To edit a le in the emacs editor move to the subdirectory con taining the le and enter emacs lename where lename is the name of the le to be edited 4 The Linux workstation has the most recent version of emacs on it which brings up a nice GUI widget The buttons and pull down menus are obvious in their operation as such we won t describe the keyboard commands that one could also use inside of emacs H Additional Information about the Linum Workstation 1 Most of you are not used to talking to the operating system by issuing commands at the system prompt in a terminal window though you should practice getting used to it a Many Unix operating systems have a front end window ing system based on the X Window protocol called CDE Common Desktop Environment b Most Linux operating systems 69 Red Hat use the gnome front end which is also based on X Windows c As such you can following the same protocol that you do on a Window s machine by double clicking on icons and through the use of pop up menus from the toolbar located at the bottom of the main console screen 2 One nal comment about using the Linux workstation node fes ter a Since you can only log into this workstation from the con sole in Brown Hall 260 you will have to develop a protocol Donald G Luttermoser E TSU b C d 3 Also note that once you are supplied with a password to this machine change it to something that is relatively cryptic from the console and remember your new password If you forget for using the computer Limit your time for each session to 1 hour per student per day Feel free to go past this limit if there is no one waiting to use the machine Let me know if you are haVing trouble getting time on the machine If this happens I set up a users schedule to allocate time Another suggestion to avoid congestion on the worksta tion is to write the rst draft of your codes on a PC then mail it to your Z accounts You can then log into the mail server machine from the web browser software on fester and clip and past the text displayed on the web browser to an emacs window you can start from the Linux prompt make sure the lename started with emacs matches what you are pasting from the web browser your password see me and I will issue you a new one 4 Since the printers in Brown Hall 260 are not behaVing well I have set the default printer to be physlabprt3etsuedu in the lecture room e Brown Hall 264 PHYS40075007 Computational Physics Course Lecture Notes Section IX Dr Donald G Luttermoser East Tennessee State University Version 41 Abstract These Class notes are designed for use of the instructor and students of the course PHYS 40075007 Computational Physics taught by Dr Donald Lutternioser at East Tennessee State University IX Numerical Solution of Ordinary Differential Equations ODE A Introduction 1 Many important problems in engineering and the physical sci ences require the determination of a function satisfying an equa tion containing derivatives of the unknown function Perhaps the most familiar example is Newton s Second Law of Motion d2rt drt m it dt IX 1 trt a The position 7 05 of a particle of mass m is acted upon by a force F which may be a function of time t position 705 and the velocity d7 tdt b To determine the motion of the particle acted upon by a given force F it is necessary to nd a function 7 satisfying Eq IX 1 c It also is important to set up a coordinate system rst with respect to the motion of the mass before setting up the equations to be solved d If the force is due to gravity then d2rt dtg mg IX 2 Integrating Eq IX 2 we have drt t lX 3 dt 91 61 Mt gt2 0115 Cg IX 4 where 01 and Cg are constants of integration PH Y5 4007 5007 Computational Physics e To determine rt completely it is necessary to specify two additional conditions such as the position and velocity of the particle of some instant of time gt these initial conditions then give the value of 01 and Cg These are often referred to as time dependent problems f There are some problems where it is more convenient to give conditions at the boundaries of the integration path gt boundary conditions These are often referred to as time independent problems g In order to completely solve a differential equation one needs either initial or boundary conditions If a differential equation DE depends on a single independent variable and only an ordinary derivative appears in the DE then such a DE is called an ordinary differential equation ODE a For the following de nitions let s assume that the variable 06 represents the independent variable similar to what is often done in a mathematics course and y Mac is the dependent variable y is given by some function u of the independent variable b The order of an ODE is the order of the highest derivative that appears in the equation Eqs IX l and IX 2 are second order differential equations and Eq IX 3 is a rst order DE c It is convenient to follow the usual notation in the theory of DEs to represent duxdx d2uxdx2 d uxdx with the notation 7 y yn Thus Eq IX 1 is written as F ltxyyl y gt 0 IX 5 Donald G Luttermoser ETSU IX 3 d Occasionally other letters will be used instead of x for the dependent variable and y for the dependent variable ie function the meaning will be clear from the context e We shall assume that it is always possible to solve a given ODE for the highest derivative obtaining if f 0621213 y y 1 IX6 f A second important classi cation of ODEs is according to whether they are linear or nonlinear The DE in Eq IX 5 is said to be linear if F is a linear function of the variables yy y Thus the general linear ODE of order n is ao06y a106y 1 an06y We IX7 g An equation which is Lot of the form in Eq IV 7 is a nonlinear equation For example the equation for the motion of a pendulum is nonlinear 12639 g W Z s1n639 0 IX 8 3 If a DE depends upon several independent variables it is called a partial differential equation PDE gt the DE contains par tial derivatives eg 8815 instead of ordinary derivatives eg ddt The wave equation is a good example of a PDE 82u06 t 821406 If 2 7 7 a 8962 Bit 39 IX 9 B Numerical Methods 1 Terminology a Shooting or marching methods The solution is cal culated step by step by starting at one boundary and in tegrating toward the other The step size is the change b C d e f g PH Y5 4007 5007 Computational Physics 6 9 A06 in independent variable used in a shooting or marching scheme Iterative method A repetitive process by which suc cessively more accurate approximations to the solution are obtained An iteration is one cycle of the repetitive process Difference equation An approximation to a DE where a derivative is replaced by a quotient Relaxation method A solution is calculated every where at once by solving a set of difference equations in an iterative fashion Computational mesh or grid The independent vari able is represented by a set of discrete values 6 9 a set layers of given thickness in a stellar atmosphere called grid points zones or cells The A06 or A7 AMi etc is called the grid spacing or mesh size or inter val Ax E gal1 06 A model then becomes the set of physical properties spec i ed at all the grid points 69 PlTiFZpZ at zones 06 IX ll An evolution is a sequence of models at different times tn Each model is a generation or cycle Each successive advance forward in time is called the time step and At E th tn IX 12 is called the time step size Donald G Luttermoser ETSU h k Truncation error TE is the per step for shooting schemes or per mesh for relaxation schemes error in the calculation of the dependent variables Cumulative truncation error CE is the total such error across the grid Roundoff error is error introduced by the nite number of digits carried by the computer The higher precision you use the smaller the roundoff error The order of a numerical scheme is the power of the mesh size or step size in the highest order terms which are ac curately represented by a numerical scheme Some books use the TE some the CE to de ne this and so there is often some confusion Explicit schemes One where values at the next step are obtained by a direct algebraic computation involving only values from the previous step Implicit schemes One where new values must be solved for interatively 2 Expansions a Often there are times especially at boundaries where singularities appear in an equation e divide by zero During these times one must use expansions to avoid these formal singularities For instance in the hydrostatic equilibrium equation dP 7 GMT d T p as r gt 0 M gt 0 so the RHS gt 00 2 IX13 m a PH Y5 4007 5007 Computational Physics One must use the boundary conditions to develop Taylor expansions away from the singularities i The Taylor expansion takes the form 1 f06h f06hfl06 h2f 06 m IX14 where the symbol means higher order terms that are usually dropped from the derivation ii An alternative equivalent form of the Taylor se ries used in numerical analysis is fwm mme M axm where C is a value between 06 and 06 h iii We have not dropped any terms in Eq IX 15 this expansion has a nite number of terms Tay lor s theorem guarantees that there exists some value C for which Eq IX 15 is true but it doesn t tell us what that value is With Taylor expansions in mind we can transform the derivative formula from WWW x fx f DC 16 to we W h o wherex C xh i This equation is known as the right derivative IX17 formula ii The last term on the right is the truncation er ror gt it is the error introduced by the truncation of the Taylor series Donald G Luttermoser ETSU IX 7 iii Since we don t know the value of C a priori the f term truncate is usually dropped and we say the error was we make by neglecting this term is the truncation error and as such Eq IX 17 is often written as mf2 w where the truncation error term is now just speci ed by its order in h 0h IX 18 iv Note that the truncation error is linear in h the smaller we make h the smaller our TE however the more calculations we have to make which re sults in larger roundoff errors v Keep in mind that the truncation error depends on the approximation used in the algorithm whereas the roundo error depends on the hardware 1 Note that we also can reduce the size of the TE by intro ducing a centered formula for the derivative equation fxh fx h 7 f i fling 2h IX 19 i This formula is said to be centered in 06 ii While this formula looks very similar to Eq IX 17 there is a big difference when h is nite since a Taylor expansion of this form of the derivative takes on the following form f06hf06 h for 2h Wf IX20 IX 8 PHYS 40075007 Computational Physics where f 3 is the 3rd derivative of f and 06 h S C S 56 h iii The TE is now quadratic in h a big improvement over Eq IX 17 iv From this formalism it can be shown that the Taylor expansion of the second derivative is 2 1 fx l2f4c7 h 12 IX 21 where 06 h S C S xh Again the TE is quadratic in h e Selecting values of h What do you pick for h i First de ne the absolute error as 5 true value computed value IX 22 ii Neglecting roundoff error to make the TE term in Eq IX 20 small with respect to the other term in this equation then choose 65 h lt IX23 iii Generally we don t know f3 but often we can set a bound For example if sinx then f3C S 1 so if we want an absolute error of 5 m 10 6 a typical value then we should take h m 2 x 10 iv If we cannot estimate an upper bound then arbi trarily pick a value of h use it try a smaller value Donald G Luttermoser ETSU 3 of h compare the two answers and if they are close enough assume everything is ne and your answer is converged If not keep on choosing smaller val ues of h e iterate until the above is true This is Lot universally true however Shooting Methods 3 b Assume we have a DE given by dydx f x y as shown in Figure IX l Now given 06yj for allj g 239 obtain xi1yi1 This is the standard technique in the shoot ing or marching method The are a variety of ways to carry out such a calculation Euler s Method Assume we wish to follow the motion of a mass m using Newton s 2nd Law of Motion The equation of motion that we want to solve numerically is IX 24 where t IX 25 and d is the acceleration Euler s method uses the right derivative formula see Eq lX l8 where we replace the grid step h with the time step 739 The equations of motion become 07 5W m T IX26 w 07 170 IX27 170 r 270 Tth 1705 072 IX 28 Ht T Ht 71705 072 IX 29 IX 10 PH Y5 4007 5007 Computational Physics Y117 Yi Xi Xi1 X Figure IX 1 Data points given by the DE dydx Notice that 707 C72 i We introduce the notation fn E fn Dr n 12 IX 30 so f1 ft 0 Our equations for the Euler method dropping the error term now take the form 17M Un75n IX31 m1 Elma IX32 ii The calculation of trajectory would proceed as follows Specify the initial conditions F2 and 172 Choose a time step 739 Calculate the acceleration given the current Fand U Use the Euler method to compute the new Fand U CPIPLUM Go to step 3 until enough trajectory points have been computed Donald G Luttermoser ETSU IX 11 iii The truncation error in such a scheme is 1 2d27 TE WW IX33 239 so that Euler s method is a 1st order accurate method iv The cumulative error is 1 N d2 CE m 7 27 N 2 Tl dt2 av tow 7 239 2 N dt2 tN to d2 T N WW 2 Exth t0d r W IX34 NleleH or CE N c EulerCromer and Midpoint Methods A modi ed version of Euler s method is UnH Um Tdn IX 35 Fn1 F7 TUn1 i The updated velocity is used in the second equa tion This formula is called the Euler Cromer method ii Notice that the TE is still of order 72 iii We can modify this further to come up with the midpoint method UnH Um Tdn IX 37 FM Fn 7 1 T n IX 38 IX 12 d 8 PH Y5 4007 5007 Computational Physics iv Plugging the velocity equation into the position equation we see that 1 Fnl Fm 717 id g IX 39 v The TE is still of order 72 in the velocity equa Unfortu nately this midpoint method gives good results for tion but for position the TE is now 73 relatively few physical systems projectile motion is one of them SecondOrder Runge Kutta Method Runge Kutta RK schemes establish higher derivatives by calculating intermediate or provisional values of y in the interval 062306241 i Second order RK involves 2 substeps per step 1 9amp1 9239 Axi xz a 1112 1X40 A 239 2 yi1 9239 706 lfltxiayi f56i1 yf1l a IX 41 where the superscript p represents provisional ii Since RK schemes are widely used in computa tional physics especially 4th order RK we will de vote an entire subsection to it see below PredictorCorrector Methods These use 06yj for j g 239 to establish higher derivatives dnydx We will not cover these methods in detail in this course Relaxation Henyey Methods In these methods one solves all equations at all grid points at once In our example here let s assume we have a spherical distribution of gas Donald G Luttermoser ETSU a Grid First one needs to set up a grid of independent variable values from one boundary to the other Let s assume that 7 2 corresponds to this grid then the grid is represented as 73239 01N or N 1 grid points b Model The model is de ned as the set of dependent variables at each grid point 69 PlTiMnLn at 7 2 for 4N 4 unknowns c Then the Differential Equations take on the form 17 f1PTMTLTr IX42 17 f2PTMTLTr IX43 1 f3PTMTLTr IX44 dd f4PTMTLTr IX45 d Differencing Schemes i Forward differencing Pi1 Pi Ari f1 E7 39 39 ii Centered differencing M i f M AH i 1 2 l 2 7 IX 47 iii Backward differencing B B L f1 3241 EH 39 39 DC 48 An IX 13 IX 14 e 0 PH Y5 4007 5007 Computational Physics Boundary Conditions i Center to boundary 1 01P07T07M7 07L7 0 02PoaT0MT0LT0 IX 49 ii Surface 16 boundary 2 81PN7TNJMTNJL7 N 0 82PN7TN7MTN7LTN 0 IX50 Difference Equations i Choose one of the differencing schemes and apply to Eqs IX 42 through IX 45 ii Plug in known 7 2 Values which results in 4N dif ference equations which be hnear algebraic equations for 4N 4 unknowns iii These equations can be written formally as 91 sz T13 Mm Lm Pi17 T1417 Mri17 LM1gt 7 0 IX51 92 sz T13 Mm Lm Pi17 T1417 Mri17 LM1gt 7 0 IX52 93 B7 T13 Mm Lm P1417 T1417 Mri17 Lri1 7 0 IX53 94 sz T13 Mm Lm Pi17 T1417 Mri17 LM1gt 0 IX54 iv We have such equations for 239 0 l N l IX 51 to IX 54 plus IX 49 and IX 50 4N 4 equations gt in 4N 4 unknowns v Great The only problem is that these equations are not linear They can t be solved directly A gen eralized Newton Raphson iterative scheme is used to overcome this di iculty Donald G Luttermoser ETSU g Method Solve linearized equations for corrections 7 gt IX 15 Guess to solution lt 7 Are corrections small enough yes Solution Add correction to the old guess T One iterates until the guesses converge to some prescribed accuracy determined by the size of the corrections PU TU i Consider one of the difference equations In gen jth iteration 39th j correction eral 2 138 gt7 0 ii Say we know PZj Tim and want to calculate the jth correction To generate linear equations we can solve for the jth correction plug PZj SBj Tim 6110 into gm Taylor expand keep terms up to rst order only in the jth corrections and then set it equal to zero IX 16 PH Y5 4007 5007 Computational Physics gmPWgt6PWNWP l6PXLN 0 gnur utf tuy j 3 j DQQQ algal 63 BB 62 an 1X55 39 3998132 6P l4aei iii We have now linearized the set of difference equa tions The partial derivatives can be evaluated be cause they depend only on the known quantities 1351 Tim iv By casting all of Eqs IX 5 1 to IX 54 and IX 49 amp IX 50 we get 4N 4 linear for 4N 4 unknown equations corrections v Formally 71 6 E x T aXamp matrix of correction error partials vector vector Donald G Luttermoser ETSU IX 17 m 01 P gt 6P0 I 39 J 02 P0 39 139 6M5 91 P0 J39 j 5L 92 Po 6P 139 1039 93 P0 6T1 g4 Pom 39 6X 3 E j j 91 P1 SALiv1 7 6L 1 j 6PM 931 gt 6T 131911 gt 139 SAXV 31 Pg 6L 139 7 TN 32 PN 739 7 0 0 0 370300 0 0 0 0 0 0 0 0 X X X 0 0 0 0 0 0 0 0 o 3sz 31 o o o o 0 Sign X X X X X X X 0 0 0 0 0 X X X X X X X X 0 0 0 0 0 H X X X X X X X X 0 0 0 0 0 A a a a a a a a a 0 0 0 0 31 2 3 2 22 2 352 2 371 2 2 3152 2 352 2 0 0 0 0 0 X X X X X X X X 0 0 0 0 0 X X X X X X X X 0 0 0 0 0 X X X X X X X X 0 0 0 0 0 0 0 0 0 3793123 X X X X 0 0 0 0 0 0 0 0 X X X X X IX 18 PH Y5 4007 5007 Computational Physics vi The 2 matrix has a banded structure and is mostly empty vii To solve for the jth correction this matrix must be inverted Numerically the number of operations involved is Full 2 4N4 gtlt 4N4 N WV Actual N 824N 4 Banded Matrix h Gaussian Elimination Many methods now exist for inverting banded matrices In the last section of the notes we were introduced to the Gaussian elimination scheme for solving sets of linear equations 1 here describe this technique for ODEs Donald G Luttermoser ETSU i Start at one corner and solve one block of equa tions for some variables in terms of the others 69 for upper left 2 x 4 block 6P5 apo Mg b1306ng CPO IX57 6T0 aTo Mg bT06L ggt 6T0 ii Store the a bc s iii Substitute Eq IX 57 in the next block and re peat and above 69 aM39r O M Sl39 l l CM39I O 6L3 atom mm cm IX58 6P1 apl Mfi39 bpl Li Cpl 6T1 aTl Mfi39 le Ll l CT become the new set of Eqs IX 57 to substitute into the next block iv Solving the last set of equations at the other end of the matrix gives 6MT SM and 512334 6L9 v Now go backwards and back substitute into the Eqs IX 58 using the stored a b c s to nd all the jth corrections IX 19 IX 20 PH Y5 4007 5007 Computational Physics i Convergence Criterion The iterative process is said to converge when the corrections are as small as desired 69 all lt C const lt5Pil B c N 102 to 106 Most codes at least 1 D converge after N3 to 5 iterations C FourthOrder Runge Kutta Method 1 When numerically solving ODE one typically rewrites 2nd and higher order DEs into a set of lst order DEs a A 2nd order equation such as 1 a IX 59 b This can be rewritten as 2 equations dx 7 lX 60 v d do 7 lX 61 a d 2 As shown above Euler s method is a rst order method which is graphically represented in Figure IX 2 3 The Runge Kutta method is essentially a modi ed Euler s method a Use the derivative at one step to extrapolate the midpoint value use the midpoint derivative to extrapolate the function at the next step see Figure lX B b Evaluates the derivative function twice at each step 739 Cu mulative error is of order C72 a second order method Donald G Luttermoser ETSU IX 21 1 Poor Guessf T Figure IX 2 Euler s method is a rstorder method which is not necessarily accurate XD initial value A A m X Better Guess A2 m i h X gm X I A39ctual X XR1 X hquot A2 X A ftgXa 0 Figure IX 3 The RungeKutta method is a secondiorder method which is more accurate than Euler s method IX 22 6 7 PH Y5 4007 5007 Computational Physics Runge Kutta methods achieve better results than Euler by using intermediate computations at intermediate grid steps The fourthorder rule is the favorite method as it achieves good accuracy with modest computational complexity 3 b C d e Use derivative of rst step to get trial midpoint Use its derivative at rst step to get second trial midpoint Use its derivative to get a trail end point lntegrate by Simpson s Rule using average of two mid point estimates The cumulative error is of fourth order and truncation error is C75 The fourth order RK scheme mathematically is where 1 a a A 1 m T m 67111 2F2 2F3 F4 IX62 F1 fli 15 IX63 a a 1 F2 f ETFl t ET a a 1 1 F3 f ETFQ t ET 134 fr TF3 t 739 IX66 Compared with Euler 4th order RK has 4 times more calcula tions per step but uses the fourth root as many steps to achieve convergence see Figure lX 4 Donald C Lutterrnoser ETSU IX 23 Midpoint A Initial Value I I I 4 t0 1 2 Quadratic Simpson Rule Extrapolation Figure IX 4 Fourtheorder RK With the E s in Eqs D963 to D966 being represented With Al s in this gure Also fis given by X in this gure 8 You may wonder Why 4th70rder and not 8th or 23rd70rder BungeiKutta Well the higher order methods have better trune cation error but also require more computation hence more roundoff errorl The optimum for RK schemes is 4th orderl 9 Sometimes accuracy can be improved through use of an adaptive step size Adaptive methods are fairly easy to incorporate and I refer you to Numerical Recipes for a description on how to incorporate theml 10 The following program is taken from Numerical Recipes and is a Fortran 77 version of the 4theorder RK scheme IX 24 PHYS 40075007 Computational Physics 96 ltltltltltltltltltltltltltltltltltltltltltltltltltltltltlt RK4 RK4 RK4 gtgtgtgtgtgtgtgtgtgtgtgtgtgtgtgtgtgtgtgtgtgtgtgtgtgtgt SUBROUTINE RK4Y DYDX N X H YOUT DERIVS Given values for N variables Y and their derivatives DYDX known at X use the 4thorder RungeKutta method to advance the solution over an interval H and return the incremented variables as YOUT which need not be a distinct array from Y The user supplies the subroutine DERIVSX Y DYDX which returns the derivatives DYDX at X 96969696969696 IMPLICIT REAL8 AHOZ EXTERNAL DERIVS PARAMETER NMAX 10 DIMENSION YN DYDXN YOUTN YTNMAX DYTNMAX DYMNMAX HH H 05DO H6 H 6DO XH X HH gt1 First step gt1 DO 11 I 1 YTI YI 11 CONTINUE gt1 N HHDYDXI Second step gt1 CALL DERIVSXH YT DYT DO 12 I 1 N YTI YI HHDYTI 12 CONTINUE gt1 Third step gt1 CALL DERIVSXH YT DYM DO 13 I 1 N YTI YI HDYMI DYMI DYTI DYMI 13 CONTINUE gt1 Fourth step gt1 CALL DERIVSXH YT DYT DO 14 I 1 N YOUTI YI H6DYDXIDYTI2DODYMI 14 CONTINUE Donald G Luttermoser ETSU IX 25 RETURN END D The Adams Method The ShampineGordon Routine 1 In 1975 LF Shampine and MK Gordon wrote a book enti tled Computer Solution of Ordinary Di erential Equations The Initial Value Problem 3 b This textbook described numerical techniques in solving non sti initial value problems in ordinary differential equations The ODE code itself is comprised of a few subroutines and a driver program as shown in the example program below The subroutines are i DE lntegrates a system of up to 20 rst order ODEs of the form DYDT FT Y1 Y2 YNEQN Y given at T This subroutine integrates from T to TOUT On re turn the parameters in the call list are initialized for continuing the integration The user has only to de ne a new value of TOUT and call DE again ii STEP lntegrates a system of rst order ODEs over one step normally from T to T H using a modi ed divided difference form of the Adams Pece formulas Local extrapolation is used to im prove absolute stability and accuracy The code ad justs its order at step size to control the local error per unit step in a generalized sense Special devices are included to control roundoff error and to detect when the user is requesting too much accuracy see program odef below for further details IX 26 anaaaaaaaaaaaan PH Y5 4007 5007 Computational Physics INTRP The methods used in subroutine STEP approximate the solution near 06 by a polynomial Subroutine INTRP approximates the solution at mam by evaluating the polynomial there Information de ning this polynomial is passed from STEP so INTRP cannot be used alone see program odef be low iii c Note that the user of these routines also has to supply a subroutine called F which sets up the derivative equations YP to be solved by DE An example of such a subroutine is provided following the listing of the odef le below Note that these subroutines are supplied to you on the Project 2 Web Page for this course if you wish to use them If you use any of these routines please note that when ever you see a mark in the far right of the program these lines have to be modi ed to re ect the machine pre cision em for single and double precision of the machine on which you are carrying out these calculations You can determine this precision with machinf program supplied to you on the Project 2 Web Page ltltltltltlt File odef gtgtgtgtgtgt subroutine odefneqnyttoutrelerrabserriflagworkiwork implicit real8ahoz double precision subroutine ode integrates a system of neqn first order ordinary differential equations of the form dyidt fty1y2yneqn yi given at t the subroutine integrates from t to tout on return the parameters in the call list are set for continuing the integration the user has only to define a new value tout and call ode again the differential equations are actually solved by a suite of codes de step and intrp allocates virtual storage in the arrays work and iwork and calls de de is a supervisor which directs the solution it calls on the routines intrp to advance the integration and to interpolate at output points step uses a modified divided difference form of the adams pece ode step and Donald G Luttermoser ETSU aanaanaananaanaanaananaaaananananaaaaaaaaaaaana I 27 formulas and local extrapolation it adjusts the order and step size to control the local error per unit step in a generalized sense normally each call to in the direction of tout integrates beyond tout step advances the solution one step for reasons of efficiency de internally though never beyond t10toutt and calls intrp to interpolate the solution at tout an option is provided to stop the integration at tout but it should be used only if it is impossible to continue the integration beyond tout this code is completely explained and documented in the text computer solution of ordinary differential equations the initial value problem by l f shampine and m k gordon the parameters represent f double precision subroutine ftyyp to evaluate derivatives ypidyidt neqn number of equations to be integrated integer4 y solution vector at t real8 t independent variable real8 tout point at which solution is desired real8 relerrabserr relative and absolute error tolerances for local error test real8 at each step the code requires dabslocal error le dabsyrelerr abserr for each component of the local error and solution vectors iflag indicates status of integration integer4 work real8 arrays to hold information internal to iwork integer4 which is necessary for subsequent calls first call to ode the user must provide storage in his calling program for the arrays in the call list yneqn work10021neqn iwork5 declare f in an external statement supply the double precision subroutine ftyyp to evaluate dyidt ypi fty1y2 and initialize the parameters yneqn neqn number of equations to be integrated y t starting point of integration tout point at which solution is desired relerrabserr relative and absolute local error tolerances iflag 11 indicator to initialize the code normal input is 1 the user should set iflag1 only if it is impossible to continue the integration beyond tout all parameters except f vector of initial conditions neqn and tout may be altered by the IX 28 aanaanaanaananaaaaaaaaaaaanananaaaaaaaa subsequent calls to subroutine ode the integration only define a new tout reach tout and the user wants to continue he just calls again the output value of subsequent calls is to stop the integration internally at the new tout ie change output be changed by the user before continuing remain unchanged PH Y5 4007 5007 Computational Physics code on output so must be variables in the calling program output from ode neqn unchanged y solution at t t last point reached in integration normal return has t tout tout unchanged relerrabserr normal return has tolerances unchanged iflag3 signals tolerances increased iflag 2 normal return integration reached tout 3 integration did not reach tout because error tolerances too small relerr abserr increased appropriately for continuing 4 integration did not reach tout because more than 500 steps needed 5 integration did not reach tout because equations appear to be stiff 6 invalid input parameters fatal error the value of iflag is returned negative when the input value is negative and the integration does not reach tout ie 3 4 5 workiwork information generally of no interest to the user but necessary for subsequent calls ode returns with all information needed to continue if the integration reached tout the user need and call again if the integration did not iflag is the appropriate input value for the only situation in which it should be altered iflag2 iflag2 to input error tolerances may all other parameters must a 3 subroutines de and step be sure they are set before using ode contain machine dependent constants logical startphase1nornd dimension yneqnwork1iwork5 external f Donald G Luttermoser ETSU 000000000 I 29 data ialphaibetaisigiviwigiphaseipsiixihiholdistart 1 itoldidelsn113253850627576888990919293 iyy 100 iwt iyy neqn ip iwt neqn iyp ip neqn iypout iyp neqn iphi iypout neqn ifiabsiflag eq 1 go to 1 start workistart gt 00d0 phase1 workiphase gt 00d0 nornd iwork2 ne 1 1 call defneqnyttoutrelerrabserriflagworkiyy 1 workiwtworkipworkiypworkiypoutworkiphi 2 workialphaworkibetaworkisigworkivworkiwworkig 3 phase1workipsiworkixworkihworkiholdstart 4 workitoldworkidelsniwork1norndiwork3iwork4 5 iwork5 workistart 10d0 ifstart workistart 10d0 workiphase 10d0 ifphase1 workiphase 10d0 iwork2 1 ifnornd iwork2 1 return end subroutine defneqnyttoutrelerrabserriflag 1 yywtpypypoutphialphabetasigvwgphase1psixhhold 2 starttolddelsgnnsnorndkkoldisnold implicit real8ahoz ode merely allocates storage for de inconvenience of a long call list described in the comments for to relieve the user of the consequently de is used as ode this code is completely explained and documented in the text computer solution of ordinary differential equations the initial value problem by l f shampine and m k gordon logical stiffcrashstartphase1nornd dimension yneqnyyneqnwtneqnphineqn16pneqnypneqn 1 ypoutneqnpsi12alpha12beta12sig13v12w12g13 external f I 30 PH Y5 4007 5007 Computational Physics anon 0000 00000 the only machine dependent constant is based on the machine unit roundoff error u which is the smallest positive number such that 10u gt 10 u must be calculated and fouru40u in the following data statement before using de machin calculates u inserted in subroutine data fouru888d15 inserted the routine fouru and twou20u must also be step before calling de the constant maxnum is the maximum number of steps allowed in one call to de the user may change this limit by altering the following statement data maxnumSOO test for improper parameters fouru 40 d1mach4 ifneqn lt 1 go to 10 ift eq tout go to 10 ifrelerr lt 00d0 or eps dmax1relerrabserr ifeps le 00d0 go to 10 ififlag eq 0 go to 10 isn isign1iflag iflag iabsiflag ififlag eq 1 go to 20 ift ne told go to 10 ififlag ge 2 10 iflag 6 return abserr lt 00d0 go to 10 and iflag le 5 go to 20 on each call set interval of integration and counter for number of steps adjust input error tolerances to define weight vector for subroutine step 20 del tout t absdel dabsdel tend t 100d0del ifisn lt 0 tend tout nostep O kle4 O stiff false releps relerreps abseps abserreps ififlag eq 1 go to 30 Donald G Luttermoser ETSU a anon 0000 I 31 ifisnold lt 0 go to 30 ifdelsgndel gt 00d0 go to 50 on start and restart also set work variables x and yy store the direction of integration and initialize the step size 30 start true x t do 40 l 1neqn 4O yyl yl delsgn dsign10d0del h dsigndmax1dabstoutXfourudabsxtoutX if already past output point interpolate and return 50 ifdabsxt lt absdel go to 60 call intrpxyytoutyypoutneqnkoldphipsi iflag 2 t tout told t isnold isn return if cannot go past output point and sufficiently close extrapolate and return 60 ifisn gt 0 dabstout x ge fourudabsx go to 80 h tout x call fxyyyp do 70 l 1neqn y1 yyl hypl iflag 2 t tout told t isnold isn Or 70 return test for too many steps 80 ifnostep lt maxnum go to 100 iflag isn4 ifstiff iflag isn5 do 90 l 1neqn y1 yyl t x told t isnold 1 90 I 32 0000000 PH Y5 4007 5007 Computational Physics return limit step size set weight vector and take a step 100 h dsigndmin1dabshdabstendXh do 110 l 1neqn wtl relepsdabsyyl abseps call stepxyyf neqnhepswtstart 1 holdkkoldcrashphipyppsi 2 alphabetasigvwgphase1nsnornd 110 test for tolerances too small ifnotcrash go to 130 iflag isn3 relerr epsreleps abserr epsabseps do 120 l 1neqn 120 y1 yyl t x told t isnold 1 return augment counter on number of steps and test for stiffness 130 nostep nostep 1 kle4 kle4 1 ifkold gt 4 kle4 O ifkle4 ge 50 stiff go to 50 end true subroutine stepxyfneqnhepswtstart 1 holdkkoldcrashphipyppsi 1 holdkkoldcrashphipyppsi 2 alphabetasigvwgphase1nsnornd implicit real8ahoz double precision subroutine step integrates a system of first order ordinary differential equations one step normally from x to xh using a modified divided difference form of the adams pece formulas extrapolation is used to improve absolute stability and accuracy the code adjusts its order and step size to control the local error local
Are you sure you want to buy this material for
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'