Quantitive Methods in Hydro
Quantitive Methods in Hydro HYD 510
Popular in Course
Popular in Environmental Science & Engineering
This 11 page Class Notes was uploaded by Dr. Wendell Schuppe on Thursday October 15, 2015. The Class Notes belongs to HYD 510 at New Mexico Institute of Mining and Technology taught by Staff in Fall. Since its upload, it has received 31 views. For similar materials see /class/223633/hyd-510-new-mexico-institute-of-mining-and-technology in Environmental Science & Engineering at New Mexico Institute of Mining and Technology.
Reviews for Quantitive Methods in Hydro
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/15/15
New Mexico Tech Hyd 510 Hydrology Program Lab Notes 5 Linear Algebraic Quantitative Methods in Hydrology Equations and Solvers Linear Algebraic Equations What is a linear algebraic equation A linear equation is an equation in the form 611x1 azx2 anxn b where x1 x2 xn are variables and a1 a2 an and b are constants not all zero All variables are only to the rst power and the variables are not arguments in logarithms trigonometric functions etc There are also no products or square roots of variables The following equations are linear 2x 3y 4 z 5 3x 1y 2 x1 2x2 3x3 nxn 1 The following equations are not linear x2 y2 1 xy 2 sin x e l 7x1 3x2 2 x 3 As you may recall a more general de nition of a linear equation uses the notation Lx f 0 where L is an operator x is the dependent variable and f is a known in this case f0 The operator L is linear if Laxa Lx which also implies Lbyb Ly Lxy Lx L01 where a and b are arbitrary constants More compactly L is linear if Laxby a Lx b L01 Note that X can be a vector as it is in most of the algebraic examples given above eg X x1 x2 etc Systems of Linear Algebraic Equations A system of linear algebraic equations is a collection of linear equations that we solve at the same time The system is only solvable if there are as many equations as unknown variables The following system contains two equations and two unknown variables x and y 5l New Mexico Tech Hydrology Program Quantitative Methods in Hydrology 2x y 8 x y 2 Each of these equations is the equation for a line solution to this system of equations is the x y coordinate where these lines intersect The lines plotted on the following plot Notice that they at F2 F4 This is the solution to the system of equations Note that we can add the two equations obtain 3996 or subtract them to obtain x2y 10 these equations also intersect these lines at 24 Therefore we could replace either of the two equations with their sum or difference and still nd same solution ie the same intersection point we multiplied one of these equations by a constant would still obtain the same solution to the system equations It is possible that a system of equations has no solutions like the following system 2xy8 2xy4 The two lines de ned by these equations are and have no points in common It is also possible that a system of equations has in nite number of solutions 2xy8 4x2y16 In this case the two lines coincide and therefore have all points in common 52 Hyd 510 Lab Notes 5 Linear Algebraic Equations and Solvers The are intersect to Both of the Also if we of parallel an New Mexico Tech Hyd 510 Hydrology Program Lab Notes 5 Linear Algebraic Quantitative Methods in Hydrology Equations and Solvers Solving Systems of Linear Algebraic Equations Systems of equations are usually solved by rst putting them in matrix form Recall that we can use matrix multiplication to multiply a 2x2 matrix A by a 2x1 vector z to get the 2x1 vector b as follows A 2 an an z x all all y Azall allxa11xallyblb all all y al1x ally bl So if we had the 2x2 system of equations given above the one that had a solution 2x y 8 x y 2 we could rewrite it in matrix form If we use a112 a1 l and b18 then the rst equation is equivalent to the rst row in the matrix product If we use a21l an l and bz2 then the second equation is equivalent to the second row in the matrix product Therefore our system of equations is equivalent to the following matrix equation 2 l x 8 1 1 y 2 The rst way we will solve the matrix equation Azb is by matrix inversion Although this is conceptually the easiest method to understand it is computationally the most dif cult method and is seldom used Matrix Inversion If we wanted to solve the equation 4x 8 we would multiply both sides by 14 to get 14 4x 14 8 or x2 This operation is analogous to matrix inversion With the matrix equation Az b we multiply both sides by the matrix inverse Al to get A39lAz A39lb Remember that the product of a matrix and its inverse is the identity matrix so this becomes IzA391b But the product of the identity matrix and a vector is just the vector so we have z A39lb Therefore we can obtain our solution vector x by simple matrix multiplication As an example we will use the 3x3 system of equations 53 New Mexico Tech Hyd 510 Hydrology Program Lab Notes 5 Linear Algebraic Quantitative Methods in Hydrology Equations and Solvers x 4y 42 3 3y 22 7 3x 3y 82 1 which as a matrix equation Azb is 14 4x3 032y7 33821 Where the matrix A and vectors z and b are given by respectively 1 4 4 x 3 A0 3 zy b7 3 3 z 1 The inverse of this matrix A is 02308 05641 02564 A391 00769 02564 00256 01154 01154 00385 Multiplying this matrix by b371T we obtain x3 2 12 Why don t we use matrix inversion in practice The expense in this calculation is the cost number of operations and complexity of calculating a matrix inverse If the matrix is nxn the cost is on the order of n3 The direct numerical methods described below for solving Azb have a cost on the order of n2 If eg n 104 as it is on a typical groundwater model then matrix inversion costs roughly 10000 times more takes 10000 times longer than a direct solver If the direct solver takes 1 second when you hit the enter button the solution pops up almost instantaneously If instead you had used a matrix inversion you would go out for a long lunch and come back two hours later and it would still be running Suppose the solution with the direct solver the problem ran in a hour With inversion you could sail around the world solo climb Mt Everest and lead an expedition up the Amazon River When you got back it still would still be running assuming the power hadn t been interrupted by a summer Monsoon thunderstorm Don t use inversion Numerical Solution Methods In these notes we ll visit two fundamental approaches to solving linear algebraic equations like Az b that don t use inversion We call them equation solvers The rst is aDirect Solution Method in which an exact solution is obtained subject only to roundoff error The second approach involves Iterative or Indirect Solution Methods in which we guess the solution for the vector z and then see how well it balances the equations gives a result Az that matches the 54 New Mexico Tech Hyd 510 Hydrology Program Lab Notes 5 Linear Algebraic Quantitative Methods in Hydrology Equations and Solvers vector b then guess again Iterative methods provide formal strategies for guesses that can and often do converge so that the estimated vector z approaches the true vector Of course they are not perfect They have truncation error as well as roundoff error as we truncation the number of iterations We ll visit Jacobi Iteration the simplest and slowest converging of these iterative methods Direct Solution Methods There are generally two steps in direct methods The rst step is to decompose the matrix A into a more easily solved version of itself such as an upper triangular form The second step is to solve the resulting simpli ed equations by back substitution Decomposition is the most expensive step since it involves the most operations We ll review two such methods Gauss Elimination and LU Decomposition In Gauss Elimination the decomposition step is referred to as forward elimination Gauss Elimination Recall that with a system of equations we can add two equations together or multiply an equation by a constant without changing the solution Since the rows in a matrix equation represent linear equations we can add rows of the matrix equation together or multiply a row by a constant and the solution will not change With Gaussian elimination we add rows together and multiply rows by constants in such a way that we obtain an upper triangular matrix for A Then it is straightforward to solve for our unknown variables As an example of Gaussian elimination we will again use the 3x3 system of equations x4y 4z3 3y227 3x3y82l which is written in the form of a matrix equation Azb as 14 4x3 032y7 33821 With Gaussian elimination we rewrite the matrix A as an augmented matrix which means that we attach the column vector b to the right side of A 14 3 I A3903 7 33 1 Now we use row addition and scalar multiplication to get the matrix to the left of the dashed line into upper triangular form In other words we need to get zeros below the diagonal We do this 55 New Mexico Tech Hyd 510 Hydrology Program Lab Notes 5 Linear Algebraic Quantitative Methods in Hydrology Equations and Solvers rst with the leftmost column and the topmost row then we work our way down each column and then across to the remaining columns We ll start with the rst column Since all is on the diagonal it is nonzero in an upper triangular matrix so we do nothing to the first row Since an is zero already we don t do anything to the second row However an is nonzero so we must perform some operations on the third row If we multiply the rst row by 3 and add the result to the third row we can replace the existing third row with this new result and we are left with 8310 14 43 I A390 3 27 0 9 20i 8 Now we ll do the same thing with the second column We won t change an or an because they can be nonzero in an upper triangular matrix But we need to eliminate 832 If we multiply the second row by 3 and add it to the third row we can replace the existing third row with a new one and we are left with a320 14 4l3 I A390 3 27 0 0 26i13 Now that we have an upper triangular matrix we can rewrite the augmented matrix as a matrix equation l4 4x3 032y7 0026213 Now we solve for x y and 2 using back substitution Multiplying the last row of the matrix with the vector we have 26213 so obviously 212 Multiplying the second row of the matrix with the vector we have 3y2z7 and with 21 2 we get y2 Finally multiplying the rst row of the matrix with the vector we have x4y 4z3 and with 212 and y2 we get x 3 Our solution vector is bT3 2 12 The basic algorithm for Gaussian elimination is 1 To obtain all zeros below the diagonal element in column 139 replace the jth row rj with aijaii rirj for all j gt 139 2 Remember that the b vector must also be included in the row operations 3 After obtaining an upper triangular matrix use back substitution to obtain the solution vector x Start with the nth element of the solution vector xnbnam 56 New Mexico Tech Hyd 510 Hydrology Program Lab Notes 5 Linear Algebraic Quantitative Methods in Hydrology Equations and Solvers Then work backward using the following algorithm Xjbjajnxnaj n1Xn1aj j1Xj1ajj LU Decomposition Another method for solving a matrix equation or a system of linear equations written in matrix form is called LU decomposition For the system Azb we use the product of matrices PALU where L is a lower triangular matrix U is an upper triangular matrix and P is a matrix that has a 1 in each row Each ofthese matrices is the same size as A We ll use MATLAB to calculate P L and U We can multiply both sides of Azb by P to get PAzPb and then we can replace PA with LU LUzPb Since Pb is the product of a matrix and a vector the result is also a vector that we ll call b Also Uz is the product of a matrix and a vector so its result is a vector we ll call y Uzy So we ll solve two sets of equations in sequence Lyb and Uzy Since L is lower triangular we can use forward substitution same as back substitution that we used with Gaussian elimination but we start with the rst row to solve for y Once y is de ned we can solve Uzy with back substitution because U is upper triangular As an example we will use the same 3x3 system of equations as before x4y 4z3 3y227 3x3y821 which as a matrix equation Azb is 1 4 4 x 3 0 3 y 7 3 3 z 1 The LU factorization of A gives 1 0 0 3 3 8 0 0 1 L010U03 2 P010 13 1 1 0 0 263 1 0 0 Using matrix multiplication you can verify that PALU We can solve for b Pbl73T Using forward substitution on Lyb we obtain y17133T Then solving Uzy with back substitution we obtain z32l2T which is the same solution we got with Gaussian elimination Applications Note that we solve linear algebra problems in a variety of applications including chemical speciation and the nite difference methods of solving PDEs and 2quotd order BVP ODEs With PDEs involving time and space such as a timedependent diffusion equation we typically solve the lst order time part using a Backwardintime Euler Method and the 2quotd order space part using 57 New Mexico Tech Hyd 510 Hydrology Program Lab Notes 5 Linear Algebraic Quantitative Methods in Hydrology Equations and Solvers nite differences The nite difference equations result in a matrix problem of the form Azb at each time step which must be solved before proceeding to the next time step Thus you must solve the matrix problem a number of times equal to the number of time steps The cool part about direct methods is that you need decompose A only once then you can use this decomposition over and over again only having to backsubstitute at each time step This saves great expense However it also assumes that the original PDE is linear Suppose we are solving a timedependent diffusion problem What would happen if the diffusion coefficient were nonlinear depending on solute concentration Then matix A would change with the value of the dependent variable z or Az The matrix problem Az zb would also be nonlinear To use the direct methods above we need to linearize the matrix problem for example by prescribing the value of z used to determine A We do this using an iterative method usually fixed point iteration or Newton iteration For example with xed point iteration the linearized equation is Az1 zb where i is the xed point iteration number We have to rebuild A decompose it and perform backward substitution of this equation each nonlinear iteration Such calculations are common in chemical speciation codes as well as nite difference codes Indirect or Iterative Solution Methods Iterative methods sometimes called numerical relaxation methods are approximate Given the matrix equation Azb we guess the solution vector z and then see how well it balances the equations gives a result Az that matches the vector b then guess again Iterative methods provide formal strategies for guesses that can and often do converge so that the estimated vector z approaches the true vector Of course they are not perfect They have truncation error as well as roundoff error as we truncation the number of iterations We ll review Jacobi Iteration the simplest and slowest converging of these iterative methods Keep in mind that these are methods for solving linear algebraic equations Don t mix up these iterative equation solvers with iterative methods of handling nonlinearities Just as for direct methods if the problem is nonlinear you need to iterate on that as well Consequently in non linear problems there are nonlinear xed point or Newton iterations buried inside each iteration of the indirect matrix solver That is there are iterations inside of iterations nested for loops If the problem is only mildly nonlinear you can try to handle both in one loop but be careful Jacobi iteration The Jacobi iteration method is based on the following manipulation of the original matrix problem Rewriting Azb as a summation we have Az b 3 2611ij b1 1 Suppose that we guess the values of z and then see if they satisfy this equation If the guess is good then we are done If not we can make another guess and try again and so on iteratively If we make these new guesses intelligently then we should get closer and closer to balancing the equation each time 58 New Mexico Tech Hydrology Program Lab Notes 5 Linear Algebraic Quantitative Methods in Hydrology Equations and Solvers Jacobi iteration is the most basic iterative scheme At each iteration k 1 it solves for zH1 given all the other 2 j 139 from the previous iteration k n n 7 HI k 7 2L1sz 7a 2sz 7bx or 11 11 Jsr 2f L b 7 Z avz fori 12n ampk 01 g z n j 1 To start this process you need an initial guess 2 for each value ofz Ifthe problem is suf ciently well behaved and for decent initial guesses the method converges We need a stopping criteria to stop the process when the solution is suf ciently close We typically use a comparison of old and new values of z If for any i the absolute value of the difference l2 1 7 zfl exeeds the roundoff error of the machinelanguage eps in MATLAB we will continue to iterate there are many other possible criteria We also need to specify the maximum number of iterations we will allow in order to insure that the method doesn t run forever Its typical to call this number maxit Other iterative methods The next most advanced iterative method is call the GaussSeidel GS method It simply notes that in performing the calculation above for 21 you ve already calculated new values for 2 j lt i Why use the old values of these when new ones are available 5m mm Gaussseldelukesalrmslhaltme m Thus 1 1 K zH1 Lbx 7 2amp1sz 7 2asz fori 12n ampk 01 an 11 Jxl GS usually converges just a little bit faster than Jacobi There are large families of iterative methods Some combine direct methods and indirect methods by breaking the domain into blocks Iteration is used between blocks and in each iterative step direct methods are used with in a block In nite difference methods the block is usually taken to be a line of nite difference node points Then there are iterative methods that look a lot like Newton s method of nonlinear iteration They search for a vector 2 that minimizes some measure of the balance of Az b These sophisticated equation solving algorithms are now common in applied codes An example is the Preconditioned Cory ugate Gradient method the gradient is a derivative used to drive the minimization 59 New Mexico Tech Hyd 510 Hydrology Program Lab Notes 5 Linear Algebraic Quantitative Methods in Hydrology Equations and Solvers Linear Equations Exercises in building solvers 1 Write a code to perform Gaussian Elimination Given a square nxn nonsymmetric coefficient matrix A and nxl load vector b write a Gaussian Elimination code to solve the system of linear algebraic equations1 Axb for unknown vector x A suggested pseudocode generic algorithm outline is Input n an nxn matrix A and an nxl vector b Forward Elimination for j l to n l for i jl to n m1 alla1 an 0 b1 b1 7 mutbj for k j 1 to n all an m1 ajk Backsubstitution for i n downto 1 x1 b1 for j il to n x1 x1 7 autx x1 xlan Output Solution vector x 2 Write a code to perform Jacobi iterations Given a square nxn nonsymmetric coefficient matrix A and nxl load vector b write an code to solve the system of linear algebraic equations Axb for unknown nxl vector x using Jacobi iteration l k xl1 b1 61 u 261ny foril2nampk0l 11 1 To start this process you need an initial guess x10 for each value of x If the problem is sufficiently well behaved and for decent initial guesses the method converges We need a stopping criteria to stop the process when the solution is sufficiently close Let s use a comparison of old and new values of x If for any 139 the absolute value 1 Note the shift in the notation for the vector of unknowns from z to X It is traditional to use x for this purpose as in the text Chapter 20 so I ve made the shift in this assignment 5 10 New Mexico Tech Hyd 510 Hydrology Program Lab Notes 5 Linear Algebraic Quantitative Methods in Hydrology Equations and Solvers of the difference lxlk 1 xlkl exeeds the roundoff error of the machinelanguage eps in MATLAB we will continue to iterate there are many other possible criteria We also need to specify the maximum number of iterations we will allow in order to insure that the method doesn t run forever Call this number maxit A suggested pseudocode is Input Input n A b x0 eps maxit Initialization CheckMatrix Condition zeros on diagonal for i 1 ton if lanl lt eps then print lt eps algorithm fails because lanl and exit print b x0 eps maxit kO x x0 Perform Iterations 10 while k S maxit k kl for i l to n xoldl x1 for i l to n sum O for j l to n if j ii then sum sum autxoldj x1 b1 sumam Check convergence for i l to n k1 k if Ix xl Zepsx1 then return to 10 Backsubstitution see Gaussian Elimination Example Output Print k x and exit inside k loop if k exceeds maxit and exit Print Algorithm failed to converge 5ll
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'