Numerical Analysis

PREFACE

About the Author

Harold A. Toomey, M.S. in Electrical and Computer Engineering, is currently a Software Engineer for Novell in Provo, Utah. While minoring in mathematics at Brigham Young University, he tutored students in calculus, then tutored C programming at BYU's Electrical Engineering Department. Not content with the provided FORTRAN algorithms while taking several numerical methods courses, he began coding numerical algorithms in C. The introductory text used for these numerical analysis courses was "Numerical Analysis."

History of "Numerical Analysis Algorithms in C"

BYU's mathematics department expressed an interest in having all of the algorithms found in the "Numerical Analysis" text programmed in C, along with a few of their favorites still in FORTRAN. Version 3.0 was finally completed in December 1988. BYU was the first university to purchase a university site license. This software is being used for their numerical methods courses today and has been tested by hundreds of students. Their input has resulted in several other versions, culminating into version 4.2. Version 4.0 became necessary for the fourth edition of the "Numerical Analysis" text. Since 1988, several universities and scores of students have purchased these programs to be used in college course work and on the job. See the file "revhist.doc" (for revision history) for a complete overview of the history of "Numerical Analysis Algorithms in C."

Acknowledgements

The author would like to express his appreciation to the many individuals who made suggestions for improvement on the previous versions of these algorithms. These include the professor who gave directions for the first version: G. S. Gill, Brigham Young University (also a reviewer for the third edition of the text), and Bruce Cardwell who supervises the Numerical Analysis Laboratory also at Brigham Young University. Special thanks also go to Jay Lawlor, M.S. Electrical Engineering, for giving timely feedback while using the algorithms for a numerical methods class at BYU. In particular, thanks also goes to Holly Z. Toomey for typesetting previous versions of the Examples Book.

CONTENTS

PREFACE ‑iv‑

1. Introduction 1‑1

1.1 Getting Started 1‑1

1.2 Purpose of the Programs 1‑2

1.3 For Instructors 1‑2

1.3.1 "Numerical Analysis" Authors' Recommendations 1‑2

1.3.2 Homework Helpers 1‑2

1.3.3 Modifying Programs 1‑3

1.3.4 Intentionally Introducing Errors 1‑3

1.4 Product Support 1‑3

2. Installation 2‑1

2.1 Basic Installation Procedures 2‑1

2.2 Uploading to Mainframe Computers 2‑2

3. "Numerical Analysis Algorithms in C" Files 3‑1

3.1 Algorithm Files 3‑1

3.3 Supporting C Source Code 3‑5

3.4 Documentation Files 3‑6

3.5 Utility Files 3‑6

3.6 Batch, Script and Command Files 3‑7

3.7 File Structure Chart 3‑8

3.8 File Name Translation Table from 3rd to 4th Edition 3‑8

3.9 4th Edition Differences 3‑9

4. Step-By-Step Examples on Various Computers 4‑1

4.1 Need List 4‑1

4.2 Customizing Naautil.c 4‑1

4.3 Example Using MS-DOS, Microsoft C and the P-Edit Editor 4‑2

4.4 Example Using UNIX, cc and the vi Editor 4‑5

4.5 Example Using a Macintosh and THINK C 4‑8

4.6 Example Using VAX/VMS, CC and the EDIT/EDT Editor 4‑12

5. For Those New to C 5‑1

5.1 Mathematical Operators 5‑1

5.2 Mathematical Functions 5‑2

5.3 General Language Hints 5‑5

5.4 Language Transition Kit 5‑6

6. Helps and Hints 6‑1

6.1 Generally Nice To Know 6‑1

6.1.1 Professor's Favorites, Must Have, Algorithms 6‑1

6.1.2 Homework Helper Algorithms 6‑1

6.1.3 Optional Title 6‑1

6.1.4 Optional File Saving 6‑2

6.1.5 Finding Functions 6‑2

6.1.6 Using Default Inputs 6‑2

6.1.7 Changing Arithmetic Precision 6‑2

6.1.8 Using Floating-Point Numbers in Functions 6‑3

6.1.9 The Pow() Function 6‑4

6.1.10 Implementing SIG-Digit Rounding/Truncation 6‑4

6.1.11 Floating-Point Output Alignment 6‑5

6.2 Converting Programs into Functions 6‑5

6.2.1 An Example Using Simpson's Rule 6‑7

6.3 Using Input Files (*.IN) 6‑8

6.4 Using Output Files (*.OUT) 6‑10

6.5 Explanation of the Naautil.c File 6‑10

6.5.1 #Define Flags 6‑10

6.5.2 Flag Default Settings 6‑11

6.5.3 Description of the Routines 6‑12

6.6 Using Naautil.c as Object Code 6‑14

6.6.1 MS-DOS 6‑15

6.6.2 UNIX 6‑15

6.6.3 Macintosh 6‑15

6.6.4 VAX/VMS 6‑16

6.7 Supporting C Source Code Usage List 6‑16

6.8 "Numerical Analysis" Text Errors and Corrections 6‑17

6.8.1 3rd Edition Errors 6‑17

6.8.2 4th Edition Errors 6‑18

6.9 Watch for These Run-Time Errors 6‑20

6.9.1 Stack Space 6‑20

6.9.2 Division By Zero 6‑20

6.9.3 Null Pointer Assignments 6‑20

6.9.4 No Disk Space 6‑21

6.9.5 Floating-Point Accuracy 6‑21

6.9.6 Program Stuck in an Infinite Loop 6‑21

7. Useful Utilities 7‑1

7.1 Convert.c - Converting Files from Extended ASCII to Standard ASCII 7‑1

7.1.1 Why Convert.c is Needed 7‑1

7.1.2 How to Use Convert.c 7‑2

7.2 List.com - A Better TYPE Command 7‑3

7.3 Time-Saving Batch, Script and Command Files 7‑3

7.3.1 CC.BAT 7‑3

7.3.2 CCC 7‑5

7.3.3 VAXCC.COM 7‑6

8. The Equation Evaluator Routines 8‑1

8.1 What the Routines Do 8‑1

8.2 How to Insert the Routines into a Program 8‑1

8.3 An Example Using Simpson's Rule 8‑2

8.4 Using Eqeval.c As Pre-Compiled Object Code 8‑2

8.5 Valid Math Operators and Functions 8‑3

8.6 Sample Equations 8‑4

8.7 Possible Error Messages 8‑4

8.8 List of Algorithms Using the Equation Evaluator Routines 8‑5

8.9 Limitations 8‑6

8.10 Trade-Offs 8‑6

9. Portability 9‑1

9.1 C vs ANSI C 9‑2

9.2 IBM PCs and MS-DOS 9‑3

9.3 UNIX Workstations 9‑3

9.4 Macintosh Computers 9‑4

9.5 VAX Mainframes 9‑5

9.6 Tested Compilers 9‑5

10. Sample License Agreements 10‑1

10.1 Individual License Sample 10‑1

10.2 University/Corporation Site License Sample 10‑3

11. Packaging Information 11‑1

11.1 MS-DOS Diskettes 11‑1

11.1.1 5¼" 1.2M High Density Diskettes 11‑1

11.1.2 5¼" 360K Low Density Diskettes 11‑2

11.1.3 3½" 1.44M High Density Diskettes 11‑2

11.1.4 3½" 720K Low Density Diskettes 11‑2

11.2 Macintosh Diskettes 11‑2

11.2.1 3½" 800K Macintosh Diskettes 11‑3

12. Purchasing Information 12‑1

12.1 $20.00 Club 12‑1

12.2 Order Form 12‑1

References 12‑2

Appendix A: C Source Code for 041.C A‑1

Appendix B: C Source Code for NAAUTIL.C B‑1

Appendix C: Language Comparison Charts C‑1

C.1 C vs Ada C‑2

C.2 C vs BASIC C‑8

C.3 C vs C++ C‑13

C.4 C vs FORTRAN 77 C‑14

C.5 C vs Pascal C‑20

Appendix D: Sample Programs in Other Languages D‑1

D.1 Ada D‑2

D.1.1 SIMPSON.ADA D‑2

D.1.2 NAAUTIL.ADA D‑4

D.1.3 SIMPSON.IN D‑6

D.1.4 SIMPSON.OUT D‑6

D.2 BASIC D‑7

D.2.1 SIMPSON.BAS D‑7

D.2.2 SIMPSON.IN D‑8

D.2.3 SIMPSON.OUT D‑9

D.3 C D‑10

D.3.1 SIMPSON.C D‑10

D.3.2 NAAUTIL.H D‑11

D.3.3 SIMPSON.IN D‑14

D.3.4 SIMPSON.OUT D‑14

D.4 C++ D‑15

D.4.1 SIMPSON.CPP D‑15

D.4.2 NAAUTIL.HPP D‑16

D.4.3 SIMPSON.IN D‑18

D.4.4 SIMPSON.OUT D‑18

D.5 FORTRAN 77 D‑19

D.5.1 SIMPSON.FOR D‑19

D.5.2 SIMPSON.IN D‑21

D.5.3 SIMPSON.OUT D‑21

D.6 Pascal D‑22

D.6.1 SIMPSON.PAS D‑22

D.6.2 NAAUTIL.INC D‑24

D.6.3 NAAMATH.INC D‑25

D.6.4 SIMPSON.IN D‑25

D.6.5 SIMPSON.OUT D‑25

1. Introduction

"Numerical Analysis Algorithms in C" contains 116 stand-alone programs implementing the algorithms found in the texts:

"Numerical Analysis", third and fourth edition,

Richard L. Burden & J. Douglas Faires, 1988.

Each program is written in ANSI C to make them more portable to other computer systems. They should run on any computer with a reasonable C compiler, such as IBM PCs, UNIX workstations, VAXes, and Macintoshes.

The "Numerical Analysis" text, hereafter referred to as "the text", covers the following numerical topics:

Chapter 1 - Mathematical Preliminaries

Chapter 2 - Solutions of equations in one variable

Chapter 3 - Interpolation and polynomial approximation

Chapter 4 - Numerical differentiation and integration

Chapter 5 - Initial-value problems for ordinary differential equations

Chapter 6 - Direct methods for solving linear systems

Chapter 7 - Iterative techniques in matrix algebra

Chapter 8 - Approximation theory

Chapter 9 - Approximating eigenvalues

Chapter 10 - Numerical solutions of nonlinear systems of equations

Chapter 11 - Boundary-value problems for ordinary differential equations

Chapter 12 - Numerical solutions to partial differential equations

From these topics, "Numerical Analysis Algorithms in C" has programmed routines for: vector and matrix manipulation, linear equations (LU decomposition/backsolving, matrix inversion, etc.), matrix/vector norms, eigenvalue/vectors, complex number and polynomial manipulation, least-square polynomial approximation, FFTs, numerical integration, root finding, solution of nonlinear equations, Taylor polynomial approximation, cubic splines, derivatives, ordinary and partial differentiation.

This User's Manual will help you to use these programs to their fullest potential. It will walk you through an example, tutor you if you are unfamiliar with the C language, introduce you to several useful utilities, and assist you when running these programs on different computer systems.

1.1 Getting Started

To install "Numerical Analysis Algorithms in C" onto your computer system, see Chapter 2 - "Installation." If you are new to the C programming language, you may wish to read through Chapter 5 - "For Those New to C." If you want a detailed example using various C compilers and operating systems, see Chapter 4 - "Step-By-Step Examples on Various Computer Systems."

This software package contains about 1.5M bytes of files. If disk space is limited, then just copy the eight supporting ".c" files ("complex.c", "eqeval.c", "gaussj.c", "naautil.c", "naautil2.c", "naautil3.c", "round.c" and "trunc.c") and the desired algorithms onto your disk. The eight supporting files require about 100K of disk space. If you are running these algorithms from a floppy disk, be sure to leave the write protect tab off so the programs can save their output to a file. If this is undesirable, see Sub-Section 6.1.4 - "Optional File Saving."

If you feel comfortable with C, go ahead and compile and run an algorithm. The source code is very readable and user friendly. To see what the algorithm numbers correspond to, see Section 3.1 - "Algorithm Files." This is the most important list in this manual and should be printed out for frequent reference. Section 3.1 is also given in the file "readme.doc" for your convenience.

1.2 Purpose of the Programs

These programs are fast, but are not optimized for speed. As stated by the authors in the text's preface:

"Although the algorithms will lead to correct programs for the examples and exercises in the text, it must be emphasized that there has been no attempt to write general-purpose software. In particular, the algorithms have not always been listed in the form that leads to the most efficient program in terms of either time or storage requirements."

The purpose of these programs is to teach students numerical methods, not programming and optimization skills. For a good book of general-purpose mathematical software, see the book "Numerical Recipes in C" listed in the references. These programs can also be used as a tool for building other programs. Once the algorithms are understood, they can be more easily enhanced for general-purpose applications.

1.3 For Instructors

This software package is intended to be used by instructors of numerical methods/analysis courses. The best way to learn numerical methods is to program the algorithms from scratch and have them run on a computer. This is a time consuming process and may take a "good" programmer from 1 to 5 hours per program. Students can best benefit from these programs AFTER taking the appropriate numerical analysis courses.

1.3.1 "Numerical Analysis" Authors' Recommendations

The authors of the text "Numerical Analysis" mention in the preface that:

"Actual programs are not included because, in our experience, this encourages some students to generate results without fully understanding the method involved."

In other words, as an instructor, you may consider giving your students only selected main algorithms, and definitely not the "Homework Helpers" algorithms as discussed below.

1.3.2 Homework Helpers

Roughly half of the included programs are labeled as "Homework Helpers." Most of these programs modify the given text algorithms to satisfy the homework exercises in the text. An example of this is turning Algorithm 2.4 - Secant Method ("024.c") into the Method of False Position ("024B.c"). Use these "homework helpers" to correct homework assignments. Do NOT just give these out to your students. Most modifications will take only a short time to implement, once the algorithm is understood.

1.3.3 Modifying Programs

These algorithms are given as a learning tool. Modifying them is part of the learning process. These algorithms may be modified by the instructor or by the students, even though this package is copyrighted. They may not, however, be altered to be resold for profit without prior written consent from the programmer. See the sample licensing agreements in Chapter 10 for more details.

1.3.4 Intentionally Introducing Errors

As an alternative to withholding these programs from your students, you may wish to give them a copy with intentionally introduced errors. This would cause them to search the entire program over for correctness, bridging the gap between giving too little or too much information.

1.4 Product Support

If questions arise, ranging from getting these algorithms to work with your compiler to adapting a particular algorithm to a specific application, just call CARE-FREE SOFTWARE at 1-801-785-0464. The programmer will answer your questions at no charge other than the normal phone charges on your monthly phone statement. Enhancements, recommendations and bug reports are always welcomed.

2. Installation

2.1 Basic Installation Procedures

The "Numerical Analysis Algorithms in C" programs do not come with an installation program. To install these algorithms onto your computer, do the following steps:

1. Make a set of backup diskettes. See your operating system manual for specifics.

2. Make another set of "working" diskettes or copy the diskettes onto your hard disk. All 500+ files combined require less than 1.5M bytes of disk space. Only a couple of the files are required at a time to get the algorithms to work properly, making them useful even on systems without a hard disk.

3. You may want to convert each file on the "working" disk from extended ASCII to standard ASCII. This is usually required for Macintoshes, most UNIX computers, and VAXes. Failure to do so may result in scrambled looking output characters. Use "convert.exe", as explained in Section 7.1, to do this task relatively easily. Macintosh disks ordered from Care-Free Software have had this step done already.

4. It is recommended that the algorithms be placed in their own sub-directory (or Macintosh folder), such as "naa42." This sub-directory can be created and entered by typing one of the following sets of commands:

MS-DOS:

C:\> MD NAA42 - make directory

C:\> CD NAA42 - change directory

C:\NAA42> DIR /P - show directory contents

UNIX:

% mkdir naa42 - make directory

% chdir naa42 - change directory

% pwd - show current directory

% ls -alF - show directory contents

VAX/VMS:

$ CREATE/DIR [SMITH.NAA42] - make directory

$ SET DEFAULT [.NAA42] - change directory

$ SHOW DEFAULT - show current directory

$ DIR/SIZE/DATE - show directory contents

5. To be able to run every program from a floppy diskette, eight support files are required:

complex.c naautil.c round.c

eqeval.c naautil2.c trunc.c

gaussj.c naautil3.c

These files require about 100K bytes of disk space. The desired algorithm files such as "041.c" are also needed. The majority of the algorithms need only "naautil.c" which is about 20K bytes large.

6. If the programs do not compile correctly, you may need to change some flags inside the "naautil.c" file. Use your text editor to modify this file. The contents of "naautil.c" should be self-documenting. These flags are defined near the top of the file. See Section 6.5 - "Explanation of the Naautil.c File" if more detailed information is desired.

In the event that nothing seems to be working, you can set both the EQ_EVAL and the FILE_SAVE flags to FALSE. This will disable the options to save the output to a file and to use the Equation Evaluator routines, but the algorithms will usually work. These two options use variable length argument lists, which may not work on older compilers.

7. If all else fails, ask another C programmer for help or call CARE-FREE SOFTWARE for free technical support.

2.2 Uploading to Mainframe Computers

To get these programs onto many workstations or mainframe computers, communications software is usually required. A well-supported communications protocol is known as Kermit. An example using Kermit looks something like this:

NOTE: This example uses CALL/ProComm to transfer files onto a VAX/UNIX workstation.

1. Log onto the mainframe using CALL, ProComm or your favorite communications package. Select kermit as the transfer protocol. Use binary mode to send files containing extended ASCII characters. Use ASCII mode if the files have been converted to standard ASCII by the "convert.exe" program. Binary mode is slower than ASCII mode. Remember that C files are case sensitive.

2. On the mainframe, change to an appropriate directory and type:

For a VAX, type:

$ use kermit (Do NOT type "$ kermit")

Kermit-32> set file_type binary (or: set file_type ascii)

Kermit-32> receive

For a UNIX workstation, type:

% kermit

Kermit-32> set binary (or: set ascii)

Kermit-32> receive

3. On your PC, immediately issue the file sending commands.

For CALL, type:

[F9] File Send Kermit

File to transfer: filename

For ProComm, type:

[ALT] K

2) Send

Please enter filespec: filename

4. Patiently wait as the file(s) are transferred to the mainframe. The use of wild cards is recommended (ie - *.C instead of filename).

5. Exit kermit on the mainframe.

Kermit-32> exit

$ logout

A host full of other issues have been left to the user, such as baud rate, parity, stop bits, duplex, use of wild cards, etc. These are unique to each computer system and communications software package.

You may want to convert the files from extended ASCII to standard ASCII (using "convert.c") before uploading them to a mainframe computer. If you plan to view and print your work on an IBM PC but compile and run the algorithms on a mainframe, you may want to keep the files in extended ASCII.

Test your preferences using Algorithm 4.1 ("041.c"). It uses three different extended ASCII characters to form an integral sign: '!', '#' and '"'. "Convert.c" changes these three characters into standard ASCII: '[', '|' and ']'.

3. "Numerical Analysis Algorithms in C" Files

This software package contains 116 algorithms. Each algorithm has been coded as a stand-alone program. Each program prompts for input, executes the algorithm as described in the text "Numerical Analysis", and prints the results. Other math packages provide only subroutines, requiring a programmer to insert them inside a program and either hard code or prompt for the inputs and print the outputs.

The files are catagorized as follows, where "nnn" represent algorithm numbers like "041" for Algorithm 4.1:

a. nnn.C Algorithms from the text "Numerical Analysis" fourth edition. (57 total)

b. nnnA.C Algorithms not found in the text. Included as "Professor Favorites, Must Have" as recommended by mathematics professors at Brigham Young University. (6 total)

c. nnnB.C, nnnC.C, and nnnD.C

Algorithms included as "Homework Helpers." Some are asked for in the homework exercises while others are for helping with important concepts covered in the text. These can save hours of coding on the homework exercises. (53 total)

d. *.C NAA supporting files containing 57 functions. (8 total)

e. *.IN Input files used to test each algorithm. They match the inputs to the example problems presented after each algorithm in the text. (116 total)

f. *.OUT Output files used to test each algorithm. They match the outputs to the example problems presented after each algorithm in the text. (116 total)

g. *.EXE Executable programs for each algorithm. The default functions (like f(x)) are the same as those used in the example problems presented after each algorithm in the text. These programs must be purchased separately and are currently available only for MS-DOS and Macintosh computers. (116 total)

h. *.DOC Documentation in simple text file format. Includes "readme.doc", "revhist.doc" and "usersman.doc."

Each program was tested on the sample problems given in the text just after the algorithm description. These sample solutions are found in the OUT sub-directory in files named with a ".out" extension. Their inputs are found in the IN sub-directory in files named with a ".in" extension.

Over two-thirds of the algorithms need to be compiled only once. They are marked with an asterisk (*) on the table below. Of these algorithms, nearly half are able to prompt you for an equation during run-time. See Chapter 8 - "The Equation Evaluator Routines" for more details.

3.1 Algorithm Files

CHAPTER 1 Mathematical Preliminaries

COMPLEX.C - "Numerical Recipes in C" Complex Number Routines

EQEVAL.C - Equation Evaluator Routines

GAUSSJ.C - "Numerical Recipes in C" Gauss-Jordan Matrix Solver

NAAUTIL.C - "Numerical Analysis Algorithms in C" Utilities I (standard)

NAAUTIL2.C - "Numerical Analysis Algorithms in C" Utilities II (extended)

NAAUTIL3.C - "Numerical Analysis Algorithms in C" Utilities III (complex)

ROUND.C - Rounds a floating-point value to SIG significant digits

TRUNC.C - Truncates a floating-point value to SIG significant digits

011B.C* - Taylor Polynomial Approximation Algorithm 1.1B

CHAPTER 2 Solutions of Equations in One Variable

021.C* - Bisection (or Binary-Search) Algorithm 2.1

022.C* - Fixed-Point Algorithm 2.2

023.C - Newton-Raphson Algorithm 2.3

024.C* - Secant Algorithm 2.4

024B.C* - Method of False Position (or Regula Falsi) Algorithm 2.4B

024C.C - Modified Newton-Raphson Method Algorithm 2.4C

025.C* - Steffensen Algorithm 2.5

026.C* - Horner Algorithm 2.6

027.C* - Müller Algorithm 2.7

028A.C* + Complex Polynomial Solver (CPOLY) Algorithm 2.8A

CHAPTER 3 Interpolation and Polynomial Approximation

031.C* - Neville's Iterated Interpolation Algorithm 3.1

031B.C* - Neville's Iterated Interpolation (with rounding) Algorithm 3.1B

031C.C* - Aitken's Iterated Interpolation Algorithm 3.1C

032.C* - Newton's Interpolatory Divided-Difference Formula Algorithm 3.2

033.C* - Hermite Interpolation Algorithm 3.3

034.C* - Natural Cubic Spline Algorithm 3.4

035.C* - Clamped Cubic Spline Algorithm 3.5

CHAPTER 4 Numerical Differentiation and Integration

040B1.C - 1st Derivative Approximation (for functions) Algorithm 4.0B1

040B2.C* - 1st Derivative Approximation (for tabulated data) Algorithm 4.0B2

040B3.C - 1st Derivative Approximation (for functions w/TOL) Algorithm 4.0B3

040C1.C - 2nd Derivative Approximation (for functions) Algorithm 4.0C1

040C2.C* - 2nd Derivative Approximation (for tabulated data) Algorithm 4.0C2

040D1.C* - Richardson's Extrapolation Algorithm 4.0D1

040D2.C* - Richardson's Extrapolation (with rounding) Algorithm 4.0D2

041.C* - Composite Simpson's Rule Algorithm 4.1

041B.C* - Composite Trapezoidal Rule Algorithm 4.1B

041C.C* - Composite Midpoint Rule Algorithm 4.1C

041D.C* - Newton-Cotes Formulas for Integrals (8 total) Algorithm 4.1D

042.C* - Adaptive Quadrature Algorithm 4.2

043.C* - Romberg Algorithm 4.3

043B.C* - Gaussian Quadrature Algorithm 4.3B

044.C - Composite Simpson's Rule for Double Integrals Algorithm 4.4

044B.C - Composite Trapezoid Rule for Double Integrals Algorithm 4.4B

044C.C - Gaussian Quadrature for Double Integrals Algorithm 4.4C

045.C - Composite Simpson's Rule for Triple Integrals Algorithm 4.5

045B.C - Composite Trapezoid Rule for Triple Integrals Algorithm 4.5B

045C.C - Gaussian Quadrature for Triple Integrals Algorithm 4.5C

CHAPTER 5 Initial-Value Problems for Ordinary Differential Equations

051.C* - Euler Algorithm 5.1

051B.C* - Midpoint, Modified Euler, and Heun's Methods Algorithm 5.1B

052.C* - Runge-Kutta (Order Four) Algorithm 5.2

053.C - Runge-Kutta-Fehlberg Algorithm 5.3

054.C* - Adam's Fourth-Order Predictor-Corrector Algorithm 5.4

054B.C* - Adams-Bashforth (all four) and Milne's Methods Algorithm 5.4B

054C.C* - Milne-Simpson Predictor-Corrector Algorithm 5.4C

055.C* - Adam's Variable Step-Size Predictor-Corrector Algorithm 5.5

056.C* + Extrapolation Algorithm 5.6

057.C - Runge-Kutta for Systems of Differential Equations Algorithm 5.7

057B.C - Euler's Variable Step-Size for Systems Algorithm 5.7B

058.C - Trapezoidal with Newton Iteration Algorithm 5.8

CHAPTER 6 Direct Methods for Solving Linear Systems

060B.C* - Matrix Inverter Algorithm 6.0B

060C.C* - Determinant of a Matrix Algorithm 6.0C

060D.C* - Matrix Multiplier Algorithm 6.0D

061.C* - Gaussian Elimination with Backward Substitution Algorithm 6.1

061B.C* - Gaussian Elimination with Backward Substitution Algorithm 6.1B

(with rounding)

061C1.C* - Gauss-Jordan Method Algorithm 6.1C1

061C2.C* - Gauss-Jordan Method (with rounding) Algorithm 6.1C2

061D1.C* - Gaussian-Elimination - Gauss-Jordan Hybrid Method Algorithm 6.1D1

061D2.C* - Gaussian-Elimination - Gauss-Jordan Hybrid Method Algorithm 6.1D2

(with rounding)

062.C* - Gaussian Elimination with Maximal Column Pivoting Algorithm 6.2

062B.C* - Gaussian Elimination with Maximal Column Pivoting Algorithm 6.2B

(with rounding)

063.C* - Gaussian Elimination with Scaled Column Pivoting Algorithm 6.3

063B.C* - Gaussian Elimination with Scaled Column Pivoting Algorithm 6.3B

(with rounding)

064.C* - Direct Factorization Algorithm 6.4

064B.C* - Direct Factorization which solves AX=B Algorithm 6.4B

064C.C* - Direct Factorization with Maximal Column Pivoting Algorithm 6.4C

(3rd edition)

065.C* - LDL^t Factorization Algorithm 6.5

065B.C* - LDL^t Factorization which solves AX=B Algorithm 6.5B

066.C* - Choleski Algorithm 6.6

066B.C* - Choleski which solves AX=B Algorithm 6.6B

067.C* - Crout Reduction for Tridiagonal Linear Systems Algorithm 6.7

CHAPTER 7 Iterative Techniques in Matrix Algebra

070B.C* - Vector and Matrix Norms Algorithm 7.0B

071.C* - Jacobi Iterative Algorithm 7.1

072.C* - Gauss-Seidel Iterative Algorithm 7.2

073.C* - Successive Over Relaxation (SOR) Algorithm 7.3

074.C* - Iterative Refinement (with rounding) Algorithm 7.4

074B.C* - Iterative Refinement (single-precision) Algorithm 7.4B

CHAPTER 8 Approximation Theory

080B.C* - Least-Squares Polynomial Approximation Algorithm 8.0B

081.C* + Fast Fourier Transformation Algorithm 8.1

CHAPTER 9 Approximating Eigenvalues

091.C* - Power Method Algorithm 9.1

091B.C* - Power Method with Aitken's Delta² Method Algorithm 9.1B

092.C* - Symmetric Power Method Algorithm 9.2

093.C* - Inverse Power Method Algorithm 9.3

094.C* - Wielandt's Deflation Algorithm 9.4

094B.C* - Wielandt's Deflation using Power Method for lambda1 Algorithm 9.4B

O095.C* - Householder Method Algorithm 9.5

095B.C* - Householder Method (3rd edition) Algorithm 9.5B

095C.C* - Householder Method for Non-Symmetric Matrices Algorithm 9.5C

(Upper Hessenberg)

095D.C* - Householder Method (with rounding) Algorithm 9.5D

096.C* - QR Algorithm Algorithm 9.6

096B.C* - QL Algorithm (3rd edition) Algorithm 9.6B

CHAPTER 10 Numerical Solutions of Nonlinear Systems of Equations

101.C - Newton's Method for Systems Algorithm 10.1

101A.C - Steffensen's Method for Systems Algorithm 10.1A

102.C - Broyden's Method for Systems Algorithm 10.2

103.C - Steepest Descent Method (with F(x) and J(x)) Algorithm 10.3

103B.C - Steepest Descent Method (with G(x) and gradG(x)) Algorithm 10.3B

CHAPTER 11 Boundary-Value Problems for Ordinary Differential Equations

111.C - Linear Shooting Algorithm 11.1

112.C - Nonlinear Shooting with Newton's Method Algorithm 11.2

112B.C - Nonlinear Shooting with Secant Method Algorithm 11.2B

113.C - Linear Finite Difference Algorithm 11.3

113B.C - Linear Finite Difference (Richardson's Extrapolation) Algorithm 11.3B

114.C - Nonlinear Finite Difference Algorithm 11.4

114B.C - Nonlinear Finite Difference (Richardson's Extrapolation) Algorithm 11.4B

115.C - Piecewise Linear Rayleigh-Ritz Algorithm 11.5

116.C - Cubic Spline Rayleigh-Ritz Algorithm 11.6

CHAPTER 12 Numerical Solutions to Partial-Differential Equations

121.C - Poisson Equation Finite-Difference (Elliptic) Algorithm 12.1

122.C* - Heat Equation Backward-Difference (Parabolic) Algorithm 12.2

122B.C* - Heat Equation Forward-Difference (Parabolic) Algorithm 12.2B

122C.C* - Heat Equation Richardson's Method (Parabolic) Algorithm 12.2C

123.C* - Crank-Nicolson (Parabolic) Algorithm 12.3

124.C - Wave Equation Finite-Difference (Hyperbolic) Algorithm 12.4

125.C - Finite-Element Algorithm 12.5

126A.C - Parabolic Equations With Newton Iteration in 1-D Algorithm 12.6A

127A.C - Parabolic Equations With Newton Iteration in 2-D Algorithm 12.7A

128A.C - Elliptic Equations With Newton Iteration in 2-D Algorithm 12.8A

129A.C - Biharmonic Equation Using Gauss-Jordan Method Algorithm 12.9A

The '+'s above mean the program may need a larger stack when compiled and linked.

The '*'s above mean the program needs to be compiled only once.

3.3 Supporting C Source Code

The eight files below are needed to compile each and every program. Most algorithms require only one or two of them at a time.

COMPLEX.C

"Complex.c" contain several routines for operating on complex numbers. It originated from the book "Numerical Recipes in C" and is only used in "naautil3.c."

EQEVAL.C

"Eqeval.c" contains the Equation Evaluator routines. These routines enable a program to enter and evaluate an equation during run-time. It is useful within algorithms that need to evaluate a single function such as f(x) or f(y,t). It is used by 34 algorithms. See Chapter 8 - "The Equation Evaluator Routines" for more details on this file.

GAUSSJ.C

"Gaussj.c" is a Gauss-Jordan matrix solver routine. It originated from the book "Numerical Recipes in C." It is used by only 9 of the algorithms.

NAAUTIL.C

"Naautil.c" contain important routines used by all of the algorithms. Most are for dynamically allocating memory for arrays. Some of the routines originated from the book "Numerical Recipes in C." See Section 6.5 - "Explanation of the Naautil.c File."

NAAUTIL2.C

"Naautil2.c" contains more dynamically allocated memory routines for less-used data types. it is used only 2 times.

NAAUTIL3.C

"Naautil3.c" contains more dynamically allocated memory routines for complex data types. It is used only 3 times.

ROUND.C

"Round.c" rounds a floating-point value to SIG significant digits. Only 9 algorithms currently use this function. See Sub-Section 6.1.10 to see how this file is used.

TRUNC.C

"Trunc.c" truncates, or chops, a floating-point value to SIG significant digits. None of the algorithms use this function, but it can easily replace "round.c."

3.4 Documentation Files

Previous versions of "Numerical Analysis Algorithms in C" consisted of only two document files: "readme.doc" and "math.doc." With version 4.2, these documents have been consolidated and greatly expanded into this User's Manual ("usersman.doc"). Three document files are included as listed below.

README.DOC

"Readme.doc" gives a list of all the algorithms as well as an order form. This information can also be found inside the User's Manual.

REVHIST.DOC

"Revhist.doc" gives a detailed list of all changes made to each version of "Numerical Analysis Algorithms in C". It lists the additions, corrections, and changes made to each algorithm, to the supporting files, and to the documentation.

USERSMAN.DOC

"Usersman.doc" is this User's Manual in DOS text format. This format is readable by all text editors and word processors. It can be read using MS-DOS's "type" command or the "list.com" utility included with the diskettes.

3.5 Utility Files

041EE.C

"041ee.c" is an example of how to integrate the equation evaluator routines into an algorithm.

041FUN.C

"041fun.c" is an example of Algorithm 4.1 turned into a stand-alone function.

CONVERT.C

"Convert.c" is the C source code for a utility which translates text files into standard seven-bit ASCII files. It is useful before placing these algorithms on non-MS-DOS computers, such as UNIX and VAX computers. See Section 7.1 - "Convert.c - Converting Files from Extended ASCII to Standard ASCII."

CONVERT.EXE

"Convert.exe" is the MS-DOS executable of "convert.c."

LISTALL

"Listall" is a text file listing all source code files on the root directory of the distribution disks. It can be used with "convert.exe" to convert all the programs at once.

LISTOUT

"Listout" is a text file listing all output files in the OUT sub-directory of the distribution disks. It can be used with "convert.exe" to convert all of the output files at once.

LIST.COM

"List.com" is an MS-DOS program which acts as a better "TYPE" command. It uses the arrow keys and other editing keys to view text files. "List.com" does not allow you to edit files, just view them. It is a public domain program. See Section 7.2 - "List.com - A better TYPE Command" for instructions on how to use it.

3.6 Batch, Script and Command Files

Three commands text files are included to simplify the task of compiling and running the algorithms on different computer systems.

CC.BAT

"Cc.bat" is an MS-DOS batch file used for compiling, running and viewing a Microsoft C 5.0 program. It can be easily altered to allow for linking to "naautil.c" and "eqeval.c" object files, speeding up the compile time. It can also be altered to increase the stack size of a program.

CCC

"Ccc" is a UNIX script file used for compiling, running, and viewing a C program. It can be easily altered to allow for linking to "naautil.c" object code, speeding up the compile time.

VAXCC.COM

"Vaxcc.com" is a VAX/VMS command file used for compiling and linking a mathematical VAX C program. It can be easily altered to allow for linking to "naautil.c" object code, speeding up the compile time.

3.7 File Structure Chart

The chart below describes how the files are organized on the distribution diskettes.

/ (root)

+))))))))0))))))))0))))))))0))2)))))0))))))))0)))))))),

* * * * * * * *.C *.DOC UTIL LANGS IN OUT EXE * * * * *

* * * * *

*.* * *.IN *.OUT *.EXE

* (OPTIONAL)

+)))))))))))0)))))))))))0)))))2)))))0)))))))))))0))))))))))),

* * * * * *

ADA BASIC C CPP FORTRAN PASCAL

* * * * * *

SIMPSON.ADA SIMPSON.BAS SIMPSON.C SIMPSON.CPP SIMPSON.FOR SIMPSON.PAS

NAAUTIL.ADA SIMPSON.IN SIMPSON.H SIMPSON.HPP SIMPSON.IN NAAUTIL.INC

SIMPSON.IN SIMPSON.OUT SIMPSON.IN SIMPSON.IN SIMPSON.OUT NAAMATH.INC

SIMPSON.OUT SIMPSON.OUT SIMPSON.OUT SIMPSON.IN

SIMPSON.OUT

3.8 File Name Translation Table from 3rd to 4th Edition

This translation table correlates the third edition text algorithms with the fourth edition text algorithms. The B and C extensions indicate algorithms that were changed or replaced from the third edition and retained with the fourth edition algorithms.

Edition * Edition Edition * Edition Edition * Edition

3rd * 4th 3rd * 4th 3rd * 4th

))))))))3))))))))) ))))))))3))))))))) ))))))))3)))))))))

2.1 * 2.1 5.3 * 5.3 8.6 * 9.2

2.2 * 2.2 5.4 * 5.4 8.7 * 9.3

2.3 * 2.3 5.5 * 5.5 8.8 * 9.4

2.4 * 2.4 5.6 * 5.6 8.9 * 9.5

2.5 * 2.5 5.7 * 5.7 8.10 * 9.6B

2.6 * 2.6 5.8 * 5.8 9.1 * 10.1

2.7 * 2.7 6.1 * 6.1 9.2 * 10.2

3.1 * 3.1 6.2 * 6.2 9.3 * 10.3

3.2 * 3.2 6.3 * 6.3 10.1 * 11.1

3.3 * 3.3 6.4 * 6.4 10.2 * 11.2

3.4 * 3.4 6.5 * 6.4C 10.3 * 11.3

3.5 * 3.5 6.6 * 6.6 10.4 * 11.4

4.1 * 4.1 6.7 * 6.7 10.5 * 11.5

4.2 * 4.2 8.1 * 7.1 10.6 * 11.6

4.3 * 4.3 8.1 * 7.1 11.1 * 12.1

4.4 * 4.4 8.2 * 7.2 11.2 * 12.2

5.1 * 5.1 8.3 * 7.3 11.3 * 12.3

5.2 * 5.2 8.4 * 7.4 11.4 * 12.4

8.5 * 9.1 11.5 * 12.5

3.9 4th Edition Differences

In the fourth edition's PREFACE, pages vii-viii list the "CHANGES IN THE FOURTH EDITION". The specifics of these changes are listed below.

Renamed Algorithms: 4.1, 4.4, 7.1, 7.2, 9.2, 10.1, 11.2

New to 4th Edition: 4.5, 6.5, 9.6

Modified in 4th Edition: 9.5B

Discontinued in 4th Edition: 6.4C, 9.6B

4. Step-By-Step Examples on Various Computers

This chapter gives four step-by-step examples on several different computer systems. The example will use Algorithm 4.1 - Composite Simpson's Rule for Integration ("041.c") and will compute the integral of f(x) = 2*cos(x) from 1 to 2 using 20 intervals.

Eight steps are typical every time an algorithm is used. These steps are:

Step #1 - Change to Correct Directory (operating system)

Step #2 - Retrieve Algorithm (editor)

Step #3 - Edit Algorithm (editor)

Step #4 - Save Modifications (editor)

Step #5 - Compile Algorithm (compiler)

Step #6 - Run Program (operating system)

Step #7 - View Output (operating system)

Step #8 - Print Output (operating system)

For two-thirds of the algorithms, Steps 2-4 are unnecessary and Step 5 needs to be done only once. These files are marked with an asterisk ('*') in the table in Section 3.1.

The examples below will cover these eight steps on four different computer systems: MS-DOS PCs, UNIX, Macintoshes, and VAXes. Before following any of these examples, first check the need list below and configure your "naautil.c" file.

4.1 Need List

For this example the files "naautil.c" and "041.c" are needed. "Naautil.c" and "041.c" are listed in Appendices A and B to be conveniently referred to during this example. A simple text editor and a C compiler are also required. The C compiler should be ANSI compatible if at all possible. This will save you from possible incompatibility problems.

It is recommended that you try this example out on your computer system as you read this section. Be sure to modify only COPIES of the original algorithms so the algorithms can be used over and over again without problems.

4.2 Customizing Naautil.c

The first decisions to be made are what options and flags you would like to use or set inside the "naautil.c" file. These flags are usually set only once. An explanation of each flag is given below.

ANSI:

If your compiler supports the ANSI C standard, then set ANSI to TRUE. Set ANSI to FALSE only if the program will not compile with it set to TRUE. This flag mostly effects function prototype styles.

ANSI_FUNCT:

Set this flag to TRUE to use the ANSI style for declaring functions over the K&R style. This flag must be set to TRUE if using THINK C 4.0 on a Macintosh.

FILE_SAVE:

If you would like to save the output to a file, then set FILE_SAVE to TRUE. The output is still printed to the screen as you run the program. Set it to FALSE if you do not want to save the output to a file.

TITLE_PROMPT:

If you would like to be prompted for an optional title at the start of each program, then set TITLE_PROMPT to TRUE. This is useful when the output is to be handed in as homework, allowing the user's name or the problem number to be entered. No title is printed to the output file if the [ENTER] key is hit by itself. Set it to FALSE if you do not want to be bothered with entering a title every time you run an algorithm.

EQ_EVAL:

Several of the algorithms require a single function to be evaluated. Set EQ_EVAL to TRUE if you wish to enter the function during run-time instead of at compile time. A couple of simple modifications MUST be made to your algorithm BEFORE this option will be effective. See Chapter 8 - "The Equation Evaluator Routines" for instructions on using this option.

NAAUTIL_OBJ:

This option is useful for users who wish to speed up the compilation process. See Section 6.6 - "Using Naautil.c as Object Code" for more details.

These examples assume the following default settings:

FLAG SETTING

ANSI TRUE

ANSI_FUNCT FALSE

FILE_SAVE TRUE

TITLE_PROMPT TRUE

EQ_EVAL FALSE (Is set to TRUE in "041ee.c")

NAAUTIL_OBJ FALSE

The ANSI, ANSI_FUNCT and OLD_UNIX_OS flags may need to be changed if your compiler varies from the ANSI standard. See Section 6.5 - "Explanation of the Naautil.c File" for a more thorough explanation of the "naautil.c" flags.

4.3 Example Using MS-DOS, Microsoft C and the P-Edit Editor

This example uses the following software:

Operating System: MS-DOS on an IBM PC

Compiler: Microsoft C 5.0

Editor: WordPerfect's P-Edit Editor

No special "naautil.c" flags need to be set.

This example assumes the files were installed onto the "C" drive in the "\NAA42" sub-directory. The DOS prompt will be represented by "C:\NAA42> ".

Step #1 - Change to Correct Directory

Assuming the "Numerical Analysis Algorithm in C" files are located in the "\NAA42" sub-directory of the "C" drive, go there by typing:

C:\> CD \NAA42 - changes directories

C:\NAA42> DIR /P - shows directory's contents

Step #2 - Retrieve Algorithm

Invoke your text editor and retrieve the algorithm file:

C:\NAA42> PE 041.C

The file "041.c" is now loaded and is ready for editing. A text editor is preferred over a word processor. If you plan to use a word processor as your editor, be sure to retrieve and save all files as text-only files.

Step #3 - Edit Algorithm

You must now modify the function f(x). F(x) is listed twice - once as text and once as the actual function call. All functions are defined at the top of each program. To quickly find where modifications are necessary, search for the '$' character. This character is used exclusively for locating lines of code that need updating in all "Numerical Analysis Algorithms in C" files.

Search for the first '$':

[F2] $ [F2] - search

The first '$' should be found at line 22 of "041.c."

Change line 22 from: char *eq_text_f = "f(x) = sin(x)";

to: char *eq_text_f = "f(x) = 2*cos(x)";

This string of text will be printed as output exactly as it appears inside the quotations when the program is run.

Now search for the second '$':

[F2] $ [F2] - search

The second '$' should find the function itself on line 31 of "041.c."

Change line 31 from: return (sin(x));

to: return (2.0 * cos(x));

Step #4 - Save Modifications

Now save the file "041.c" with the above changes and exit the editor:

[F7] Y [ENTER] Y Y - save and exit

Step #5 - Compile Algorithm

Now compile and link "041.c" into the executable file "041.exe." At the prompt type:

C:\NAA42> CL 041.C

The batch file "cc.bat" can also be used in place of the "CL" command. See Sub-Section 7.3.1 on using "cc.bat." If the program requires a larger stack than the default size, using "CL 041.C /link /ST:4096" will increase the stack from 2K bytes to 4K bytes in Microsoft C 5.0.

Step #6 - Run Program

To run "041.exe", at the DOS prompt type:

C:\NAA42> 041

The ".exe" extension can be left off. Answer the prompts with the predetermined inputs. The screen should look something like this:

64444444444444444444444444444444444444444444444444444444444447

5 ‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑ 5

5 "Numerical Analysis Algorithms in C" v4.2 5

5 ‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑‑ 5

5 5

5 Enter an optional title [ie ‑ Set 2.1, Problem 2 a) ]. 5

5 ‑‑‑‑> User's Manual Example 5

5 5

5 Composite Simpson's Rule ‑ Algorithm 4.1 5

5 5

5 f(x) = 2*cos(x) 5

5 5

5 Enter endpoint a: 1 5

5 Enter endpoint b: 2 5

5 Enter number of intervals on [a,b], n: 20 5

5 Interval number h = 0.05 5

5 5

5 !2 5

5 XI = * f(x) dx = 0.13565288875 5

5 "1 5

5 5

5 Required 21 functional evaluations. 5

5 5

5 Output saved into file "041.out". 5

94444444444444444444444444444444444444444444444444444444444448

As indicated by the output, a file named "041.out" is created which contains the output results in a ready-to-print format.

Step #7 - View Output

To view the contents of the output file "041.out", use either the DOS "type" command, the "Numerical Analysis Algorithms in C" utility "list.com", or your text editor. See Section 7.2 for instructions on the usage of the "list.com" utility.

C:\NAA42> TYPE 041.OUT - Using DOS's "type"

C:\NAA42> UTIL\LIST 041.OUT - Using "list.com"

If the file's contents are accurate, then you are ready to print out a copy to be turned in as homework.

Step #8 - Print Output

To print out the output file from DOS, type:

C:\NAA42> PRINT 041.OUT

This step can also be done from within most text editors. WARNING: Be careful not to print the executable file "041.exe". It will waste reams of paper.

4.4 Example Using UNIX, cc and the vi Editor

This example uses the following software:

Operating System: UNIX

Compiler: cc

Editor: vi

You may need to set the OLD_UNIX_OS flag to TRUE if your C compiler requires the include file <varargs.h> instead of <stdarg.h> for variable length argument lists. See your system's "/usr/include" sub-directory to determine which include file will be used.

The percent ('%') character will be used to represent the UNIX shell prompt.

Step #1 - Change to Correct Directory

Assuming the "Numerical Analysis Algorithm in C" files are located in the "naa42" sub-directory, go there by typing:

% cd naa42 - changes directories

% pwd - shows current directory

% ls -alF - shows directory's contents

Step #2 - Retrieve Algorithm

Invoke the vi editor and retrieve the algorithm file:

% vi 041.c

The file "041.c" is now loaded and is ready for editing.

Step #3 - Edit Algorithm

Search for the first '$':

/$ - search

The first '$' should be found at line 22 of "041.c."

Change line 22 from: char *eq_text_f = "f(x) = sin(x)";

to: char *eq_text_f = "f(x) = 2*cos(x)";

This string of text will be printed as output exactly as it appears inside the quotations when the program is run.

Now search for the second '$':

n - search (next occurrence)

The second '$' should find the function itself on line 31 of "041.c."

Change line 31 from: return (sin(x));

to: return (2.0 * cos(x));

Here are a few vi editing commands you should know for future reference:

i Enters insert mode (Exit this mode using [ESC])

R Enters typeover mode (Exit this mode using [ESC])

r Replace character

w Moves forward one word

b Moves backward one word

x Deletes a character

dw Deletes a word

dd Deletes a line

cw Changes a word (follow text by an [ESC] key)

:# Go to line number #

:w Saves (writes) editor contents to a file

:q Quits (exits) the editor

ZZ Exits the editor saving all changes (Same as ":wq")

[ESC] Exits insert, typeover, and other editing modes

/string Searches forward for string

?string Searches backwards for string

n Continue search for string

Arrow keys, ^g, ^h, ^j, ^k, or [SPACE] move the cursor

Step #4 - Save Modifications

Now save the file "041.c" with the above changes and exit the editor:

:wq - write and quit

ZZ - save and exit (faster to type than ":wq")