PETE NUMERICAL METH
PETE NUMERICAL METH PETE 301
Popular in Course
Popular in Petroleum Engineering
This 84 page Study Guide was uploaded by Lenore Medhurst on Wednesday October 21, 2015. The Study Guide belongs to PETE 301 at Texas A&M University taught by Staff in Fall. Since its upload, it has received 75 views. For similar materials see /class/225888/pete-301-texas-a-m-university in Petroleum Engineering at Texas A&M University.
Reviews for PETE NUMERICAL METH
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/21/15
PETE 301 PETROLEUM ENGINEERING NUMERICAL METHODS LECTURE NOTES AND STUDY GUIDE Peter P Valk August 187 2004 Contents I Introduction 01 Course Objectives 02 Resources 03 Course structure II Tools H N Introduction to Excel and VBA lll Starting Excel saving a le l l l l l l l l l l l l l l l l l l l l l l l l l Workbook worksheet views printing 121 Review Spreadsheet basics l l l l l l l l l l l l l l l l l l l l l l l l l l l l l l Graphing With Excel l l l l l l l l l l l l l l l l l l l l l l l l l l l l l Numerical methods inside Excel l l l l l l l l l l l l l l l l l l l l l l l Database Management With Excel Excel and the WWW Keyboard macros and What they can7t do l l l l l l l l l l l l l l l l l l VBA variable sub function module l l l l l l l l l l l l l l l l l l l l l 110 Structured Programming Subs Decisions and Loops l l l l l l l l l l l l H to HHHHHHH 00N 01gtgtCA3 Mathemat ica III Engineering models and numerical methods 3 g Mathematical Modeling and Engineering Problem Solving Conservation Laws and Engineering Accuracy and Precision Approximations and RoundOff Errors l l l l l Error propagation Existence and uniqueness Convergence criterion Order of the method rate of convergence stability sensitivity l l l l l l Classi cation of problems and metho s 31 3 2 3 3 3 4 3 5 3 6 i i i i i i i i i i i i i i i i i i i i Methods related to single variable problems 41 Roots of equations and extrema of functions l l l l l l l l l l l l l l l l 411 Bisection 412 False position l l l l l l l l l l l l l l l l l l l l l l l l l l l l l Gama 01 G 413 Newton7s method 414 Optimization Golden ratio search for a minimum 42 Fixedpoint iteration as a general framework for iterative processes 43 Representing7 manipulating functions 431 Numerical Differentiation of functions that can be evaluated every where 432 Numerical lntegration 433 Simpsonls 13 rule 44 lnterpolation7 Smoothing7 Curve fit7 Least squares 4 41 Least Squares and its Numerical Aspects 45 Numerical solution of ordinary differential equations 431 Basic methods for solving ODE Methods related to multivariable problems 51 Matrices7 vectors 52 System of linear equations 53 LU Decomposition and Backsubstitution 54 System of nonlinear equations7 NewtonRaphson method 55 Multivariable extrema 56 Solution of system of ordinary differential equations RK4 57 Partial differential equations and reservoir simulation 571 Reservoir Simulation Integral Transforms7 Spectral methods7 inversion of the Laplace trans form 61 Transformations 62 Laplace transform 621 Numerical inversion of the Laplace transform IV Appendix 7 00 Study guide 71 Be able to Assignments7 test problems Error propagation borehole volume 82 Error propagation onegallon cube 83 Error propagation barrel 84 Error propagation hydrocarbon volume Root nding Newton iteration z factor Root nding solution of the PR equation of state Root nding heat balance Root nding friction factor 39 Well deliverability isoterm Optimization Maximizing Net Present Value Integration of a function Numerical integration of discrete data pore volume 85 86 87 88 89 810 811 812 813 814 Straightline formation volume factor model 1 74 815 Straightline formation volume factor model 2 74 816 Straightline gas in place 75 817 Straightline FlowAfter Flow Test IPR 75 818 Nonlinear least squares oil viscosity as a function of pressure and temper ature 76 819 ODE Production rate decline from a differential equation 76 820 ODE pressure versus depth in a shutin well 77 821 Meaning of matrices and vectors chemical component balance 77 822 System of linear equations mixing 77 823 System of linear equations With many righthand sides log interpretation 78 824 Nonlinear System of Equations Chemical Equilibrium of Methan Combus tion 825 Units Darcy s law 80 9 Tables 81 Part I Introduction 01 Course Objectives 0 Learn the underlying principles behind some commonly used numerical methods 7 Understand the route from arithmetic and calculus to numerical methods 7 Understand the basic concepts including uniqueness approximation error con vergence and stabi ityi 0 Learn how to program a numerical method using modern programming concepts 7 Problem solving steps 7 Structured programming 7 Use of subroutines modules 7 Learn how to make use of Visual Basic for Applications 7 Learn how to use an integrated program development system providing debug ging facilities 0 Recognize the main features of a mathematical problem arising in a technical con text 7 Single or multivariable 7 Linear or nonlinear 7 Involves algebraic differential or integral equations 7 Direct or iterative technique should be applied 7 Explicit or implicit scheme is appropriate 0 Develop skills to formulate a problem starting from physical principles 7 Start With conservation principles 7 Include material properties 7 Add problem speci c conditions 7 Learn how to use various resources in solving a complex problem 7 Use of software systems Subroutine libraries modules 7 Internet resources 0 Improve ability to critically analyze numerical results 0 Improve documentation skills 02 Resources The concepts and methods of technical computing are summarized in the excellent book of Chapra and Canale Numerical Methods for Engineers CC In order to make ef cient use of our time we request you to read the material before the lecturer Any similar book can be used however The Visual Basic for Applications Excel program development environment will be available in the lab Whenever you are asked to solve a problem by computer we request an Excel Workbook with Visual Basic source code and results Excel and VBA for Excel has a vast literature We will give you the opportunity to experiment with many of the numerical methods using a Living Textbook developed particularly for this course To m e full use of the Living Textbook you will need a software system called Mathematical This course however does not require you to learn how to program in Mathematical In fact to run any portion of the Living Textbook you will need minimum knowledge which can be acquired in 5 minutes Many of the speci c petroleum engineering problems also will be provided in the form of Living Textbooki You may experiment with them and may use the results for checking your Visual Basic for Applications VBA program results Basic References 1 SC Chapra and Rf P Canale Numerical Methods for Engineers McGraw Hill any edition any year or any similar title 2 Either of the following or any similar title a Larsen Engineering with Excel Prentice Hall 2002 b Liengme A guide to Microsoft Excel 2002 for scientists and engineers Butterworth Heinemann 2002 c Bloch Excel for Engineers and Scientists 2002 3 All of your Petroleum Engineering Books Special thanks go to L Piper Bi Maggard Rf Archer Rf Wattenbarger i 03 Course structure H to 9 F 01 533 Engineering problem solving tools and programming a Engineering problem solving approaches software tools b Excel c VBA d Structured and Modular Programming Basics of numerical methods a Taylor polynomial Numerical errors Error propagation b lteration Convergence Order Stability Sensitivity c Classi cation of problems and methods Methods related to single variable problems a Roots of equations and extrema of functions Bisection False position Newton golden ratio search b Numerical differentiation and integration Finite difference Trapezoid rule impson s ru c Interpolation Smoothing Curve t Least squares d Numerical solution of ordinary differential equations Methods related to multivariable problems a Matrices vectors b System of linear equations c Gauss GaussJordan LU decomposition Special cases lterative methods d Multivariable nonlinear equations NewtonRaphson e Multivariable extrema NelderMead simplex Fully nonlinear least squares Partial differential equations and reservoir simulation a Classi cation and Numerical Solution of PDEs b Transient solution of the diffusivity equation c Reservoir simulation Transformation methods a Spectral methods b Laplace transform inversion Part II Tools Chapter 1 Introduction to Excel and VBA 11 Starting Excel saving a le To create a new workbook On the File menu click New and then click Blank Workbook on the New Workbook task panel To save a workbook On the File menu click Save As In the File name box type a new name for the workbook In the Save as type list click a le format if you want other than ixls In general it is useful to click the right mouse button Most of the time it will offer a menu item you really want 12 Workbook worksheet views printing 121 Review A Microsoft Excel workbook is a le that contains one or more worksheets which you can use to organize various kinds of related information You can enter and edit data on several worksheets simultaneously and perform calculations based on data from more than one worksheeti When you create a chart you can place the chart on the same worksheet as its related data or on a separate chart sheeti To View and scroll independently in different parts of a worksheet you can split a worksheet horizontally and vertically into separate panesi Splitting a worksheet into panes allows you to View different parts of the same worksheet side by side and is useful for example when you want to paste data between different areas of a large worksheeti To see how it will look like in print On the View menu click Page Break Previewi Before you do any work make sure that you macros are allowed security level is low and you can save an excel spreadsheet to your own diski 13 Spreadsheet basics The difference between relative and absolute references is the key concept to understand Relative references A relative cell reference in a formula such as Al is based on the relative position of the cell that contains the formula and the cell the reference refers to If the position of the cell that contains the formula changes the reference is changed If you copy the formula across rows or down columns the reference automatically adjusts By default new formulas use relative referencesr For example if you copy a relative reference in cell B2 to cell B3 it automatically adjusts from Al to A2r Absolute references An absolute cell reference in a formula such as Al always refer to a cell in a speci c location If the position of the cell that contains the formula changes the absolute reference remains the same If you copy the formula across rows or down columns the absolute reference does not adjust By default new formulas use relative references and you need to switch them to absolute references For example if you copy a absolute reference in cell B2 to cell B3 it stays the same in both cells Alr 14 Graphing with Excel Import a text le Click the cell where you want to put the data from the text le To ensure that the external data does not replace existing data make sure that the worksheet has no data below or to the right of the cell you click On the Data menu point to lmport External Data and then click lmport Data In the Files of type box click Text Filesr Check out also On the Data menu Text to Columnsr Creating an XY Scatter Plot You can create a chart on its own sheet or as an embedded object on a worksheet You can also publish a chart on a Web page To create a chart you must rst enter the data for the chart on the worksheet Arrange your data in columns with x values in the rst column and corresponding y values andor bubble size values in adjacent columnsr Then select that data and use the Chart Wizard to step through the process of choosing the chart type XY and the various chart options or use the Chart toolbar to create a basic chart that you can format laterr Best t equations may be calculated for chart data via the Add Trendline dialogr Trendline equations are displayed by checking the correct box under the Options tab Errors bars can be added to a chart by left clicking on a data series and accessing the Format Data Series dialog box Editing an Existing Graph Check out Creating Graphs with Multiple Curves Check out More graphs types Check out 15 Numerical methods inside Excel Goal Seek The Goal Seek dialog box can be used to solve for roots of a single equationr An initial guess is needed Goal Seek performs calculations which improve the solution accuracy Solver The Solver can be used to solve both linear and nonlinear optimization problems ie maximize or minimize an objective function subject to known constraints The Solver Parameters dialog box refers to the following Cell for the objective function value which is to be maximized or minimized Cells for the the values of the indepen dent decision variables lf the values of the independent variables are subject to the constraints then cells constrained 16 Database Management with Excel Check out Search Filter and Sort in the Data menu 17 Excel and the WWW Excel les may be converted to HTML web documents The HTML Wizard guides you through the process 18 Keyboard macros and what they can t do If you perform a task repeatedly in Microsoft Excel you can automate the task with a macro A macro is a series of commands and functions that are stored in a Microsoft Visual Basic module and can be run whenever you need to perform the task For example if you often enter long text strings in cells you can create a macro to format those cells so that the text wraps Recording keyboard macros When you record a macro Excel stores information about each step you take as you perform a series of commands You then run the macro to repeat or 77play back77 the commands If you make a mistake when you record the macro corrections you make are also recorded Visual Basic stores each macro in a new module attached to a workboo Making a macro easy to run You can run a macro by choosing it from a list in the Macro dialog box To make a macro run whenever you click a particular button or press a particular key combination you can assign the macro to a toolbar button a keyboard shortcut or a graphic object on a worksheet Keyboard macros are to replay your keystrokes If you want to program your own algorithms you need VBA programming 19 VBA variable sub function module Visual Basic for Applications The programming language of this course is Visual Basic for Applications VBA which is included with the Excel spreadsheet application This combination of programming language and spreadsheet has many powerful features which make it an excellent tool for problem solving 1 A modern structured programming language with many intrinsic procedures a Mathematics b Data type conversion c String manipulation d File system 2 Ability to control all of Excells features through the Microsoft Excel Objects in cluding a Data management b Graphics c Solver for optimization and root nding 3 Ability to add custom functions which can be used like 77builtin77 Excel functions 4 Capabilities allowing the development of event driven user friendly programs through a variety of interface too s a Dialog and input boxes b Custom user forms with spinnerscheckboxes etci c Custom menus and toolbars 5 The ability to extend the language through object oriented programming with the use of Class Modules Example VBA011 program After Hahn If a stone is thrown vertically upward with an initial speed u its vertical displacement s after a time t has elapsed is given by a simple formula In this formula g is the acceleration due to gravity Air resistance has been ignored We would like to compute the value of s given u and ti Note that we are not concerned here with how to derive the formula rather we want to compute its value The logical preparation of this program is as follows 1 Place the values of g u and t into variables which can be used for calculations 2 Compute the value of s according to the formula 3 Output the value of s which is the result 4 Stop This plan may seem trivial to you and a waste of time writing down Yet you would be surprised how many beginners preferring to rush straight to the computer try to program step 2 before step 1 It is well worth developing the mental discipline of planning your program rstif pen and paper turns you off why not use your word processor You could even enter the plan as comment lines in the programi Option Explicit Const g As Double 981 gravitational acceleration ms 2 Sub VBA011 vertical motion under gravity Dim s As Double vertical displacement m Dim t As Double time 5 Dim 11 As Double initial velocity ms u30 t 6 sutg2t 2 MsgBox quots quot amp CStrs ampquot metersquot s is converted to a string End Sub To make this program run do the following 1 From a new Excel workbook use the Tools Macros Visual Basic Editor menu choice to start the VBA integrated development environment includes editor compiler help and debuggeri 2 In 77Project77 window of VBA right click on your new Excel workbook usually 77VBAProject Bookl and insert a new module 3 Cut and paste the text of VBAOll shown in the above text box into the newly created modulei Note that the 77 Option Explicit77 is not part of the program but is a module optioni Note the colors of comments green language elements blue and names black 4 Left click to place the cursor anywhere inside the subroutine and use the menu to choose Run Run SubUser Formi 5 If you have questions on the language elements shown in blue place the cursor on the word and hit the F1 key Some remarks 1 VBA performs automatic syntax checking and formatting of code 2 The single quote character 7 denotes the beginning of a comment 3 ln PETE 301 every module should require explicit variable typing by using 77Option Explicit77 module option 4 The strange way of declaring g makes it a named constant its value cannot be changed The fact that it is placed outside the subroutine makes it valid in the whole module scope module 5 There is no input All needed data are hard wired into this simple subroutinei Later we will learn how to read from the worksheetsi 6 Output is done using a message box Later we learn how to put results into the worksheetsi 7 This module contains only one subroutine macroi In fact many of them can be placed into one modulei Builtin or intrinsic functions Visual Basic has many builtin functions You will often need the math functions Abs Atn Cos Exp Fix lnt Rnd Sgn Sin Sqr Tani Log is naturallogi If you wish tenbased logarithm use LogxLog10i Example Sin function Returns a Double specifying the sine of an angle Syntax Sinnumber The required number argument is a Double or any valid numeric expression that expresses an angle in radians Remarks The Sin function takes an angle and returns the ratio of two sides of a right triangle The ratio is the length of the side opposite the angle divided by the length of the hypotenuse The result lies in the range 1 to 1 To convert degrees to radians multiply degrees by pilSO To convert radians to degrees multiply radians by lSOpi Those functions are useful but they are not customizable Even if you nd a Visual Basic function that is 77this close77 to what you need you cant worm your way into the innards of Visual Basic to change the way it works You can however create a function of your own A great number of other functions are also available For instance both Excel and Visual Basic have functions that return a random number between 0 and l The Excel function is named Rand and the Visual Basic function is named Rnd You can use the Excel function in a worksheet cell and with some tricks also in your VBA macro but you can use the Visual Basic function only in a macro Data Types The following table shows the supported data types Byte 1 byte 0 to 255 Boolean 2 bytes True or False lnteger 2 bytes 32768 to 32767 Longlong integer 4 bytes 2l47483648 to 2147483647 Single single floatingpoint 4 bytes 3402823E38 to 3402823E38 Double double floatingpoint 8 bytes 1E308 to 18E308 String variablelength 10length 0 to appr 2 billion Variant with numbers 16 bytes see range of a Double Arrays of any data type require 20 bytes of memory plus 4 bytes for each array dimension plus the number of bytes occupied by the data itself Dim Statement Declares variables and allocates storage space Dim NumberUfEmployees As Integer Dim x As Double y as Double When Option Explicit appears in a module as it should in PETESOI you must ex plicitly declare all variables us1ng the D1m Private Public ReDim or Static statements If you attempt to use an undeclared variable name an error occurs at compile time Note Dim ab As Integer will declare b as an lnteger but a will be the default Variant You have to write Dim a As Integer b As Integer if you wish to declare both as integeri Numerics VBA makes a clear distinction between integer and real Integers are used for counting thingsi Real numbers include a fractional parti Real numbers stored as integers are rounded to the nearest integeri Dim As Integer In As Double 4 3 1t n resu 1s 1 n 53 result is 2 n 72 result is 2 n 52 result is 4 m 52 result is 25 You should be aware of rules regarding priorities For instance a frequent error is to rite z pVRT instead of the correct 2 pwIMT Strings To convert any other value to a string use the Cstr intrinsic function Objectsinput and output The structure 77VVithEnd With77 and the 77operator77 point can be understood if you are familiar with the concept of object 77ThisVVorkbook77 is a hierarchical object It has several Workseets a whole collection of them One of them is called 77Sheet277i The following example shows how to read a number into variable 77a77 from cell 77A4 of Sheet2 with ThisWorkbookWorksheetsquotSheet2quot End With This will output date and time into cells 77A4 and 77A577 of 77Sheet2 Sub VBA014 Dim d As String t As String d Date Intrinsic function t Time Intrinsic function Display the result with ThisWorkbookWorksheetsquotSheet2quot Cells1 4 d Cells1 5 t End With End Sub Objects have properties and methods associated with themi For instance when you are formatting a cell you change some of its properties using some methods available Advanced topics Variable scope lifetime visibility The following key words are important Static Public Private The scope of a variable is from where it can be reachedaltered Variables declared in a procedure can be used only within the procedure local variable for the procedure 15 Variables declared on the module level can be used in the entire modul in any of the procedures within the module global variable for the module i Static variables keep their value between two occasions of entering into their scope In other words their lifetime is not limited to one call In general a procedurelevel variable looses its value between two calls unless it is declared Statici Both modulelevel variables and all procedures are Public by default meaning that they are visible even from other modules or even other applications You can make them Private and then they are visible only within the module We will often use module level variables globals to simplify engineering problem solving when some data are needed in various parts of the programi Advanced topics Argument passing by Value or by Reference All arguments are passed to procedures by reference unless you specify otherwise This is efficient because all arguments passed by reference take the same amount of time to pass and the same amount of space 4 bytes within a procedure regardless of the arguments data type You can pass an argument by value if you include the ByVal keyword in the pro cedures declaration Arguments passed by value consume from 2 16 bytes within the procedure depending on the arguments data type Larger data types take slightly longer to pass by value than smaller ones Errors7 Debugging and Testing Compilation errors are errors in syntax and construc tion like spelling mistakes that are picked up by the compiler during compilation the process whereby your program is translated into machine code They are the most fre quent type of error The compiler prints messages which may or may not be helpful when it encounters such an error Generally there are three sorts of compiler errors 4 i Ordinary errorsthe compiler will attempt to continue compilation after one or more of these errors has occurred eigi missing End If statement 53 i Fatal errorsthe compiler will not attempt further compilation after detecting a fatal error eigi Program too complicatedtoo many strings 9 Warningsthese are not strictly errors but are intended to inform you that you ave done something unusual which might cause problems later eigi Expression in If construct is constant or variable has not been assigned a value but used or variable has not been declared if implicit typing is use Some typical errors g 981 comma instead of decimal point x 1 2 t 3 unpaired parentheses Cell215 Wrong function name Cell instead of Cells If ab Missing then x3 Else x4 End If If a program compiles successfully it will run Errors occurring at this stage are called runtime errors and are invariably fatal ie the program crashes i An error message such as 1 Floating point division by zero 2i Over ow 3 Negative argument in the function sqr is generated Over ow is quite common It occurs for example when an attempt is made to compute a real expression which is too large or when you divide by a previously not assigned variable which might be almost zero because it contains some trashi Negative argument to Sqr or nonpositive argument to Log is often the result of a previous logical errori 39 es a program will give numerical answers to a problem which appear inexplicably different from w at we know to e t e correct mathematical solution This can be due to rounding error which results from the nite precision available on the computer eigi two or four bytes per variable instead of an in nite num err Run the following program tion Explicit Sub VBA022 Dim sumx As Double Sllmx Do sumx sumx 0001 DebugPrint sumx If sumx 02 Then Exit Do Loop End Sub Watch the Immediate Windowl You will nd that you need to crash the program to stop eg with ctrl brea i The reason is x never has he value 02 exactly because of rounding error In fact x misses the value of 02 by about 10 i Now try this If abssu1nx 02 lt 1ei6 Then Exit Do Programming is not about avoiding errors from the beginning In fact that would be almost impossible Real programming is about localizing your errors by debugging and correcting them Check your work 1001 times Use test problems Print out intermediate results Investigate limiting cases Message Boxes7 Custom Dialogue Boxes 110 Structured Programming Subs Decisions and Loops Using a Subroutine Subroutines are an important component of 77top down77 program design Subroutines allow the programmer to build and test program procedures that perform a single well de ned task This can make complex programs easier to design test and debug All executable code must be contained in a procedure Procedures can t be nested within other procedures Example Let in your Excel workbook 77Sheet77 contain the following Cells Al and B1 should contain the X and Y coordinates of a point Cells A2 and B2 should contain the X and Y coordinates of a second point This program will read the coordinates of the two points and calculate the slope and intercept of the unique straight line de ned by the two points and output the results in a message box These two subroutines should be contained in a brand new module within your VBA projecti Because the MyLine subroutine has arguments it cannot be executed directly like the previously discussed subroutines by using the VBA Run menui Option Explicit Private Intercept As Double Sub VBA015 Dim X1 As Double X2 As Double Y1 As Double Y2 As Double the two points from the worksheet quotSheetsquot With ThisWorkbook Worksheets quot Sheets quot x1 Cell Calculate the Slope and Intercept Call MyiLineX1 Y1 X2 Y2 Slope Display the result With ThisWorkbook Worksheets quot Sheets quot Cel quot Slope quot Cells3 2 ope Cells4 1 quotInterceptquot Cells4 2 Intercept d With End Sub A subroutine used inside Sub MyLinex As Doubley As Doublexx As Doubleyy As Doublem As Double m yyiy xxix Intercept y i m t x End Sub Note Watch for scope of the various variables lntercept is a module level variable but Slope is not lntercept cannot be seen from other modu esi The IfThenElse Decision Structure xample after Hahn Most banks offer differential interest rates more for the rich less for the poor The Else part may be missing and you may create more complex structure with Else If If balance lt 3000 Then apr 0 0 Else apr 006 Loops the heart of programming math Countable Do the block for i 1 then i 2 etc For i 1 To 10 u C u Next 1 Count down Do it forj 100 thenj 90 etc For i 100 To 50 Step 10 block Next 1 Forever Loop forever unless 77H stops it Do If quotlogicaliexpressionquot Then Exit Do quotblockquot Condition controlled Do while logical expression is true While quotlogicaliexpressionquot quotblock quot Wen 77Select Case structure makes your code readable Use the Select Case statement as an alternative to using Elself in lfwThenwElse statements when comparing one expres sion to several different values While lfwThenwElse statements can evaluate a different expression for each Elself statement the Select Case statement evaluates an expression only once at the top of the control structure Select Case performance Case 1 bonus salary t 01 Case 2 3 bonus salary t 009 Case Is gt 3 bonus 100 Case Else bonu 0 End Select Arrays Arrays are declared the same way as other variables using the Dim Static Private or Public statements The difference between scalar variables those that arenlt arrays and array variables is that you generally must specify the size of the array An array whose size is speci ed is a xedsize arrayl An array whose size can be changed while a program is running is a dynamic arrayl hether an array is indexed from 0 or 1 depends on the setting of the Option Base statement If Option Base 1 is not speci ed all array indexes begin at zero Declaring a xed array In the following line of code a xedsize array is declared as an lnteger array having 10 rows and 10 co umns 19 Option Base 1 Dim MyArray10 10 As Integer The first argument represents the rows the second argument represents the columns As With any other variable declaration unless you specify a data type for the array the data type of the elements in a declared array is Varianti To Write code according to PETESOl conventions explicitly declare your arrays to be of a data type other than Varianti Arrays can also be dimensioned 77onthe fly using the ReDim commandi It is used for dynamic allocation of arrays on the procedure leveli useful related function is UBoundOi It returns the largest available subscript for the indicated dimension of an array Chapter 2 Mathematica Many of the concepts methods and applications of this course are now programmed as 77Living Textbook Using the software Mathematica any section notebook of the 77Living Textbook77 can be read on the screen andor run as a program The particular examples can be modi ed and rerun providing a unique possibility for experimenting with numerical methods and petroleum engineering applications Mathematica is a complete 39 t for 39 and 39 39 in science and technology delivering computing power to over one million people around the world Breakthrough new features such as an innovative typesetting system that can do math now make Mathematica even easier to use Known for delivering quick accurate numeric and symbolic solutions Mathematica is ideal for creating interactive scienti c papers technical reports presentations and courseware that include text active formulas two and threedimensional graphics and customizable buttons and palettes A growing library of application packages provides specialized capabilities in areas such as engineering nance statistics data analysis optics astronomy and fuzzy logic The associated web site is at httpwwwwolframcom While it is very useful and therefore strongly recommended that you should become familiar with the basics of the fully integrated environment for technical computing Math ematica this course is not intended to teach programming in Mathematica In fact to make use of the Living Textbook you need to know only a very few things The basic concept is the notebook A Mathematica notebook is something between a word processortype document and a Fortrantype program When a notebook is running the results of the computations are appearing immediately after the input lines To start a Mathematica notebook in Windows Double click on the notebook To see details A Mathematica notebook consists of a hierarchical structure of cells You can think of the main cells as Chapters Every Chapter contains sections every section contains subsections etc At start usually you will see only the main cell names that is the Chapter titles Then you can open or close any Chapter by clicking the little triangle at the left 0 run the whole notebook Select Kernel Evaluate notebook from the menu If the system asks a question about evaluating initialization cells respond clicking 77yes To see what is happening use the mouse and the scroll bar or the PgUp and PgDn keys 0 make a small change and rerun an example Go to the item you want to change with the mouse cursor Use the Del or Backspace keys and type what you want Once you are ready select Kernel Evaluate Cell For those interested in shortcuts you can use the Enter key on the numeric keypad as well Mathematica notebooks may contain a large amount of graphics Since printing graph ics may overload the system and especially the printer we ask you not to print the note books unless it is necessary Before printing make sure you deleted unnecessary parts Living Textbook Content LTSOZl LTSOSl LTSOS2 LTSOSS LTSO41 LTSO61 LTSO71 LTSO72 LTSOSl LTSOQl LTSOQ2 LTSlOl LTSlO2 LTSlZl LTSlSl Taylor Series Expansion Remainder Term and Truncation Error RootFinding More about iterations Extrema of single variable function Numerical Integration of functions Ordinary Differential Equations Basic Concepts Numerical Solution of the ODE lnitial Value Problem Reservoir MaterialBalance Vectors Matrices Linear system of Equations Direct Solution of Systems of Linear Equation Systems of Linear Equations GaussSiedel lterative Method Solution of System of NonLinear Equations Minimization of a Multi Variable Function OilWater 2D IMPES Spectral Methods Discrete Fourier Transform Part III Engineering models and numerical methods Chapter 3 Mathematical Modeling and Engineering Problem Solving 31 Conservation Laws and Engineering 32 Accuracy and Precision Approximations and Round Off Errors Precision is related to repeatability or to the number of digits obtained Accuracy is related to how near we are to the 77true77 value Obviously accuracy is more important than precision ln analysis of the numerical methods we often start with the Taylor series expansion fzh fzhfzllh2f z2l Several terms on the right hand side will enter the 77formula77 of the method The rst term left out is the leading error term The additional terms are usually neglected in the analysis Truncation error difference between true result for actual input and the approximate result produced by given algorithm using exact arithmetic in general it is estimated by the leading error term Rounding error difference between result produced by given algorithm using exact arith metic and result produced by same algorithm using limited precision arithmetic Computational error is sum of truncation error and rounding error but one of these usu ally dominates Since usually input data are in error and we do a lot of 77steps the most important issue is How does the error propagate Measures of the error in approximating a number with true value I by an approxima tion 5 True Error El z 7 E Relative Error 75 Percentage Relative Error 75 X 100 Of course in general the true value will not be known so these error measurements cannot be carried out rigorously However in many cases upper bounds for the values of these errors can be found which give us enough information to put the result into perspective 33 Error propagation A new perspective The result of a computation is a multi variable function of the input values For simplicity consider two variables The rst order Taylor approximation 8 8 1A1 iAy BI IALyAwfw 8y Rearranging and introducing Af Any Ay 7 we obtain the basic relatiuon of 77 error calculus 5f W Af llEHAer HaiyllAerm In this error calculus AI means always a positive quantity it represents the maximum possible value of the actual error Derive what are the two partial derivatives of fzy fziy fzw fIy From there we obtain that for summationsubstraction the absolute errors are added For multiplication and division the relative errors are added Example Find the error in a laboratory determination of the gas zfactor given that p 3245 psi Ap 3p8i meaning i3 psi 1 1977 fig Av 0003 fig n ll 7 mole An 0001 lb 7 mole T 7396701 AT 003 0R Differentiating 2 pv nRT we obtain AZ llnETllAPJr ll llAH ll llAnJr NW HAT A simpler form is A2 HzH 2 810 3 Does the result have any dimension Why did we not use more decimal gures Why do we use absolute value Why do we not use absolute value in the p term Does it matter that T is in the denominator Note lb mole is a rather obsolete unit it is 453592 mol or sometimes the mass or weight of it R 1073 psi ftSlb 7 mole 0R is the universal gas constant in the English or eld system of units Error propagation in a numerical method What happens with the error inherited The stability of a computational algorithm errors originated from truncation and round off are not ampli ed from step to step but are rather kept under control Illconditioned problem A problem is illconditioned if results are very sensitive to input data regardless of the method used 34 Existence and uniqueness Convergence criterion Review the following concepts Existence and uniqueness of solution Graphical versus numerical solution Analytical versus numerical method to get the results numerically Direct vs iterative explicit vs implicit Convergence Mathematically a series converges to a limit if the difference of the nth term and the limit can be made smaller than any selected positive epsilon by selecting a large enough n Typically in many numerical problems a sequence of approximate values will be ob tained which hopefully converge towards the desired solution zoxlxgzgz4 iiii If they appear to be converging then a measure something like the percentage relative error is used 5a X 100 i1 If however the I can be nearzero an absolute error criterion is better Why If this error approximation 5a drops below some userspeci ed tolerance then the process is stopped The results should be checked if possible in some other independent way for the required accuracy 35 Order of the method rate of convergence stabil ity sensitivity Since most methods are derived from some Taylor expansion considerations the theoret ical order of a method is related to the truncation of the Taylor series Examples will be given in the numerical differentiation sectioni Broadly speaking the order of the method is the exponent of the variable stepsize in the leading error termi Total error truncation and roundoff Error propagation in a numerical method is as important as the magnitude of the error Practical rate of convergence is often determined by numerical experimentation For this purpose we need simple examples with known exact solutions If we plot estimates of errors on a loglog plot we will see that the error is decreasing with stepsize and the relation appears to be a straight line lf for instance the for every log cycle of the stepsize the error changes two log cycles the rate is quadratici Often instead of the effect of stepsize we are more interested in the effect of number of iterations TradeOff Regarding 7 stepsize The point is that smaller stepsize provides larger accuracy but after a while the increased number of calculations bring more roundoff errorsi Therefore there exists always an optimumi TradeOff regarding complexity is very similar It is better to use a higher order method then a lower order method up to a certain point Too much complexity however may back re the same way as to many small steps An algorithm is stable if errors originated from truncation and roundoff are not ampli ed from step to step but are rather kept under control The opposite is unstable A problem is wellconditioned if results are not very sensitive to input data at least if an appropriate method is used The opposite is called illconditione 36 Classi cation of problems and methods Single variable multi variable our main classi cation Other can be linear vsi nonlinear direct vsi iterative explicit vsi impliciti 0 Single variable methods Root of a single equation Minimum of function Manipula tion of functions differentiation integration functions given in form of expression or algorithm discrete data point smoothing curve tting diff and int Ordinary Differential Equation ODE 0 Multi variable problems Linear Algebra Matrices vectors Systems of linear equa tions direct special iterative Nonlinear System of nonlinear equations Minimum of multi variable function general and least squares System of ODE Partial Differen tial Equations Reservoir simulation 0 Other Transformation methods Chapter 4 Methods related to single variable problems 41 Roots of equations and extrema of functions Single and multiple roots bracketing versus open methods Nonlinear equation may have multiple root where both function and derivative are zero 411 Bisection In the bisection method first the user estimates approximately the position of the root providing a lower and upperbound bracketing the root If the interval a b is a bracket then fafb S 0 A first estimate of the root is then computed as the midpoint between the two bounds z a b2 We can improve on this estimate by determining which side of X the root lies and then moving the bracket accordingly To do this we have to overwrite either a or b by the current midpoint 1 depending on which previous function value has the same sign as At this point we have a new bracket and can continue with the next iteration Each iteration halves bisects the bracket and therefore halves the maximum possible error hence the term 77bisection We can calculate apriori the necessary number of steps to reach a certain fraction of the initial bracket 412 False position The search starts with an interval ab It is assumed that the function changes sign in the interval fafb S 0 as in the case of the bisection method Writing the equation of the line passing through the two points a u and b b 7 a yew ff Web a and requiring that y should be zero we obtain for the next approximation of the root 7 f b C b fbfa bia It is easily seen that c is inside the interval a In Now we calculate and overwrite either a or b with 5 depending on which function value has the same sign as After one step again we have a bracket but it is smaller than the initial one Note 7 f a C a fafb ail is the same During the process the size of the bracket does not tend to zero The criterion to stop the iteration is the following the deviation of the new location 5 from any of the old locations a or I should be less than a speci ed tolerance 413 Newt0n7s method The method also called NewtonRaphson method may be used to solve a general equa tion 0 provided the lefthand side can be differentiated 77analytically If an initial guess is known we can calculate the function and its derivative at 10 Writing the equation of the tangent line 9 1610 fIOWEO 1 Substituting y 0 we nd that the line crosses the z axis at lo HID This is the iteration formula of the Newton7s method The obtained 11 can be then renamed to 10 and hence we can continue the iterations 1IO Example Suppose that 13 z 7 3 Then fz 312 1 The following program starts the iteration with z 2 It uses two internal functions F for and DF for its derivative f It stops either when the absolute value of the step is less than lE 6 or after 20 iterations Note that there are two conditions that will stop the While loop either 77convergence or the completion of 20 iterations Otherwise the program could run inde nitely Option Explic it Globals Group 1 Coefficients Dim c3 As Double c2 As Double c1 As Double c0 As Double Globals Group 2 Iteration control parameters Dim MaxIts As Integer Eps As Double x0 As Double leftihandiside function Function Fx As Double As Double K c1 x c0 End Funct i on derivative Function DFx As Double As Double F3c3x 22c2xc1 End Function Main sub Sub Newtonx As Double fx As Double Converged As Boolean Dim Its As Integer d As Double Initialize Its 0 Converged False x x0 Loo While Not Converged And Its lt MaxIts d 7Fx DFx X X d Its It s Converged 1 Absd lt Eps End Sub Driver program Sub VBA041 Values Group 1 Coefficients c31 c20 c11 c03 Values Group 2 Iteration control MaxIts 0 Eps 0000001 x0 Declaration of Group 339 Arguments Dim x As Double fx As Double Converged As Boolean Call Newtonx fx Converged With Worksheets1 CellS2 1 quotfxquot CellS2 2 fx End With If Converged Then MsgBox quotConvergence criterion satisfied quot Else MsgBox quotConvergence criterion NUT satisfied quot nd If End Sub To solve another problem you do not have to rewrite the Whole program It is enough to Change the functions F and DF7 and some of the input Failures in the Newton Method Though it has quadratic convergence there are likely to be problems in some cases One problem can occur if f is very large in which case the graph y f is curving very rapid y A second problem occurs if f is very small which can happen if the function is almost parallel to the 1 axis near the root or if there is a second root very close by in which case 0 somewhere between the roots or if the graph of y f just touches the z axis at the root but does not cross it Multiple rootl It is often observed that the method converges from a 77near enough77 starting point but diverges if the starting point is selected without carer 414 Optimization Golden ratio search for a minimum For minimizing a function of one variable we need a bracket for the solution analogous to sign change for nonlinear equation A realvalued function is unimodal on interval a b if there is unique location 10 such that fis strictly decreasing for z lt 10 that is on the left and strictly increasing for z gt 10 that is on the right Suppose fis unimodal on a b and let 11 and 12 be two points within interval with 11 lt 12 i Evaluating and comparing and fzg we can discard either 1212 or a 11 with minimum secured to lie in the remaining subinterv i To repeat the iteration we need to carry out only textbfone new function evaluation To reduce the length of interval by a xed fraction at each iteration each new pair of points must have the same relationship with respect to the new interval that the previous pair had with respect to previous interval To accomplish this choose the relative positions of the two points to satisfy the golden ratio 739 7 l2 m 0618 Whichever subinterval is retained its length will be 0618 relative to previous interval and one interior point retained will be at position either at 0618 or 0382 42 Fixedpoint iteration as a general framework for iterative processes Many iterative methods for solving nonlinear equations can be viewed as Ik1 5 y where 91 is constructed from the function If you understand the above formalism you can easily create the 77g function77 for example for the Newton iteration gz z 7 awHz The scheme is called xedpoint iteration or direct iterationl There are four basic cases monotone convergence 0 oscillatory convergence o monotone divergence o oscillatory divergence Almost all iterative processes manifest one of the above behavior types even if the derivation of the method has nothing to do with direct iteration 43 Representing manipulating functions 431 Numerical Differentiation of functions that can be evalu ated everywhere When an analytical solution to the first derivative of a given function is difficult or in convenient then a numerical method can be used to provide an approximate though often highly accurate computation Many methods exist with increasing complexity they perform in general computations to a higher accuracy A basic understanding of the meaning of Truncation Error and Round off Error is crucial Intuitive error analysis of numerical differentiation Take the function Calculate the derivative at x 7r 3 using the simplest forward divided difference Use step size lE 5 lE lO lE 15 lE 20 until your calculator stops giving reasonable result Use radian switch if you need Forward Difference Approximation rst derivative The Forwarddifference approximation method for the numerical differentiation of a func tion can be derived by considering differentiation from first principles or by considering Taylor7s expansion The differentiation from first principles derivation is intuitive 1 fzh 7 M f I g h As a computer cannot divide by zero the computed nite version of this expression is f I h i f 1 Hz This is the Forwarddifference approximation for the numerical derivative of a function it has the most basic form for a numerical derivative and is the least accurate The derivation through Taylorls expansion provides more information on the method frhfzf zf rm Rearrange and represent all remaining terms in one remainder term fz flt1hL g I fZ We do not know z and hence neglect the term and hence obtain the forward difference approximation previously derived With this derivation however we see clearly that the error is approximately ft2 f z ie proportional to h This method is of first order 33 A plot of the error on a loglog scale would result in a 450 straight line To 77minimize77 the error choose a small value of h but h should not be too small as 7roundoff errors7 in the machine arithmetic increase as h decreases Backward difference approximation of the rst derivative is very similar f1 fI h 4 Hz g f z Again we neglect the term containing 2 Therefore f z 7 f z 7 h m w lt gt hlt gt and the method is of rst order The centraldifference approximation of the rst derivative gives reasonably accurate results and is often implemented 1 fltzhgt emeh 2h Interestingly it does not make use of the function value at the base point at all It is a second order method and the loglog plot of the error on an example with known exact result would give a straight line with slope two The centraldifference approximation of the second derivative is 1mm m f1 h 7 2 21f1h This is a second order method 432 Numerical Integration Review geometric meaning Review physical engineering meaning Trapezoid rule The aim is to calculate numerically the integral of a function from a to b The basic idea is to evaluate the function at speci c locations between the limits summing the function evaluations will give an approximation to the integral area under the curve There are many formulae for numerical integration also known as quadrature with increasing complexity these formulae provide a greater accuracy The 7Trapezoidal Rule7 is the simplest Consider integrating a known function over the interval a b The Trapezoidal Rule gives the following expression for the exact integral 1 If h f fb hg f 7 2 a 12 2 where h is the interval a b f z is the second derivative of the function evaluated at some unknown point 2 between a and b The value of f z is generally not known 2 is unknown though bounded between a and b and so this term is omitted from the 34 solution when we perform the numerical calculation This leads in general to an error a 7truncation7 error in the numerical evaluation of the integral The accuracy of the calculated result can be increased by multiple application We will cut the interval ab into n 77panels 1 h fa fz1 fz2 frn1 1 where I0a zlah zga2h In71anilh Inb hb7an 433 Simpson7s 13 rule The simple Simpsonls rule uses a parabolic approximation to the curve The multiple application of Simpson s rule is usually the rst choice of engineers because of its reliability and simplicity h I W g u 4f11 2f12 41617171 f0 where n is even and h b 7 an General idea of Gauss quadratures A quadrature rule is weighted sum of nite num ber of sample values of the integrand function Both the trapezoid and Simpsonls rules are in this general class Closed formulas Containing left and right points Open for mulas NOT containing left and right points Advantage of open formulas eg improper integrals Richardson extrapolation applied to numerical integration An error estimate can be formulated by making use of the fact that the error has a known order The plan is the following 0 Do calculation with h and with h2 0 Express remainder term as a power of h o The difference between the two results can be used to improve our nal approxima tion Try to extrapolate to 77infinitely small77 h Example repeated Simpsons rule has a leading error term of 4th order Assume we did the calculations with h and with h2 Therefore we have two computed values 11 and 12 However what we really want to know is I that is still unknown 11 1 Ch4 12 IC In the above equations C is also unknown To eliminate it multiply by 24 and fl 2412 241 Ch4 7 1 I Ch4 Then subtract 2412 7 11 24 71 Finally rearrange and obtain estimate of the unknown 1 I w 1612 7 11 15 Note that we do not really need C and we have not calculated it 44 Interpolation Smoothing Curve t Least squares Review Purposes for lnterpolation Plotting smooth curve through discrete data points Reading between lines of table etc If the data is noise corrupted is of nite accuracy we rather use approximation smooth ing instead of strict interpolation Interpolating by simple function families 0 Polynomials o Piecewise polynomials splines o Trigonometric functions will be considered among 77spectral77 methods 0 Exponential functions 0 Rational functions An interpolating polynomial is one which goes exactly through a set of data points In general a set of n 1 data points can be interpolated by a polynomial of degree n Such interpolation leaves no degree of freedom the interpolating polynomial is unique Unfortunately it is only useful for fairly small sets of data points An interpolating polynomial of large degree is required for a large set of data points and such high degree polynomials tend to have many oscillations and will approximate the data set very badly Cubic spline is piecewise 3rd order polynomial that is twice differentiable It is used extensively to interpolate a large number of data points Rational function is a fraction both the numerator and the denominator are polynomials 36 Differentiation and integration of noisecorrupted data 0 Differentiation is more sensitive to error in input data You might want to use quite large step size why 0 Integration is 77smoothing out Stepsize is the smaller the better in almost all cases Approximation by leastsquares The meaning of the criterion Sum of Squares as a measure of 77how good the t is Other requirements might include o Smoothness 0 Least number of parameters 0 Extrapolating power etc Two basic cases 0 Model is based on 77 rst principles77 0 Model is just a convenient vehicle 441 Least Squares and its Numerical Aspects When we want smooth approximation we need some degrees of freedom For instance to tting a straight line 2 parameters to ten points leaves us with 8 degrees of freedom There are in nitely many straight lines and we need a criterion to select a unique one The least squares criterion is the most popular one It minimizes the sum of squared deviations between observations and values calculated from the model The LSQ straight line t is also called linear regression line A word of caution Least squares methods are very adversely affected by 77 bad data For example ifa data set is basically linear but one faulty data point is a long way outside of the linear trend of the rest of the data then this one point can drastically change the regression line This happens because we square the deviation of the data point from t e ine Hence one point with a extra large deviation contributes a huge amount to the objective function The regression line will therefore be moved a long way towards the one deviant point and away from the main linear trend in order to minimize the deviation of this faulty point This is a big factor in many engineering applications and special methods are used to try to detect faulty data points and eliminate them Or to use other possible measures of the 77 goodness of t77 such as sum of absolute deviations Straight line t linear regression Suppose that some the data points should t a straight line y mm b if it were not for experimental andAmeasurement errors The problem is to nd the straight line characterized by 771 and b which best ts the data points in the least squares sense ymzb 37 Observations 11791 12792m1n7yn The derivation of formulae for the slope and intercept should be familiar by nowi The nal formulae are n n n n 2211 7 2213 Useful subroutines for straightline t Option Explicit Option Base 1 Sub StraightLineFitx As Double y As Double In As Double b As Double Dim n As Integer n UBoundX m n scalprodx y s sumx sumy n scalprodx x 7 sumx 2 b sumy s m sumx n End Sub Function sumx As Double As Double Dim n As Integer i As Integer n UBoundX 0 For i1 to n sumsum xi Next End Funct i on Function scalprodx As Double y As Double As Double Dim n As Integer i As Integer n UBoundX scalprodscalprod xi yi Next End Funct i on Transformations to straight line form Models of physical phenomena are in few cases in straightline formi Many 2parameter nonlinear models can be7 however7 easily transformed into straightline formi In the loglog paper does the same for a powerlaw model For other nonlinear models you should be able to gure out the transformation and make the correspondence between the straight line parameters slope and intercept and the real model parameters We will see many examples Excel services trendline option on graphs with the Show Trendline Option 77slope77 and 77intercept77 spreadsheet functions linear regression in data analysis packager 45 Numerical solution of ordinary differential equa tions Meaning Ordinary differential equation ODE all derivatives are with respect to single independent variable often representing timer Equations with higher derivatives can be transformed into a system of rst order equations The ODE f ft y is nonautonomous The ODE f is autonomousi Solution embeds into the directional eld The ODE f fty does not by itself determine a unique solution To single out a particular solution one must specify a value yo of the solution function at some point to This is called initial value Order of ODE determined by highestorder of derivative involvedi Higher order ODEs may be narrowed down77 also by using boundary values 451 Basic methods for solving ODE Euler Eulerls method uses the forward difference approximation for and can be sum marized by the marching formu a Ii 1121 h 9139 yiil h flxiilyyiill RK4 This method for solving y f 1 starting at 10 yo is one of the most popular ones because it is quite accurate and stable The formulae are Ii 1121 h 9139 9121 g 39 k1 27 21 k4 where 64 flxiil h M71 h 39 k3 Note that the 16 values are various approximations to the slope of the unknown function Heun Heun7s method improves on Eulerls method by a 77better77 choice of the slope of the line followed from the point In yn to the next point zn1yn1i Instead of choosing this slope as dygquot Imyn as in the Euler method choose it as the mean between this slope and the slope at the destination point reached by the Euler methodi Simple Heun Ii M71 h y yFl h 39 f1i717 will M71 39f1iilyyiill flIz39yygl 39 Iterated Heun simplest predictorcorrector 71 M71 fl1i717 M71 flIz39yyi continue until 77convergence77i This is a predictor corrector implicit iterative method Needs a stopping criterion for inner iterationi It is a onestep second order method and is stable In practice we use higher order variantsi Singlemulti step7 predictorcorrector7 explicitimplicit So far we7ve considered selfstarting one step methodsi They use only yi Multiplestep methods use more than one previously calculated value that is to use yi and yFl etc The traditional predictorcorrector solution methods use one equation to predict the value of yn1 using the previously computed yn and other previously computed yi values This is followed by a second equation which is used to correct the value given by the rst equation They provide some improvement especially for stiff problems but need a jumpstarter for instance RK4 In choosing step size for advancing numerical solution of ODE we want to take large steps to reduce computational cost but must also take into account both stability and accuracy Adaptive step size is adapted automatically from error estimate A method is called implicit if some kind of iteration is involved Predictorcorrector methods are impliciti A 77stifl 7 ODE corresponds to physical process whose components have disparate time scales or whose time scale is small compared to interval over the interval which it is studied on Chapter 5 Methods related to multivariable problems 51 Matrices vectors What s an Array What s a Matrix Whats the Difference An array is a collection of objects each occupying a position relative to the other objects in the collection An array of antennas for example usually refers to a group of antennas laidout in either a straight line or in a planet The objects of the array can be real physical objects like antennas or spark plugs or intangible objects like numbers When an engineer talks about a matrix she is generally talking about a numerical arrayi There is the further implication that the laws of Matrix Algebra dictate how an array can be combined With andor operated on by other numerical arraysi Although the matrices encountered by engineering students are usually onedimensional numbers positioned in a straight line or twodimensional numbers positioned in a plane they can be higherdimensional as we A subroutine to get the scalar product of two vectors Function scalprodx As Double y As Double As Double Dim n As Integer i As nteger n scalprod scalprod xi yi ex End Function A subroutine to multiply two compatible matrices Sub matmulta As Double b As Double c As Double multiplies matrices a anb o e Dim row As Integer cola As Integer Dim i As Int Dim s As Don C rowb As Integer col As Integer eger j As Integer k As Integer le cola UBounda 2 rowb r UBoundb 1 col 7 UBo d If N t cola rowb Then MsgBox quotmatprod a and b are not compatiblequot Stop End If If UBound c 1 lt row 0r UBoundc 2 lt col Then MsgBox quotmatprod c is not large enough to contain the resul quot 1 To row For j 1 To col 0 cola s s ai k t bk j End Sub A subroutine to multiply matrix A by a column vector b Sub matvecproda As Double b As Double CO As Double Dim row As Integer col As Integer Dim i As Integer j A Dim s As Double row UBounda 1 col UBounda 2 If row lt Uboundb 0R row lt Uboundc Then s Integer MsgBox quotmatvecprod Vector b or c is not compatible with matrix aquot Stop d I For i 1 To row 5 For J 1 To col s s ai j t bj Next j Ci s Next i End Sub Matrix is lower triangular if all entries above main diagonal are zero aij 0 for i lt j Matrix is upper triangular if all entries below main diagonal are zero aij 0 for i gt ji Review linear dependence and rank 52 System of linear equations Given m X n matrix A and mvector b and unknown nvector x satisfying Ax b System of equations asks Can b be expressed as linear combination of columns of A if so coef cients of linear combination given by components of solution vector x Example 1 Consider the system 2173y8 4175yl6 The determinant is nonzero Two lines crossing at one point Example 2 Consider the system 2173y8 4176yl6 Thus the determinant is 2 6 7 3 4 0 And indeed the second equation is 2 times the rst Geometrically this corresponds to two coincident lines in terms of rank the coef cient matrix has rank 1 the augmented matrix has rank 1 Example 3 Consider the system 2173y8 4zi6yl7 Thus the determinant is 2 6 7 3 4 0 There is no solution Geometrically this corresponds to two parallel lines never crossing in terms of rank the coef cient matrix has rank 1 the augmented matrix has rank 2 Existence and Uniqueness Statement with Rank The rank of a matrix tells us how many independent rows or how many independent columns the coef cient matrix has 7 it may be less than the total number of equations lf rank of the coef cient matrix is less than the number of equations two cases can happen 0 Consistent if rank of coeff matrix rank of augmented matrix in nite number of solutions 0 NonConsistent if rank of coeff matrix lt rank of augmented matrix No solution We can we express the solution with the inverse of the coef cient matrix as x A lb but it is usually not a good idea to get the solution via the inverse The reason for this is simple to calculate the inverse of the matrix A is more computational work than to solve the original system of equations Naive Gauss elimination A system of linear equations can be represented in matrix forms of various types eg zl2zg 7313 5 1 2 73 11 5 211 7 312 13 3 becomes 2 73 1 12 3 711127213 0 71 1 72 13 0 To solve the system a method algoiithm called Gaussian Elimination is used starting with the augmented matrix containing the coefficient matrix on the left with one extra column containing the right hand side values a similar method is used to nd A 1i The augmented matrix is changed by operations which do not change the solutions of the associated system of equations The allowed changes are the three elementary row operations a lnterchanging two rows b multiplying or dividing a row by a nonzero constant c Replacing a row by the sum of that row and a multiple of another row adding a negative multiple of a row 7 eg adding 72 X row 2 to row 1 is the same as subtracting the positive multiple subtracting 2gtlt row 2 from row 1 7 so that this row operation also includes subtractioni The objective is to change the coefficient matrix to a special form of matrix called upper t7iangular in which all entries in the bottom left corner below the diagonal are zero Sometimes the goal is to create an echelon form which is upper triangular in which the first nonzero in each row is a 11 The equations represented by this upper triangular form have exactly the same solu tions but are much easier to solve eg Solving the example above a label like Pt2 R2 72 R1 means replace row 2 by the sum of row 2 p us 72 times row 1 At each stage one entry is a special entry called the pivot A pivot shown in bold face in the next example is the value in the pivot Tow which is used to modify the other rows and is in the column where zeros are being create 1 1 2 7315 R R 72R1 a 2 73 1 3 2 2 71 1 72 0 R R3R1 1 2 73 5 0 77 7 77 0 0 72 2 This has created an upper triangular matrix the first three columns To further simplify this divide the second row by 77 and the third row by 72 creating an echelon 44 form The corresponding equations can easily be solved by starting with the last equation which gives the 13 value then nding 12 from the next equation up and so on back substitution 1 2 11 212 7 313 5 11 5 313 7 212 2 Solve third 0 1 71 1 12 7 13 1 12 1 13 0 Solve second 0 0 13 71 13 71 Solve rst This version is called naive elimination because we assumed that we can always divide by the diagonal element in the next row eg 7 in iteration No two Naive GaussJordan elimination In this version of the elimination we reduce the system to an equivalent system with the identity matrix and hence the solution can be read off from the last column Option Explicit Option Base 1 Sub GaussJordanA As Double Augmented matrix Dim Piv lt As Double TarElt As Double Dim n As In teger Dim PivRow As Integ ber of equations er TarRow As Integer pivot row target row Dim i As Integer j ntege n UBound 1 For PivRow 1 To n process every row PivElt APivRow PivRow choose pivot element If PiVElt 0 Then MsgBox quotZero pivot element encounteredquot Stop End If or 1 To n 1 APivRow j APivRow j PivElt divide whole row ex For TarRow 1 To n now replace all other rows If Not TarRow PivRow Then TarElt ATarRow PivRow For j 1 To n 1 ATarRow j ATarRow j APivRow j t TarElt Next j End If Next T Next PivRow End Sub Answer the questions Where is the coef cient matrix at input Where is the 77right hand side At output Where is the solution What happened to the coef cient matrix For both the naive Gauss and GaussJordan methods in some cases the diagonal entry is simply zero and so it is not an option to use it In general solving linear systems of equations it is necessary to use different pivot choices than the diagonal values used in the rst example above 53 LU Decomposition and Backsubstitution An advanced form of the Gauss elimination is called LU decomposition substitutioni The procedure is divided into two parts the rst one operates on the coef cient matrix and does the computationally demanding half of the job the second one called substitution uses the results of the rst part to obtain a solution to a given righthand side Note that the decomposition destroys the original coef cient matrix overwrites it by the L and U matrix elements An additional vector lndx is used to remember the row changes found necessary during 77partial pivoting The lndx vector then enters the substitution routinei Note that it is declared private but global for the modul containing the LUDecompose and LUSubstitute routines If in spite of partial pivoting the decomposition can not be done7 then the matrix is singular and the procedure stops With a message In other words7 partial pivoting is enough7 full pivoting does usually not improve the method The solution comes out from the substitution routine as x Option Explicit Option Base 1 Sub VBAOQZO LU decompositiontsubstitution driver im n As Integer With WorksheetsquotSheet1quot Cells1 1 quotNumber of equations nquot End With With WorksheetsquotSheet1quot End with Dim i As Integer j As Integer ReDim An n As Double bn As Double xn As Double ReDim Indxn As Integer With WorksheetsquotSheet1quot or n For j 1 To n Ai j Cellsi 1 j Nextj bi Cellsi 1 n 1 Next 1 End With Call LUDecomposeA Indx Call LUSubstituteA Indx b x With WorksheetsquotSheet1quot Cells1 h 3 quotSolutionquot For i 1 To n Cellsi 1 h s xi Next i End With End Sub It is a good practice to put the LUDecomposition and LUSubstitute into a separate module Option Explic it Option Base 1 Sub LUDecomposeA As Double Indx As Integer input Ann output Ann and indxn Const Tiny As Double 1E16 Dim n As Integer n UBoundA 1 Dim i As Int Dim dum As D For k 1 To n i 1 eger j As Integer k As Integer m As Integer ouble m k For i 1 T0 n If AbSAi k gt AbsAm k Then m i Next 1 Indxk m dum Am k End If If Absdum gt Tiny Then 1 dum Else MsgBox quotSingular coefficient matrixquot top S End If For 1 1 To n M1 k 41 k dum ext 1 For j k 1 To n dumi m j Am J Ak J AkJ dum Forik1Ton Ai jAij Aik tdum xti If AbSAn n lt Tiny Then MsgBox quotSingular coefficient matrixquot Stop End If End Sub The LUDecomposition overwrites the coef cient matrix by the decomposed LU formi You have to be careful not to call it more than once Notice that the backsubstitution can be called as many times as many righthand sides we have as far as the coef cient matrix of the system is the same Sub LUSubstituteA As Double Indx As Integer b As Double x As Double input Ann containing the LU decomposition output Xn solution Dim n As Integer n UBoundA 1 Dim i As Integer j As Integer k As Integer Dim dum As Double For i r 1 To xi bi Next i For k 1 To n i 1 1 Indxk um xi x1 xk xk dum For J k 1 To n xj xj Aj k t dum Next Next k For k n To 1 Step 1 xk xk Ak k For j 1 To k i 1 xj xj 7 Aj k t xk End Sub Sparse systems Tridiagonal Thomas algorithm Sparse methods can be applied to all kinds matrices of special form and even to a general sparse matrix of no particular type but the methods are quite complicated in general A special type of matrix called tridiagonal in which the only non zeros are on the main diagonal and immediately below it and a ove it ie on the diagonals immediately below and above the main diagonal a very effective algorithm is due to Thomas eg Given the tridiagonal matrix shown store the three diagonals in a three column array shown 21000 12100 a01212 A02400gtb22434 00132 511020 00024 The storage saving in the above example is only from 25 down to 151 However in larger matrices the reduction is very signi cant For an n X n matrix with n2 entries the storage will be reduced to 37L Hence if n 1000 for a tridiagonal matrix then the full storage is 1 000 000 and the reduced storage is 30001 The algorithm not only saves storage but it enormously reduces the number of computations This is due to the fact that we skip multiplications and additions because we know that the result will be zero T e Thomas algorithm is extensively used in reservoir simulation 48 Option Explicit 0 tion Base 1 Sub Thomasa As Double b As Tridiagonal system of equations a is subdiagonal b is diagonal c is super diagonal a1 and Cn are not used d is the right hand side x is the result Dim n As Integer i As Integer n 7 UBoundb ReDim wn As Double gn As Double Double c As Double d As Double x As Double w1 b 1 EU dl W1 For i To n Next xn 7 g For 1 n 7 7 1 To 1 Step 71 xi 7 gi 7 61 xi 1 wi i Next End Sub Note that matrix operations and solution of linear systems are also provided by excel as spreadsheet functions7 macrosi Direct vs Iterative Gauss elimination7 GaussJordan quot 39 39 LU J r quot39 are DIRECT methods In direct methods we can predict the number of multiplications and additions if the number of equations n is known lterative methods are generalization of the direct iteration we learned for one variable I 91 equationsi lterative methods are used for very large but sparse systems reservoir simulationi You can not predict the number of operations for an indirect method and you need a starting guess and a convergence criterion The iterative method we are going to learn is the GaussSeidel method The method is illustrated for a general 3 X 3 systemi The strategy is to express 11 from the rst row7 12 from the second row7 etc Then start with a guess and calculate a new 11 from the rst equation Then with the best available approximants calculate 12 and so forth In other words7 as soon as a variable is updated in one equation its new value is used on the right hand side of the next equation a1111 a1212 a1313 b1 a2111 a2212 a2313 b2 a3111 a3212 3313 113 11111 b1 7 11212 7 aiszs 12212 b2 7 12111 7 12313 assets b3 7 13111 7 13212 k l k k a1111 b1 7 a1212 7 1313 k 1 1 k a2212 b2 7 a2111 7 2313 k 1 k 1 k 1 331ng 113 7 3190ler 329 r It works ne on diagonally dominant matricesi Diagonally dominant matrices the absolute value of the diagonal element is greater than the sum of the absolute values of all the other elements in a given rowi ldentity matrix is diagonally dominant 54 System of nonlinear equations NewtonRaphson method The method may be used to solve a general equation fx 0 provided the lefthand side can be differentiated 77analytically Note that boldface denotes matrix and vector We have n equations to solve and 11 components of the unknown to nd If an initial guess is known and at that location we can calculate the function fx0 and its partial derivatives collected into the Jacobian matrix Jx0 The rst element in the rst row of the Jacobian matrix will contain the partial derivative of f1 with respect to 11 the second element in the rst row will contain the partial derivative of f1 with respect to 12 and so on Writing the equation of the tangent plane fx0 Jx0x 7 x0 0 Substituting y 0 we nd that the line crosses the x axis at x1 where x1 7 x0 Ax Jx01 fx0 This is the iteration formula of the Newton7s method The obtained 11 can be then renamed to 10 and hence we can continue the iterations In practice we do not use the matrix inverse we just solve Jx0Ax 7fx0 for Ax using LU Decomposition Cautious Approach We calculate the Ax and multiply it by a less than one a if we encounter convergence problems Alpha can be adjusted during the iteration process Questions How many subroutines does the user have to provide for a NewtonRaphson method Two subroutines one to evaluate f and the other for the Jacobian matrix J NR Pros and Contras Second order convergence Needs 77analytical77 Jacobian Convergence as in the case of single variable Good if you start from the vicinity of the solution but might diverge Basically it is included in most sophisticated engineering codes Other issues Convergence criterion Sometimes it is formulated in terms of function values 77errors77 That is dangerous because often the functions are math constructions without clear physical meaning and the unit is blurred In practice x always has a well de ned unit and 77eps77 has to correspond to it If the elements of x are of different dimensions or vary several orders of magnitude use normalization NR variants If the Jacobian is not available analytically we can estimate it by nite difference or Quasi Newton 1 We can continuously improve its estimate with every new function evaluation or Quasi Newton 2 Even better continuously improve the estimate of the lNVERSE OF THE JACOBIAN with every new function evaluation 50 NewtonRaphson method example Solve f1 11 3 3 f2 7 2 We Will use a VBA programi The input is from the spreadsheet o kmax from Cells1 2 maximum number of iterations 20 0 eps from Cells2 2 tolerance LOB6 0 alpha from Cells3 2 a 1 0 starting point 11 and 12 from Cells4 2 and Cells5 2 eg 1 and l The program contains a main driving subroutine and the two user subroutinesi Note that it also needs some auxiliary subroutines linear solver LU decomposition and backsubstitution presented previously The linear solver pair can be placed into another modulei Option Explicit tion Base 1 Sub VBA101 NewtontRaphson Method Example Dim n As Integer n 2 ReDim Xn As Double fn As Double Jn n As Double Do u e Dim alpha As Double eps As Double 5 As Double er Dim k As Integer kmax As Integer i As Integer With WorksheetsquotSheet1quot ell quot quot End With With WorksheetsquotSheet1quot With WorksheetsquotSheet1quot kmax el End With With WorksheetsquotSheet1quot with WorksheetsquotSheet1quot Cells1 4 quotIterationquot Cells1 5 quotX1quot Cells1 6 quotX2quot Cells1 7 quotf1quot Cells1 8 quotf2quot End With iteration loop For k 1 To kmax Next 1 Call Jacobianx J Call LUDecomposeU Indx Call LUSubstituteJ Indx b dx 5 For i 1 To n s s dxi 2 xi xi alpha it dxi Next 1 With WorksheetsquotSheet1quot k Cellsk 1 8 f2 End With If Sqrs lt eps Then Exit For Next k End Sub Sub flmX As Double f As Double f1 x1 x2 2 r 3 f2 x1 x2 2 End Sub Sub Jacobianx As Double J As Double J 1 J122x2 J2 1 1 x2 J2 2 x1 x2 2 End Sub 52 55 Multivariable extrema fx A min where f is scalar x is vector There are two basic types unconstrained and constrained Visualization Surface plot and contur plot Most frequent application Nonlinear loeastsquares t 1quot A 39 for 39 and quotE 39 77 in principle77 equivalent fz7 gt min and gradient 0 Just to solve 77gradient 077 is somewhat dangerous it is better watch fx Direct Methods Not using any more information only calculated function values at previous points and at new point A popular direct method is Method of polytopes due to Nelder Mead A 77simplex77 or 77polytope77 is a cluster of n1 points for instance for twovariable functions a simplex consists of three points forming a triangle The Nelder Mead simplex method method of polytopes uses the idea of a 77wan dering77 simplex with possible deformation shrinkage for instance Given a starting simplex the program calculates the function values at the nodes and depending on the results tries to improve the simplex by deleting the worst point and bringing in a better one If the attempts are not successful the simplex is shrinks In case of success it becomes larger Geometrically this means the simplex is wandering and deforming in each 77iteration77 step adopting itself to the shape of the contour lines while systematically approaching the best location Other methods Gradient based steepest descent and second derivative based Newton Excel service Solver Fully nonlinear least squares LS objective function is a sum to be minimized Linear least squares straight line t Generalized Linear Least Squares after transform of variables straight line t max 2 parameters y mm b more than two variables but still linear y 10111 10213 psxixs because the parameters 101102 p3 go into the model linearly Nonlinear least squares several parameters 101102103 the dependence of the response is nonlinear on them P111 9 1P211 P312 Observations 1117 1127 91 1217 1227 y2 l 17117 In27 yn The Objective function is Q lty1 7 P1I11 gtlty2 7 P1I21 gt7 lty 7 10117 gt 1P2111P3172 1P2I21P3122 n 1P2In1 P31n2 The minimization problem is nd 101102103 which minimizes the objective function Q7 is usually solved by a special method called GaussNewton Marquardt Example X1 X2 y l 2 1154 Given the observations 3 4 2045 5 4 3125 7 2 5526 The estimated parameters are 151 3267 132 02257 1 0777 Excels allpurpose 77Solver77 is also an efficient tool to minimize the objective function even when additional constraints have to be taken into account 56 Solution of system of ordinary differential equa tions RK4 A system of first order equations has one function variable the variable being dif ferentiated for each equation in the system and all of these function variables are functions of one independent variable such as I or t The Initial Value Problem lVP is of the form 1 i f1lt1791792 31 f2 1791792 with initial conditons7 for some known a In 162 91W 9 92 a 93 where the right hand side functions fhfg are any functions of the independent variable I and the two function variables yl yg but not of the derivatives of these A two equation system is shown above7 but it can be extended in the obvious way to more differential equations The algorithm now involves vectors watch for boldface Ii 121 h M yH g k1 2k2 2k3 191 where k3 f 1121 7Yi71 5 39kzl k4 flziil hd iq h 39 ksl The vector operations can be done by simple loops 54 57 Partial differential equations and reservoir sim ulation Review Physical meaning Often from conservation principles Initial and boundary conditions Transient versus steadyState With initialvalue problem solution is obtained by starting with initial values and marching forward in time step by step generating successive rows in the solution table Coordinate systems Cartesian or Radial Classi cation 7 Heat equation u um parabolic 7 Wave equation um um hyperbolic 7 Laplace equation um uyy 0 elliptic Numerical approaches Finite difference nite element Semianalytic transfor mation Finite difference approximation of derivatives The diffusivity equation includes derivatives of pressure in space and in time To solve the diffusivity equation numerically we must nd ways to represent these derivatives Getting the derivatives approximated correctly is an important part of getting an accurate numerical solution Transient solution of the diffusivity equation Finite difference implicit Method Linear reservoir constant pressure boundaries 0210 910 77 E where 7 k n 7 M 45 porosity M viscosity ct compressibility total k permea iit We want to solve the transient problem The evolution of the pressure eld if the initial condition is constant pressure pi and the pressures at the two ends are kept constant pleft and gingm respectively First discretize the region on interest over both space 1112 7117 1 and timet0t17 tm tm Replace the analytical derivatives by numerical approximations Approximation of the second partial derivative in the xdirection central Pi1 210139 10121 AI2 Approximation of the rst partial derivative in time between time levels n and n1 So far we have not stated what timestep n or n1 the terms on the left are being evaluate An 77implicit77 solution scheme means that the pressures at the nodes will be evalu ated at the new timestep n1 in the second derivative space term 77101431 2P 10511 Pin1 7 pi A352 At It is useful to introduce a new notation 2 a 77 At Then the system of equations become 0 giVen Pleft 102 7 2 ap1 1025 api for i 2 37 nl p2 given Fright This is a tridiagonal system with the unknowns 1 1 1 pl 7 103 pupil P2 In matrix representation for n 57 as an example 1 cl 0 0 0 p1 d1 a2 22 CQ 0 0 p2 d2 0 as 1 3 03 0 P3 d3 0 0 a4 114 64 P4 d4 0 0 0 a5 25 p5 d5 For instance7 b1 17 Cl 0 and d1 is nothing else7 but the given pleft The system can be solved by the Thomas algorithm Example data porosity 0 2 Viscosity 0 8 cp compressibility ct l 10 5 lpsi permeability k 5 mcl reservoir length L 3000 ft initial pressure pi 2800 psi pressure at X 0 ft pleft 1000 psi pressure at X 3000 ft gingm 2800 psi time t 0100 days With the units indicated7 00063316 77 4W5 and hence 7 45 A902 0 00633k At Choose appropriate mm and ntl Program Option Explicit Option Base 1 Sub VBA111 Linear reservoir constant pressures at the two end points Dim nx As Integer nt As Integer ipr As Integer Dim it As Integer ix As Integer Dim phi As Double mu As ouble ct As Double k As Double Dim xlen As Double pleft As Double pright As Double pini As Double tend As Double Dim alpha As Double dx As Double dt As Double t As Double Dim i As Integer ip As Integer 1 With Worksheetsquotsheet1quot Cells 1 1 Z quotPorosity fraction quot phi Cellsi 2 i Z i 1 61151 1 Z quotViscosit GP mu Z Cellsi 2 1 Cellsi 1 Z quotCompressibility lpsi quot ct Z Cellsi 2 1 Z Cellsi 1 Z quotPermeability Ind quot k Z Cellsi 2 i Z Cells1 1 Z quotLength of reservoir ft quot xlen Z Cells1 2 i Z Cells1 1 Z quotInitial press es p quot pl 1 Cells1 2 1 Z Cells1 1 Z quotLeft pressure ps quot pleft Z Cel s1 2 i Z i 1 Cellsi 1 Z quotRight pressure psi quot pright Z Cellsi 2 i Cellsi 1 Z quotNumber of grid points in x nx Z Cellsi 2 i Z Cellsi 1 Z quotEnd time day quot tend Z Cellsi 2 i Z Cellsi 1 Z quotNumber of time steps quot nt Z Cellsi 2 i Cellsi 1 Z quotPrint at every quot ipr Z Cellsi 2 ReDim anx As Double bnx As Double CDX As Double ReDim dnx As Double pnx As Double b1 Z 1 For i Z 2 To nx Z 1 bi Z 2 alpha Next 1 bnx Z 1 6 31 Z 0 For i Z 2 To nx Z 1 Ci Z Z1 Next i d d1 Z pleft dnx Z pright Initialize t Z 0 For i Z 1 To nx pi Z pini Next i ip Z 1 Cells1 6 Z quotLoc ftZgtquot ForiZ1Tonx ells16i ZiZ1 dx Nexti Cellsip 4 Z quottime day quot ip Z ip 1 ing r it Z 1 To nt Cellsip 6 Z quotp ps1Zgtquot Z it t dt For i Z 2 To nx Z 1 di Z alpha i pi Next i Call Thomasa b c d p If it Mod ipr Z 0 Then Cellsip 4 Z t 1 To nx Cellsip 6 i Z Pi Next i 1 1p Z ip 1 End f Next it End With End Sub Note that this 77simulator77 uses the Thomas procedure It is a good practice to put the Thomas solver into a separate module Explicit Solution Scheme Rewrite the above program to obtain an explicit solu tion algorithm Experiment with the time step Verify the following statement Stability lmplicit stable Explicit at a less than 2 it is unstable 571 Reservoir Simulation Gridblock Finite Volume approach Two dimensional x and y direction Oil Material Balance Equation Oil in place in block i V2750 N 7 B0 where Vp Aszhqb Oil in place changes from V1750 B0 to V1750 n1 B0 during timestep n1 Thus7 the rate of accumulation is i V1750 n1 7 V1750 n At B0 B0 Flow from block il to block i EAST direction 7 000633kkmhAy pi1 7 pi 7 T kr I 414 i MDBD 7A1 7 E MB 0E pm 101 A similar equation can be written from the WEST direction The net in ow is 1 k r TWltMBgt0WP1gt1 PiTE BgtDEP11 Pz Including a possible well term go and the previously determined accumulation term we obtain the gridblock cell balance 59 i V1750 n1 7 V1750 n 7 i V1750 n1 7 V1750 n At B0 B0 4quot 7 At B0 B0 Flow term can be implicit n 1 or explicit We have similar equations for water and gas The most important unknowns are the cell pressuresi Often we solve an implicit scheme for the pressures then update everything else saturations by an explicit scheme IMPE S Special issues Well models Fully implicit schemes various solvers ADI method of alternate directions Chapter 6 Integral Transforms Spectral methods inversion of the Laplace transform 61 Transformations Function to function Laplace Data set to data set Discrete Fourier Transform Spectral methods Seismics7 UltraSound etc FFT algorithm 62 Laplace transform Forward and Inverse Transform Some Rules Analytical integration or numerical solution 621 Numerical inversion of the Laplace transform Put the Laplace Inversion Stehfest Module and the program into Excel Run it The program calculates the Inverse of the Laplace transform Fslsi The result should be the unit step function one fore every positive time Check if you get back the value 10 for the given time points Run iti Modify the program to calculate the integral from zero to t of the unit step function7 by dividing the Laplace transform by s 0 Now change the Laplace transform to another function and run again for times between zero and 10 Select one of the following pairs and check the result using the known inverse 0 Now calculate the integral of the function using the inverse Laplace transform methodi Check the value of the integral F s M 1 6705 05 1 m arctan The main program rst 77initializes77 the algorithm with Stlniti The parameter is the number of terms n in the Stehfest algorithmi In this case we select 16 The actual calculation is invoked by Call StCalC t ft The whole driving program is only a couple of lines Option Explicit Sub Funs As Double F As Double F 1 5 End Sub Sub VBAleO Dim i As Integer np As Integer Dim t As Double ft As Double Call StInit16 With Worksheetsquotsheet1quot Cells1 1 quotnumber of time points quot np Cells1 2 Cells1 3 quottquot Cells1 4 quotftquot For i 1 To np t Cellsi 1 3 Call StCalCt ft Cellsi 1 4 ft Next i End With End Sub It is a good programming practice to put the following Stehfest procedures initial ization and calculation into a separate module Opt ion Explic it Private Stn As Integer Private Stc As Double Function Min11 As Integer 12 As Integer If 11 lt 12 Then M111 11 e M111 12 End Function Function Max11 As Integer 12 As Integer If 11 gt 12 Then Max 11 E se Max 12 End Function Sub StInitn As Integer Stn Max2 Min16 2 Intn 2 ReDim StcStn As Double ReDim Fact Stn As Double Dim n2 As Integer i As Integer k As Integer Dim ck As Double sk As Double n2 Stn 2 Fact0 1 For i 1 To Stn Facti i it Facti 1 Next For 1 1 To Stn ck 0 For k Inti 1 2 To Mini n2 5 k ck ck sk n2 Fact2k Fa u 1Facti 1 H ext Stci 1 i n2 ck Me i End Sub Sub StCalct As Double ft As Double Dim ln2pt As Double 5 As Double F As Double Dim i As Integer ln2pt Log2 t ft 0 For i 1 To Stn 5 ln Call Funs F ft ft Stci F Next 1 ft ln2pt ft d Sub Part IV Appendix Chapter 7 Study guide 71 Be able to if Compare scope of module level and procedure level variables 2 Use arrays 3 Write meaningful comments to coder 4 Subdivide a programming project into smaller units 5 Use defensive programming strategies to check for obvious errors in input pa rameters 6 Use the interactive debugger to trace logic and numerical errors in code 7 Give a simple explanation of the term oating point number77i 8 Give a de nition of roundoff error i 9 Give a simple example of over owi 10 Identify the truncation error from a Taylor series expansion 11 ldentify at least two important differences between symbolic and numeric computationsi 12 Estimate the error of a simple calculation if errors associated With the input are given 13 Write any scalar equation in the form 0 14 De ne multiple rooti 15 Explain the role of bracketingi 16 Write a simple logical expression that expresses the condition for an interval to contain bracket the root7 17 Manually perform a feW steps of the bisection method 18 ldentify situations that cause Newton7s method to fail 19 Calculate the number of steps for the bisection method to reduce the uncer tainty interval xtimes to O tom OH wwwmmmmmmm mworooo ozmpbw CA3 CA3 wwwww oo ozmgtgt CAD to gtgtgtgtgtgt GHQ gtgtgtgtgtgtgtgt mgtgtw g N i Qualitatively compare the convergence rates of bisection falseposition and Newton7s method i Plot any as a means of graphically identifying the location of roots i Describe the role of global variables in nding the roots of fz a b i i 0 where a b i i i are parameters and the method returns the value of I that gives f 0 for fixed values ofa b i Compare effects of noise on numerical differentiation and integration Derive the approximation order of a simple nite difference formula i Manually perform multiple application of Simpsonls rule r Manually perform Richardson extrapolation i Recognize if a least squares problem is linear can be linearized is nonlinear i Be able to determine degrees of freedom i Identify connection between linearized and original model parameters i Identify the criteria for obtaining a least squares t What is minimized i Manually evaluate the formulas for the slope and intercept of a line t to data i Derive the transformations for tting data to nonlinear twoparameter func tions of various type i Make the connection between slope and intercept of the linearized model with the original model parameters i Recognize initial value versus boundary value problem i Manually calculate a few steps of Eulerls and Heun7s methods i Program the RK4 method i Describe concepts singlemulti step self starting or not predictorcorrector i Express conditions of linear dependence of vectors row or column in terms of matrix rank and viceversa i Explain the condition of singularity of an n n matrix in terms of linear inde pendencei i Express matrix rank as a measure of linear independence i Write the formal solution to AI 12 i Explain why it is not a good idea to use the formal solution as a computational procedure for solving AI i i Name the solution algorithm most commonly used for solving AI 12 i Describe the reason for pivoting i Answer the question Is pivoting a remedy for illconditioned systems i Write the preferred expression for solving LI b or Uz b when L is lower triangular and U is upper triangular What algorithm makes use of the solution for these systems i Convert a higher order ODE to an equivalent system of coupled first order DEsi mmmmmmmmmm u b prom ozmtb39wmwrorooo i Use NewtonRaphson method for a system of nonlinear equations i Recognize the basic PDE types Relate transient to steady state i Name methods to solve PDEs numerically i Solve the transient diffusivity equation by implicit nite difference Be able to use the Thomas algorithm Describe basic approaches to Reservoir simulation i Discuss explicit versus implicit approach i Work with units involved Discuss relation between well scheduling and boundary conditions i Be able to use a simulator and interpret results i Have an overview of spectral methods Be able to use an available module in this case for numerical Laplace inversion Chapter 8 Assignments test problems 81 Error propagation borehole volume A 6000 ft i l deep borehole has a diameter 5 inch i 004 inch What is the volume in bbl Give your answer with error estimate 82 Error propagation onegallon cube What is the side length of a onegallon i5 cube Indicate the unit and error estimate of your answer 83 Error propagation barrel Consider a right circular cylinder that has a volume 1 BBL British barrel i006 The height is exactly twice the diameter Estimate the diameter indicating the unit and the uncertainty 84 Error propagation hydrocarbon volume The volume occupied by hydrocarbons in a certain reservoir can be calculated from VH0 Vq5l 7 510 where V is the reservoir volume 45 is the porosity and 51 is the water saturation Assuming V415 25105 ft3 i15 and 51 065 i 012 calculate the hydro carbon volume and give both the absolute and the relative errors of the calculated value Note that the water saturation error is given as an absolute errorl 85 Root nding Newton iteration 2 factor An equation of state written in 77zfactor form77 results in the following nonlinear equation to be so ved 24 25002 71300 0 You use Newtonls method to nd the root starting from z 1000 we will call it the zeroth approximation of the root Carry out the rst iteration step and give the rst approximation of the root Do not continue the iteration process 86 Root nding solution of the PR equation of state he Peng Robinson equation is one of the cubic equations of state lt1 VMltVMbgtLbltVMebgtgt VM b RT where R is the universal gas constant p is the pressure T is the temperature VM is the molar volume and at and b are de ned as follows b 007780 PC 2 2 ac 045724 PC 2 a ac 1037464 1542260 7 026992w217 TTc It is somewhat easier to work with dimensionless quantities such as the deviation factor 2 Then the equation takes the form hp at 212 3122172 all 122 123103 23lt 1gt22lpltW m W 21 WW W or simply 230222512co0 where the coef cients can be read of from the previous expression A substance is characterized by three properties critical temperature To critical pressure p0 and acentric factor w For nbutane To 4252 K p0 375 MPa w 0193 What is the deviation factor 2 and the corresponding molar volume of butane at T 37315 K and p 15 MPa Note there are three roots The smallest one might be the right value if the substance is liquid The largest one is the solution if the substance is in gas state The middle root has no physical meaning So which is the right one We can turn to the saturation curve for help The tension of butane at 100 0C is 153 MPa so the actual pressure given is less than the gasliquid equilibrium pressure Therefore at 15 MPa the butane is gas Advanced Another way is to calculate Gibbs free energies from the equation of state for the two extreme roots and select the one which has less free energy 87 Root nding heat balance The enthalpy of 1 kg methane can be calculated from the following polynomial h 746621983T70001031T2424310 5T37294610 9T4721810 13T5 where h is the speci c enthalpy in kJkg and T is the absolute temperature in K We wish to determine the temperature of 5 kg methane in the nal state if the initial state is 300 K and the enthalpy is increased by 12000 kJ 88 Root nding friction factor The DarcyMoody friction factor f is a function of the Reynolds number Re and the relative roughness e D If the Reynolds number is less than or equal to 2100 then the ow regime is laminar and f 64 Re If the Reynolds number is greater than 2100 then the friction factor can be obtained solving the following nonlinear t i 7 e D 251 equatlon called Colebrook equatlon 7 2010910 lt SJ Ra Write a complete VBA function that has two inputs the values of Re and eD It should calculate the friction factor Hint For bracketing methods choose f as unknown for open methods choose Why 89 Root nding Well deliverability The ln ow Performance Relationship IPR gives the production rate depending on the reservoir pressure and wellbore owing pressure It characterizes the ability of the reservoir to deliver For instance 0 QIPR 4 38 2 10 where the production rate 4 is in MSCFD and both the average reservoir pressure p and wellbore owing pressure pwf are in psi The Tubing Performance Relationship TPR gives the production rate depending on wellbore owing pressure It characterizes the ability of the wellbore to provide vertical lift For instance 05 quR 9368 pal 7105 From the TPR it can be seen that in this case the minimum well bore owing pressure to provide natural lift is 103 psi ln steadystate production 4119 quR Find the production rate at various reservoir average pressures Hint rst rearrange the equation into 77fx077 form 810 Root nding isoterm ash n the simplest isoterm ash calculation we know the inlet composition of the mixture 21 22 2 mole fractions and the distribution factors 1 2 n We wis to nd the mole fraction of vapor V and the composition of the vapor phase and the composition of the liquid phase The equlibrium relations are 1 amp Ii The material balance relations are My Ii1 V 2i Rachford combined all these equations into one nonlinear equation with one un known ii Kl702i 70 i1 K71V1 Once solved for V the composition of the individual phases can be calculated from Ki 2139 yi Kl71V1 Ii KrmH 71 Solve the problem for the following input data Component K 2139 C1 610 03396 C2 90 00646 C5 0035 06068 811 Optimization Maximizing Net Present Value One MSCF additionally produced gas brings 35 US revenuei Since the revenue is realized at the end of the year it has to be discounted For optimum sizing a well stimulation treatment we need to know the costs and bene ts Let us consider matrix acidizingi The size of the treatment is best characterized by the volume of acid spenti All costs depend on this extensive variable The bene ts can be measured by the revenue from the additional oil production in the next three years However the revenue we get back in the future is 77 less valuable77 than the dollar we spend today hence it has to be discounted depending on the year it is realized in The discount ratio is 12 NPV is the difference of discounted revenue and cost Our goal is to maximize the NPV Cost of one gallon of organic acid mixture is 34 US Pumping charge 9 US gali Fixed cost is 10000 US for a single treatment Total cost in US depends on the acid spent v in gal cost 10 000 34 9v Bene t comes from increased production rate The production rate in MSCFD is given by 30 000 7 63 s where v is the acid spent and s is the skin factor modi ed by the acid treatment The pretreatment skin factor is 7 Depending on the injected volume the skin factor decreases according to the relation S 10i1vl3 but never becomes negative The discounted bene t for a 3year horizon is therefore expressed in US as 3 3 5 30 000 30 000 benefit Z W365 W m i1 101v13 Find the optimum amount of acid to spend on this well 72 812 Integration of a function PV The production rate q STBD from an oil well declines according to the Hyperbolic Decline Law n 7 4139 q 7 1 Km where K 000207 a is the decline constant 11 03 is the Arp exponent 11 100 STBD is the initial production rate t is the number of days elapsed from the start of the production It can be shown that for the given decline the time when the production rate reaches the economic limit of 20 STBD is ta 1000 days If we wish to calculate the discounted revenue Present Value of the production we have to multiply produced quantities by the oil price 0 STBi A 77 dollar tomorrow is however less than a dollar today77 so we have to discount the revenue using a yearly discount rate i Year and a discounting method which is 77continuous77 in this case c 25 STB i 12 Year Using the continuous discounting method the Present Value of one barrel oil pro duced in the future it days from now can be expressed by 55 255 The Present Value of the cumulative production from 77now77 to abandonment can be calculated from the integral a 71 PV ce q Zldt o 1 Km I 813 Numerical integration of discrete data pore volume A reservoir of area A 3208 acre 113974 108 ft2 has a pay thickness of 40 ft The pore space is completely lled with oili The porosity was measured as at various depths from the top of the formation Depth ft porosity 45 0 2916 10 03265 20 03187 30 03193 40 02854 Calculate the oil in place in reservoir cubic feet using a multiple application of the trapezoid rule and b multiple application of Simpsonls 13 rule 814 Straightline formation volume factor model 1 Given p5 2012 psi bubble point pressure Data observed p psi B0 resBBLSTB 1 1500 262 1600 1279 1800 1298 Determine the parameters of the nonlinear model describing the B0 1 a IxPb 7 10 What is the best estimate of the B0 at the bubble point B0 815 Straightline formation volume factor model 2 Consider the following model of Formation Volume Factor B0 as a function of pressure p BO aeMpipb where B0 is in resBBLSTB p in psi and p5 is the known bubble point pressure p5 3007 psi The model parameters a and b are to be nd The following laboratory data are available p psi B0 resBBLSTB 1 500 070 1500 1175 2500 1301 Determine the Formation Volume Factor at the bubble point p5 using the above mo e 816 Straightline gas in place Production and static eld pressure data for a gas eld is given below Craft and Haw ins pz psia Gp MMM SCF 6553 0393 6468 1642 6393 3226 6329 4260 6246 5504 6136 7539 6080 8749 The material balance for a volumetric gas reservoir is P Pi Pi i 7 G 2 z p Where p psia is static eld pressure z is deviation factor for corresponding to the pressure Gi MSCF original initial gas in place Gp produced gas MSCF and the subscript 239 refers to initial Find original gas in place from the above data using straightline form and leastsquares t Hint from the straight line form estimate slope and intercept then interpret them to obtain rst pizi then Ci GPGltLGZ39gtE Pi 239 817 Straightline FlowAfter Flow Test IPR Alternative linear form A frequently used lPR equation 2 n P 4 q39max 1 7 P Where 4 production rate BOPD p average reservoir pressure psi pwf welbore owing pressure psi Data observed pwf psi q BOPD 2500 1099 2000 1925 1500 2587 1000 3092 Find the Absolute Open Flow Potential Hint ll out the following table rst Original model Linearized model Parameter l gmax m Parameter 2 n b lndependent variable p X Dependent variable 4 y 818 Nonlinear least squares oil viscosity as a func tion of pressure and temperature Consider the following model of oil viscosity 0 for a certain eld as a function of pore pressure7 p and layer temperature T 0 1112 1Clp15 where 0 is in up p in psia7 T in 0R The model parameters a7 b and c are to be determined by the method of nonlinear least squares using a general purpose minimization program eg7 Solver The available data are p psi T7 0R 0 up 500 460 0764 1500 480 124 2500 500 219 1500 520 102 500 540 0833 Write module containing a VBA function that calculates the objective function to be minimized The rst line of the function should read as function SSQa As Double b As Double c As Double As Double Use global modulelevel variables for all data except for the three model param eters a 127 c which are arguments of the function 819 ODE Production rate decline from a differ ential equation Consider the following differential equation describing the production rate decline of a well 000025 L5 fracdqodt 7 1 01 Ito1 q where go is the production rate in BOPD andt is time in days The initial production rate is 498 BOPD The factor describes the decrease of ow ef ciency with time that is a consequence of evolving damage Use the RK4 method with stepsize h 1007 1007 1007 65 days to calculate the production rate at t 365 days Advanced derive and solve another ODE for the cumulative production 820 ODE pressure versus depth in a shutin well A gas well is shutin at the wellhead no ow Write a differential equation for the pressure assuming the deviation factor is a known function of pressure The temperature can be considered constant Density p M kg m3 Molar mass of gas M 027 kgmol Temperature T 375 K Deviation factor 2 117 cl 3 lO F Pressure at the surface po 32 106 Pa The solution of the differential equation should be a function providing pressure p7 Pa as a function of depth D7 In What is the order of the ODE What is the initial value condition in this case What method can be used to obtain the solution numerically Repeat the derivation for a uid with density p pg 1 Cfp where p0 235 Icym3 and Cf 0310 5 Pa 1 and the pressure at the surface is as before p0 32106 Pa 821 Meaning of matrices and vectors chemical component balance Right a conservation equation in matrixvector form7 for instance for a distillation column 822 System of linear equations mixing Three petroleum products Pl7 P2 and P3 are available The Octane number deter mined by two independent methods research R0 and motor MO are shown in the following table P1 P2 P3 R0 835 912 961 MO 843 887 979 You wish to create a mixture which is of Octane number ROMO2 89 deter mined by any of the methods The question is what mass fractions to use if there is a solution at all It is rea sonable to assume that the resulting R0 and MO numbers are the weighted sums of those of the constituents with the weighting factors being the mass fractions Note Since the density of the uids is the same corresponding mass fractions and volume fractions are also the same Select a reasonable notation for the mass fractions Write a system of equations for the unknown fractions If there are less equations than variables think about one using the de nition of mass fractions If there are more equations than variables think about dependence and delete the unnecessary equation Is there a unique solution Are the fractions all positive Do they sum up yielding one 823 System of linear equations with many right handsides log interpretation A basic approach to determine lithology of simple formations is neutrondensity crossploti Let V1 V2 and V3 denote the volume fractions of sand clay and pore space which is assumed to be lled by water The sum of the volumes fractions is of course unity Let N denote the neutron value dimensionless and p the density in gcm3 Note the of cial Sl density unit is kgm3 i As we all know 1 gcm3 1000 kgm3 The log evaluation is based on the simultaneous application of the two relations N 72V1 37V2 100V3 p 2 65V1 2 41V2 1i0V3 expressing that both the neutron value and the density is a linear combination of the volume fractions For instance if V1 076 V2 0 05 and V3 019 then N 1933 and p 23245 gcmgi If we could measure N and p then we could nd out the volume fraction solving the system 72V1 37V2 100V3 N 2 65V1 2 41V2 1i0V3 p V1 V2 V3 1 Since logs are taken densely say meter by meter we need to solve for the unknown V1 V2 and V3 with many righthand sidesi 78 Write a program to do the calculation calling the LUDecomposition once and the LUsubstitute many times Depth7 ft N p 96013 1000 193 232 1001 1656 237 1002 1589 238 1003 1054 247 1004 1995 233 1005 1974 232 1006 183 234 1007 276 219 1008 2306 225 1009 1449 240 1010 30 215 Additional challenge Using the LU solver pair7 calculate the inverse of the coef cient matrix Check if the original matrix multiplied by the obtained inverse is the identity matrix 824 Nonlinear System of Equations Chemical Equi librium of Methan Combustion Consider the following two chemical reactions CH4 H20 C0 3H2 anyl 3235 COHQO COngHg anM 03726 Notation y means mole fraction The sum of mole fractions is 1 Ky is the equi librium constant expressed with mole fractions7 dimensionless 111191 1nyCH4 1nyH20 1nyco 31I1yH2 1 Kyz i 111 900 111 szo 111 9002 111 sz The an values are given for 1000 K ame temperature and atmospheric pressure The following is the molar composition of the initial mixture Component initial mole frac 1 CH4 2 H2 06 3 C O 0 4 C 02 0 5 H2 0 What is the equilibrium composition Hint Let us start from 1 mole mixture and let 11 and 12 denote the reaction extents expressed in mole The reaction extent tells us how many moles of the lefthand side of a chemical reaction has been converted into the righthandside The following table shows the mole fractions as the function of the two reaction extents notice that the rst reaction causes mole number change7 while the second does not 1 CH4 721 2 H2O 013111227112 3 Co 4 CO 3105 43212 5 H2 12x1 We can reduce the problem to two equations and two unknowns using the reaction extentsi 825 Units Darcy s law If the permeability is k 10 md 9 87 10 15 m2 and the uid Viscosity is M 5 cp 0005 Pa s7 what is the Darcy velocity caused by 05 psift 11310 Pam pressure gradient Assume the cross sectional area perpendicular to the ow is 1 ft2 00929 2 What is the ow rate in resftSs and in resBBLD in 72133 and in mSday Chapter 9 Tables Quantity Length L Time t Mass m Chem substance n Temperature T Current 1 Velocity V Area A Volume V Speci c volume Molar volume Density Acceleration a Force F Pressure p also stress Work Heat lnti Energy Enthalpy Power Speci c heats mass based Electri Potential Voltage Electri Resistance Electri charge Flow rate Shear rate Viscosity Permeability Compressibility Work t Heat m A T PowerCurrent PotentialCurrent Current Time A Volume At A Velocity A Length Stress shear rate Pressure Length Ve ocity Viscosity l Pressure Sl Unit m kg Sl Other name N Pa Nm2 Jm3 J N m Pa m3 J J W Js Pa mSs J kg K V WA Q VA C Multiplier Pre x name Symbol 10g giga G 106 mega M 103 kilo k 10 3 milli m 10 6 micro M 011 deci d 0101 centi c Circumference of circle d 7r rea of circle d2 7r4 Volume of sphere d3 7r6 Surface area of sphere d2 7r Volume of right circular cylinder th 7r4 Pyramid area of base h3 Kinetic energy AKE mA Q Potential energy APE mgAh Hydrostatic pressure Ap pgAh Darcy s Law QA VDWCy 7 2 Pipe Mech Energ Bal Ah imW 16 Pipe Mech Energ Bal Ap gpAh gA Q f ump power ApVA ApQ Pipe Reynolds number de nition Re w p Friction factor def Darcy7Moody7Stanton f Apfrictf72 54 H for laminar ow HagenPoissule E for turbulent ow f functionRe7 relative roughness Acceleration of graVity 918066 ms2 Max Water Density 9991973 kgm3 Liquid water speci heat 1 atm 4180 Jkg K Heat of vapi of water 1 atm 21257 106 Jkg Air Density 273 K7 1 atm 1293 kgm3 Air 51701 114 Air relative molecular mass 012891 kgmol Temperature Abs zero Water freezing pl 1 atm Water boiling pl 1 atm 0 K 27315 K 373 K 273115 0C 0 DC 100 0C 0 R 49167 R 67167 R 45967 0F 32 0F 212 0F pV nRT R up 7 cu 7 GP lcleal gas H i E J P 3 N R 87314m0 5 87314m 2lgfjbf87314m0lf pillfig R1I986m 1545m 10773m Other units Other unit Expressed With S1 unit inc 070254 In ft 073048 In mile 160973 In acre 40469 In2 gal 37854 X 10 3 m3 bbl 0158987 In3 Liter 0001 m3 lbm 0745359 kg lbf 474482 N psi 689478 Pa bar 1 X 105 Pa Atmosphere 101325 Pa mmof Hg 13332 Pa degR 07555556 K hp 74570 W BTU 10551 J min 60 s hour 3600 s clay 86400 s md 9869 X 10 16 m2 CP 0001 Pa s Inside the English system 1 ft 12 in 1 mi 5280 ft 1 ac 43560 ft2 1 sq mi 640 ac 1 BBL 42 gal 576146 1 lbf 327174 43 1 psi 1 fig 144 1 atm 14696 1 BTU 778717 4 25037 1 hp 42741 JET 1 8p graVI liq 8734 3 6274 spi graVI liq lbmole 45359 mol 379 M multiply by 1000
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'