Sci Prog Atmos Sci
Sci Prog Atmos Sci MSC 321
Popular in Course
Popular in Marine Science
This 53 page Class Notes was uploaded by Rowan Ward on Thursday September 17, 2015. The Class Notes belongs to MSC 321 at University of Miami taught by Staff in Fall. Since its upload, it has received 14 views. For similar materials see /class/205715/msc-321-university-of-miami in Marine Science at University of Miami.
Reviews for Sci Prog Atmos Sci
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: 09/17/15
Chapter 11 Finite Difference Approximation of Derivatives 1 11 Introduction The standard de nition of derivative in elementary calculus is the following 7 ux Ax 7 u r 7 111130 Ax 111 Computers however cannot deal with the limit of Ax 7 07 and hence a discrete analogue of the continuous case need to be adopted In a discretization step7 the set of points on which the function is de ned is nite7 and the function value is available on a discrete set of points Approximations to the derivative will have to come from this discrete table of the function Figure 111 shows the discrete set of points xi where the function is known We will use the notation ui to denote the value of the function at the 24th node of the computational grid The nodes divide the axis into a set of intervals of width Axl x111 7 xi When the grid spacing is xed7 ie all intervals are of equal size7 we will refer to the grid spacing as Ax There are de nite advantages to a constant grid spacing as we will see later 1 12 Finite Difference Approximation The de nition of the derivative in the continuum can be used to approximate the derivative in the discrete case N Ax 7 11111 7 ul 1 Ax Ax where now Ax is nite and small but not necessarily in nitesimally small7 ie This is known as a forward Euler approximation since it uses forward differencing 112 77 78 CHAPTER 1 1 FINI TE DIFFERENCE APPR OXIMATI ON OF DERIVATIVES Figure 111 Computational grid and example of backward forward and central approximation to the derivative at point x The dash dot line shows the centered parabolic interpolation while the dashed line show the backward blue forward red and centered magenta linear interpolation to the function lntuitively the approximation will improve ie the error will be smaller as Ax is made smaller The above is not the only approximation possible two equally valid approximations are backward Euler 11950 M 113 Ax Ax Centered Dz erenee Ax 7 7 Ax i qu 7 ui1 uml 2Ax 2Ax 114 All these de nitions are equivalent in the continuum but lead to different approx imations in the discrete case The question becomes which one is better and is 113 TAYLOR SERIES 79 there a way to quantify the error committed The answer lies in the application of Taylor series analysis We brie y describe Taylor series in the next section7 before applying them to investigate the approximation errors of nite difference formulae 1 13 Taylor series Starting with the identity uxux lv u15ds 115 Since is arbitrary7 the formula should hold with replaced by u 7 ie7 umwmww mm Replacing this expression in the original formula and carrying out the integration since is constant we get x 7 m m u s d5 d5 117 The process can be repeated with Mmwmfww ma 7 I 7 I I i iyn I m m m m i uzl x 011 1 2 u x1 u 5 d5 d5 d5 119 This process can be repeated under the assumption that is suf ciently differentiable7 and we nd 7 2 7 I n W um z e muzi WWW RM 39 39 1110 where the remainder is given by RM mmu 1 d5 1 1111 Equation 1110 is known as the Taylor series of the function about the point xi Notice that the series is a polynomial in x 7 the signed distance of z to m and the coef cients are the scaled derivatives of the function evaluated at xi 80 CHAPTER 1 1 FINI TE DIFFERENCE APPR OXIMATI ON OF DERIVATIVES If the n 1 th derivative of the function u has minimum m and maximum M over the interval then we can write fmltdsn1 R1Mltdsgt 1 1112 n1 n1 x 7 71 1 x 7 lt lt m n1 R39Lfli 1113 which shows that the remainder is bounded by the values of the derivative and the distance of the point z to the expansion point xi raised to the power 71 1 If we further assume that um is continuous then it must take all values between m and M that is x 7 man Rn1 Um 5 1114 for some g in the interval 1131 Taylor series and nite differences Taylor series have been widely used to study the behavior of numerical approxi mation to differential equations Let us investigate the forward Euler with Taylor series To do so7 we expand the function u at n1 about the point xi 3U Ax 3211 Ax 3311 An Am a mi 2 g m 3 g m 1115 The Taylor series can be rearranged to read as follows Am 7 311 7 Am 3211 Ax 3311 11 16 Am 3s m T 2 3x2 m 3 M3 m 39 Truncation Error where it is now clear that the forward Euler formula 112 corresponds to truncat ing the Taylor series after the second term The right hand side of equation 1116 is the error committed in terminating the series and is referred to as the trun cation error The tuncation error can be de ned as the difference between the partial derivative and its nite difference representation For suf ciently smooth functions7 ie ones that possess continuous higher order derivatives7 and suf ciently small Am the rst term in the series can be used to characterize the order of magnitude of the error The rst term in the truncation error is the product of the second derivative evaluated at 1 and the grid spacing Am the former is a property of the function itself while the latter is a numerical parameter which can be changed Thus7 for nite the numerical approximation depends lineraly on the parameter Azi If we were to half Am we ought to expect a linear decrease 113 TAYLOR SERIES 81 in the error for suf ciently small Ax We will use the big Oh77 notation to refer to this behavior so that TE 0Axi In general if Am is not constant we pick a representative value of the grid spacing7 either the average of the largest grid spacing Note that in general the exact truncation error is not known7 and all we can do is characterize the behavior of the error as Ax a 0 So now we can write 311 7 7 ui1 W as mi 0Az 1117 mi The taylor series expansion can be used to get an expression for the truncation error of the backward difference formula u Azil 6211 Asgil 3311 3 3x3 1118 1i 1i 1i where AMA 7xi1 We can now get an expression for the error corresponding to backward difference approximation of the rst derivative A il 3311 3 6x3 7 7 Axi1 7 311 A iil 3x 7 Aiiil 6211 2 6x2 1119 Truncation Error It is now clear that the truncation error of the backward difference7 while not the same as the forward difference7 behave similarly in terms of order of magnitude analysis7 and is linear in As 6711 as ui uiil 7 0A 1120 AM lt z lt gt 1239 Notice that in both cases we have used the information provided at just two points to derive the approxirnation7 and the error behaves linearly in both instances Higher order approximation of the rst derivative can be obtained by combining the two Taylor series equation 1115 and 1118 Notice rst that the high order derivatives of the function u are all evaluated at the same point xi and are the same in both expansions We can now form a linear combination of the equations whereby the leading error term is made to vanish In the present case this can be done by inspection of equations 1116 and 1119 Multiplying the rst by AQA and the second by An and adding both equations we get 7611 7 i AiiilAii 6311 as r 3 6x3 39 12 2 1121 There are several points to note about the preceding expression First the approx imation uses information about the functions u at three points MA and n1 1 Ui1 W A if A441 A i A iil 1 A iil 82 CHAPTER 11 FIN TE DIFFERENCE APPROXIMATI ON OF DERIVATIVES Second the truncation error is TE 0Am1Ax and is second order that is if the grid spacing is decreased by 12 the TE error decreases by factor of 22 Thirdly the previous point can be made clearer by focussing on the important case where the grid spacing is constant Am1 Am Am the expression simpli es to 7 Azz 3 ui1 U171 7 5M 31 3x3 2A7 1122 m M Hence for an equally spaced grid the centered difference approximation converges quadratically as Ax 7 0 311 7 7 U11 U171 3m 2 2M 0Az 1123 1239 Note that like the forward and backward Euler difference formula the centered dif ference uses information at only two points but delivers twice the order ofthe other two methods This property will hold in general whenever the grid spacing is con stant and the computational stencil ie the set of points used in approximating the derivative is symmetric 1132 Higher order approximation The Taylor expansion provides a very useful tool for the derivation of higher or der approximation to derivatives of any order There are several approaches to achieve this We will rst look at an expendient one before elaborating on the more systematic one In most of the following we will assume the grid spacing to be constant as is usually the case in most applications Equation 1122 provides us with the simplest way to derive a fourth order ap proximation An important property of this centered formula is that its truncation error contains only odd derivative terms 1111 7 114 311 sz 3311 Az 3 AH 3711 szm 62m11u 7777777m 7m 2Az 3m 31 3x3 51 3x5 71 M7 2m 116z2m1 1124 The above formula can be applied with As replace by 2Ax and 3Ax respectively to get 1112 7 114 7 311 2A22 3311 2A24 3511 2A26 3711 4m amp 1 31 6x3 1 51 6x5 1 71 6x7 1 0Af 3925 1113 7 ulng 7 311 3A22 3311 3A24 3511 3A26 3711 6Az 1 a 31 6x3 1 51 6x5 1 71 6x7 1 0Af 3926 It is now clear how to combine the different estimates to obtain a fourth order approximation to the rst derivative Multiplying equation 1124 by 22 and 113 TAYLOR SERIES 83 substracting it from equation 11257 we cancel the second order error term to get 8ui1 7 ill1 7 11139er 7 ill2 7 3U i 4Az4 3511 20Ax6 3711 12m amp 5x 62657 7 a oltm8gt 1127 Repeating the process for equation but using the factor 32 and substracting it from equation 11267 we get 27ui1 7 ui1 7 ui3 7 uiig i 6711 7 9A4 6511 7 6711 48Az 7 3m 5 3x5 7 3x7 Although both equations 1127 and 1128 are valid7 the latter is not used in practice since it does not make sense to disregard neighboring points while using more distant ones However7 the expression is useful to derive a sixth order ap proximation to the rst derivative multiply equation 1128 by 9 and equation 1128 by 4 and substract to get 1 0Az8 1128 45uH1 7 ui1 7 9ui2 7 ui2 ui3 7 uiig i 9U 6711 8 60Az 7 Q M 1129 The process can be repeated to derive higher order approximations 1133 Remarks The validity of the Taylor series analysis of the truncation error depends on the existence of higher order derivatives If these derivatives do not exist7 then the higher order approximations cannot be expected to hold To demonstrate the issue more clearly we will look at speci c examples Example 1 The function sin 7 is in nitely smooth and differentiable7 and its rst derivative is given by um 7TCOS 7 Given the smoothness of the function we expect the Taylor series analysis of the truncation error to hold We set about verifying this claim in a practical calculation We lay down a computational grid on the interval 71 S x S 1 of constant grid spacing Ax 2M The approximation points are then xi Z39Az 717 239 01 7 M Let E be the error between the nite difference approximation to the rst derivative7 111 and its analytical derivative um 61 111 7 um 1130 The numerical approximation 111 will be computed using the forward difference7 equation 11177 the backward difference7 equation 11207 and the centered dif ference approximations of order 27 4 and 67 equations 11227 11277 and 1129 We will use two measures to characterize the error 6139 and to measure its rate of decrease as the number of grid points is increased One is a bulk measure and 84 CHAPTER 1 1 FINI TE DIFFERENCE APPR OXIMATI ON OF DERIVATIVES 1 001 05 0005 0 0 705 70005 711 70 5 0 0 5 1 70011 70 5 0 0 5 1 Figure 112 Finite difference approximation to the derivative ofthe function sin 7 The top left panel shows the function as a function of x The top right panel shows the spatial distribution of the error using the Forward difference black line7 the backward difference red line7 and the centered differences of various order magenta lines for the case M 1024 the centered difference curves lie atop each other because their errors are much smaller then those of the rst order schemes The lower panels are convergence curves showing the rate of decrease of the rms and maximum errors as the number of grid cells increases 113 TAYLOR SERIES 85 consists of the root mean square error7 and the other one consists of the maximum error magnitude We will use the following notations for the rms and max errors M llEllz M 26 1131 i0 Helloo giggled 1132 The right panel of gure 112 shows the variations of E as a function of z for the case M 1024 for several nite difference approximations to um For the rst order schemes the errors peak at i12 and reaches 001 The error is much smaller for the higher order centered difference scheme The lower panels of gure 112 show the decrease of the rms error on the left7 and maximum error on the right as a function of the number of cells M It is seen that the convergence rate increases with an increase in the order of the approximation as predicted by the Taylor series analysis The slopes on this log log plot are 1 for forward and backward difference7 and 27 4 and 6 for the centered difference schemes of order 27 4 and 67 respectively Notice that the maximum error decreases at the same rate as the rms error even though it reports a higher error Finally7 if one were to gauge the ef ciency of using information most accurately7 it is evident that for a given M7 the high order methods achieve the lowest error Example 2 We now investigate the numerical approximation to a function with nite differentiability7 more precisely7 one that has a discontinuous third derivative This function is de ned as follows lt 7T 7TCOS7T 0lt 7mg 7 x 7r Notice that the function and its rst two derivatives are continuous at z 07 but the third derivative is discontinuous An examination ofthe graph of the function in gure 113 shows a curve7 at least visually the so called eye ball norm The error distribution is shown in the top right panel of gure 113 for the case M 1024 and the fourth order centered difference scheme Notice that the error is very small except for the spike near the discontinuity The error curves in the lower panels show that the second order centered difference converges faster then the forward and backward Euler scheme7 but that the convergence rates of the fourth and sixth order centered schemes are no better then that of the second order one This is a direct consequence ofthe discontinuity in the third derivative whereby the Taylor expansion is valid only up to the third term The effects of the discontinuity are more clearly seen in the maximum error plot lower right panel then in the mean error one lower left panel The main message of this example is that for 86 CHAPTER 1 1 FINI TE DIFFERENCE APPR OXIMATI ON OF DERIVATIVES x 10 15 1 1 05 05 g w 0 0 705 705 71 71 71 70 5 0 0 5 1 71 705 0 0 5 1 X X 100 100 N 7 E g 10 5 g 10 5 E 10710 1 2 3 4 10710 0 l 2 3 4 10 10 10 10 10 10 Figure 113 Finite difference approximation to the derivative of a function with discontinuous third derivative The top left panel shows the function which7 to the eyeball norm7 appears to be quite smooth The top right panel shows the spatial distribution of the error M 1024 using the fourth order centered difference notice the spike at the discontinuity in the derivative The lower panels are convergence curves showing the rate of decrease of the rms and maximum errors as the number of grid cells increases 113 TAYLOR SERIES 87 functions with a nite number of derivatives7 the Taylor series prediction for the high order schemes does not hold Notice that the error for the fourth and sixth order schemes are lower then the other 37 but their rate of convergence is the same as the second order scheme This is largely coincidental and would change according to the function 1134 Systematic Derivation of higher order derivative The Taylor series expansion provides a systematic way of deriving approximation to higher order derivatives of any order provided of course that the function is smooth enough Here we assume that the grid spacing is uniform for simplicity Suppose that the stencil chosen includes the points mj such that 239 7 l S j S 7 There are thus l points to the left and 7 points to the right of the point 239 where the derivative is desired for a total of r l 1 points The Taylor expansion is mAs mA2M 3 um l mm A 1 2x 3 szm mA4M mAx5 4 m 5 uim 111 1 1133 for m fl 7 Multiplying each of these expansions by a constant am and summing them up we obtain the following equation T T T Ax 3U Z amuimvlt 2 am ul lt Z mam m7lm70 m7lm70 m7lm70 T sz 3211 2 7 7 07 a 2x 6x2 7 Axs 3311 3 7 7 07 a 3x 6x3 T m4a if m7lm70 m 4 64 T Ax5 3511 m7lm70m am l MS 139 1134 It is clear that the coef cient ofthe k th derivative is given by bk 2ilm0 mkam Equation 1134 allows us to determine the r l coef cients am according to the derivative desired and the order desired Hence ifthe rst order derivative is needed at fourth order accuracy7 we would set b1 to 1 and b2 0 This would provide us with four equations7 and hence we need at least four points in order to deter mine its solution uniquely More generally7 if we need the k th derivative then the 88 CHAPTER 1 1 FINI TE DIFFERENCE APPR OXIMATI ON OF DERIVATIVES highest derivative to be neglected must be of order k 10 7 17 and hence k 10 7 1 points are needed The equations will then have the form T bq Z 771an qk7q12kpi1 1135 m7lm70 where qu is the Kronecker delta qk 1 is q k and 0 otherwise For the solution to exit and be unique we must have l r k p Once the solution is obtained we can determine the leading order truncation term by calculating the coef cient multiplying the next higher derivative in the truncation error series T bk1 Z mkpam 1136 m7lm70 Example 3 As an example of the application of the previous procedure7 let us x the stencil to r 1 and l 73 Notice that this is an off centered scheme The system of equation then reads as follows in matrix form 3 72 71 1 LS b1 73 2 fly 12 072 b2 73 23 713 13 01 b3 1137 73 72v Th4 14 01 b4 If the rst derivative is desired to fourth order accuracy7 we would set b1 1 and b2 07 while if the second derivative is required to third order accuracy we would set b134 0 and b2 1 The coef cients for the rst example would be 13 1 12 7 1 12 ail i E 718 1138 11 3 1 135 Discrete Operator Operators are often used to describe the discrete transformations needed in ap proximating derivatives This reduces the lengths of formulae and can be used to derive new approximations We will limit ourselves to the case of the centered difference operator qu 7 meg smut W 1139 uil 7 uiil 6w 2m 2 um 0Az2 1140 62W um 0Az2 1141 2Az 114 POLYNOMIAL FITTING 89 The second order derivative can be computed by noticing that 611 61 6111 61mm 0Az2 1142 6 11 7 ui i 0A 2 1143 m A96 7 1461 96 1 6 6 i 0 A 2 1144 Elt 7 i M gt 39 uH 7 2u u1 W umOAx2 1145 The truncation error can be veri ed by going through the formal Taylor series analysis Another application of operator notation is the derivation of higher order for mula For example7 we know from the Taylor series that A 2 6216 um 0464 1146 lf 1 can estimate the third order derivative to second order then I can substitute this estimate in the above formula to get a fourth order estimate Applying the 6 operator to both sides of the above equation we get 2 63621111 631 ATqu 0Az4 um 0Az2 1147 Thus we have sz 2 2 Rearranging the equation we have 3 umlm lt1 i if 62W 0Az4 1149 114 Polynomial Fitting Taylor series expansion are not the only means to develop nite difference approxi mation An another approach is to rely on polynomial tting such as splines which we will not discuss here7 and Lagrange interpolation We will concentrate on the later in the following section Lagrange interpolation consists of tting a polynomial of a speci ed defree to a given set of pairs The slope at the point 1 is approximated by taking the derivative of the polynomial at the point The approach is best illustrate by looking at speci c examples 90 CHAPTER 1 1 FINI TE DIFFERENCE APPR OXIMATI ON OF DERIVATIVES 1141 Linear Fit The linear polynomial i 7 ii i 7 144 W 7 A95 1117 The derivative of this function yields the forward difference formula 7 7 11111 7 ul uw 7 3m Ax 1239 1151 A Taylor series analysis will show this approximation to be linear Likewise if a linear interpolation is used to interpolate the function in 111 3 x 3 xi we get the backward difference formula 1 142 Quadratic Fit It is easily veri ed that the following quadratic interpolation will t the function values at the points 1 and mil 95 7 95095 7 95m 95 7 95H95 7 95m 95 7 95H95 7 95 LZW quot 15171 7 quotuquot 11 2Az2 Ax 2Az2 1152 Differentiating the functions and evaluating it at we can get expressions for the rst and second derivatives 3L2 Ui1 7 71 7 1153 h m 2Az 3sz Ui1 7 2W U171 1154 3x2 B Ax Notice that these expression are identical to the formulae obtained earlier A Taylor series analysis would con rm that both expression are second order accurate 1143 Higher order formula Higher order fomula can be develop by Lagrange polynomials of increasing degree A word of caution is that high order Lagrange interpolation is practical when the evaluation point is in the middle of the stencil High order Lagrange interpola tion is notoriously noisy near the end of the stencil when equal grid spacing is used7 and leads to the well known problem of Runge oscillations Spectral methods that do not use periodic Fourier functions the usual sin and cos functions rely on unevenly spaced points 114 POLYNOMIAL FITTING 91 Equa rSpaczd Lagmngz Immolan39on EquallrSpaczd Lagmngz Immolan39on xquot CX 05 n v x n u Iquot quotx quot7 I r H u 39 v 705 l r I l I l L I 1 I 399 1 5 I39 r K I 72 39 39 725 39 39 75 l I I II II 755 u r 705 o 05 1 x r 705 o 05 GWL obamrSpaczd Lagmngz Immolan39vm Figure 114 Illustration of the Runge phenomenon for equally spaced Lagrangian interpolation upper gures The right upper gure illustrate the worsening am plitude of the oscillations as the degree is increased The Runge oscillations are suppressed if an unequally spaced set of interpolation point is used lower panel here one based on Gauss Lobatto roots of Chebyshev polynomials The solution black line refers to the exact solution and the dashed lines to the Lagrangian in terpolants The location of the interpolation points can be guessed by the crossing of the dashed lines and the solid black line 92 CHAPTER 11 FIN TE DIFFERENCE APPROXIMATI ON OF DERIVATIVES To illustrate the Runge phenomenon well take the simple example of interpo lating the function 1 WW in the interval S 1 The Lagrange interpolation using an equally spaced grid is shown in the upper panel of gure 114 the solid line refers to the exact function f while the dashed colored lines to the Lagrange interpolants of different orders In the center of the interval near z 0 the difference between the dashed lines and the solid black line decreases quickly as the polynomial order is increased However near the edges of the interval the Lagrangian interpolants oscillates between the interpolation points At a xed point near the boundary the oscillations7 arnplitude becornes bigger as the polynomial degree is increased the amplitude of the 16 order polynomial reaches of value of 17 and has to be plotted separately for clarity of presentation This is not the case when a non uniform grid is used for the interpolation as shown in the lower left panel of gure 114 The interpolants approach the true function in the center and at the edges of the interval The points used in this case are the Gauss Lobatto roots of the Chebyshev polynomial of degree N 7 1 where N is the number of points 1155 Chapter 5 Fortran Procedures Subprograms and Functions Numerical programming often requires repeating the same calculations with differ ent data and variables Examples include having to integrate different functions or the same function over different intervals repeatedly Having to recode the same chunk of fortran statement is wasteful of time and error prone as typos and bugs can be introduced in the process7 moreover7 the time we have invested in validating a speci c piece of code would have to be duplicated again and again Most pro gramming languages provide a mechanism to avoid these repetitions and pitfalls ln FORTRAN this is done simply and elegantly by turning algorithmically related code segments into subprograms or procedures There are many kinds of sub programs in FORTRAN subroutines7 functions7 and BLOCK DATA subprograms Moreover7 subroutines and functions can come in two flavors7 internal and exter nal The func function already encountered in the trapezoidal rule integration is one example of an internal procedure Subprograms are essential to build complex software as it allows programmers to think of their code in terms of algorithmic units and permits some degree of code abstraction The bene ts is that the de velopment and validation of code can be done exibly and in manageable units Validated software can also be reused without having to worry repeatedly about its correctness 51 Modular Design The example of trapezoidal integration will be used again to illustrate the use of external procedures the code for an adaptive trapezoidal integration is listed in T Clearly one can see a division of tasks among the different portions ofthe code More speci cally7 one can distinguish the following phases in the program 0 An initialization phase that primarily reads in the input data lines 8 9 33 34CHAPTER5FORTRANPROCEDURESSUBPROGRAMSA lFUNCUYONS 1 2 3 4 5 6 7 A A program trap implicit none integer N real 2 ab integer maxiteration integer i iteration counters real dx Ans xi number of intervals interval extrema maximum number of iterations allowed error Ansold tolerance print Enter a btolerance initial N maxiteration read ab tolerance N maxiteration do iteration 1maxiteration dx b arealN Ans O5 myfunca myfunCb do i 1N 1 xi idx a Ans Ans myfuncxi enddo Ans Ans dx if iteration gt 1 then error absAns Ansold if error lt tolerance then exit else N 2N Ansold Ans endif else error 00 endif print N Ans error enddo print N Ans error stop contains real function myfuncx implicit none real intentin x myfunc x5 return end function myfunc end program trap iteration count grid size Evaluate function at a and b Integral estimate with N points one iteration before check absolute error converged exit loop if not converged do following double number of intervals store old answer print final answer Figure 51 Code for adaptive trapezoidal integration 52 FUNCTIONS AND SUBROUTINES 35 o The calculation of the trapezoidal estimate of the integral lines 11 17 0 The convergence check to see if the result caculcated is reliable lines 1828 The logic of the program is pretty simple the user input the parameters that in clude an initial estimate for the number of intervals and the desired accuracy The program then iterates maxiteration times calculating the integral with different values of N If ther error in two consecutives estimates is less then the required tolerance the loop is exited line 21 else the number of intervals is doubled line 23 the previous result is stored line 24 and the calculation is repeated starting from line 10 One drawback of the present design is that we have to check at line 18 if have already an estimate of the answer before performing the convergence check 52 Functions and subroutines The code listed in gure already uses a function call myfunc to calculate the function It was used to avoid having to explicitly code the formula for the ini tialization of the sum line 12 where it is used twice and its repeated evaluation in line 15 Now that the code has been packaged in the function myfunc Any changes to the integrand need to be performed once by modifying that function the rest of the code would remain untouched It is important to stress at this point that myfunc is a self contained procedure that receives an input X and returns an output named myfunc At run time the function is called twice at line 12 once for a and b and upon each iteration of the loop at line 15 Each time the dummy argument of the function X is associated in turn with the variables a b and Xi One important feature of functions is that they can accept multiple arguments as input but can only return one argument as output Another functional unit amenable for packaging77 in a function is the code that actually evaluates the integral It needs three input variables and produces exactly one output The code to achieve this is shown in gure T The procedure Trapezoid is a real function that receives the interval parameters and the num ber of cells and that returns in its name the trapzoidal estimate to the integral The list of variables in the function statement are called dummy arguments because they are only place holders for the actual arguments that will be passed to the procedure These dummy arguments must have their data type declared it is also good practice to declare their intent whether they are input argument intent in output arguments intentout or both intent inout Finally it is possible to use with some stretch a subroutine call to package the input statements and another one to package the iterative convergence check on the trapezoidal estimate The extra code needed to achieve these results is shown in gure 36CHAPTER5FORTRANPROCEDURESSUBPROGRAMSA lFUNCUYONS 1 2 3 4 5 6 7 8 9 11 28 29 30 31 32 33 34 35 36 37 38 39 4O 41 42 43 44 program trap implicit none integer N number of intervals real 2 ab interval extrema integer maxiteration maximum number of iterations allowed integer i iteration counters real dx Ans xi error Ansold tolerance print Enter a btolerance initial N maxiteration read ab tolerance N maxiteration Ansold TrapezoidabN first estimate of integral do iteration 2maxiteration iteration count N 2N double number of intervals Ans TrapezoidabN new estimate of integral error absAns Ansold absolute error if error lt tolerance then exit converged exit loop else Ansold Ans save old answer if not converged endif print N Ans error monitor adaptation progress enddo print N Ans error print final answer stop contains Numerical integration via trapezoidal rule real function TrapezoidabN implicit none integer intentin N real intentin ab integer i xdx res dx b arealN res 05 myfuncamyfunCb do i 1N 1 x idx a res res myfuncx enddo Trapezoid dx res return 45 end function Trapezoid AAAlz A AVALHAN All mnwAnAzAn mAwmnAm NAL A 52 PTUVCTYONSAANISUBIUDUTYNES 37 program numinteg implicit none integer NmaXiteration real abainteg tolerance call InputabtoleranceNmaXiteration call AdaptTrapezoidaintegabtolerlanceNmaXiteration print lntegral 7 ainteg stop contains the functions Trapezoid and myfunc should go below here I Adaptive integration I subroutine AdaptTrapezoidaintegabtoleranceNmaXiteration implicit none integer intentinout N maXiteration real intentin ab tolerance real intentout ainteg integer iteration real Ans Ansolderr Ansold TrapezoidabN first estimate of integral do iteration 2maXiteration iteration count N 2N double number of intervals Ans TrapezoidabN new estimate of integral error absAns Ansold absolute error if error lt tolerance then exit converged exit loop else Ansold Ans save old answer if not converged endif print N Ans error monitor adaptation progress enddo return end subroutine AdaptTrapezoid end program numinteg Input parameters subroutine InputabtoleranceNmaXiteration implicit none integer intentout N maxiteration number of subinterval real intentout ab tolerance interval print Enter a btolerance initial N maxiteration read ab tolerance N maxiteration return end subroutine Input 4 All 14 4 m n w A n Abln wmAnmnm 38 CHAPTER 5 FORTRAN PROCEDURES SUBPROGRAMS AND FUNCTIONS 521 Remarks about the main program 0 The main code has been shortened because the work has been relegated to subroutines The logical sequence of execution is however neatly encapsulated by the sub program units by giving subroutines informative names Different subrou tines in a single program must have different names A subroutine is invoked by the call statement A call is a kind of branch transferring execution to the subroutine When a subroutine is nished con trol is returned to the next statement after the invoking call The variables that are used in each subroutine invocation are passed to the routine in a parenthesized list usually referred to as an argument list more speci cally they are called actual arguments in the call statement 522 Remarks about the routines o The beginning and end of a subroutine are de ned by the pairs of subroutine and end subroutine statements The pair de ne the scope of a subrou tine The end statement must be present to tell the compiler there are no more source line in the routine Classical fortran does not support the end subroutine statement only end The name given in the subroutine must match the one in the call statement for the subroutine to be invoked When Input is invoked control is transfered to the rst executable statement inside the routine which is print Enter ab N After the return control is transfered back to the rst statement after call Input Every subprogram must execute a return or a STOP which returns control to the Operating System 0 The arguments listed in the subroutine statement must match in number order and type the argument appearing in the call statement that invoked the subroutine However the names of the vaiiables need not be the same as the name in the calling routine This is essential if the same routine need to be called repeatedly with different actual parameters The variables listed in the subroutine statement are called dummy argument because their names are unrelated to those in other routines including the calling routine For example the routine Trapezoid returns its result in variable res which is then stored in variable ainteg in subroutine AdaptTrapezoid The argument list is declared in each subroutine You can see that an at tribute is attached to each argument that informs the compiler how we are going to use these arguments Some of these are solely input arguments 53 CALL B Y REFERENCE 39 and will not be changed by the routine these are given the intentin at tribute Others have no value upon entry to the subroutine but their value is assigned upon exist they are given the intentout attribute Finally some arguments have a value upon entry that is modi ed in the course of the subroutine execution these are given the intentinout attribute The intent attribute is not part of classical FORTRAN but is extremely useful in documenting the code by immediately telling the programmer what vari ables are modi ed Furthermore declaring the attribute allows the compiler to rst check and alert the programmer if a value is changed inadvertantly and second to perform some optimization The variables X and i are local to subroutine Trapezoid and they have nothing to do with variables having the same name in other routines includ ing the calling routine The scope of these variables is local to subroutine Trapezoid Local variables are not necessarily saved across subroutine invocation even though some compiler do so by default If we want to insure that local vari ables keep their values across we would use the attribute save in declaring the variable The main program and the subroutines it calls are each compiled separetely one at a time This separate compilation is the reason that the names of variables in dummy argument lists and local variables having nothing to do with the names of variables in other routines The compiler in translating each routine into machine language has no way of knowing about any other routines to be used in the program Only the loader knows about all the pieces that must be linked together to create the executable by that time all variable names would have been translated into memory addresses Separate compilation is one ofthe great features of FORTRAN as it allows the building and use of libraries Subprograms can invoke other subprograms to any depth The only exception is that subprogram can invoke themselves in classical FORTRAN Recursion is actually a new feature added to FORTRAN 90 53 Call by reference How are arguments transferred from the invoking subroutine to the called sub routine One could imagine that a subroutine copies actual arguments into its dummy arguments at invocation performs the calculations and then the output arguments are copied back after the return That would actually be a waste of memory and CPU both precious commodities in the early days of computing 40 CHAPTER 5 FORTRAN PROCEDURES SUBPROGRAMS AND FUNCTIONS FORTRAN uses instead what is known as call by reference or call by name because what is passed to the subroutine is not the value of the actual argument but its memory address Subprograms manipulates the actual parameters directly by addressing the memory location where they are stored Address calculation takes place entirely behind the scenes and memory locations are not directly avail able to the FORTRAN programmer So if a subrprogram changes the value of a dummy argument7 the value of the actual argument gets changed immediately For this reason7 passing arguments that have the same name can lead to surprising results if their values are changed and can be potentially dangerous This is called aliasing of arguments 54 Function Subprograms Function subprograms are very similar to subprograms but there are some impor tant differences Functions return a single variable in the name of the function itself7 no sep arate variable need be devoted to the return value A subroutine does not have a type and does not need to have one But a function7s name returns a value and hence must have a type The compiler knows about its built in functions like sqrt but it has to be told the type of user de ned functions There are speci c functions for every variable type7 dsqrt for example returns a double precision variable Fortran7 however7 allows some functions to have generic names7 and the function invoked depends on the data type of the input argument For example sqrt 4 O returns a 32 bit real result whereas sqrt 4d0 returns a 64 bit real since 4dO is a double precision constant A subroutine is invoked with call whereas a function is invoked simply by mentioning its name in an expression A subroutine is allowed to not have any parameters7 but classical FORTRAN insists that a function must have at least one parameter The scope of a function are delimited by the function and end function statements The functions name must be assigned a result and the compiler should com plain if it is not 0 Functions and subroutines can invoke other functions Special considerations should be taken to handle recursive functions 55 INTERFACE BLOCKS AND MODULES 41 55 Interface Blocks and Modules The use of functions and subroutines hinges on effectively passing the data be tween the calling and called procedures correctly There are many pitfalls that expose programmers to errors7 such as missing arguments7 wrong data types cor respondence7 The language introduces a number of constructs to minimize the chances of introducing bugs The rst one consists of using interface blocks that de ne the type of arguments passed7 and the other relies on using modules lnterface checking is but one of the bene ts of using modules and we will defer their presentation to the next section 551 Interface block An interface block is used to alert the calling program to the list and type of arguments a subroutine expects It is basically7 the subroutine or function but stripped of all its execution statements and declaration of local variables An interface block would look like the following interface function myfunc X real intentin X end function myfunc end interface and would be inserted in the program unit that calls the function myfunc Another example is interface subroutine InputabN real intentout ab integer intentout N end subroutine Input end interface 552 Modules and procedures The above is an example of a user coded interface When the program unit called is contained within another program unit or within a module7 the interface block is built automatically by the compiler and the programmer need not bother with cutting and pasting these pieces of codes Modules and procedures are powerful tools that allow program units to share data and procedures For the moment we are going to limit ourselves to using these constructs solely for the purpose of building these interfaces automatically An example module that gathers the functionality required in estimating the integral and verifying it is shown in gure 420HAPTER5FORTRANPROCEDURESSUBPROGRAMSA lFUNCUYONS module trapez contains I Adaptive integration subroutine AdaptTrapezoidaintegabtoleranceNmaXiteration implicit none integer intentinout N maXiteration real intentin ab tolerance real intentout ainteg integer iteration real Ans Ansolderr Ansold TrapezoidabN do iteration 2maXiteration N 2N Ans TrapezoidabN error absAns Ansold if error lt tolerance then exit else Ansold Ans endif print N Ans error monitor adaptation progress enddo return end subroutine AdaptTrapezoid I first estimate of integral iteration count double number of intervals new estimate of integral absolute error converged exit loop save old answer if not converged Numerical integration via trapezoidal rule I real function TrapezoidabN implicit none integer intentin N real intentin ab integer i real XdX res dX b arealN res 05 myfuncamyfuncb do i 1N 1 X idX a res res myfuncx enddo Trapezoid dX res return end function Trapezoid end module trapez E br mlc E AI AAAJ JA Avnmwl Ln A m n n w N 1n4n 14 F mnmn1 m N m A doch m 5 6 RECURSIVE PROCEDURES 43 program numinteg use trapez implicit none integer Nmaxiteration real abainteg tolerance call InputabtoleranceNmaxiteration call AdaptTrapezoidaintegabtolerlanceNmaxiteration print lntegral 7 ainteg stop contains the functions Trapezoid and myfunc should go below here Figure 55 Usage of a module T the corresponding code that uses the module is shown in T Notice that the main program has the statement use trapez to instruct the compiler to associate the subroutine AdaptTrapezoid with the one declared within the module The module can be in a separate le or within the same le as the main program In either cases7 the use requires that the module trapez be compiled prior to the main program Hence if it is included in the le containing the main program7 it must precede it The compiler can hence process and produce the information that the module need provide the main program with 56 Recursive Procedures Fortran by default does not permit functions of subroutines to call themselves either directly or indirectly However7 complex algorithms can often be cast elegantly in terms of recursive calls7 that is the algorithm divides the steps into subtasks that are identical by require different data In that sense recursive procedures resemble do loops by are much more powerful They are also more expensive as a procedure call7 and as local variable must be created for each instance the procedure is called Let us revisit the example of trapezoidal integration Our adaptive procedure doubled the number of points uniformely7 regardless of the shape of the function This can be wasteful in places where the function behaves mostly like a straight line since the trapezoidal error is already small there The error analysis reveals that the errors are largest in regions where the curvature of the function is large A more astute adaptation strategy would resolve the high curvature regions7 only sparse sampling of the function would be necessary in regions where it behaves 44 CHAPTER 5 FORTRAN PROCEDURES SUBPROGRAMS AND FUNCTIONS linearly One strategy to achieve this objective is via the following pseudo code 1 1 trapezoid estimate of the integral is fa i fb Slbia 2 2 2 trapezoid estimate of the integral is bia fa ab fb Sl bia S i i i m 2 2 2flt2gt2 22f where zm a b2 is the mid point of the interval 3 Clearly if SQ 51 then the curve is a straight line in the interval a S x 3 b7 and we are done The answer SQ might still be acceptable if 52 7 51 lt 6 a speci ed tolerance 4 If the error is unacceptable7 then I would need to divide the interval a b into two subintervals7 1 xml and xm b and calculate the trapezoidal estimate separately in each one But these last two steps are exactly identical to the step started at 1 We merely need to repeat steps 1 37 with the tolerance speci ed at 62 because 1 divided the interval to calculate these integrals The code listed in implements the algorithm described above Notice that the statement recursive is required ahead of the function declaration to make the function recursive A critical issue in any recursive procedure is to make sure that it terminates at some points In the present instance the stopping critertion is embedded within the error bounds which does not call the function if it drops below the tolerance 5 6 RECURSIVE PROCEDURES 45 Numerical integration via trapezoidal rule recursive real function TrapezoidRecabtolerance implicit none real intentin abtolerance real xmdx estl est2 error myfunc dX b a interval width xm O5ab interval midpoint estl 05 myfuncamyfuncb dX one trapeze estimate est2 05 dX myfuncxm estl two trapeze estimate error absest2 estl if error gt 30tolerance then estl TrapezoidRecaXmO5tolerance print xm est2 TrapezoidRecxmb05tolerance TrapezoidRec estl est2 else TrapezoidRec est2 endif return end function TrapezoidRec Figure 56 Recursive trapezoidal integration Cell numbers Figure 57 Quadrature points for the recursive trapezoidal rule7 the grid spacing for these boxes is shown in the right hand gure Chapter 4 Flow Control The programs presented so far are example of straightline code in which the executable statements are performed exactly once in the order they appear in Oc casionally we need the computer to process the statement out of order7 depending on the value of the input or of intermediate calculations7 or more then once The process of controling the order of code execution is called ow control We will visit two types of ow controls in the present lecturebranching and looping 41 Conditionals Branching and Logical 411 Quadratic Roots To make our discussion more concrete we will focus on the simple example of calculating the roots of a quadratic equation a2bc0 41 The two roots of the equations are given by 7b i xbz 7 4ac 20 42 i Clearly the nature of the two roots depend on the sign of the discriminant d b2 7 4 For d gt 0 there are two distinct real roots7 for d 0 there is a double root7 and for d lt 0 there are pairs of complex roots Before we launch on how we code a solution of the quadratic equation we il lustrate the process of turning the mathematical steps into FORTRAN code The original charting relied on visual diagrams Here we will use the simple alternative of pseudo code 1 input abc 22 CHAPTER 4 FLOW CONTROL 2 form discriminant d b2 7 4ac 3 a if d lt 0 print a message that no real exist and exit b if d 0 print xi ib2a c ifdgt0printxi gra 412 IF statements Clearly the execution ow depends on the value of d and a test is required to determine which FORTRAN statements should be executed The if statement is designed to do just that It takes the form if logicaleXpression then conditional execution statement 1 conditional execution statement 2 endif The logicaleXpression is the test to be performed If the test returns a value of true the statements between the if and endif are executed7 otherwise they are skipped and execution continues after the endif statement Hence one code that can perform the above calculation is as follows program quadraticroots implicit none real abcdX1X2 read abc print a a 7 b b c c d b2 40ac if d gt 00 then d sqrtd X bd20a X2 b d20a print X1 X1 print X2 X2 endif if d 00 then X1 bsqrtd20a X2 X1 print X1X2 Xl endif if d lt 00 then 41 CONDITIONALS BRANCHING AND LOGICAL 23 print No real root exist endif stop end program quadraticroots The if statement allows multiple forms to be used to build more complicated tests A rst alternative form is the simple if elseif else form if logicalexpressionA then StatAl StatA2 elseif logicalexpressionB then StatBl StatB2 elseif logicalexpressionC then StatCl StatC2 else StatDl StatD2 endif The expressions logicaLeXpressionA7 logicaLeXpressionB7 logicaLeXpressionB7 are tested in order7 and the rst one to test true will have its enclosed statement executed7 all others will be skipped If none of the expressions evaluates to true the statements after the else are executed The following relational operations are built into fortran mathematics F77 or F90 F90 altb a1tb altb agtb agtb agtb agb a1eb altb aZb ageb agtb a b aeqb a ay b aneb ab 24 CHAPTER 4 FLOW CONTROL program quadraticroots implicit none real abcdX1X2 read abc print a a 7 b b C C d b2 40ac if d gt 00 then d sqrtd X bd20a X2 b d20a print X1 X1 print X2 X2 elseif d 00 then X1 bsqrtd20a X2 X1 print X1X2 X1 else print No real root eXist endif stop end program quadraticroots This form allows us to skip one unnessary test Finally it is important to note that if statements can be nested Each level however7 should have its own well de ned endif A third version of the quadratic roots program would read if d lt 00 then print No real root eXist else if d 00 then X1 bsqrtd20a X2 X1 print X1X2 X1 else d sqrtd X1 bd20a X2 b d20a print X1 X1 print X2 X2 endif endif 41 CONDITIONALS BRANCHING AND LOGICAL P Q notP PorQ PandQ PeqvQ PneqvQ true true false true true true false true false false true false false true false true true true false false true false false true false false true false 413 logical data type and associated operators We mentioned above that the test expression must evaluate to true for the enclosed statements to be executed To go along with these logical operations of true or false FORTRAN has an intrinsic data type called logical to express these ideas The logical data types can only have one of two values true or false An example of their use is program quadlogic implicit none real abcdx1x2 logical positive read abc d b2 40ac if d lt 00 then positive false else positive true endif if positive then d sqrtd X1 bd20a X2 b d20a print X1X2 else print No real roots endif stop end program quadlogic Certain operations are allowed on logical statements they are The not operator ips the state The PandQ is true if and only P and Q are true it is false otherwise The PandQ is true if P and07quot Q is true it is false otherwise The PeqvQ is true if P and Q have the same value it is false otherwise The PneqvQ is true if P and Q have different values it is false otherwise 26 CHAPTER 4 FLOW CONTROL 414 select case statement FORTRAN 90 allows a slightly neater form of branching called the case statement It takes the form select case expression case caseselectora fortran statementsa case caseselectorb fortran statementsb default this is optional defaults fortran statements end select The major difference between the if and case statements is that the caseselector must consists of constants or named constants of type integers7 logical or character 42 Do Loops 421 Numerical Integration Occasionally the need arises to go over some statements multiple times7 most often when implementing iterative type procedures We will take the simple example of calculating a de nite integral of the form A b mm 43 and recalling that the integral is the area under the curve The simple trape zoidal rule can be used for this calculation whereby the interval a S x S b is divided into small segments denoted by xill S x S zi where m i171 239 017 7N Note that the width of the interval is constant and is Ax JXJ V The area under the trapezoids is tabulated below We will use the simple shorthand notation of fl We now need merely to add up the small trapezoids AAZ to get an approximation AN to the area under the curve This addition yields N AN 2AA 44 i1 AW f02f12f22f3quot 2fi 2fN722fN71fN45 2 A gf1 l f2 i f3 i quotl fi fN71 46 42 DO LOOPS 1 5t SAW xi m m u m as as m uz u Figure 41 Graphical illustration of the trapezoidal rule integration The graph of the continuous approximation is shown in red and the discrete trapezoidal area are shown in black The error between the two area estimate for a single discrete interval is shown in gray Ax Ax Ax Ax Ax Ax Ax Ax Ax Ax Table 41 Trapezoidal integration 28 CHAPTER 4 FLOW CONTROL It is known that the error in approximating the integral can be bounded by b C1 12N2 3322 dzf EA7Alt 7 lN 7 dz 47 This is an example of an error bound7 and reading it is somewhat tricky for non mathematician The theorem asserts that the error depends on the width of the interval and the maximum curvature of the function in the interval7 both intrinsic properties of the function f and not of the numerical parameter These will not change if the numerical parameters are changed The only numerical parameter is N7 and the error decreases like the square of the number of interval This is an example of quadratic convergence since doubling the number of interval to 2N will decrease the error by a factor of 4 422 Indexed Do loop Clearly the trapezoidal rule is repetitive in nature and what we need to do is initialize the estimate7 evaluate the function and sum it up Here is a code to do just that using the do construct program trapezoid implicit none integer iN real abdXX read abN dX b arealN aint O5funafunb do i 1N X idX a aint aint funx enddo print Integral 7 aint stop end program trapezoid The statements between the do and enddo statement will be repeated N times The variable i is the loop index7 and ranges from 1 to N in the present instance In the rst run through the loop i will be initialized to 17 and will be incremented by 1 after passing the enddo statement Execution then resumes at the beginning of the loop7 and the loop index i is compared to the ending index N7 if iltN the loop is repeated7 each time increment the loop index i by 17 otherwise it is skipped and execution resumes at the rst executable statement after the enddo statement The general form of the looping statement is 42 DO LOOPS do i ibeginiendistride fortranstatements enddo Here ibegin and iend are the lower and upper limit7 and istride is the the increment When the last parameter is absent7 the increment defaults to l the increment must be speci ed for all other values All these must be integer variables All the following are valid loops do i 11o2 print i print the first 4 odd numbers in ascending order enddo do i 2o2 2 print i print the first 10 even numbers in descending order enddo The loop index must never be modi ed in the body of the loop The following is not permissible do i 1100 i 2i 1 PROBLEM loop index modified print i enddo Although it is impossible to branch into a do loop7 it is possible to branch out of it in one of two ways either by using the exit or the cycle statement The exit transfers control out of the loop to the statement after the enddo and interrupts the count while cycle just skip the count Here are two examples do i 1200 if modi20 then cycle skip to the end of the loop and restart endif if igt100 then exit skip to the end of the loop and breakout endif print i prints only odd i enddo 423 Do While loops Not all repetitive calculations have their iteration count known7 sometimes it de pends on the result of intermediate calculations An example is to decide on the number of intervals needed to compute the trapezoidal rule to a given tolerance 30 CHAPTER 4 FLOW CONTROL One strategy is to calculate the integral rst using N interval to obtain a rst approximation AN the number of interval is then doubled and the results would then be AZN If the difference between AN and AZN is small we can be assured that we have reached an error level we can live with otherwise the number of interval is doubled yet again The FORTRAN construct needed to carry out the logic is a do while expression expression must involve variables that are modified during the loop enddo Thus the loop keeps being repeated as long as the test expression evaluates to true the loop terminates when evaluates to false For the integration example we would have7 for example iter O convergencefalsegoontrue aintoldO 0 do while goon insert code for trapezoidal calculation iter iter 1 err absaint aintold if err lt tolerance then convergence true goon false else if iter gt maxiter then goonfalse endif endif enddo 424 In nite do loops A slightly different version of the loop uses triggers an ini te do loop with the exit mechanism embedded within the body of the loop The format is do if loopexpressiontest then exit endif enddo 43 THE INFAMOUS GOTO STATEMENT 31 Clearly the loopexpressiontest must evaluate to true at some point for the loop to exit This form is useful when the programer need to go through the body of the loop at least once before triggering the test 43 The infamous GoTo statement In the middle ages of computing7 programmers frequently used the infamous GoTo statement to skip executable code The format is usually of the form goto labelnb where labelnb is the numeric label of a statement label For example the exe cutable statements in the quadratic root program can be written as if d lt 00 goto 30 if d 00 goto 20 10 d sqrtd X1 bd20a X2 b d20a print X1 X1 X2 X2 goto 40 20 X1 b20a print X1X2 X1 30 print no real root 40 Stop The above example is a relatively easy read as the code is short and the listing ows mostly downward A slightly more contorted code can use the goto to effect loops as in program looping impli C it none integer iis is O i O 10 if i lt 10 goto 40 i i 1 is is i go to 10 40 print is is stop end program looping This code is already spaghetti code in that the downward ow is interrupted by an upward jolt with the goto 10 statement The logic of the algorithm is obfuscated The goto statement can be easily abused and the resulting code ends up looking 32 CHAPTER 4 FLOW CONTROL like a pasta dish with entertwined coding noodles everywhere It is to be avoided as much as possible Chapter 9 Arrays All variables we have used so far are scalar single values meaning that the storage space is reserved to store only real or one double precision Often time however we need to work with a collection of these variables Many engineering and scienti c calculations also involve vectors and matrices which are rectangulars arrays of numbers Fortran provides a special data structure called an array to represent vectors and matrices An array is essentially an ordered list of scalar elements each of which can be referred to by its location in the array 91 Vectors The program in 91 declares two arrays sets their value and compute the mean of one of them The new feature introduced here is in line 4 ofthe program where two vectors are declared X11 and f 11 whose dimension is 11 this is done simply by appending the size in parenthesis 11 to the name of the array There are various ways of declaring arrays in FORTRAN depending on where in the code they are declared ie a main program a subroutine or in a module However array declarations with known dimensions at compile time are allowed anywhere in the code ie the dimensions should be either constants such as the number 11 in the example or named constants declared in an integer parameter statement The latter is preferred in large codes as it simpli es changing the array dimension later by speci cally naming the array dimension The following are valid array declarations in classical FORTRAN integer parameter n2OO real dimensionn Xf real y dimension y 11 By default the index of a FORTRAN array starts at 1 it starts at 0 in O The default can be overridden by specifying explicitly the bounds of an array as in 57 58 CHAPTER 9 ARRAYS program vecl implicit none real parameter pi3141592654 real meandXX11f11 integer i dx 1010 sum 00 do i 111 Xi i 1dX fi Xi sinpiXi mean mean fi print Xi fXi enddo mean mean11 print mean mean stop end program vecl integer xcord 220 In this case the array index starts at 2 and ends at 20 so that the array has 20 2 1 23 elements In classical FORTRAN the elements of an array are stored in contiguous mem ory locations7 so that X occupies a total of 11 adjacent words An element outside the dimension of an array7 such as XO or x127 is a memory location that is probably occupied by some other data such as other variables in the program or by some machine instruction belonging the present executable or possibly the op erating system Unpredictable behavior ensues if your code tries to use an array element outside the dimensioned size ofthe array The particular behavior depends on the data that happens to be in that memory location In the actually happy event that you try to access operating system memory a segmentation fault occurs the program dies an ignominous death and you know for sure something in the code is wrong In other circumstances the erroneous value used does not cause the code to crash but produces wrong or unexpected results but without any warning that something went wrong Such out of bounds array indexing must be zealously gaurded against by carefully checking the array indexing logic in the code never refer to an array indece outside the dimensioned size It is a very good idea to turn array bounds checking in the compiler options when you are testing your code This is usually C in most compilers it causes the compiler to issue a detailed warning of where the bound trespass occured Array bounds checking slows execution and should be turned o once you are ready for production runs 92 MATRICES AND M ULTIDIMENSIONAL ARRAYS 59 92 Matrices and Multidimensional Arrays Matrices in fortran are represented by arrays with two dimensions For example a two and three dimensional array would be declared as follows integer i2d43 i3d456 The matrix i2d has 4 rows the left most dimension and 6 columns the right most dimension and is easy to visualize as shown below i2d11 i2d12 i2d13 i2d21 i2d22 i2d23 i2d31 i2d32 i2d33 i2d4 1 i2d42 i2d43 It is of course harder to visualize three dimensional arrays unless one uses a depth77 for the third right most dimension one can then imagine 6 two dimensional slices of i3d45 z More array dimension can be declared simply by making the parenthesized list of dimensions longer Standard fortran supports up to 7 dimensional arrays The way the arrays are arranged in memory is an important topic that has bearing on how to pass arrays to subroutines But before we get to that it is informative to present some jargon associated with arrays 0 rank is the number of dimensions of an array A vector has rank 1 i2d has rank 2 and i3d has rank 3 extent is the number of elements in a dimension For example i3d has an extent of 4 along dimension 1 5 along dimension 2 and 6 along dimension size is the total number of elements in an array and is simply equal to the product of the extents i2d has size 12 and i3d has size 120 0 Shape is an ordered list of the extents so that i3d hsa shape 456 Do loops and arrays go hand in hand Accessing the members of a multi dimensional array requires nested do loops doj13 doi14 doi14 doj12 print i2dij print i2dij enddo enddo enddo enddo The code on the left access the array column wise whereas the code on the right accesses it row wise The index notation allows easy access to all the array members in any order however access by column wise order is the most ef cient since it takes advantage of of memory traf c rules in modern computers For arrays with large dimensions the left hand code would execute substantially faster then the right hand code 60 CHAPTER 9 ARRAYS 93 Layout of arrays in memory Computer memory is strictly linear with one byte following another like the box cars in a railroad train So there is no natural way of storing 2 dimensional let alone multi dimensional data So multi dimensional arrays must be folded into one dimensional sequences of memory locations ie the Visual is a line rather then a rectangle FORTRAN stores arrays in columnmajor order the rt col umn is listed rst from top to bottom followed by the second column etc Thus the left most index is changing fastest unlike an odometer for example The matrix 11 12 13 14 M 21 22 23 24 31 32 33 34 would be arranged as follows in memory M11 M21 M31 M12 M22 M32 M13 M23 M33 M14 M24 M34 931 Ef cient accent of arrays We are now in a position to explain the ef ciency of accessing memory by column wise order If an array cannot t in the fast memory inside a CPU portions of it must reside in the computers main memory where access time is slow access to the entire array would then require shuttling the data piecewise to the CPU registers where they are acted upon The CPU brings contiguous arrays chunks into memory to save time instead of element by element since every time the CPU reaches into main memory a time penalty is exacted Accessing the array in its contiguous order the way it is laid out in memory ensures that this penalty is minimized If the array ts in registers then these considerations do not apply Recently computer manufacturers have been adding di ferent of types of memory to computers so called L1 and L2 cache memory to mitigate the slowness of main memory the term non uniform memory access NUMA refers to the implementation of different memory hiearachies with di ferent access speeds 94 Passing arrays to subroutines The passing of array arguments particularly multi dimensional arrays to subrou tine is somewhat of a delicate issue that relates directly to the layout of arrays in memory When an array is passed to a subroutine what is actually passed is the memory address of the rst array element So for a subroutine to manage the array elements correctly it has to be informed about the rank and the extent of 94 IMSSHKARRAYSTYSUBROLUYNES 61 the dimensions A subroutine that transposes a matrix would have the following code program examp integer parameter m3n4 integer intentin Amn ATnm Bmn BTnm Cm1n CTnm integer ijij ij 0 do j 1n do i 1m ij ij 1 Aij nj1 1 Cij Aij Eij Aij enddo Em1j O enddo call Transpose1AATmn UK since the entire array is passed call Transpose1BBTmn UK since B has the right size call Transpose1CCTmn wrong since only a portion of C is passed call Transpose2CCTmnm1n UK since leading dimension of C is passed stop end program examp subroutine transpose1PQmn implicit none integer intentin mn integer intentin Pmn integer intentout Qnm integer ij do j 1n enddo return end subroutine transposel subroutine transpose2PQmnlP1Q implicit none integer intentin mnlPlQ integer intentin PlPn integer intentout QlQm integer ij do j 1n do i 1m Bji Aij 62 CHAPTER 9 ARRAYS enddo enddo return end subroutine transpose2 subroutine PrettyPrintPldPmn implicit none integer intentin ldPmn real intentin PldPn integer ij do i 1m do j 1n write6fmt i4 Pij enddo write6 enddo return end subroutine PrettyPrint The rst call to transpose would accomplish the desired result since the entire matrices are passed to the subroutine Sometimes it is useful to declare a large vector in the program and use it as a matrix inside the subroutine This is possible if care is taken in declaring enough element for the vector 95 Allocatable Arrays The new FORTRAN standard 90 and newer allows the dynamic allocation of memory via the allocate statements FORTRAN is a little strict about how it allocates arrays7 particularly as it seeks to preserve the optimization of memory access An array has to be declared explicitly allocatable before one can allocate memory for it at run time This is accomplished with the allocatable attribute The following example allocates a two dimensional array program aralloc implicit none integer mn ij integer allocatable ind do read mn if mn ne 0 and ngtO then ifallocatedind deallocateind deallocate before reallocating allocateindmn allocate array 95 ALLOCATABLE ARRAYS 63 else exit endif do j 1n do i 1m indij j 1mi enddo enddo enddo stop end program aralloc The example also include a number of functions that can be used to inquire the status of allocatable objects Allocatable arrays in FORTRAN must have their rank declared7 you cannot allocate a two dimensional array declared as a one dimensional array The restriction is meant to allow the compiler to carry out the same type of optimization performed on allocatable arrays as regular arrays o The allocatable attribute indicates that memory for the array will be de cided at run time7 and the compiler should defer the memory allocation The rank of the array is declared by the comma separated list of colons The extent along each dimension and the size are deferred until execution time The allocate statement makes available the requested storage and returns the address of the rst element The allocate statement accepts an op tion argument that lets the user know whether the memory allocation was successful allocate indmn statis The stat keyword refers to the name of the optional argument whose value is returned in integer is a non zero value indicates an error It is an error to writeread from an un allocated array The deallocate statement returns the memory to the Operating System The array cannot be addressed and the values stored in it are lost 0 The allocated function returns a logical true if the array is allocated and a false otherwise 64 CHAPTER 9 ARRAYS 96 Array Sections FORTRAN 90 allows the manipulation of portions of the array by using the colon notation For example the following code assign the value 1 to the rst three elements of the array and the value two to the rest program arsec implicit none integer id6 ja42 ka55 id13 1 id362 id1 52 1 odd index elements id1 id3 id5 id2 62 6 even index elements id2 id4 id6 ja ka1 41 2 stop end program arsec The information triplets represent a sort of implied do loop and has the form expressionl expression2 expressionS Each of the expressions must evaluate to an integer7 the rst expression gives a lower bound7 the second an upper bound7 and the third the stride between the elements If the lower bound is omitted it defaults to 17 if the upper bound is omitted it defaults to the maximum bound allocated or declared7 and nally if the third expression is omitted it defaults to 1 The sections can be used effectively in multi dimensional array even though the margin of error increases with the array rank The following example copies a 4 gtlt 2 section of a 5 gtlt 5 The most common useful of this construct is to pass rows of the matrix to subroutines or functions do i 1nrows a ddotrai1ncols vec ncols enddo However this practice is discouraged as it involves copying the data to a temporary array at the entrace and exit to the subroutine the penalty in CPU time increases with the size of the array Chapter 14 Finite Volume Method and Time Discretization 141 Introduction The previous chapter focussed primarily on the issues of turning the partial dif ferential equation into a nite volume formulation7 and on the various ways of reconstructing the function at cell edges from cell values The result is a spatially discretized ordinary di erential equation that we need to integrate in time Here we focus on identifying suitable time discretization that are both accurate and sta ble The stability analysis is however7 more complicated now since we are dealing with systems of ordinary di erential equations Hence instead of a scalar ampli cation factor we have an ampli cation matrix whose impact on the growth of the numerical solution must be investigated This requires the analysis of the structure of the matrix and its eigenvalues The Von Neumann stability analysis provide a simpler approach to study the time stability of the nite di erence scheme 142 The Forward Euler ODE system For the sake of simplicity and later on necessity we assume that the ow eld is constant and positive u gt 0 Furthermore we need an initial condition which we assume to be T705 0 T2 and an inlet boundary condition which take as as T where formulation becomes then and are known quantities The piecewise constant FV ClTj i 7 dt 7 Act 7 3 1 141 ClTj 7777771 77 23N 142 dt u m 7 J 7 7 7 14ZCHAPTER 14 FINITE VOLUME METHOD AND TIME DISCRETIZATION Note thatj 1 is special because ofthe need to apply the inlet boundary condition The above is a system of ODE that must be advanced in time with some suitable time stepping scheme The Forward Euler scheme7 even though it is only rst order in time7 is one of the simplest scheme which allow the unknowns to be marched in time explicitly So the fully discrete system takes the form TT 7 TT n1 n 7 T T i uAt A CE I 143 i M Tj Tj71gt7 N E 144 M711 1 7 MT 145 where u is a dimensionaless quantity known as the Courant number it represents the ratio of the distance covered by the ow in time At over the grid spacing The above equation can also be written in matrix form as T1 1 1 7 ILL T1 T2 In 1 7 u T2 T121 M 1 i M T121 7T7 M 1 i M 7T7 Tj1 M 1 i M Tj1 T3171 1 1 7 M T3171 TN M 1 i M TN 146 The left hand side vector represent the unknown cell averages at time tn1 while the right hand side consists of the evolution matrix M7 the cell average vector at time tn7 and the so called load vector which contain the forcing on the system7 in the present case it contains only the inlet boundary condition Although the matrix M here is relatively simply it is a bi diagonal matrix consisting of 1 7 u on the main diagonal7 u on the rst subdiagonal7 and zero everywhere else7 on is still required to analyze the eigenvalues of this matrix 143 The VonNeumann stability analysis The Von Neumann stability analysis allows us to determine the suitability of the time stepping scheme without complicated matrix algebra Instead the problem is modi ed slightly by assuming the boundary conditions are periodic in space This precludes studying the impact of applying the boundary on the schemes stability7 but however7 this is rarely a problem in practice7 particular for such an essential boudary condition as the inlet condition Under these circumstance7 all rows of the u 143 THE VONNEUMANN STABILITY ANALYSIS 143 matrix are similar consisting of 17 on the diagonal and u on the subdiagonal The second idea of the Von Neumann analysis is to apply a discete Fourier transform7 ie to decompose the signals into waves of arbitrary wavenumbers km T Z mama 147 m1 where is the amplitude of the wave with wavenumber km at time tn Since the system is linear7 it is enough to study the impact of a single wave7 and apply the principle of superposition to calculate the evolution of the original signal The overall stability will be guaranteed if stability is satis ed for all wavenumbers The subscript m will henceforth be dropped and we set Ty TWPMW which we insert in equation 142 to get Tn1 ikzj MTnez kij71 1 7 Mj m ikzj Dividing through by Tnemw7 and using Am 3 7 771 we get the ampli cation factor for wavenumber k Tn1 kA 7 7 7 72 z A 7 Tn 7 1 M 6 149 For stability we require that lAl S 1 which yields the following condition W M lt17 u 6in lt17 u new 17 M M1 xemm 67mm M2 17 217 u17 cos kAx kA 174M17Msm271 Since lAlZ must be less or equal to 1 for stability we need to impose only that M1 7 M Z 0 since all other terms are positive This product is positive only if both terms are positive or both are negative The negative possibility is clearly not possible since it would require u lt 0 and u gt 1 The only possibility is hence for 0 S u S 1 Hence if the Courant number is less then 17 the scheme will be stable for all wavenumbers This is called conditional stability7 and it puts a limit on the maximum time step one can take to advance in time u Am 7 lt E 7 mm 1 At lt u 1414 The maximum time step allowed cannot exceed the time it takes the ow to cross one grid box Notice that if the wind were to blow in the negative direction u lt 0 and we in sisted on using downstream values as opposed to upstream values in reconstructing the uxes7 then u lt 0 and the scheme will be unconditionally unstable 144CHAPTER 14 FINITE VOLUME METHOD AND TIME DISCRETIZATION W T l l l l l l l l 1 n1 uz n3 n4 n5 U7 I8 us I5 kAxn Figure 141 Ampli cation factor for the FE in time7 Donor cell scheme for various Courant numbers A detailed analysis of the ampli cation factor for this scheme will reveal its behavior for all waves supported on the computational grid Notice for example that if A is the wavelength of wavenumber k then km 27139 1415 and hence kA is proportional to the ratio of grid cells per wavelength If a lot of grid cells are used per wavelength so that kA 27TA ltlt 17 the ampli cation factor is close to 1 and the waves are not damped Notice also that on a discrete grid7 the wavelength cannot be arbitrarily small since the smallest wave that can be represented on the grid has wavelength 2A so that kmaxAm 7139 This shortest wave will be damped at a rate determined by the Courant number A 1 7 4M1 7 u Figure 141 shows the dependence of the ampli cation factor on kA and u Notice that on the left side of the spectrum kA ltlt 17 where waves are well resolved7 the ampli cation factor is close to 17 and the waves are not damped as one would expect from solving the PDE exactly The unresolved waves on the other hand with kAx close to 7r are damped The maximum damping occurs for a Courant number of 12 The damping is less for all other Courant numbers
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'