• Author / Uploaded
• Shamit Rathi

Citation preview

PETE 301 PETROLEUM ENGINEERING NUMERICAL METHODS LECTURE NOTES AND STUDY GUIDE Peter P. Valk´o August 18, 2004

Contents I

Introduction 0.1 0.2 0.3

II

4

Course Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Resources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Course structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Tools

5 6 7

8

1 Introduction to Excel and VBA 1.1 Starting Excel, saving a file . . . . . . . . . . . . . . 1.2 Workbook, worksheet, views, printing . . . . . . . . 1.2.1 Review . . . . . . . . . . . . . . . . . . . . . 1.3 Spreadsheet basics . . . . . . . . . . . . . . . . . . . 1.4 Graphing with Excel . . . . . . . . . . . . . . . . . . 1.5 Numerical methods inside Excel . . . . . . . . . . . . 1.6 Database Management with Excel . . . . . . . . . . 1.7 Excel and the WWW . . . . . . . . . . . . . . . . . 1.8 Keyboard macros and what they can’t do . . . . . . 1.9 VBA, variable, sub, function, module . . . . . . . . . 1.10 Structured Programming: Subs, Decisions and Loops

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

9 9 9 9 9 10 10 11 11 11 11 17

2 Mathematica

21

III

23

Engineering models and numerical methods

3 Mathematical Modeling and Engineering Problem Solving 3.1 Conservation Laws and Engineering . . . . . . . . . . . . . . . 3.2 Accuracy and Precision, Approximations and Round-Off Errors 3.3 Error propagation . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Existence and uniqueness, Convergence criterion . . . . . . . . 3.5 Order of the method, rate of convergence, stability, sensitivity . 3.6 Classification of problems and methods . . . . . . . . . . . . . .

. . . . . .

24 24 24 25 26 27 27

4 Methods related to single variable problems 4.1 Roots of equations and extrema of functions . . . . . . . . . . . . . . . . . 4.1.1 Bisection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.2 False position . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29 29 29 29

1

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

4.2 4.3

4.4 4.5

4.1.3 Newton’s method . . . . . . . . . . . . . . . . . . . . . . . . 4.1.4 Optimization: Golden ratio search for a minimum . . . . . Fixed-point iteration as a general framework for iterative processes Representing, manipulating functions . . . . . . . . . . . . . . . . . 4.3.1 Numerical Differentiation of functions that can be evaluated where . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Numerical Integration . . . . . . . . . . . . . . . . . . . . . 4.3.3 Simpson’s 1/3 rule . . . . . . . . . . . . . . . . . . . . . . . Interpolation, Smoothing, Curve fit, Least squares . . . . . . . . . 4.4.1 Least Squares and its Numerical Aspects . . . . . . . . . . Numerical solution of ordinary differential equations . . . . . . . . 4.5.1 Basic methods for solving ODE . . . . . . . . . . . . . . . .

5 Methods related to multi-variable problems 5.1 Matrices, vectors . . . . . . . . . . . . . . . . . . . . . . . 5.2 System of linear equations . . . . . . . . . . . . . . . . . . 5.3 LU Decomposition and Backsubstitution . . . . . . . . . . 5.4 System of nonlinear equations, Newton-Raphson method . 5.5 Multivariable extrema . . . . . . . . . . . . . . . . . . . . 5.6 Solution of system of ordinary differential equations: RK4 5.7 Partial differential equations and reservoir simulation . . . 5.7.1 Reservoir Simulation . . . . . . . . . . . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . . . . . . . . . . every. . . . . . . . . . . . . . . . . . . . . . . . . . . .

33 34 35 36 37 39 39

. . . . . . . .

41 41 43 46 50 53 54 55 59

. . . . . . . .

. . . . . . . .

. . . . . . . .

30 32 32 33

6 Integral Transforms, Spectral methods, inversion of the Laplace transform 6.1 Transformations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Laplace transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2.1 Numerical inversion of the Laplace transform . . . . . . . . . . . .

61 61 61 61

IV

64

Appendix

7 Study guide 7.1 Be able to ... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

65 65

8 Assignments, test problems 8.1 Error propagation: borehole volume . . . . . . . . . 8.2 Error propagation: one-gallon cube . . . . . . . . . . 8.3 Error propagation: barrel . . . . . . . . . . . . . . . 8.4 Error propagation: hydrocarbon volume . . . . . . . 8.5 Root finding; Newton iteration: z factor . . . . . . . 8.6 Root finding: solution of the PR equation of state . 8.7 Root finding: heat balance . . . . . . . . . . . . . . . 8.8 Root finding: friction factor . . . . . . . . . . . . . . 8.9 Root finding: Well deliverability . . . . . . . . . . . 8.10 Root finding: isoterm flash . . . . . . . . . . . . . . 8.11 Optimization: Maximizing Net Present Value . . . . 8.12 Integration of a function: PV . . . . . . . . . . . . . 8.13 Numerical integration of discrete data: pore volume

68 68 68 68 68 69 69 70 70 71 71 72 73 73

2

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

8.14 8.15 8.16 8.17 8.18 8.19 8.20 8.21 8.22 8.23 8.24 8.25

Straight-line: formation volume factor model 1 . . . . . . . . . . . . . . . Straight-line: formation volume factor model 2 . . . . . . . . . . . . . . . Straight-line: gas in place . . . . . . . . . . . . . . . . . . . . . . . . . . . Straight-line: Flow-After-Flow Test (IPR) . . . . . . . . . . . . . . . . . . Nonlinear least squares: oil viscosity as a function of pressure and temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ODE: Production rate decline from a differential equation . . . . . . . . . ODE: pressure versus depth in a shut-in well . . . . . . . . . . . . . . . . Meaning of matrices and vectors: chemical component balance . . . . . . System of linear equations: mixing . . . . . . . . . . . . . . . . . . . . . . System of linear equations with many right-hand-sides: log interpretation Nonlinear System of Equations: Chemical Equilibrium of Methan Combustion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Units: Darcy’s law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9 Tables

74 74 75 75 76 76 77 77 77 78 79 80 81

3

Part I

Introduction

4

0.1

Course Objectives

• Learn the underlying principles behind some commonly used numerical methods. – Understand the route from arithmetic and calculus to numerical methods. – Understand the basic concepts including uniqueness, approximation, error, convergence, and stability. • Learn how to program a numerical method using modern programming concepts – Problem solving steps – Structured programming – Use of subroutines, modules – Learn how to make use of Visual Basic for Applications. – Learn how to use an integrated program development system providing debugging facilities. • Recognize the main features of a mathematical problem arising in a technical context. – Single- or multi-variable – Linear or non-linear – Involves algebraic, differential or integral equations – Direct or iterative technique should be applied – Explicit or implicit scheme is appropriate • Develop skills to formulate a problem starting from physical principles. – Start with conservation principles – Include material properties – Add problem specific conditions – Learn how to use various resources in solving a complex problem – Use of software systems; Subroutine libraries, modules – Internet resources • Improve ability to critically analyze numerical results • Improve documentation skills

5

0.2

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 efficient use of our time, we request you to read the material before the lecture. 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 make full use of the Living Textbook you will need a software system called Mathematica. This course, however, does not require you to learn how to program in Mathematica. 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 specific petroleum engineering problems also will be provided in the form of Living Textbook. You may experiment with them and may use the results for checking your Visual Basic for Applications (VBA) program results.

Basic References 1) S.C. Chapra and R. 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, B. Maggard, R. Archer, R. Wattenbarger .

6

0.3

Course structure

1. Engineering problem solving tools and programming (a) Engineering problem solving approaches, software tools (b) Excel (c) VBA (d) Structured and Modular Programming 2. Basics of numerical methods (a) Taylor polynomial, Numerical errors, Error propagation (b) Iteration, Convergence, Order, Stability, Sensitivity (c) Classification of problems and methods 3. 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; Simpson’s rule (c) Interpolation, Smoothing, Curve fit, Least squares (d) Numerical solution of ordinary differential equations 4. Methods related to multi-variable problems (a) Matrices, vectors (b) System of linear equations (c) Gauss, Gauss-Jordan, LU decomposition Special cases, Iterative methods (d) Multivariable nonlinear equations: Newton-Raphson (e) Multivariable extrema: Nelder-Mead simplex, Fully nonlinear least squares 5. Partial differential equations and reservoir simulation (a) Classification and Numerical Solution of PDEs (b) Transient solution of the diffusivity equation (c) Reservoir simulation 6. Transformation methods (a) Spectral methods (b) Laplace transform inversion

7

Part II

Tools

8

Chapter 1

Introduction to Excel and VBA 1.1

Starting Excel, saving a file

To create a new workbook: On the File menu, click New, and then click Blank Workbook on the New Workbook task pane. 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 file format (if you want other than .xls) (In general, it is useful to click the right mouse button. Most of the time it will offer a menu item you really want.)

1.2 1.2.1

Workbook, worksheet, views, printing Review

A Microsoft Excel workbook is a file 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 worksheet. When you create a chart, you can place the chart on the same worksheet as its related data or on a separate chart sheet. To view and scroll independently in different parts of a worksheet, you can split a worksheet horizontally and vertically into separate panes. 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 worksheet. To see how it will look like in print: On the View menu, click Page Break Preview. 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 disk.

1.3

Spreadsheet basics

The difference between relative and absolute references is the key concept to understand:

9

Relative references A relative cell reference in a formula, such as A1, 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 references. For example, if you copy a relative reference in cell B2 to cell B3, it automatically adjusts from =A1 to =A2. Absolute references An absolute cell reference in a formula, such as $A$1, always refer to a cell in a specific 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 $A$1.

1.4

Graphing with Excel

Import a text file Click the cell where you want to put the data from the text file. 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 Import External Data, and then click Import Data. In the Files of type box, click Text Files. Check out also: On the Data menu: Text to Columns. 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 first enter the data for the chart on the worksheet. Arrange your data in columns, with x values in the first column and corresponding y values and/or bubble size values in adjacent columns. 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 later. Best-fit equations may be calculated for chart data via the Add Trendline dialog. 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.

1.5

Numerical methods inside Excel

Goal Seek The Goal Seek dialog box can be used to solve for roots of a single equation. An initial guess is needed. Goal Seek performs calculations which improve the solution accuracy. 10

Solver The Solver can be used to solve both linear and nonlinear optimization problems, i.e. 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 independent (decision) variables If the values of the independent variables are subject to the constraints, then cells constrained.

1.6

Database Management with Excel

Check out Search, Filter and Sort in the Data menu.

1.7

Excel and the WWW

Excel files may be converted to HTML (web) documents. The HTML Wizard guides you through the process.

1.8

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 ”play back,” 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 workbook. 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.

1.9

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 11

(c) String manipulation (d) File system 2. Ability to control all of Excel’s features through the Microsoft Excel Objects, including (a) Data management (b) Graphics (c) Solver (for optimization and root finding) 3. Ability to add custom functions which can be used like ”built-in” Excel functions 4. Capabilities allowing the development of event driven, user friendly programs through a variety of interface tools (a) Dialog and input boxes (b) Custom user forms with spinners,checkboxes, etc. (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 t. 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 first-if pen and paper turns you off why not use your word processor? You could even enter the plan as comment lines in the program. Option Explicit Const g As Double = 9.81 ’ gravitational acceleration, m/s^2 Sub VBA011() ’vertical motion under gravity Dim s As Double ’ vertical displacement, m Dim t As Double ’ time, s Dim u As Double ’ initial velocity, m/s

12

u = 30 t = 6 s = u * t - (g / 2) * t ^ 2 MsgBox "s = " & CStr(s) &" meters." End Sub

’ s is converted to a string

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 debugger). 2. In ”Project” window of VBA, right click on your new Excel workbook (usually ”VBAProject (Book1)”) and insert a new module. 3. Cut and paste the text of VBA011, shown in the above text box into the newly created module. Note that the ”Option Explicit” is not part of the program, but is a module option. 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 Sub/User Form. 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 (’) denotes the beginning of a comment. 3. In PETE 301 every module should require explicit variable typing by using ”Option Explicit” 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 subroutine. (Later we will learn how to read from the worksheets.) 6. Output is done using a message box. Later we learn how to put results into the worksheets. 7. This module contains only one subroutine (macro). In fact many of them can be placed into one module. Built-in (or intrinsic) functions Visual Basic has many built-in functions. You will often need the math functions: Abs, Atn, Cos, Exp, Fix, Int, Rnd, Sgn, Sin, Sqr, Tan. Log is natural-log. If you wish ten-based logarithm, use Log(x)/Log(10).

13

Example: Sin function Returns a Double specifying the sine of an angle. Syntax: Sin(number) 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 pi/180. To convert radians to degrees, multiply radians by 180/pi. Those functions are useful, but they are not customizable. Even if you find a Visual Basic function that is ”this close” 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 1. 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 Boolean Integer Long(long integer) Single (single floating-point) Double (double floating-point) String (variable-length) Variant (with numbers)

1 byte 2 bytes 2 bytes 4 bytes 4 bytes 8 bytes 10+length 16 bytes

0 to 255 True or False -32,768 to 32,767 -2,147,483,648 to 2,147,483,647 -3.402823E38 to 3.402823E38 -1.E308 to 1.8E308 0 to appr 2 billion 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 NumberOfEmployees As Integer Dim x As Double, y as Double

When Option Explicit appears in a module, (as it should in PETE301) you must explicitly declare all variables using the Dim, Private, Public, ReDim, or Static statements. If you attempt to use an undeclared variable name, an error occurs at compile time. Note: Dim a,b As Integer

will declare b as an Integer, but a will be the default (Variant). You have to write Dim a As Integer, b As Integer

14

if you wish to declare both as integer. Numerics VBA makes a clear distinction between integer and real. Integers are used for counting things. Real numbers include a fractional part. Real numbers stored as integers are rounded to the nearest integer. Dim n = n = n = n = m =

n As Integer, m As Double 4/3 !result is 1 5/3 !result is 2 7/2 !result is 2 5/2 !result is 4 5/2 !result is 2.5

You should be aware of rules regarding priorities. For instance a frequent error is to write z = p*v/R*T

instead of the correct: z = p*v/(R*T)

Strings To convert any other value to a string, use the Cstr intrinsic function. Objects,input and output The structure ”With/End With” and the ”operator” point can be understood if you are familiar with the concept of object. ”ThisWorkbook” is a hierarchical object. It has several Workseets (a whole collection of them). One of them is called ”Sheet2”. The following example shows how to read a number into variable ”a” from cell ”A4” of Sheet2: With ThisWorkbook.Worksheets("Sheet2") a =.Cells(1, 4) End With

This will output date and time into cells ”A4” and ”A5” of ”Sheet2”: Sub VBA014() Dim d As String, t As String d = Date() ’ Intrinsic function t = Time() ’ Intrinsic function ’ Display the result With ThisWorkbook.Worksheets("Sheet2") .Cells(1, 4)= d .Cells(1, 5)= t End With End Sub

(Objects have properties and methods associated with them. 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 reached/altered. 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). 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 procedure-level variable looses its value between two calls, unless it is declared Static. Both module-level 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 program. 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 argument’s data type. (You can pass an argument by value if you include the ByVal keyword in the procedure’s declaration. Arguments passed by value consume from 2 16 bytes within the procedure, depending on the argument’s data type. Larger data types take slightly longer to pass by value than smaller ones.) Errors, Debugging and Testing Compilation errors are errors in syntax and construction, 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 frequent 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: 1. Ordinary errors-the compiler will attempt to continue compilation after one or more of these errors has occurred, e.g. ”missing End If statement” 2. Fatal errors-the compiler will not attempt further compilation after detecting a fatal error, e.g. Program too complicated-too many strings 3. Warnings-these are not strictly errors, but are intended to inform you that you have done something unusual which might cause problems later, e.g. ”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 used) Some typical errors: g = 9,81 x = 1 + (2 * 3 .Cell(2,1)=5 If (a=b) x=3 Else x=4 End If

’ ’ ’ ’

comma instead of decimal point unpaired parentheses Wrong function name: Cell instead of Cells Missing then

If a program compiles successfully, it will run. Errors occurring at this stage are called run-time errors, and are invariably fatal, i.e. the program ”crashes”. An error message, such as 1. Floating point division by zero 16

2. Overflow 3. Negative argument in the function sqr() is generated. Overflow 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 trash. Negative argument to Sqr or non-positive argument to Log is often the result of a previous logical error. At times a program will give numerical answers to a problem which appear inexplicably different from what we know to be the correct mathematical solution. This can be due to rounding error, which results from the finite precision available on the computer. , e.g. two or four bytes per variable, instead of an infinite number. Run the following program : Option Explicit Sub VBA022() Dim sumx As Double sumx = 0.1 Do sumx = sumx + 0.001 Debug.Print sumx If (sumx = 0.2) Then Exit Do Loop End Sub

Watch the Immediate Window! You will find that you need to crash the program to stop, e.g. with ctrl-break. The reason is: x never has the value 0.2 exactly, because of rounding error. In fact, x misses the value of 0.2 by about 10−9 . Now try this: If (abs(sumx - 0.2) < 1e-6) 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 100+1 times. Use test problems. Print out intermediate results. Investigate limiting cases. Message Boxes, Custom Dialogue Boxes

1.10

Structured Programming: Subs, Decisions and Loops

Using a Subroutine Subroutines are an important component of ”top down” program design. Subroutines allow the programmer to build and test program procedures that perform a single, well defined 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.

17

Example Let in your Excel workbook, ”Sheet3” contain the following: Cells A1 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 defined 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 project. Because the My Line subroutine has arguments, it cannot be executed directly, like the previously discussed subroutines, by using the VBA Run menu. Option Explicit Private Intercept As Double Sub VBA015() Dim X1 As Double, X2 As Double, Y1 As Double, Y2 As Double Dim Slope As Double ’ Get the two points from the worksheet, "Sheet3": With ThisWorkbook.Worksheets("Sheet3") X1 = .Cells(1, 1) Y1 = .Cells(1, 2) X2 = .Cells(2, 1) Y2 = .Cells(2, 2) End With ’ Calculate the Slope and Intercept Call My_Line(X1, Y1, X2, Y2, Slope) ’ Display the result With ThisWorkbook.Worksheets("Sheet3") .Cells(3, 1)=" Slope:" .Cells(3, 2)= Slope .Cells(4, 1)= "Intercept:" .Cells(4, 2)= Intercept End With End Sub ’A subroutine used inside: Sub My_Line(x As Double,y As Double,xx As Double,yy As Double,m As Double) m = (yy - y) / (xx - x) Intercept = y - m * x End Sub

Note: Watch for scope of the various variables! Intercept is a module level variable, but Slope is not. Intercept cannot be seen from other modules. The If/Then/Else Decision Structure Example (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 < 3000 Then apr = 0.03 Else apr = 0.06 End If

18

Loops: the heart of programming math Countable: Do the block for i = 1, then i = 2 , etc.: For i = 1 To 10 "block..." Next i

Count down: Do it for j = 100, then j = 90 etc.: For i = 100 To 50 Step -10 block Next i

Forever: Loop forever unless ”If” stops it: Do If "logical_expression..." Then Exit Do "block..." Loop

Condition controlled: Do while logical expression is true: While "logical_expression..." "block..." Wend

”Select Case” structure makes your code readable Use the Select Case statement as an alternative to using ElseIf in If...Then...Else statements when comparing one expression to several different values. While If...Then...Else statements can evaluate a different expression for each ElseIf statement, the Select Case statement evaluates an expression only once, at the top of the control structure. Select Case performance Case 1 bonus = salary * 0.1 Case 2, 3 bonus = salary * 0.09 Case Is > 3 bonus = 100 Case Else bonus = 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 aren’t arrays) and array variables is that you generally must specify the size of the array. An array whose size is specified is a fixed-size array. An array whose size can be changed while a program is running is a dynamic array. Whether an array is indexed from 0 or 1 depends on the setting of the Option Base statement. If Option Base 1 is not specified, all array indexes begin at zero. Declaring a fixed array In the following line of code, a fixed-size array is declared as an Integer array having 10 rows and 10 columns: 19

Option Base 1 Dim MyArray(10, 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 Variant. To write code according to PETE301 conventions, explicitly declare your arrays to be of a data type other than Variant. Arrays can also be dimensioned ”on-the fly”, using the ReDim command. It is used for dynamic allocation of arrays on the procedure level. A useful related function is UBound(). It returns the largest available subscript for the indicated dimension of an array.

20

Chapter 2

Mathematica Many of the concepts, methods and applications of this course are now programmed as ”Living Textbook”. Using the software Mathematica any section (notebook) of the ”Living Textbook” can be read on the screen and/or run as a program. The particular examples can be modified and rerun providing a unique possibility for experimenting with numerical methods and petroleum engineering applications. Mathematica is a complete environment for calculating and communicating 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 scientific papers, technical reports, presentations, and courseware that include text, active formulas, twoand three-dimensional graphics, and customizable buttons and palettes. A growing library of application packages provides specialized capabilities in areas such as engineering, finance, statistics, data analysis, optics, astronomy, and fuzzy logic. The associated web site is at http://www.wolfram.com. 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, Mathematica, 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 processor-type document and a Fortran-type 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. To run the whole notebook: Select Kernel, Evaluate notebook from the menu. If the system asks a question about evaluating initialization cells, respond clicking ”yes”. To see what is happening use the mouse and the scroll bar or the PgUp and PgDn keys. To 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

21

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 graphics may overload the system and especially the printer, we ask you not to print the notebooks, unless it is necessary. Before printing make sure you deleted unnecessary parts. Living Textbook LTS021 LTS031 LTS032 LTS033 LTS041 LTS061 LTS071 LTS072 LTS081 LTS091 LTS092 LTS101 LTS102 LTS121 LTS131

Content Taylor Series Expansion, Remainder Term and Truncation Error Root-Finding More about iterations Extrema of single variable function Numerical Integration of functions Ordinary Differential Equations Basic Concepts Numerical Solution of the ODE Initial Value Problem Reservoir Material-Balance Vectors, Matrices, Linear system of Equations Direct Solution of Systems of Linear Equation Systems of Linear Equations: Gauss-Siedel Iterative Method Solution of System of Non-Linear Equations Minimization of a Multi-Variable Function Oil-Water 2D IMPES Spectral Methods: Discrete Fourier Transform

22

Part III

Engineering models and numerical methods

23

Chapter 3

Mathematical Modeling and Engineering Problem Solving 3.1

Conservation Laws and Engineering

3.2

Accuracy and Precision, Approximations and RoundOff Errors

Precision is related to repeatability or to the number of digits obtained. Accuracy is related to how near we are to the ”true” value. Obviously, accuracy is more important than precision. In analysis of the numerical methods we often start with the Taylor series expansion: f (x + h) = f (x) + h · f 0 (x)/1! + h2 · f 00 (x)/2! + . . . Several terms on the right hand side will enter the ”formula” of the method. The first 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 arithmetic and result produced by same algorithm using limited precision arithmetic. Computational error is sum of truncation error and rounding error, but one of these usually dominates. Since usually input data are in error and we do a lot of ”steps”, the most important issue is: How does the error propagate? Measures of the error in approximating a number with true value x by an approximation x: 24

True Error, Et = x − x Relative Error,

x−x x

Percentage Relative Error,

x−x 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.

3.3

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 first order Taylor approximation: f (x + ∆x, y + ∆y) = f (x, y) +

∂f ∂f ∆x + ∆y + . . . ∂x ∂y

Rearranging and introducing ∆f = f (x + ∆x, y + ∆y) − f (x, y) we obtain the basic relatiuon of ”error calculus”: ∆f = k

∂f ∂f k∆x + k k∆y + . . . ∂x ∂y

In this error calculus ∆x means always a positive quantity, it represents the maximum possible value of the actual error. Derive what are the two partial derivatives of f f f f

=x+y =x−y =x·y = x/y

From there we obtain, that for summation/substraction 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 z-factor, given that: p = 3245 psi v = 1.977 f t3 n = 1 lb − mole T = 739.67 o R

∆p = 3 psi ∆v = 0.003 f t3 ∆n = 0.001 lb − mole ∆T = 0.03 o R

Differentiating z = pv/(nRT ) we obtain 25

(meaning: ±3 psi)

p −pv v ∆z = k nRT k∆p + k nRT k∆v + k n−pv 2 RT k∆n + k nRT 2 k∆T

A simpler form is ∆v ∆z = kzk ∆p p + v +

∆n n

+

∆T T



= 2.8 · 10−3

(Does the result have any dimension? Why did we not use more decimal figures? 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 453.592 mol, or sometimes the mass or weight of it; R = 10.73 psi f t3 /(lb − mole o R ) is the universal gas constant in the ”English” or ”field” 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 roundoff) are not amplified from step to step but are rather kept under control Ill-conditioned problem A problem is ill-conditioned if results are very sensitive to input data (regardless of the method used)

3.4

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 n-th 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 obtained which, hopefully, converge towards the desired solution: x0 , x1 , x2 , x3 , x4 , ...... If they appear to be converging, then a measure something like the percentage relative error is used: xi+1 − xi εa = × 100 xi+1 (If, however, the x can be near-zero, an absolute error criterion is better. Why?) If this error approximation εa drops below some user-specified tolerance, then the process is stopped. (The results should be checked if possible, in some other independent way for the required accuracy.)

26

3.5

Order of the method, rate of convergence, stability, sensitivity

Since most methods are derived from some Taylor expansion considerations, the theoretical order of a method is related to the truncation of the Taylor series. Examples will be given in the numerical differentiation section. Broadly speaking: the order of the method is the exponent of the variable step-size in the leading error term. Total error: truncation and round-off 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 log-log plot, we will see that the error is decreasing with step-size, and the relation appears to be a straight line. If, for instance, the for every log cycle of the step-size the error changes two log cycles, the rate is quadratic. Often instead of the effect of step-size we are more interested in the effect of number of iterations. Trade-Off Regarding ”step-size” The point is, that smaller step-size provides larger accuracy, but after a while the increased number of calculations bring more round-off errors. Therefore there exists always an optimum. Trade-Off 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 backfire the same way as to many small steps. An algorithm is stable if errors (originated from truncation and round-off) are not amplified from step to step but are rather kept under control. The opposite is unstable. A problem is well-conditioned, if results are not very sensitive to input data (at least if an appropriate method is used). The opposite is called ill-conditioned.

3.6

Classification of problems and methods

Single variable, multi variable (our main classification) Other can be linear vs. nonlinear, direct vs. iterative, explicit vs. implicit. • Single variable methods: Root of a single equation Minimum of function Manipulation of functions (differentiation, integration) functions given in form of expression or algorithm discrete data point (smoothing, curve fitting + diff and int) Ordinary Differential Equation (ODE)

27

• Multi variable problems: Linear Algebra, Matrices, vectors Systems of linear equations direct, special, iterative Nonlinear System of nonlinear equations Minimum of multi variable function (general and least squares) System of ODE Partial Differential Equations Reservoir simulation • Other: Transformation methods

28

Chapter 4

Methods related to single variable problems 4.1

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.

4.1.1

Bisection

In the bisection method first the user estimates, approximately, the position of the root providing a lower- and upper-bound bracketing the root. If the interval [a, b] is a bracket, then f (a)f (b) ≤ 0. A first estimate of the root is then computed as the mid-point between the two bounds. x = (a + b)/2 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 mid-point, x, depending on which previous function value has the same sign as f (x). 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 ”bisection”. We can calculate a-priori the necessary number of steps to reach a certain fraction of the initial bracket.

4.1.2

False position

The search starts with an interval [a, b]. It is assumed that the function changes sign in the interval, f (a)f (b) ≤ 0 as in the case of the bisection method. Writing the equation of the line passing through the two points (a, f (a)) and (b, f (b)) y − f (b) =

f (b) − f (a) (x − b) b−a

29

and requiring that y should be zero, we obtain for the next approximation of the root: c=b−

f (b) f (b)−f (a) b−a

It is easily seen, that c is inside the interval [a; b]. Now we calculate f (c) and overwrite either a or b with c, depending on which function value has the same sign as f (c). After one step, again we have a bracket, but it is smaller than the initial one. (Note: f (a) c = a − f (a)−f (b) a−b

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 c from any of the old locations a or b should be less than a specified tolerance.

4.1.3

Newton’s method

The method (also called Newton-Raphson method) may be used to solve a general equation f (x) = 0, provided the left-hand-side can be differentiated ”analytically”. If an initial guess is known, we can calculate the function and its derivative at x0 . Writing the equation of the tangent line: y − f (x0 ) = f 0 (x0 )(x0 − x) Substituting y = 0 we find that the line crosses the x axis at x1 = x0 −

f (x0 ) f 0 (x0 )

This is the iteration formula of the Newton’s method. The obtained x1 can be then renamed to x0 and hence we can continue the iterations. Example. Suppose that f (x) = x3 + x − 3 Then f 0 (x) = 3x2 + 1 The following program starts the iteration with x = 2. It uses two internal functions: F for f (x) and DF for its derivative: f 0 (x). It stops either when the absolute value of the step is less than 1E-6, or after 20 iterations. Note that there are two conditions that will stop the While loop: either ”convergence”, or the completion of 20 iterations. Otherwise the program could run indefinitely.

30

Option Explicit ’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 ’left-hand-side function Function F(x As Double) As Double F = c3 * x ^ 3 + c2 * x ^ 2 + c1 * x + c0 End Function ’derivative Function DF(x As Double) As Double DF = 3 * c3 * x ^ 2 + 2 * c2 * x + c1 End Function ’Main sub Sub Newton(x As Double, fx As Double, Converged As Boolean) Dim Its As Integer, d As Double ’Initialize Its = 0 Converged = False x = x0 ’Loop While (Not Converged And Its < MaxIts) d = -F(x) / DF(x) x = x + d Its = Its + 1 Converged = Abs(d) < Eps Wend ’Clean-up fx = F(x) End Sub ’ Driver program Sub VBA041() ’Values Group 1: Coefficients c3 = 1: c2 = 0: c1 = 1: c0 = -3 ’Values Group 2: Iteration control MaxIts = 20: Eps = 0.000001: x0 = 2 ’Declaration of Group 3: Arguments Dim x As Double, fx As Double, Converged As Boolean Call Newton(x, fx, Converged) With Worksheets(1) .Cells(1, 1) = "x" .Cells(1, 2) = x .Cells(2, 1) = "fx" .Cells(2, 2) = fx End With If Converged Then MsgBox ("Convergence criterion satisfied ") Else MsgBox ("Convergence criterion NOT satisfied ") End 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 DF, and some of the input.

31

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 00 (xi ) is very large in which case the graph y = f (x) is curving very rapidly. A second problem occurs if f 0 (xi ) is very small which can happen if the function is almost parallel to the x− axis near the root, or if there is a second root very close by (in which case f 0 (x) = 0 somewhere between the roots), or if the graph of y = f (x) just touches the x - axis at the root but does not cross it. (Multiple root.) It is often observed, that the method converges from a ”near enough” starting point, but diverges if the starting point is selected without care.

4.1.4

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 real-valued function f () is unimodal on interval [a; b] if there is unique location x0 such that f is strictly decreasing for x < x0 (that is on the left), and strictly increasing for x > x0 (that is on the right). Suppose f is unimodal on [a; b], and let x1 and x2 be two points within interval, with x1 < x2 . Evaluating and comparing f (x1 ) and f (x2 ), we can discard either [x2 ; b] or [a; x1 ], with minimum secured to lie in the remaining subinterval. To repeat the iteration, we need to carry out only textbfone new function evaluation. To reduce the length of interval by a fixed 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 p this, choose the relative positions of the two points to satisfy the golden ratio: τ = ( (5) − 1)/2 ≈ 0.618. Whichever subinterval is retained, its length will be 0.618 relative to previous interval, and one interior point retained will be at position either at 0.618 or 0.382.

4.2

Fixed-point iteration as a general framework for iterative processes

Many iterative methods for solving nonlinear equations can be viewed as xk+1 ⇐ g(xk ) where g(x) is constructed from the f (x) function. (If you understand the above formalism, you can easily create the ”g-function” for example, for the Newton iteration: g(x) = x − f (x)/f 0 (x).) The scheme is called fixed-point iteration or direct iteration. There are four basic cases: • monotone convergence • oscillatory convergence 32

• monotone divergence • 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”.

4.3 4.3.1

Representing, manipulating functions Numerical Differentiation of functions that can be evaluated everywhere

When an analytical solution to the first derivative of a given function is difficult or inconvenient 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 f (x) = sin(x) Calculate the derivative at x = π /3 using the simplest forward divided difference Use step size: 1E-5 , 1E-10 , 1E-15 , 1E-20 , . . . (until your calculator stops giving reasonable result ) (Use ”radian” switch if you need.) Forward Difference Approximation (first derivative) The Forward-difference approximation method for the numerical differentiation of a function can be derived by considering ”differentiation from first principles” or by considering ”Taylor’s expansion”. The ”differentiation from first principles” derivation is intuitive: f (x + h) − f (x) h→0 h

f 0 (x) = lim

As a computer cannot divide by zero the computed (finite) version of this expression is: f 0 (x) ≈

f (x + h) − f (x) h

This is the Forward-difference 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 Taylor’s expansion provides more information on the method: f (x + h) = f (x) +

h h · f 0 (x) + · f 00 (x) + . . . 1! 2!

Rearrange and represent all remaining terms in one ”remainder” term” f 0 (x) =

f (x + h) − f (x) h 00 + · f (z) h 2

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 (h/2) · f 00 (z) , i.e. proportional to h. This method is of first order. 33

A plot of the error on a log-log scale would result in a 45o straight line. To ”minimize” the error choose a small value of h, but h should not be too small as ’round-off errors’ in the machine arithmetic increase as h decreases. Backward difference approximation (of the first derivative) is very similar: f 0 (x) =

f (x) − f (x − h) h 00 + · f (z) h 2

Again, we neglect the term containing z, Therefore f 0 (x) ≈

f (x) − f (x − h) h

and the method is of first order. The central-difference approximation of the first derivative gives reasonably accurate results and is often implemented: f 0 (x) ≈

f (x + h) − f (x − h) 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 log-log plot of the error (on an example with known exact result) would give a straight line with slope two. The central-difference approximation of the second derivative is: f 00 (x) ≈

f (x − h) − 2f (x) + f (x + h) h2

This is a second order method.

4.3.2

Numerical Integration

Review geometric meaning Review physical (engineering) meaning Trapezoid rule The aim is to calculate, numerically, the integral of a function f (x) from a to b. The basic idea is to evaluate the function at specific 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 ’Trapezoidal Rule’ is the simplest. Consider integrating a known function f (x) over the interval [a; b]. The Trapezoidal Rule gives the following expression for the exact integral, I: I=

h h3 00 · (f (a) + f (b)) − · f (z) 2 12

where h is the interval [a; b] , f 00 (z) is the second derivative of the function evaluated at some unknown point z, between a and b. The value of f 00 (z) is generally not known (z 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 ’truncation’ 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 [a, b] into n ”panels”:   1 1 · f (a) + f (x1 ) + f (x2 ) + ... + f (xn−1 ) + · f (b) I=h 2 2 where x0 = a x1 = a + h x2 = a + 2 · h xn−1 = a + (n − 1) · h xn = b h = (b − a)/n

4.3.3

Simpson’s 1/3 rule

The simple Simpson’s rule uses a parabolic approximation to the curve. The multiple application of Simpson’s rule is usually the first choice of engineers, because of its reliability and simplicity: h (f (a) + 4f (x1 ) + 2f (x2 ) + ... + 4f (xn−1 ) + f (b)) 3 where n is even and h = (b − a)/n. I≈

General idea of Gauss quadratures A quadrature rule is weighted sum of finite number of sample values of the integrand function. (Both the trapezoid and Simpson’s 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 e.g. 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: • Do calculation with h and with h/2; • Express remainder term as a power of h; • The difference between the two results can be used to improve our final approximation; Try to extrapolate to ”infinitely small” h.

Example: repeated Simpsons rule has a leading error term of 4-th order. Assume we did the calculations with h and with h/2. Therefore we have two computed values: I1

35

and I2 . However, what we really want to know is I, that is still unknown. I1 = I + Ch4 4 I2 = I + C h24 In the above equations C is also unknown. To eliminate it, multiply by 24 and −1: 24 I2 = 24 I + Ch
4 − I1 = I + Ch4 Then subtract 24 I2 − I1 = (24 − 1)I Finally rearrange and obtain estimate of the unknown I: 16I2 − I1 15 Note that we do not really need C and we have not calculated it. I≈

4.4

Interpolation, Smoothing, Curve fit, Least squares

Review Purposes for Interpolation: Plotting smooth curve through discrete data points Reading between lines of table, etc. If the data is noise corrupted (is of finite accuracy) we rather use approximation (smoothing) instead of strict interpolation. Interpolating by simple function families • Polynomials • Piecewise polynomials (splines) • Trigonometric functions (will be considered among ”spectral” methods) • Exponential functions • 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 noise-corrupted data • Differentiation is more sensitive to error in input data: You might want to use quite large step size (why?) • Integration is ”smoothing out”: Step-size is the smaller the better in almost all cases.

Approximation by least-squares: The meaning of the criterion Sum of Squares as a measure of ”how good the fit is”. Other requirements might include: • Smoothness, • Least number of parameters, • Extrapolating power, etc. Two basic cases: • Model is based on ”first principles” • Model is just a convenient vehicle

4.4.1

Least Squares and its Numerical Aspects

When we want smooth approximation, we need some degrees of freedom. For instance, to fitting a straight line (2 parameters) to ten points leaves us with 8 degrees of freedom. There are infinitely 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 fit is also called linear regression line. A word of caution: Least squares methods are very adversely affected by ”bad data”. For example, if a 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 the line. 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 ”goodness of fit” such as sum of absolute deviations.) Straight line fit - linear regression Suppose that some the data points should fit a straight line, y = mx + b, if it were not for experimental and measurement errors. The problem is to find the straight line characterized by m b and bb which best fits the data points in the least squares sense. Model: y = mx + b 37

Observations: (x1 , y1 ) (x2 , y2 ) . . . (xn , yn ) The derivation of formulae for the slope and intercept should be familiar by now. The final formulae are Pn Pn Pn n ( i=1 xi yi ) − ( i=1 xi ) ( i=1 yi ) m= Pn Pn 2 n ( i=1 x2i ) − ( i=1 xi ) b=

Pn Pn ( i=1 yi ) − m ( i=1 xi ) n

Useful subroutines for straight-line fit: Option Explicit Option Base 1 Sub Dim n m b End

StraightLineFit(x() As Double, y() As Double, m As Double, b As Double) n As Integer = UBound(x) = (n * scalprod(x, y) - sum(x) * sum(y)) / (n * scalprod(x, x) - sum(x) ^ 2) = (sum(y) - m * sum(x)) / n Sub

Function sum(x() As Double) As Double Dim n As Integer, i As Integer n = UBound(x) sum=0 For i=1 to n sum=sum + x(i) Next End Function Function scalprod(x() As Double, y() As Double) As Double Dim n As Integer, i As Integer n = UBound(x) scalprod=0 For i=1 to n scalprod=scalprod + x(i)*y(i) Next End Function

Transformations to straight line form Models of physical phenomena are in few cases in straight-line form. Many 2-parameter nonlinear models can be, however, easily transformed into straight-line form. In the log-log paper does the same for a power-law model. For other nonlinear models you should be able to figure 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); ”slope” and ”intercept” spreadsheet functions; linear regression in data analysis package.

38

4.5

Numerical solution of ordinary differential equations

Meaning Ordinary differential equation (ODE): all derivatives are with respect to single independent variable, often representing time. Equations with higher derivatives can be transformed into a system of first order equations. The ODE f 0 = f (t, y) is non-autonomous. The ODE f 0 = f (y) is autonomous. Solution embeds into the directional field. The ODE f 0 = f (t, y) does not by itself determine a unique solution. To single out a particular solution, one must specify a value y0 of the solution function at some point t0 . This is called initial value. Order of ODE determined by highest-order of derivative involved. Higher order ODEs may be ”narrowed down” also by using boundary values.

4.5.1

Basic methods for solving ODE

Euler Euler’s method uses the forward difference approximation for marized by the marching formula:

dy dx and

can be sum-

xi = xi−1 + h yi = yi−1 + h · f [xi−1 , yi−1 ] RK-4 This method for solving y 0 = f (x, ) , starting at (x0 , y0 ) is one of the most popular ones because it is quite accurate and stable. The formulae are: xi = xi−1 + h yi = yi−1 + h6 · (k1 + 2k2 + 2k3 + k4 ) where k1 k2 k3 k4

= f [xi−1 , yi−1 ] = f [xi−1 + h2 , yi−1 + h2 · k1 ] = f [xi−1 + h2 , yi−1 + h2 · k2 ] = f [xi−1 + h, yi−1 + h · k3 ]

Note that the ki values are various approximations to the slope of the unknown function. Heun Heun’s method improves on Euler’s method by a ”better” choice of the slope of the line followed from the point (xn , yn ) to the next point (xn+1 , yn+1 ). Instead of choosing n) this slope as dy(x = f (xn , yn ) as in the Euler method, choose it as the mean between dx this slope and the slope at the destination point reached by the Euler method. Simple Heun xi = xi−1 + h yi0 = yi−1 + h · f [xi−1 , yi−1 ]
yi1 = yi−1 + h2 · f [xi−1 , yi−1 ] + f [xi , yi0 ] 39

Iterated Heun (simplest predictor-corrector) yik = yi−1 +

h 2


· f [xi−1 , yi−1 ] + f [xi , yik−1 ]

continue until ”convergence”. This is a predictor - corrector (implicit, iterative) method. Needs a stopping criterion for inner iteration. It is a one-step, second order method and is stable. In practice we use higher order variants. Single-multi step, predictor-corrector, explicit-implicit So far we’ve considered self-starting (one step) methods. They use only yi . Multiple-step methods use more than one previously calculated value: that is to use yi and yi−1 etc. The traditional predictor-corrector solution methods use one equation to predict the value of yn+1 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 first equation. They provide some improvement especially for stiff problems, but need a jump-starter (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. Predictor-corrector methods are implicit. A ”stiff” 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.

40

Chapter 5

Methods related to multi-variable problems 5.1

Matrices, vectors

What’s an Array? What’s a Matrix? What’s 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 laid-out in either a straight line, or in a plane. 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 array. There is the further implication that the laws of Matrix Algebra dictate how an array can be combined with and/or operated on by other numerical arrays. Although the matrices encountered by engineering students are usually one-dimensional (numbers positioned in a straight line) or two-dimensional (numbers positioned in a plane), they can be higher-dimensional, as well. A subroutine to get the scalar product of two vectors Function scalprod(x() As Double, y() As Double) As Double Dim n As Integer, i As Integer n = UBound(x) scalprod = 0 For i = 1 To n scalprod = scalprod + x(i) * y(i) Next End Function

41

A subroutine to multiply two (compatible) matrices Sub matmult(a() As Double, b() As Double, c() As Double) ’multiplies matrices a anb b to get c Dim row As Integer, cola As Integer, rowb As Integer, col As Integer Dim i As Integer, j As Integer, k As Integer Dim s As Double row = UBound(a, 1) cola = UBound(a, 2) rowb = UBound(b, 1) col = UBound(b, 2) If Not (cola = rowb) Then MsgBox ("matprod: a and b are not compatible") Stop End If If (UBound(c, 1) < row Or UBound(c, 2) < col) Then MsgBox ("matprod: c is not large enough to contain the result") Stop End If For i = 1 To row For j = 1 To col s = 0 For k = 1 To cola s = s + a(i, k) * b(k, j) Next k c(i, j) = s Next j Next i End Sub

A subroutine to multiply matrix A by a (column) vector b: Sub matvecprod(a() As Double, b() As Double, c() As Double) Dim row As Integer, col As Integer Dim i As Integer, j As Integer Dim s As Double row = UBound(a, 1) col = UBound(a, 2) If ( row < Ubound(b) OR row < Ubound(c) ) Then MsgBox ("matvecprod: vector b or c is not compatible with matrix a") Stop End If For i = 1 To row s = 0 For j = 1 To col s = s + a(i, j) * b(j) Next j c(i) = s Next i End Sub

Matrix is lower triangular if all entries above main diagonal are zero: aij = 0 for i < j. Matrix is upper triangular if all entries below main diagonal are zero: aij = 0 for i > j.

42

Review linear dependence and rank

5.2

System of linear equations

Given m × n matrix A and m-vector b, and unknown n-vector x satisfying Ax = b. System of equations asks: Can b be expressed as linear combination of columns of A? If so, coefficients of linear combination given by components of solution vector x. Example 1 Consider the system 2x − 3y = 8 4x − 5y = 16 The determinant is non-zero. Two lines crossing at one point. Example 2 Consider the system 2x − 3y = 8 4x − 6y = 16 Thus the determinant is 2 · 6 − 3 · 4 = 0. And indeed the second equation is 2 times the first. Geometrically this corresponds to two coincident lines. In terms of rank: the coefficient matrix has rank 1, the augmented matrix has rank 1. Example 3 Consider the system 2x − 3y = 8 4x − 6y = 17 Thus the determinant is 2 · 6 − 3 · 4 = 0. There is no solution. Geometrically this corresponds to two parallel lines never crossing. In terms of rank: the coefficient 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 coefficient matrix has – it may be less than the total number of equations. If rank of the coefficient matrix is less than the number of equations, two cases can happen: • Consistent (if rank of coeff matrix = rank of augmented matrix): Infinite number of solutions • Non-Consistent (if rank of coeff matrix < rank of augmented matrix): No solution We can we express the solution with the inverse of the coefficient matrix, as x = A−1 b 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. 43

Naive Gauss elimination A system of linear equations can be represented in matrix forms of various types: e.g.       x1 + 2x2 − 3x3 = 5  1 2 −3 x1 5 2x1 − 3x2 + x3 = 3 becomes  2 −3 1   x2  =  3   −x1 + x2 − 2x3 = 0 −1 1 −2 x3 0 To solve the system a method (algorithm) 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 find A−1 ). 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) Interchanging two rows (b) multiplying or dividing a row by a non-zero constant (c) Replacing a row by the sum of that row and a multiple of another row (adding a negative multiple of a row – e.g. adding (−2) × row 2 to row 1) is the same as subtracting the positive multiple (subtracting 2× row 2 from row 1 – so that this row operation also includes subtraction). The objective is to change the coefficient matrix to a special form of matrix called upper triangular 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 non-zero in each row is a 1. The equations represented by this upper triangular form have exactly the same solutions but are much easier to solve. e.g. Solving the example above (a label like “R20 = R2 + (−2) R1 ” means replace row 2 by the sum of row 2 plus –2 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 row which is used to modify the other rows and is in the column where zeros are being created.   1 2 –3 5   0    2 –3 1 3  R2 =0 R2 + (−2) R1 → R 3 = R 3 + R1 –1 1 –2 0 

1  0   0

2 –7 3

–3 7 –5 

1  0 0

 5 –7    0 R 3 = R3 + 5 2 –7 0

–3 7 –2

→ −3 7




R2

 5 –7  2

This has created an upper triangular matrix (the first three columns). To further simplify this, divide the second row by –7 and the third row by –2, creating an echelon 44

form. The corresponding equations can easily be solved by starting with the last equation which gives the x3 value, then finding x2 from the next equation up and so on (back substitution):   x1 + 2x2 − 3x3 = 5 x1 = 5 + 3x3 − 2x2 = 2 Solve third 1 2 –3 5  0 1 –1 1  ⇒ x2 − x3 = 1 ⇒ x2 = 1 + x3 = 0 Solve second x3 = −1 x3 = −1 Solve first 0 0 1 –1 This version is called naive elimination, because we assumed that we can always divide by the diagonal element in the next row (e.g. 7 in iteration No two.) Naive Gauss-Jordan 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 GaussJordan(A() As Double) ’Augmented matrix Dim PivElt As Double, TarElt As Double Dim n As Integer ’number of equations Dim PivRow As Integer, TarRow As Integer ’pivot row; target row Dim i As Integer, j As Integer n = UBound(A, 1) For PivRow = 1 To n ’process every row PivElt = A(PivRow, PivRow) ’choose pivot element If PivElt = 0 Then MsgBox ("Zero pivot element encountered") Stop End If For j = 1 To n + 1 A(PivRow, j) = A(PivRow, j) / PivElt ’divide whole row Next For TarRow = 1 To n ’now replace all other rows If Not (TarRow = PivRow) Then TarElt = A(TarRow, PivRow) For j = 1 To n + 1 A(TarRow, j) = A(TarRow, j) - A(PivRow, j) * TarElt Next j End If Next TarRow Next PivRow End Sub

Answer the questions: Where is the coefficient matrix at input? Where is the ”right hand side”? At output: Where is the solution? What happened to the coefficient matrix? For both the naive Gauss and Gauss-Jordan 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 first example above.

45

5.3

LU Decomposition and Backsubstitution

An advanced form of the Gauss elimination is called LU decomposition-substitution. The procedure is divided into two parts: the first one operates on the coefficient matrix and does the computationally demanding half of the job; the second one (called substitution) uses the results of the first part to obtain a solution to a given right-hand side. Note that the decomposition destroys the original coefficient matrix (overwrites it by the L and U matrix elements.) An additional vector (Indx) is used to remember the row changes found necessary during ”partial pivoting”. The Indx vector then enters the substitution routine. (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 done, then the matrix is singular and the procedure stops with a message. (In other words, partial pivoting is enough, full pivoting does usually not improve the method.) The solution comes out from the substitution routine as x. Option Explicit Option Base 1 Sub VBA092() ’LU decomposition-substitution driver Dim n As Integer With Worksheets("Sheet1") .Cells(1, 1) = "Number of equations, n" End With With Worksheets("Sheet1") n = .Cells(1, 2) End With Dim i As Integer, j As Integer ReDim A(n, n) As Double, b(n) As Double, x(n) As Double ReDim Indx(n) As Integer With Worksheets("Sheet1") For i = 1 To n For j = 1 To n A(i, j) = .Cells(i + 1, j) Next j b(i) = .Cells(i + 1, n + 1) Next i End With Call LUDecompose(A, Indx) Call LUSubstitute(A, Indx, b, x) With Worksheets("Sheet1") .Cells(1, n + 3) = "Solution" For i = 1 To n .Cells(i + 1, n + 3) = x(i) Next i End With End Sub

It is a good practice to put the LUDecomposition and LUSubstitute into a separate module.

46

Option Explicit Option Base 1 Sub LUDecompose(A() As Double, Indx() As Integer) ’ input A(n,n) ’ output A(n,n) and indx(n) Const Tiny As Double = 1E-16 Dim n As Integer n = UBound(A, 1) Dim i As Integer, j As Integer, k As Integer, m As Integer Dim dum As Double For k = 1 To n - 1 m = k For i = k + 1 To n If (Abs(A(i, k)) > Abs(A(m, k))) Then m = i Next i Indx(k) = m dum = A(m, k) If m > k Then A(m, k) = A(k, k) A(k, k) = dum End If If Abs(dum) > Tiny Then dum = 1 / dum Else MsgBox ("Singular coefficient matrix") Stop End If For i = k + 1 To n A(i, k) = -A(i, k) * dum Next i For j = k + 1 To n dum = A(m, j) A(m, j) = A(k, j) A(k, j) = dum For i = k + 1 To n A(i, j) = A(i, j) + A(i, k) * dum Next i Next j Next k If Abs(A(n, n)) < Tiny Then MsgBox ("Singular coefficient matrix") Stop End If End Sub

The LUDecomposition overwrites the coefficient matrix by the decomposed (LU) form. You have to be careful not to call it more than once. Notice that the back-substitution can be called as many times as many right-handsides we have (as far as the coefficient matrix of the system is the same.)

47

Sub LUSubstitute(A() As Double, Indx() As Integer, b() As Double, x() As Double) ’ input A(n,n) containing the LU decomposition ’ input Indx(n) index vector ’ input b(n) right hand side ’ output x(n) solution Dim n As Integer n =UBound(A, 1) Dim i As Integer, j As Integer, k As Integer Dim dum As Double For i = 1 To n x(i) = b(i) Next i For k = 1 To n - 1 i = Indx(k) dum = x(i) x(i) = x(k) x(k) = dum For j = k + 1 To n x(j) = x(j) + A(j, k) * dum Next j Next k For k = n To 1 Step -1 x(k) = x(k) / A(k, k) For j = 1 To k - 1 x(j) = x(j) - A(j, k) * x(k) Next j Next k 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 above it (i.e. on the diagonals immediately below and above the main diagonal) a very effective algorithm is due to Thomas. e.g. Given the tridiagonal matrix shown store the three diagonals in a three column array as shown:   2 1 0 0 0  1 2 1 0 0  a: 0 1 2 1 2   → b: 2 2 4 3 4 0 2 4 0 0 A=    0 0 1 3 2  c: 1 1 0 2 0 0 0 0 2 4 The storage saving in the above example is only from 25 down to 15. However in larger matrices the reduction is very significant. For an n × n matrix with n2 entries the storage will be reduced to 3n. Hence if n = 1000 for a tridiagonal matrix, then the full storage is 1, 000, 000 and the reduced storage is 3000. 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.) The Thomas algorithm is extensively used in reservoir simulation: 48

Option Explicit Option Base 1 Sub Thomas(a() As Double, b() As Double, c() As Double d() As Double, x() As Double) ’Tridiagonal system of equations: ’a is subdiagonal, b is diagonal, c is super diagonal ’a(1) and c(n) are not used ’d is the right hand side, x is the result Dim n As Integer, i As Integer n = UBound(b) ReDim w(n) As Double, g(n) As Double w(1) = b(1) g(1) = d(1) / w(1) For i = 2 To n w(i) = b(i) - a(i) * c(i - 1) / w(i - 1) g(i) = (d(i) - a(i) * g(i - 1)) / w(i) Next i x(n) = g(n) For i = n - 1 To 1 Step -1 x(i) = g(i) - c(i) * x(i + 1) / w(i) Next i End Sub

Note that matrix operations and solution of linear systems are also provided by excel as spreadsheet functions, macros.

Direct vs Iterative Gauss elimination, Gauss-Jordan elimination, LU-decomposition/substitution are DIRECT methods. In direct methods we can predict the number of multiplications and additions if the number of equations (n) is known. Iterative methods are generalization of the direct iteration we learned for one variable x = g(x) equations. Iterative methods are used for very large but sparse systems (reservoir simulation). 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 Gauss-Seidel method. The method is illustrated for a general 3 × 3 system. The strategy is to express x1 from the first row, x2 from the second row, etc. Then start with a guess and calculate a new x1 from the first equation. Then with the best available approximants calculate x2 and so forth. In other words, as soon as a variable is updated in one equation its new value is used on the right hand side of the next equation:   a11 x1 + a12 x2 + a13 x3 = b1  a11 x1 = b1 − a12 x2 − a13 x3  a21 x1 + a22 x2 + a23 x3 = b2 ⇒ a22 x2 = b2 − a21 x1 − a23 x3   a31 x1 + a32 x2 + a33 x3 = b3 a33 x3 = b3 − a31 x1 − a32 x2  a11 xk+1 = b1 − a12 xk2 − a13 xk3  1 k+1 k ⇒ a22 xk+1 = b − a x − a x 2 21 1 23 3 2  a33 xk+1 = b3 − a31 xk+1 − a32 xk+1 3 1 2 It works fine on diagonally dominant matrices. 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 row. Identity matrix is diagonally dominant 49

5.4

System of nonlinear equations, Newton-Raphson method

The method may be used to solve a general equation f (x) = 0, provided the left-hand-side can be differentiated ”analytically”. Note that bold-face denotes matrix and vector. We have n equations to solve and n components of the unknown to find. If an initial guess is known, and at that location we can calculate the function f (x0 ) and its partial derivatives collected into the Jacobian matrix J(x0 ). The first element in the first row of the Jacobian matrix will contain the partial derivative of f1 with respect to x1 , the second element in the first row will contain the partial derivative of f1 with respect to x2 , and so on. Writing the equation of the tangent plane: f (x0 ) + J(x0 )(x − x0 ) = 0 Substituting y = 0 we find that the line crosses the x axis at x1 where −1

(x1 − x0 ) = ∆x = [J(x0 )]

f (x0 )

This is the iteration formula of the Newton’s method. The obtained x1 can be then renamed to x0 and hence we can continue the iterations. In practice we do not use the matrix inverse, we just solve J(x0 )∆x = −f (x0 ) for ∆x using LU Decomposition. Cautious Approach: We calculate the ∆x and multiply it by a less than one α 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 Newton-Raphson method? Two subroutines, one to evaluate f and the other for the Jacobian matrix, J. N/R Pros and Contras: Second order convergence Needs ”analytical” 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 (”errors”). 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 defined unit, and ”eps” has to correspond to it. If the elements of x are of different dimensions or vary several orders of magnitude, use normalization N/R variants: If the Jacobian is not available analytically, we can estimate it by finite 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 INVERSE OF THE JACOBIAN with every new function evaluation 50

Newton-Raphson method example: Solve f1 = x1 + x22 − 3 f2 = xx12 + 2 We will use a VBA program. The input is from the spreadsheet: • kmax from Cells(1, 2), maximum number of iterations, (20) • eps from Cells(2, 2), tolerance, (1.0E-6) • alpha from Cells(3, 2), α, ( 1 ) • starting point x1 and x2 , from Cells(4, 2) and Cells(5, 2), e.g. (1 and 1) The program contains a main (driving) subroutine and the two user subroutines. Note that it also needs some auxiliary subroutines (linear solver: LU decomposition and back-substitution) presented previously. The linear solver pair can be placed into another module.

51

Option Explicit Option Base 1 Sub VBA101() ’Newton-Raphson Method Example Dim n As Integer n = 2 ReDim x(n) As Double, f(n) As Double, J(n, n) As Double ReDim b(n) As Double, dx(n) As Double Dim alpha As Double, eps As Double, s As Double ReDim Indx(n) As Integer Dim k As Integer, kmax As Integer, i As Integer With Worksheets("Sheet1") .Cells(1, 1) = "kmax" .Cells(2, 1) = "eps" .Cells(3, 1) = "alpha" .Cells(4, 1) = "x1 start" .Cells(5, 1) = "x2 start" End With With Worksheets("Sheet1") With Worksheets("Sheet1") kmax = .Cells(1, 2) eps = .Cells(2, 2) alpha = .Cells(3, 2) x(1) = .Cells(4, 2) x(2) = .Cells(5, 2) End With With Worksheets("Sheet1") With Worksheets("Sheet1") .Cells(1, 4) = "Iteration" .Cells(1, 5) = "x1" .Cells(1, 6) = "x2" .Cells(1, 7) = "f1" .Cells(1, 8) = "f2" End With ’iteration loop For k = 1 To kmax Call fun(x, f) For i = 1 To n b(i) = -f(i) Next i Call Jacobian(x, J) Call LUDecompose(J, Indx) Call LUSubstitute(J, Indx, b, dx) s = 0 For i = 1 To n s = s + dx(i) ^ 2 x(i) = x(i) + alpha * dx(i) Next i With Worksheets("Sheet1") .Cells(k + 1, 4) = k .Cells(k + 1, 5) = x(1) .Cells(k + 1, 6) = x(2) .Cells(k + 1, 7) = f(1) .Cells(k + 1, 8) = f(2) End With If (Sqr(s) < eps) Then Exit For Next k End Sub Sub fun(x() As Double, f() As Double) f(1) = x(1) + x(2) ^ 2 - 3 f(2) = x(1) / x(2) + 2 End Sub Sub Jacobian(x() As Double, J() As Double) J(1, 1) = 1: J(1, 2) = 2 * x(2) J(2, 1) = 1 / x(2): J(2, 2) = -x(1) / x(2) ^ 2 End Sub

52

5.5

Multivariable extrema

f (x) → 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 loeast-squares fit Equivalence (for unconstrained and differentiable) ”in principle” equivalent: f (x)− > min and gradient = 0 (Just to solve ”gradient = 0” is somewhat dangerous, it is better watch f(x) ) 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 ”simplex” or ”polytope” is a cluster of n+1 points, for instance for two-variable functions a simplex consists of three points, forming a triangle. The Nelder-Mead simplex method (method of polytopes) uses the idea of a ”wandering” 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 ”iteration” 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 fit) Generalized Linear Least Squares: after transform of variables straight line fit (max 2 parameters) y = mx + b more than two variables but still linear y = p1 x1 + p2 x3 + p3 x1 x3 because the parameters p1 , p2 , p3 go into the model linearly. Nonlinear least squares: several parameters p1 , p2 , p3 . . ., the dependence of the response is nonlinear on them y=

p1 x1 1 + p2 x1 + p3 x2

53

Observations: (x1,1 , x1,2 , y1 ) (x2,1 , x2,2 , y2 ) . . . (xn,1 , xn,2 , yn ) The Objective function is:       p1 x1,1 p1 x2,1 p1 xn,1 Q = y1 − + y2 − +. . .+ yn − 1 + p2 x1,1 + p3 1, 2 1 + p2 x2,1 + p3 x2,2 1 + p2 xn,1 + p3 xn,2 The minimization problem is: find p1 , p2 , p3 which minimizes the objective function Q, is usually solved by a special method called Gauss-Newton-Marquardt. Example

Given the observations:

x1 1 3 5 7

x2 2 4 4 2

y 1.154 2.045 3.125 5.526

The estimated parameters are: pb1 = 3.26, pb2 = 0.225, pb3 = 0.777 Excels all-purpose ”Solver” is also an efficient tool to minimize the objective function (even when additional constraints have to be taken into account.)

5.6 Solution of system of ordinary differential equations: RK4 A system of first order equations has one function variable (the variable being differentiated) for each equation in the system and all of these function variables are functions of one independent variable such as x or t. The Initial Value Problem (IVP) is of the form  dy1 dx = f1 (x, y1 , y2 ) dy2 dx = f2 (x, y1 , y2 ) with initial conditons, for some known a, k1 , k2 : y1 (a) = y10 , y2 (a) = y20 where the right hand side functions f1 , f2 are any functions of the independent variable x and the two function variables y1 , y2 (but not of the derivatives of these). A two equation system is shown above, but it can be extended in the obvious way to more differential equations. The algorithm now involves vectors (watch for boldface): xi = xi−1 + h yi = yi−1 + h6 · (k1 + 2k2 + 2k3 + k4 ) where k1 = f [xi−1 , yi−1 ] k2 = f [xi−1 + h2 , yi−1 + h2 · k1 ] k3 = f [xi−1 + h2 , yi−1 + h2 · k2 ] k4 = f [xi−1 + h, yi−1 + h · k3 ] The vector operations can be done by simple loops. 54

5.7 Partial differential equations and reservoir simulation Review Physical meaning (Often from conservation principles) Initial and boundary conditions. Transient versus steady-State With initial-value 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 Classification – Heat equation: ut = uxx (parabolic) – Wave equation: utt = uxx (hyperbolic) – Laplace equation: uxx + uyy = 0 (elliptic) Numerical approaches Finite difference, finite element, Semi-analytic (transformation) 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 find 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. η

∂2p ∂p = ∂x2 ∂t

where η= φ µ ct k

k φµct

porosity viscosity compressibility (total) permeability

We want to solve the transient problem: The evolution of the pressure field, if the initial condition is constant pressure pi and the pressures at the two ends are kept constant, plef t and pright , respectively.

55

First discretize the region on interest over both space: x1 , x2 , . . . , xi , . . . , xnx and time:t0 , t1 , . . . , tn , . . . , tnt . Replace the analytical derivatives by numerical approximations. Approximation of the second partial derivative in the x-direction: (central): pi+1 − 2pi + pi−1 (∆x)2 Approximation of the first partial derivative in time between time levels n and n+1: pi n+1 − pi n ∆t So far we have not stated what timestep (n or n+1) the terms on the left are being evaluated. An ”implicit” solution scheme means that the pressures at the nodes will be evaluated at the new timestep (n+1) in the second derivative (space) term. η

n+1 pn+1 + pn+1 pi n+1 − pi n i+1 − 2pi i−1 = 2 (∆x) ∆t

It is useful to introduce a new notation α=

1 (∆x)2 η ∆t

Then the system of equations become: pn+1 = given plef t 1 n+1 n pn+1 + pn+1 i+1 − (2 + α)pi i−1 = αpi for i = 2, 3, ... n-1 n+1 pn = given pright This is a tridiagonal system with the unknowns: n+1 pn+1 , pn+1 , . . . pn+1 1 2 n−1 , pn

In matrix representation  b1  a2   0   0 0

(for n = 5, as an example):     c1 0 0 0 p1    b2 c2 0 0    p2    p3  =  a3 b3 c3 0      0 a4 b4 c4   p4   0 0 a5 b5 p5

d1 d2 d3 d4 d5

     

For instance, b1 = 1, c1 = 0 and d1 is nothing else, but the given plef t . The system can be solved by the Thomas algorithm.

56

Example data porosity viscosity compressibility permeability reservoir length initial pressure pressure at x = 0 ft pressure at x =3000 ft time

ct k L pi plef t pright t

0.2 0.8 1 10−5 5 3000 2800 1000 2800 0-100

cp 1/psi md ft psi psi psi days

With the units indicated, η= and hence

 α=

0.00633k φµct

φµct 0.00633k

Choose appropriate nx and nt.

57



(∆x)2 ∆t

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 Double, 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 i = 1 With Worksheets("sheet1") .Cells(i, 1) = "Porosity (fraction): ": phi = .Cells(i, 2) i = i + 1: .Cells(i, 1) = "Viscosity (cp): ": mu = .Cells(i, 2) i = i + 1: .Cells(i, 1) = "Compressibility (1/psi): ": ct = .Cells(i, 2) i = i + 1: .Cells(i, 1) = "Permeability (md): ": k = .Cells(i, 2) i = i + 1: .Cells(i, 1) = "Length of reservoir (ft): ": xlen = .Cells(i, 2) i = i + 1: .Cells(i, 1) = "Initial pressures (psi): ": pini = .Cells(i, 2) i = i + 1: .Cells(i, 1) = "Left pressure (psi): ": pleft = .Cells(i, 2) i = i + 1: .Cells(i, 1) = "Right pressure (psi): ": pright = .Cells(i, 2) i i i i

= = = =

i i i i

+ + + +

1: 1: 1: 1:

.Cells(i, .Cells(i, .Cells(i, .Cells(i,

1) 1) 1) 1)

= = = =

"Number of grid points in x: ": nx = .Cells(i, 2) "End time (day): ": tend = .Cells(i, 2) "Number of time steps: ": nt = .Cells(i, 2) "Print at every: ": ipr = .Cells(i, 2)

ReDim a(nx) As Double, b(nx) As Double, c(nx) As Double ReDim d(nx) As Double, p(nx) As Double dx = xlen / (nx - 1) dt = tend / nt alpha = phi * mu * ct / (0.00633 * k) * dx ^ 2 / dt ’a For i = 2 To nx - 1: a(i) = -1: Next i: a(nx) = 0 ’b b(1) = 1: For i = 2 To nx - 1: b(i) = 2 + alpha: Next i: b(nx) = 1 ’c c(1) = 0: For i = 2 To nx - 1: c(i) = -1: Next i ’d d(1) = pleft: d(nx) = pright ’Initialize t = 0 For i = 1 To nx: p(i) = pini: Next i ip = 1 .Cells(1, 6) = "Loc (ft)->" For i = 1 To nx: .Cells(1, 6 + i) = (i - 1) * dx: Next i .Cells(ip, 4) = "time (day)": ip = ip + 1 ’ Marching... For it = 1 To nt .Cells(ip, 6) = "p (psi)->" t = it * dt For i = 2 To nx - 1: d(i) = alpha * p(i): Next i Call Thomas(a, b, c, d, p) If ((it Mod ipr) = 0) Then .Cells(ip, 4) = t For i = 1 To nx: .Cells(ip, 6 + i) = p(i): Next i ip = ip + 1 End If Next it End With End Sub

58

Note that this ”simulator” 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 solution algorithm. Experiment with the time step. Verify the following statement: Stability Implicit: stable Explicit: at α less than 2 it is unstable!

5.7.1

Reservoir Simulation

Gridblock (Finite Volume) approach (Two dimensional: x and y direction) Oil Material Balance Equation Oil in place in block i = N=

V p So Bo

where Vp = ∆x∆yhφ Oil in place changes from 

V p So Bo

to 

V p So Bo

n

n+1

during timestep n+1. Thus, the rate of accumulation is 1 ∆t

"

V p So Bo

n+1

 −

Vp So Bo

n #

Flow from block i+1 to block i (EAST direction): qi+ 12 =

0.00633kkro h∆y µo Bo



pi+1 − pi ∆x



 = TE

kr µB

 (pi+1 − pi ) oE

A similar equation can be written from the WEST direction. The net inflow is:  TW

kr µB



 (pi−1 − pi ) + TE oW

kr µB

 (pi+1 − pi ) oE

Including a possible well term qo and the previously determined accumulation term we obtain the gridblock (”cell”) balance

59

1 ∆t



V p So Bo

n+1

 −

V p So Bo

n !

1 + qo = ∆t



V p So Bo

n+1

 −

V p So Bo

n !

Flow term can be implicit (n + 1) or explicit (n). We have similar equations for water and gas. The most important unknowns are the cell pressures. Often we solve an implicit scheme for the pressures, then update everything else (saturations) by an explicit scheme (IMPES). Special issues: Well models, Fully implicit schemes, various solvers, ADI (method of alternate directions)

60

Chapter 6

Integral Transforms, Spectral methods, inversion of the Laplace transform 6.1

Transformations

Function to function (Laplace) Data set to data set (Discrete Fourier Transform) Spectral methods (Seismics, Ultra-Sound etc) FFT algorithm

6.2

Laplace transform

Forward and Inverse Transform Some Rules Analytical integration or numerical solution?

6.2.1

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 F(s)=1/s. The result should be the unit step function (one fore every positive time.) Check if you get back the value 1.0 for the given time points. • Run it. Modify the program to calculate the integral (from zero to t) of the unit step function, by dividing the Laplace transform by s.

61

• 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.) • Now calculate the integral of the function using the inverse Laplace transform method. Check the value of the integral. F(s) 1 s+0.5 1 √ s

arctan 0.1 s

f(t) e−0.5t √1 πt 1 sin 0.1t t

The main program first ”initializes” the algorithm with StInit. The parameter is the number of terms (n) in the Stehfest algorithm. 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 Fun(s As Double, F As Double) F = 1 / s End Sub Sub VBA13b() Dim i As Integer, np As Integer Dim t As Double, ft As Double Call StInit(16) With Worksheets("sheet1") .Cells(1, 1) = "number of time points: " np = .Cells(1, 2) .Cells(1, 3) = "t" .Cells(1, 4) = "f(t)" For i = 1 To np t = .Cells(i + 1, 3) Call StCalc(t, ft) .Cells(i + 1, 4) = ft Next i End With End Sub

It is a good programming practice to put the following Stehfest procedures (initialization and calculation) into a separate module:

62

Option Explicit Private Stn As Integer Private Stc() As Double Function Min(i1 As Integer, i2 As Integer) If i1 < i2 Then Min = i1 Else Min = i2 End Function Function Max(i1 As Integer, i2 As Integer) If i1 > i2 Then Max = i1 Else Max = i2 End Function Sub StInit(n As Integer) Stn = Max(2, Min(16, 2 * Int(n / 2))) ReDim Stc(Stn) 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 Fact(0) = 1 For i = 1 To Stn Fact(i) = i * Fact(i - 1) Next i For i = 1 To Stn ck = 0 For k = Int((i + 1) / 2) To Min(i, n2) sk = k ck = ck + sk^n2 *Fact(2*k)/(Fact(n2-k)*Fact(k)*Fact(k-1)*Fact(i-k)*Fact(2*k-i)) Next k Stc(i) = (-1) ^ (i + n2) * ck Next i End Sub Sub StCalc(t As Double, ft As Double) Dim ln2pt As Double, s As Double, F As Double Dim i As Integer ln2pt = Log(2#) / t ft = 0 For i = 1 To Stn s = i * ln2pt Call Fun(s, F) ft = ft + Stc(i) * F Next i ft = ln2pt * ft End Sub

63

Part IV

Appendix

64

Chapter 7

Study guide 7.1

Be able to ...

1. Compare scope of module level and procedure level variables 2. Use arrays 3. Write meaningful comments to code. 4. Subdivide a programming project into smaller units 5. Use defensive programming strategies to check for obvious errors in input parameters. 6. Use the interactive debugger to trace logic and numerical errors in code 7. Give a simple explanation of the term “floating point number”. 8. Give a definition of “roundoff error”. 9. Give a simple example of overflow. 10. Identify the truncation error from a Taylor series expansion. 11. Identify (at least) two important differences between symbolic and numeric computations. 12. Estimate the error of a simple calculation if errors associated with the input are given 13. Write any scalar equation in the form f (x) = 0 14. Define multiple root. 15. Explain the role of bracketing. 16. Write a simple logical expression that expresses the condition for an interval to contain (bracket) the root, 17. Manually perform a few steps of the bisection method. 18. Identify situations that cause Newton’s method to fail. 19. Calculate the number of steps for the bisection method to reduce the uncertainty interval x-times 65

20. Qualitatively compare the convergence rates of bisection, false-position and Newton’s method 21. Plot any f (x) as a means of graphically identifying the location of roots. 22. Describe the role of global variables in finding the roots of f (x, a, b, . . .) = 0 where a, b, . . . are parameters, and the method returns the value of x that gives f = 0 for fixed values of a, b, . . . 23. Compare effects of noise on numerical differentiation and integration 24. Derive the approximation order of a simple finite difference formula 25. Manually perform multiple application of Simpson’s rule 26. Manually perform Richardson extrapolation 27. Recognize if a least squares problem is linear, can be linearized, is nonlinear 28. Be able to determine degrees of freedom 29. Identify connection between linearized and original model parameters 30. Identify the criteria for obtaining a least squares fit. (What is minimized?) 31. Manually evaluate the formulas for the slope and intercept of a line fit to data. 32. Derive the transformations for fitting data to non-linear two-parameter functions of various type. 33. Make the connection between slope and intercept of the linearized model with the original model parameters. 34. Recognize initial value versus boundary value problem 35. Manually calculate a few steps of Euler’s and Heun’s methods 36. Program the RK-4 method 37. Describe concepts: single-multi step (self starting or not), predictor-corrector 38. Express conditions of linear dependence of vectors (row or column) in terms of matrix rank and vice-versa. 39. Explain the condition of singularity of an n · n matrix in terms of linear independence. 40. Express matrix rank as a measure of linear independence. 41. Write the formal solution to Ax = b. 42. Explain why it is not a good idea to use the formal solution as a computational procedure for solving Ax = b. 43. Name the solution algorithm most commonly used for solving Ax = b. 44. Describe the reason for pivoting. 45. Answer the question: Is pivoting a remedy for ill-conditioned systems? 46. Write the preferred expression for solving Lx = b or U x = b when L is lower triangular and U is upper triangular. What algorithm makes use of the solution for these systems? 47. Convert a higher order ODE to an equivalent system of coupled first order ODEs.

66

48. Use Newton-Raphson method for a system of nonlinear equations 49. Recognize the basic PDE types 50. Relate transient to steady state 51. Name methods to solve PDEs numerically 52. Solve the transient diffusivity equation by implicit finite difference 53. Be able to use the Thomas algorithm 54. Describe basic approaches to Reservoir simulation 55. Discuss explicit versus implicit approach 56. Work with units involved 57. Discuss relation between well scheduling and boundary conditions 58. Be able to use a simulator and interpret results 59. Have an overview of spectral methods 60. Be able to use an available module (in this case for numerical Laplace inversion)

67

Chapter 8

Assignments, test problems 8.1

Error propagation: borehole volume

A 6000 ft (± 1 %) deep borehole has a diameter: 5 inch (± 0.04 inch). What is the volume in bbl? (Give your answer with error estimate!)

8.2

Error propagation: one-gallon cube

What is the side length of a one-gallon (±5 %) cube? Indicate the unit and error estimate of your answer!

8.3

Error propagation: barrel

Consider a right circular cylinder that has a volume 1 BBL (British barrel) (±0.06 %). The height is exactly twice the diameter. Estimate the diameter, indicating the unit and the uncertainty.

8.4

Error propagation: hydrocarbon volume

The volume occupied by hydrocarbons in a certain reservoir can be calculated from VHC = V φ(1 − Sw ) where V is the reservoir volume, φ is the porosity, and Sw is the water saturation. Assuming V φ = 2.5 105 f t3 (±15%) and Sw = 0.65 (± 0.12) calculate the hydrocarbon 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 error!

68

8.5

Root finding; Newton iteration: z factor

An equation of state written in ”z-factor form” results in the following nonlinear equation to be solved: z 4 + 2.500z − 1.300 = 0 You use Newton’s method to find the root, starting from z = 1.000 (we will call it the zero-th approximation of the root.) Carry out the first iteration step and give the first approximation of the root. (Do not continue the iteration process!)

8.6 Root finding: solution of the PR equation of state he Peng-Robinson equation is one of the cubic equations of state:   at p+ (VM − b) = RT VM (VM + b) + b(VM − b) 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 defined as follows: b = 0.07780

RTC pC

aC = 0.45724

R2 TC2 pC

h i2 p at = aC 1 + (0.37464 + 1.54226ω − 0.26992ω 2 )(1 − T /TC ) It is somewhat easier to work with dimensionless quantities such as the deviation M factor z = pV RT . Then the equation takes the form         bp at 2b 3b2 p2 at b b2 b3 p3 2 z3+ − 1 z2+ p − − z+p − + + 3 3 =0 2 2 2 2 3 3 2 2 RT R T RT R T R T R T R T or simply : z 3 + c2 z 2 + c1 z + c0 = 0 where the coefficients can be read of from the previous expression. A substance is characterized by three properties: critical temperature, TC , critical pressure, pC and acentric factor, ω.

69

For n-butane: TC = 425.2 K pC = 3.75 MPa ω = 0.193 What is the deviation factor, z (and the corresponding molar volume) of butane at T = 373.15 K and p = 1.5 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 o C is 1.53 MPa, so the actual pressure given is less than the gas-liquid equilibrium pressure. Therefore, at 1.5 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.)

8.7

Root finding: heat balance

The enthalpy of 1 kg methane can be calculated from the following polynomial: h = −4662.+1.983·T −0.001031·T 2 +4.243·10−6 ·T 3 −2.946·10−9 ·T 4 +7.218·10−13 ·T 5 where h is the (specific) enthalpy in kJ/kg and T is the absolute temperature in K. We wish to determine the temperature of 5 kg methane in the final state, if the initial state is 300 K and the enthalpy is increased by 12,000 kJ.

8.8

Root finding: friction factor

The Darcy/Moody friction factor (f) is a function of the Reynolds number, Re and the relative roughness  / D . If the Reynolds number is less than (or equal to) 2100, then the flow 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    /D √ equation (called Colebrook equation): √1 = 2.0log10 3.7 + 2.51 f

Re

f

Write a complete VBA function that has two inputs: the values of Re and /D. It should calculate the friction factor, f . Hint: For bracketing methods choose f as unknown, for open methods choose ln(f ). Why?

70

8.9

Root finding: Well deliverability

The Inflow Performance Relationship (IPR) gives the production rate depending on the reservoir pressure and wellbore flowing pressure. It characterizes the ability of the reservoir to deliver. For instance:  0.62 qIP R = 4.38 p2 − p2wf where the production rate q is in MSCF/D and both the average reservoir pressure ( p ) and wellbore flowing pressure ( pwf ) are in psi. The Tubing Performance Relationship (TPR) gives the production rate depending on wellbore flowing pressure. It characterizes the ability of the wellbore to provide vertical lift. For instance:  0.5 qT P R = 93.68 p2wf − 106 From the TPR it can be seen, that (in this case) the minimum well bore flowing pressure to provide natural lift is 103 psi. In steady-state production qIP R = qT P R . Find the production rate at various reservoir average pressures. Hint: first rearrange the equation into ”f(x)=0” form.

8.10

Root finding: isoterm flash

n the simplest isoterm flash calculation we know the inlet composition of the mixture, z1 , z2 ,..., zn (mole fractions) and the distribution factors K1 , K2 ,..., Kn . We wish to find the mole fraction of vapor, V and the composition of the vapor phase ( yi ) and the composition of the liquid phase ( xi ). The equlibrium relations are: Ki =

yi xi

The material balance relations are: yi V + xi (1 − V ) = zi Rachford combined all these equations into one nonlinear equation with one unknown: n X i=1

=

(Ki − 1)zi =0 (Ki − 1)V + 1

Once solved for V the composition of the individual phases can be calculated from yi =

K i zi (Ki − 1)V + 1

xi =

zi (Ki − 1)V + 1 71

Solve the problem for the following input data: Component C1 C2 C6

8.11

Ki 61.0 9.0 0.035

zi 0.3396 0.0646 0.6068

Optimization: Maximizing Net Present Value

One MSCF additionally produced gas brings 3.5 US$ revenue. 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 benefits. Let us consider matrix acidizing. The size of the treatment is best characterized by the volume of acid spent. All costs depend on this extensive variable. The benefits 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 ”less valuable” than the dollar we spend today, hence it has to be discounted depending on the year it is realized in. The discount ratio is 1.2 . 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$ /gal. Fixed cost is 10,000 US$ for a single treatment. Total cost (in US$) depends on the acid spent: v (in gal): cost = 10, 000 + (34 + 9)v Benefit comes from increased production rate. The production rate in MSCF/D is given by 30, 000 q= 6.3 + s where v is the acid spent, and s is the skin factor modified by the acid treatment. The pre-treatment skin factor is 7. Depending on the injected volume, the skin factor decreases according to the relation 7 1 + 0.1v 1/3

s= (but never becomes negative.)

The discounted benefit for a 3-year horizon is therefore expressed (in US$) as 3 X 3.5 benef it = 365 i 1.2 i=1

30, 000 30, 000 − 7 6.3 + 7 6.3 + 1+0.1v 1/3

Find the optimum amount of acid to spend on this well! 72

!

8.12

Integration of a function: PV

The production rate q (STBD) from an oil well declines according to the Hyperbolic Decline Law: qin q= 1 (1 + Knt) n where 1 K = 0.00207 day is the decline constant n = 0.3 is the Arp exponent qi =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, c ($/STB). A ”dollar tomorrow is, however, less than a dollar today” so we have to discount the revenue using a yearly discount rate, i (% /Year) and a discounting method (which is ”continuous” in this case.) c = 25 $/STB i = 12 % /Year Using the continuous discounting method the Present Value of one barrel oil produced in the future, t (days) from now, can be expressed by i

ce− 365 t The Present Value of the cumulative production from ”now” to abandonment can be calculated from the integral Z PV =

ta

i

ce− 365 t

0

qin 1

(1 + Knt) n

dt

8.13 Numerical integration of discrete data: pore volume A reservoir of area, A = 3208 acre ( 1.3974 108 ft2 ) has a pay thickness of 40 ft. The pore space is completely filled with oil. The porosity was measured as at various depths from the top of the formation: Depth, ft 0 10 20 30 40

porosity, φ 0.2916 0.3265 0.3187 0.3193 0.2854 73

Calculate the oil in place (in reservoir cubic feet) using a) multiple application of the trapezoid rule and b) multiple application of Simpson’s 1/3 rule.

8.14 1

Straight-line: formation volume factor model

Given: pb = 2012 psi, bubble point pressure Data (observed): p psi Bo resBBL/STB 1500 1.262 1600 1.279 1800 1.298 Determine the parameters of the nonlinear model describing the Bo: Bo =

1 √ a + b pb − p

What is the best estimate of the Bo at the bubble point?

8.15 2

Straight-line: formation volume factor model

Consider the following model of Formation Volume Factor, Bo as a function of pressure, p. Bo = aeb(p−pb ) where Bo is in resBBL/STB, p in psi, and pb is the known bubble point pressure: (pb = 3007 psi). The model parameters a and b are to be find. The following laboratory data are available: p psi 500 1500 2500

Bo resBBL/STB 1.070 1.175 1.301

Determine the Formation Volume Factor at the bubble point (pb ) using the above model.

74

8.16

Straight-line: gas in place

Production and static (field) pressure data for a gas field is given below. (Craft and Hawkins) p/z, psia 6553 6468 6393 6329 6246 6136 6080

Gp (MMM SCF) 0.393 1.642 3.226 4.260 5.504 7.539 8.749

The material balance for a volumetric gas reservoir is   p pi pi = − Gp z zi zi G i where p (psia) is static field pressure, z (-) is deviation factor for corresponding to the pressure, Gi (MSCF) original (initial) gas in place, Gp produced gas (MSCF) and the subscript i refers to initial. Find original gas in place from the above data using straight-line form and least-squares fit! Hint: from the straight line form estimate slope and intercept, then interpret them to obtain first pi /zi then Gi . Alternative linear form:  Gp = G −

8.17

zi G i pi



p z

Straight-line: Flow-After-Flow Test (IPR)

A frequently used IPR equation: " q = qmax 1 − where q production rate, BOPD p average reservoir pressure, psi pwf welbore flowing pressure, psi Data (observed): pwf psi q BOPD 2500 109.9 2000 192.5 1500 258.7 1000 309.2

75



pwf p

2 # n

Find the Absolute Open Flow Potential: Hint: fill out the following table first.

Parameter 1 Parameter 2 Independent variable Dependent variable

Original model qmax n p q

Linearized model m b x y

8.18 Nonlinear least squares: oil viscosity as a function of pressure and temperature Consider the following model of oil viscosity (µo ) for a certain field as a function of pore pressure, p and layer temperature T :
µo = a · T b 1 + c · p1.5 where µo is in cp, p in psia, T in o R. The model parameters a, b and c are to be determined by the method of nonlinear least squares using a general purpose minimization program (e.g, Solver). The available data are: p, psi 500 1500 2500 1500 500

T, o R 460 480 500 520 540

µo , cp 0.764 1.24 2.19 1.02 0.833

Write module containing a VBA function that calculates the objective function to be minimized. The first line of the function should read as: function SSQ(a As Double, b As Double, c As Double) As Double Use global (module-level) variables for all data (except for the three model parameters: a, b, c which are arguments of the function.)

8.19 ODE: Production rate decline from a differential equation Consider the following differential equation describing the production rate decline of a well:

76

f racdqo dt =

0.00025 · q 1.5 1 + 0.1 · t0.1

where qo is the production rate in BOPD and t is time in days. The initial production 0.00025 rate is 498 BOPD. (The factor 1+0.1·t 0.1 describes the decrease of flow efficiency with time that is a consequence of evolving damage.) Use the RK4 method with step-size h: 100, 100, 100, 65 days to calculate the production rate at t = 365 days. (Advanced: derive and solve another ODE for the cumulative production.)

8.20

ODE: pressure versus depth in a shut-in well

A gas well is shut-in at the wellhead (no flow). Write a differential equation for the pressure assuming the deviation factor is a known function of pressure. The temperature can be considered constant. Density Molar mass of gas Temperature Deviation factor Pressure at the surface

ρ=

M V

z=

p = M RT M T 1 1+c1 ·p

p0

kg/m3 0.27 kg/mol 375 K c1 = 3 · 10−6 P1a 3.2 106 Pa

The solution of the differential equation should be a function: providing pressure (p, Pa) as a function of depth (D, m). 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 fluid with density ρ = ρ0 (1 + cf p) where ρ0 = 235 kg/m3 and cf = 0.3·10−6 P a−1 and the pressure at the surface is as before: p0 = 3.2 106 P a .

8.21 Meaning of matrices and vectors: chemical component balance (Right a conservation equation in matrix-vector form, for instance for a distillation column.)

8.22

System of linear equations: mixing

Three petroleum products P1 , P2 and P3 are available. The Octane number determined by two independent methods: research (RO) and motor (MO) are shown in the following table: 77

RO MO

P1 83.5 84.3

P2 91.2 88.7

P3 96.1 97.9

You wish to create a mixture which is of Octane number (RO+MO)/2 = 89 determined by any of the methods. The question is what mass fractions to use (if there is a solution at all)? It is reasonable to assume that the resulting RO 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 fluids 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 definition 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?

8.23 System of linear equations with many righthand-sides: log interpretation A basic approach to determine lithology of simple formations is neutron-density crossplot. Let V1 , V2 and V3 denote the volume fractions of sand, clay and pore space (which is assumed to be filled by water). The sum of the volumes fractions is, of course, unity. Let ΦN denote the neutron value (dimensionless) and ρ the density in (g/cm3 ). (Note: the official SI density unit is kg/m3 . As we all know 1 g/cm3 = 1000 kg/m3 ) The log evaluation is based on the simultaneous application of the two relations: ΦN = −2V1 + 37V2 + 100V3 ρ = 2.65V1 + 2.41V2 + 1.0V3 expressing that both the neutron value and the density is a linear combination of the volume fractions. For instance, if V1 = 0.76, V2 = 0.05 and V3 = 0.19, then ΦN = 19.33 and ρ = 2.3245 g/cm3 . If we could measure ΦN and ρ, then we could find out the volume fraction, solving the system −2V1 + 37V2 + 100V3 = ΦN 2.65V1 + 2.41V2 + 1.0V3 = ρ V 1 + V 2 + V3 = 1 Since logs are taken densely (say meter by meter) we need to solve for the unknown V1 , V2 and V3 with many right-hand-sides. 78

Write a program to do the calculation calling the LUDecomposition (once) and the LUsubstitute (many times.) Depth, ft 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010

ΦN 19.3 16.56 15.89 10.54 19.95 19.74 18.3 27.6 23.06 14.49 30.

ρ, g/cm3 2.32 2.37 2.38 2.47 2.33 2.32 2.34 2.19 2.25 2.40 2.15

(Additional challenge: Using the LU solver pair, calculate the inverse of the coefficient matrix. Check if the original matrix multiplied by the obtained inverse is the identity matrix.)

8.24 Nonlinear System of Equations: Chemical Equilibrium of Methan Combustion Consider the following two chemical reactions CH4 + H2 O = CO + 3H2 CO + H2 O = CO2 + H2

ln Ky1 = 3.235 ln Ky2 = 0.3726

Notation: y means mole fraction. The sum of mole fractions is 1. Ky is the equilibrium constant (expressed with mole fractions, dimensionless). ln Ky1 = − ln yCH4 − ln yH2 O + ln yCO + 3 ln yH2 ln Ky2 = − ln yCO − ln yH2 O + ln yCO2 + ln yH2 The lnK values are given for 1000 K flame temperature and atmospheric pressure. The following is the molar composition of the initial mixture: 1 2 3 4 5

Component CH4 H2 CO CO2 H2

initial mole frac 0.4 0.6 0 0 0

What is the equilibrium composition? Hint: Let us start from 1 mole mixture and let x1 and x2 denote the reaction extents expressed in mole. The reaction extent tells us how many moles of the left-hand-side of a chemical reaction has been converted into the right-hand-side. The following

79

table shows the mole fractions as the function of the two reaction extents (notice that the first reaction causes mole number change, while the second does not.) 1 2 3 4 5

CH4 H2 O CO CO2 H2

0.4−x1 1+2x1 0.6−x1 −x2 1+2x1 x1 −x2 1+2x1 x2 1+2x1 3x1 +x2 1+2x1

We can reduce the problem to two equations and two unknowns using the reaction extents.

8.25

Units: Darcy’s law

If the permeability is k = 10 md ( 9.87 10−15 m2 ) and the fluid viscosity is µ = 5 cp (0.005 Pa s), what is the Darcy velocity caused by 0.5 psi/ft (11310 Pa/m) pressure gradient? Assume the cross sectional area perpendicular to the flow is 1 ft2 (0.0929 m2 ). What is the flow rate in resft3 /s and in resBBL/D ? (in m3 /s and in m3 /day?)

80

Chapter 9

Tables Quantity Length, L Time, t Mass, m Chem substance, n Temperature, T Current, I Velocity, V Area, A Volume, V Specific volume Molar volume Density, Acceleration, a Force, F Pressure, p (also stress) Work Heat Int. Energy, Enthalpy Power Specific heats (mass based) Electr. Potential (Voltage) Electr. Resistance Electr. charge Flow rate Shear rate Viscosity Permeability Compressibility

Remark basic [L] basic [t] basic [M] basic [N] basic [T] basic [I] L/ t L2 L3 V/m V/n m/V ∆V /∆t ma F/A LF Work Work Work/ t Heat/(m ∆ T) Power/Current Potential/Current Current Time ∆ Volume / ∆t ∆ Velocity/ ∆ Length Stress / shear rate Pressure/ Length Velocity Viscosity 1/Pressure

81

SI Unit m s kg mol K A m s−1 m2 m3 m3 kg−1 m3 mol−1 kg m−3 m s−2 kg m s−2 kg m−1 s−2 kg m2 s−2 kg m2 s−2 kg m2 s−2 kg m2 s−3 m2 s−2 K −1 kg m2 s−3 A−1 kg m2 s−3 A−2 As m3 s−1 s−1 Pa s m2 Pa−1

SI Other name

N Pa = N/m2 = J/m3 J = N m = Pa m3 J J W = J/s = Pa m3 /s J / (kg K) V= W/A Ω= V/A C

Multiplier 109 106 103 10−3 10−6 0.1 0.01

Prefix name giga mega kilo milli micro deci centi

Symbol G M k m µ d c

Circumference of circle Area of circle Volume of sphere Surface area of sphere Volume of right circular cylinder Pyramid

d·π d2 · π/4 d3 · π/6 d2 · π hd2 · π/4 (area of base) ·h/3

Kinetic energy Potential energy Hydrostatic pressure Darcy’s Law Pipe Mech Energ Bal Pipe Mech Energ Bal Pump power Pipe Reynolds number definition: Friction factor def (Darcy,Moody,Stanton) for laminar flow (Hagen-Poissule) for turbulent flow Acceleration of gravity Max Water Density Liquid water spec. heat (1 atm) Heat of vap. of water (1 atm) Air Density (273 K, 1 atm) Air cp /cv Air relative molecular mass

9.8066 999.973 4180 2.257 106 1.293 1.4 0.2891

∆KE = 12 m∆(V 2 ) ∆P E = mg∆h ∆p = ρg∆h ∆p Q/A = VDarcy = − µk ∆x ∆p 1 L V2 2 gρ = ∆h + 2g ∆(V ) + f d 2g

m/s2 kg/m3 J/(kg K) J/kg kg/m3 kg/mol

Temperature Abs zero 0K -273.15 o C 0R -459.67 o F

Water freezing p. (1 atm) 273.15 K 0 oC 491.67 R 32 o F

82

2

∆p = gρ∆h + ρ2 ∆(V 2 ) + f Ld ρV2 ∆pV A = ∆pQ Re = dVµ ρ f = ∆pf rict Ld ρV2 2 64 f: Re f: function(Re, relative roughness)

Water boiling p. (1 atm) 373 K 100 o C 671.67 R 212 o F

Ideal gas:

pV = nRT cp − cv = R c κ = cvp J a·m3 N ·m R= 8.314 mol·K = 8.314 Pmol·K = 8.314 mol·K f t·lbf psi·f t3 BT U R= 1.986 lb−mol·o R = 1545 lb−mol·o R = 10.73 lb−mol· oR

Other units Other unit inch ft mile acre gal bbl Liter lbm lbf psi bar Atmosphere mm-of-Hg deg-R hp BTU min hour day md cP

Expressed with 0.0254 0.3048 1609.3 4046.9 3.7854 × 10−3 0.158987 0.001 0.45359 4.4482 6894.8 1 × 105 101325 133.32 0.555556 745.70 1055.1 60 3600 86400 9.869 × 10−16 0.001

SI unit m m m m2 m3 m3 m3 kg N Pa Pa Pa Pa K W J s s s m2 Pa s

Inside the English system 1 1 1 1 1 1 1 1 1 1 1

ft mi ac sq mi BBL lbf psi atm BTU hp sp. grav. (liq) sp. grav. (liq) lb-mole M

= = = = = = = = = = = : : :

12 5280 43,560 640 42 32.174 1 14.696 778.17 42.41 8.34 141.5 131.5+AP I o

453.59 multiply by 1000

83

in. ft ft2 ac gal lbm f t s2 lbf in.2 lbf in.2 lbf ft BT U min lbm gal

mol

=

5.6146

ft3

=

144

lbf f t2

=

25,037

=

62.4

:

379

lbm f t2 s2 lbm f t3

SCF