Turbo Pascal Numerical Methods Toolbox 1986 PDF
Turbo Pascal Numerical Methods Toolbox 1986 PDF
First Edition
Printed in USA
987654321
Turbo Pascal
Numerical Methods
Toolbox~
Copyright © 1986
All Rights Reserved
BORLAND INTERNATIONAL, INC.
4585 SCOTTS VALLEY DRIVE
SCOTTS VALLEY, CALIFORNIA 95066
USA
Table of Contents
Introduction ............................................................................................................... 1
Toolbox Functions ...................................................................................................... 1
About this Manual ..................................................................................................... 2
On the Distribution Disks ........................................................................................ 3
System Requirements ................................................................................................ 3
Acknowledgements .................................................................................................... 4
Table of Contents v
Example .......................................................................................................... 97
Integration Using Adaptive Quadrature and Gaussian Quadrature
(ADAPGAUS.INC) ................................................................................................ 98
Description .......................................................................................................... 98
User-Defined Function ...................................................................................... 98
Input Parameters ................................................................................................ 98
Output Parameters ............................................................................................. 99
Syntax of the Procedure Call ............................................................................. 99
Comments ........................................................................................................... 99
Sample Program ................................................................................................ 101
Example ......................................................................................................... 101
Integration Using the Romberg Algorithm (ROMBERG.INC) ........................ 102
Description ......................................................................................................... 102
User-Defined Function ..................................................................................... 102
Input Parameters ............................................................................................... 102
Output Parameters ............................................................................................ 103
Syntax of the Procedure Call ............................................................................ 103
Sample Program ................................................................................................ 103
Example ......................................................................................................... 103
Table of Contents ix
Solution to an Initial Value Problem for an nth-Order Ordinary Differential
Equation Using the Runge-Kutta Method (RUNGE_N.lNC) ........................ 178
Description ......................................................................................................... 178
User-Defined Types .......................................................................................... 180
User-Defined Function ..................................................................................... 180
Input Parameters ............................................................................................... 180
Output Parameters ............................................................................................ 181
Syntax of the Procedure Call ............................................................................ 181
Comments .......................................................................................................... 181
Sample Program ................................................................................................ 182
Example ......................................................................................................... 182
Solution to an Initial Value Problem for a System of Coupled First-Order
Ordinary Differential Equations Using the Runge-Kutta Method
(RUNGE_Sl.INC) ................................................................................................ 186
Description ............................... ;......................................................................... 186
User-Defined Types .: ........................................................................................ 188
User-Defined Functions ................................................................................... 188
Input Parameters ............................................................................................... 189
Output Parameters ............................................................................................ 189
Syntax of the Procedure Call ............................................................................ 190
Comments .......................................................................................................... 190
Sample Program ................................................................................................ 191
Example ......................................................................................................... 191
Solution to an Initial Value Problem for a System of Coupled Second-Order
Ordinary Differential Equations Using the Runge-Kutta Method
(RUNGE_S2.INC) ............................................................................................... 196
Description ......................................................................................................... 196
User-Defined Types .......................................................................................... 199
User-Defined Functions ................................................................................... 199
Input Parameters .............................................................................................. 200
Output Parameters ............................................................................................ 201
Syntax of the Procedure Call ............................................................................ 201
Comments .......................................................................................................... 201
Sample Program ............................................................................................... 203
Example ........................................................................................................ 203
Solution to Boundary Value Problem for a Second-Order Ordinary Differential
Equation Using the Shooting and Runge-Kutta Methods (SHOOT2.1NC) .. 208
Description ........................................................................................................ 208
User-Defined Types ......................................................................................... 209
User-Defined Function .................................................................................... 209
Input Parameters .............................................................................................. 209
Output Parameters ............................................................................................ 210
Syntax of the Procedure Call ............................................................................ 210
Comments .......................................................................................................... 210
Table of Contents xi
Description ................................................................................................... 242
Input Parameters ......................................................................................... 242
Output Parameters ...................................................................................... 242
Syntax of the Procedure Call ...................................................................... 242
Procedure FFT ................................................................................................. 242
Description ................................................................................................... 242
Input Parameters ......................................................................................... 243
Output Parameters ...................................................................................... 243
Syntax of the Procedure Call ...................................................................... 243
Fast Fourier Transform Applications ................................................................... 244
COMPFFT.INC ............................................................................................... 244
Description ................................................................................................... 244
Input Parameters ......................................................................................... 244
Output Parameters ...................................................................................... 244
Syntax of the Procedure Call ...................................................................... 245
Comments .................................................................................................... 245
REALFFT.INC ................................................................................................. 245
Description ................................................................................................... 245
Input Parameters ......................................................................................... 245
Output Parameters ...................................................................................... 246
Syntax of the Procedure Call ...................................................................... 246
Comments .................................................................................................... 246
COMPCNVL.INC ........................................................................................... 246
Description ................................................................................................... 246
Input Parameters ......................................................................................... 247
Output Parameters ...................................................................................... 247
Syntax of the Procedure Call ...................................................................... 247
Comments .................................................................................................... 247
REALCNVL.INC ............................................................................................. 248
Description ................................................................................................... 248
Input Parameters ......................................................................................... 248
Output Parameters ...................................................................................... 248
Syntax of the Procedure Call ...................................................................... 249
Comments .................................................................................................... 249
COMPCORR.INC ........................................................................................... 249
Description ................................................................................................... 249
Input Parameters ......................................................................................... 250
Output Parameters ...................................................................................... 250
Syntax of the Procedure Call ...................................................................... 250
Comments .................................................................................................... 250
REALCORR.INC .............................................................................................. 251
Description .................................................................................................... 251
Input Parameters .......................................................................................... 251
Output Parameters ...................................................................................... 252
The Turbo Numerical Methods Toolbox is a reference manual for both the student
of numerical analysis and the professional needing efficient routines. An elemen-
tary background in calculus and linear algebra is assumed, although many of the
algorithms use only high-school-level mathematics. A general knowledge of Turbo
Pascal® is also assumed. If you need to brush up on your knowledge of Pascal, we
suggest looking at the Turbo Pascal Reference Manual and/or the Turbo Pascal
Tutor Manual.
Before you begin using a particular routine, read through this brief introductory
chapter and then refer to the chapter that interests you.
Toolbox Functions
The major areas in numerical analysis are represented in this Toolbox, with each
chapter focusing on a particular problem. Each routine begins with a general
description of the implemented algorithm or numerical method. (References to
numerical analysis texts are provided for each numerical procedure.) User-supplied
types, functions, and input and output parameters are defined, and the syntax of
the procedure call is provided. If appropriate, a "Comments" section is also pro-
vided.
Finally, every algorithm in the Toolbox is accompanied by a general-purpose
program that handles all the necessary I/O, while allowing you to try each algo-
rithm without building any code. Handily, these sample programs will often reduce
the coding your own application may require.
As an example, let's say you want to find the roots to an equation in one variable.
First, you would read the introduction to Chapter 2, "Roots to Equations in One
Variable," and choose the numerical method best suited to your particular problem.
Second, you would run the sample program for the desired numerical method to
determine the necessary input and output. Third, you would write a Turbo Pascal
function defining your equation, using the function already coded in the sample
program as a guide. Fourth, you would run the sample program with your function
substituted for the original one. Of course, if these algorithms are to be part of a
larger program, you must build all the interfaces to the other parts of the system;
but this should only be done after you gain experience with the particular numeri-
cal method.
Several books are referred to throughout the text; complete references are listed
at the back of the book in the section entitled "References."
The body of this manual is printed in normal typeface; other typefaces serve to
illustrate the following:
Alternate This type displays program examples and procedure and function
declarations.
Italics This type emphasizes certain concepts, first-mentioned terms, and
mathematical expressions.
Boldface This type marks the reserved words of Turbo Pascal in text and in
program examples.
The routines for this Toolbox are contained on three packed disks. Their contents
and general installation instructions are covered in Chapter l.
System Requirements
All routines will run in standard Turbo Pascal version 3.0. (They will also run in
version 2.0, but you must make one change to use the sample programs; see the
section entitled "Installation" in Chapter 1.) All sample programs will run on an
IBM® PC or compatible machine using DOS 2.0 or greater.
A small portion of the Toolbox uses graphics (see Chapter 11). These programs
are for PC-DOS users only, requiring either an IBM PC Color Graphics Adapter,
an IBM Enhanced Graphics Adapter, or a Hercules Monochrome Graphics
Adapter. Recompiling these requires the Turbo Pascal Graphix Toolbox® version
1.06A or later, and Turbo Pascal 3.0. They can also be recompiled for the EGA and
several other cards using version 1.07A of the Graphix Toolbox.
We strongly recommend that anyone serious about numerical analysis invest in
hardware and software to run Turbo Pascal with support for the Intel 8087 numeri-
cal processing chip:
• Hardware: an 8087 chip plugged into the motherboard of a PC, XT® or equiva-
lent, or an 80287 chip in an AT® or equivalent.
• Software: TURBO-87.COM, the version of Turbo Pascal designed to take advan-
tage of the 8087 chip.
For machines running the Intel 8088 CPU, the increase in execution speed of
programs using real-number arithmetic is often a factor of ten or more, while for
80286 machines, the increase is only about a factor of two (Fried 1985). Perhaps
more important than speed is the increase in accuracy -16 significant figures accu-
racy for Turbo Pascal with 8087 support versus 11 significant figures for standard
Turbo. Since round-off errors are a serious concern in numerical analysis, the
increased accuracy is of great value.
All of the examples in this manual were run using Turbo-87. If you run them
without the Turbo-87, you will usually get less accurate answers.
Introduction 3
AcktWwlellgements
We refer to several products throughout the manual; following is a list of each one
and its respective company.
• Turbo Pascal, Turbo Pascal Graphix Toolbox, SideKick, SuperKey, and
Reflex: The Database Manager are registered trademarks and Turbo Pascal
Numerical Methods Toolbox is a trademark of Borland International, Inc.
• IBM, XT, and AT are registered trademarks of International Business Machines
Corp.
• MS-DOS is a registered trademark of Microsoft Corp.
This chapter provides you with everything you need to start using the routines in
this Toolbox. We'll discuss how to unpack the disks for use and list the files avail-
able once the disks are unpacked. We also briefly discuss data types and defined
constants used in the Toolbox, and the setting of compiler directives.
First, though, before we thrust you into the middle of numerical madness, let's
take a look at one way to use this Toolbox.
In late 1986 and early 1987, the America's Cup 12-meter yacht championship was
held. The 12-meter yachts are just large sailboats, but the competition is so intense
that the only way to be competitive is to use dozens of people, spend millions of
dollars, design a special boat, and spend a couple of years training for the race. The
race has become so sophisticated that many of the sailboats have on-board com-
puters and other electronic equipment.
To keep stride with other challengers, one yacht's crew used personal com-
puters, and of course, Borland software. They used Turbo Pascal to design the
boat's hull. They used Reflex: The Analyst® to maintain their databases and to dis-
play plots while the boat was sailing. And when it came time to do some mathemat-
ical modeling, again they turned to Borland for its inimitable software and chose
the Toolbox.
5
Simply speaking, the problem they had was one of "precision monitoring." It
takes a crew of very highly skilled sailors to compete in America's Cup races, but
even the best skippers cannot act with sufficient precision to win. A typical race
lasts for several hours, and the winner usually wins by only a few feet.
The electronic equipment on a boat can sense with reasonable accuracy all of the
crucial variables: boat velocity, wind velocity, boat direction, boat position, and so
on. This data must then be made available to the skipper in a coherent form, and
he/she must decide at what angle to place the rudder based on that information.
The problem is too complex to rely on intuition alone.
Even just displaying the velocity is more complex than you might think at first.
When sailing on the ocean, the waves are big enough that the velocity is in constant
flux. Fortunately, the fluctuations due to the waves represents a steadily periodic
force. By using the Fourier transforms in Chapter 10 of the Toolbox, the crew was
able to identify the periodic portion of the velocity and subtract it out. The result:
the velocity as a function of time but with the wave fluctuations eliminated. The
graph of this modified velocity is much smoother, and allows the skipper to tell
much more quickly and accurately whether the boat is accelerating or decelerating.
To measure the acceleration quantitatively, the crew used the fact that the accel-
eration is the derivative of the velocity. They were able to do this easily with the
differentiation routines in Chapter 4 of the Toolbox. They were also able to directly
measure the distance travelled by using the integration routines in Chapter 5, and
the fact that distance is the integral of the speed.
Perhaps the most difficult problem in navigating a sailboat is aiming the rudder.
You can't just aim the boat in the direction that you want to go, rather you have to
pick a direction that you can sail rapidly, depending on the wind direction. An
experienced skipper can judge this pretty well, but not well enough. Every boat is
a little different, and the best way to handle one boat is not necessarily the best way
to handle another.
So, the team ran extensive trial races with the boat to gather data on how the
boat performed in various circumstances. The data was collected automatically by
electronic instruments on board, and stored digitally on floppy disks. They then
used Reflex to manage the data and to display graphs. But they lacked the tools to
relate their data to the data they would have under actual racing conditions.
In order to predict the behavior of their boat in an actual race, the team created a
model from their collected data using the least-squares routines in Chapter 9 of the
Toolbox. With the least-squares routines, you can create a multiparameter model
and then find the values of the parameters that make the model most accurately fit
the data. With a mathematical model of the boat's behavior, the team was then able
to predict how the boat would perform under circumstances similar but not identi-
cal to its practices.
All of the Toolbox routines are contained on three disks. (Note, however, that
MS-DOS® users will receive two disks; Disk 3 is for PC-DOS users only.) Each
disk has packed files corresponding to chapters in the manual. Use the program
UNPACK.EXE to unpack the files, described in the next section, "Installation."
The files for each chapter are self-contained and do not require any files from
any other chapter, with these exceptions:
\
Routine Beginnings 7
Disk 2
UNPACK.EXE (installation program to unpack chapters)
CHAP8 (packed file with routines for Chapter 8)
CHApg (packed file with routines for Chapter 9)
CHAPIO (packed file with routines for Chapter 10)
CHAP11 (packed file with routines for Chapter 11; PC-DOS users only)
Disk 3 (PC-DOS users only)
README
README.COM (program to display README file)
LSQIBM.COM (requires IBM graphics monitor)
FFTIBM.COM (requires IBM graphics monitor)
LSQHERC.COM (requires Hercules graphics card)
FFTHERC.COM (requires Hercules graphics card)
SAMP11A.DAT (data file for LSQ*.COM)
SAMP11B.DAT (data file for FFT*.COM)
14X9.FON (from the Graphix Toolbox)
4X6.FON (from the Graphix Toolbox)
8X8.FON (from the Graphix Toolbox)
ERROR.MSG (from the Graphix Toolbox)
Installation
The files CHAP2 through CHAP11 are packed files corresponding to the chapters
in this manual. In order to use these files, you must first unpack them with
UNPACK.EXE. The syntax is as follows:
UNPACK packed-file-name target-drive
For example, the files for Chapter 2 can be extracted and put on the current
directory on drive C by placing Disk 1 in drive A, changing the logged drive to A by
typing A: at the DOS prompt, and then typing
UNPACK CHAP2 c:
Wildcards are okay to use for the packed file name, and a directory can be
specified with the target drive if it ends with a backslash (\). For example, all of the
packed files on disk can be placed in a directory C:\NUMERIC by doing the
following:
UNPACK CHAP* C:\NUMERIC\
if the directory C:\NUMERIC already exists.
Routine Beginnings 9
CHAP5 (packed file) "Numerical Integration"
ADAPGAUS.lNC ROMBERG.PAS
ADAPGAUS.PAS SIMPSON.lNC
ADAPSIMP.lNC SIMPSON.PAS
ADAPSIMP.PAS TRAPZOID.lNC
ROMBERG.INC TRAPZOID.PAS
All sample programs call the include file COMMON.INC from the disk. This
file includes procedures that are common to all sample programs. When copying
any of the sample programs to a disk, be sure to also copy the file COMMON.INC
to that disk or the sample programs will not compile.
To use the sample programs with Turbo Pascal version 2.0, rename COM-
MON2.INC (which is on Disk 1) to COMMON.INC. (You may wish to preserve a
copy of tlie original COMMON.INC file by first copying it to a file called COM-
MON3.INC). If you run the sample programs with version 2.0 and do not make
this change, the programs will compile but will handle I/O errors incorrectly.
We have made the sample programs general and easy to use. For example,
numerical input can originate from the keyboard (where improper input is trap-
ped) or from a text file; output can be sent to the printer, screen, or text file; other
refinements are also included. Since, to a beginner, the supporting code may
obscure the simplicity of calling the procedure, we have included a minimal sample
program for Newton-Raphson's method of root-finding (RAPHSON2.INC).
Routine Beginnings II
The Graphics Denws
Because graphic displays are often an essential part of numerical analysis, we have
included two demonstration programs (for PC-DOS users only) that involve dis-
play numerical results. As previously stated, graphics hardware is not necessary for
this Toolbox, but it is required for these two graphics programs. The programs are
built with subsets of the Turbo Pascal Graphix Toolbox; there are separate versions
for systems with the Hercules Monochrome Graphics card and the IBM Color
Graphics card (or good emulations of these cards).
The demonstration programs are on Disk 3. For instructions about how to run or
recompile them, see Chapter 11.
Data types that might be confused with those in the calling program have been
prefixed with the letters TN (for Turbo Numerical); for example, TNmatrix or
TNvector. You must define these variable types in your top-level program for two
reasons. First, you will probably need to use this type in your top-level program,
and the type must be defined to have the same scope as the Toolbox procedure.
Second, you will want to dimension arrays based on your particular needs. For
example, the Lagrange procedure requires the definition
type TNvector = array[O .. TNArraySize] of Real;
The identifier TNArraySize is never referred to in any of the include files. It
should be optimized by the user, although we have set a default value in each of the
sample programs. It may be replaced with an integer or byte.
TNNearlyZero is the only defined constant that must be changed when standard
Turbo Pascal is used (instead of Turbo-87); it should be changed to IE - 7. (With-
out changing this constant, you will get a syntax error when compiling with stan-
dard Turbo Pascal.) For Turbo-87, the default value of this constant is IE -15.
(There are a few exceptions to these default values; where appropriate, the "Com-
ments" section of each routine indicates the exceptions.)
Aside from the usual default values of the compiler directives in standard Turbo
Pascal, we have set the compiler directive to {$R + } in all include files that use
arrays, and to {$I-} in all sample programs. The first directive checks to see that all
array-indexing operations are within the defined bounds and all assignments to
scalar and subrange variables are within range. The latter directive disables I/O
error-checking. All the sample programs have their own I/O error-checking proce-
dures (loaded in the file COMMON.lNC), so that the {$I-} directive must remain
disabled in the sample programs. The array checker {$R +} should always be
active, since the performance penalty is slight and the advantages are significant.
Routine Beginnings 13
14 Turbo Numerical Methods Toolbox
c H A p T E R 2
Roots to Equations in One Variable
The routines in this chapter are for finding the roots of a single equation in one real
variable. A typical problem is to solve
x * exp(x) - 10 = 0 .
In general, the routines find a value of x, where x is a scalar real variable,
satisfying
f(x) = 0.0
where f is a real-valued function that you program in Pascal.
All of the methods are approximate methods, meaning that they find an approxi-
mate value of x that makesf(x) close to zero. Because of round-off error, it is usually
not possible to find the exact value of x. Furthermore, they are all iterative
methods, meaning that you specify some initial guess that is some value for x,
which you think is reasonably close to the solution. The routine repeats some calcu-
lations that replace the guess x with a more accurate guess until the required level
of accuracy is achieved.
The bisection method (BISECT.INC) returns an approximation to a root of a real
continuous function of the real variable x. This method always converges (as long as
the function changes signs at a root), but may do so relatively slowly.
The Newton-Raphson method (RAPHSON.INC) also returns an approximation
to a root of a real functionf of the real variable x. When this algorithm converges, it
is usually faster than the bisection method. If more than one root of a polynomial
equation is desired, then use Newton-Horner's method (NEWfDEFL.INC).
15
The secant method (SECANT.lNC) is similar to the Newton-Raphson method,
but doesn't require knowledge of the first derivative of the function. Consequently,
it is more flexible than the Newton-Raphson method, though somewhat slower.
Newton-Horner's method (NEWfDEFL.lNC) applies Newton's method to
real polynomials. It also uses deflation techniques to attempt to approximate all the
real roots of a real polynomial. Both the Newton-Horner and Newton-Raphson
methods are faster than the bisection and secant methods, but are undefined if
If'(x)1 < = TNNearlyZero. This is less of a problem on machines with a high-
precision math coprocessor, since TNNearlyZero is smaller.
The Newton-Horner and Newton-Raphson methods both converge around mul-
tiple roots, although convergence is slow. These algorithms depend upon an initial
approximation of the root. If the initial approximation is not sufficiently close to the
root, the Newton methods may not converge. In some instances, an initial choice
may lead to successive iterations that oscillate indefinitely about a value of x usu-
ally associated with a relative minimum or relative maximum off. In either case,
the bisection method could be used to determine the root or to determine a close
approximation to the root that can be employed as an initial approximation in the
Newton-Raphson or Newton-Horner methods.
Muller's method (MULLER.lNC) returns an approximation to a root (possibly
complex) of a complex function of the complex variable x. Although Muller's
method can approximate the roots of polynomials, we recommend that you use
Newton-Horner's method, the secant method, or (in the case of complex polyno-
mials) Laguerre's method to find the roots of polynomials.
Laguerre's method (LAGUERRE.lNC) attempts to approximate all the real and
complex roots of a real or complex polynomial. Laguerre's method is very reliable
and quick, even when converging to a multiple root. This is the best general
method to use with polynomials.
A caution when solving polynomial equations: Polynomials can be ill-
conditioned, in the sense that small changes in the coefficients may lead to large
changes in the roots.
All the root-finding routines use the function TestForRoot to determine if a root has
been found.
function TestForRoot{X, OldX, V, Tol : Real) : Boolean;
(********************************************************************)
(* Here are four stopping criteria. If you wish to *)
(* change the active criteria, simply comment off the current *)
{* criteria (including the appropriate or) and remove the comment *)
{* brackets from the criteria (including the appropriate or) you *)
(* wish to be active. *)
(********************************************************************)
begin
TestForRoot := (***************************)
(ABS{V) <= TNNearlyZero) (* V=O *)
(* *)
or (* *)
(* *)
(ABS{X - OldX) < ABS{OldX*Tol)) (* relative change in X *)
(* *)
(* *)
{* or *) (* *)
{* *) (* *)
{* (ABS{X - OldX) < Tol) *) (* absolute change in X *)
{* *) (* *)
{* or *) (* *)
{* *) (* *)
(* (ABS{V) <= Tol) *) (* absolute change in V *)
(***************************)
end; { procedure TestForRoot }
The four separate tests provided by function TestForRoot may be used in any
combination. The default criteria tests the absolute value of Y and the relative
change in X. If you wish to change the active criteria, simply comment off the
current criteria (including the appropriate or) and remove the comment brackets
from the criteria (including the appropriate or) you wish to be active.
The first criterion simply checks to see if Y is zero (TNNearlyZero is defined at
the beginning of the procedure). This criterion should usually be kept active.
The second criterion examines the relative change in X between iterations. To
avoid division by zero errors, QldX has been multiplied through the inequality.
The third criterion checks the absolute change in X between iterations.
The fourth criterion determines the absolute difference between Y and the
allowable tolerance. Note: The parameter Tol(erance) means something different in
each test. Be sure you know which tests are active when you input a value for Tol.
Description
This method (Burden and Faires 1985, 28 ff.) provides a procedure for finding a
root of a real continuous function f, specified by the user on a user-supplied real
interval [a,b]. The functionsfia) and fib) must be of opposite signs. The algorithm
successively bisects the interval and converges to the root of the function. You must
also specify the desired accuracy to which the root should be approximated.
User-Defined Function
Input Parameters
Comrnents
Sample Program
The sample program BISECT.PAS provides I/O functions that demonstrate the
bisection algorithm. To modify this program for your own function, simply change
the definition of function TNTargetF. Note that the file BISECT.INC is included
after the function TNTargetF is defined.
{---------------------------------------------}
2. Run BISECT.PAS:
Left endpoint: 0
Right endpoint: 100
Tolerance (> 0, default = 1.000E-08): 1E-6
Maximum number of iterations (>= 0, default = 100)? 100
Direct output to one of the following:
(S)creen
(P)rinter
(F) il e
Description
This example uses Newton-Raphson's algorithm (Burden and Faires 1985,42 fr.) to
find a root of a real user-specified function when the derivative of the function and
an initial guess are given. The algorithm constructs the tangent line at each iterate
approximation of the root. The intersection of the tangent line with the x-axis
provides the next iterate value of the root. You must specify the desired tolerance to
which the root should be approximated.
User-Defined Functions
Input Parameters
Comments
Newton's method involves division by the value of the derivative of the function.
Should the algorithm attempt to do any calculations at a point where the deriva-
tive is less than TNNearlyZero, the routine will stop and return an error message
(Error = 2).
Convergence is determined with the Boolean function TestForRoot described at
the beginning of this chapter.
Sample Program
The sample program RAPHSON .PAS provides I/O functions that demonstrate the
Newton-Raphson algorithm. Note that the file RAPHSON.INC is included after
the functions TNTargetF and TNDerivF are defined.
The program RAPHSON2.PAS also provides I/O functions that demonstrate the
Newton-Raphson method. It is an extremely bare-bones program and is provided
Example
Problem. Determine the solution to the equation cos(x) = x.
1. Code the following two functions into RAPHSON.PAS (or RAPHSON2.PAS):
{---------- HERE IS THE FUNCTION -------------}
function TNTargetF{x : Real) : Real;
begin
TNTargetF := Cos (x) - x;
end; { function TNTargetF }
{---------------------------------------------}
{-------- HERE IS THE DERIVATIVE -------------}
function TNDerivF{x : Real) : Real;
begin
TNDerivF := -Sin{x) - 1;
end; { function TNDerivF }
{---------------------------------------------}
2. Run RAPHSON.PAS:
Initial approximation to the root: 0
Tolerance (> 0, default = 1.000E-08): 1E-6
Maximum number of iterations (>= 0, default = 100): 100
Direct output to one of the following:
{S)creen
{P)rinter
{F)i1e
Number of iterations: 5
Calculated root: 7.39085133215161E-001
Value of the function
at the calculated root: O.OOOOOOOOOOOOOOE+OOO
Value of the derivative
of the function at the
calculated root: -1.67361202918321E+000
Description
This example uses the secant method (Gerald and Wheatley 1984, 11-13) to find a
root of a user-specified real function given two initial real approximations to the
root. The secant method constructs a secant through the two points specified by
the initial approximations. The intersection of this line and the x-axis is used as the
next best approximation to the root. The approximation to the root and its prede-
cessor are used to construct the next secant line. The process continues until a root
is approximated with specified accuracy or until a specified number of iterations
have been exceeded.
User-Defined Function
Input Parameters
Comments
The secant algorithm constructs a line through two points and finds the intersec-
tion of that line with the x-axis. If the line has a slope whose absolute values are
less than TNNearlyZero (that is, the two points have the same y-value), then it has
no intersection with the x-axis (or infinitely many if it lies on the x-axis) and the
algorithm will no longer continue. If this happens, Error 2 is returned. Error 2 will
also be returned if the absolute difference of the two initial approximations (Guessl
and Guess2) is less than TNNearlyZero.
Convergence is determined with the Boolean function TestForRoat described at
the beginning of this chapter.
Sample Program
The sample program SECANT.PAS provides I/O functions that demonstrate the
secant algorithm. Note that the file SECANT.lNC is included after the function
TNTargetF is defined.
{---------------------------------------------}
2. Run SECANT.PAS:
First initial approximation to the root: 0
Second initial approximation to the root:
Tolerance (> 0, default = 1.000E-08): 1E-8
Maximum number of iterations (>= 0, default 100): 100
Direct output to one of the following:
{S)creen
(P)rinter
(F) il e
Description
User-Defined Types
Input Parameters
Output Parameters
Degree: Integer: Degree of the deflated polynomial (> 2 if some of the roots are
not approximated).
NumRoots: Integer: Number of roots found.
Poly:TNvector; Coefficients of the deflated polynomial.
Root:TNvector: Real part of all roots found.
Imag:TNvector: Imaginary part of all roots found (nonzero for 2 at most).
Value:TNvector: Value of the polynomial at each approximate root.
Deriv:TNvector: Value of the derivative at each found root.
Iter:TNIntVector: Number of iterations required to find each root.
Error:Byte: 0: No error.
1: Maximum number of iterations exceeded.
2: The slope is zero (see "Comments").
3: Degree =:; O.
4: Tol =:; O.
5: Maxlter < O.
If a root is found, it is returned with the value of the polynomial at that root
(which should be close to zero) and with the value of the derivative at that root. If
the last two roots are complex (only two can be complex, since they are evaluated
by the quadratic formula), then the value and derivative at those points are arbi-
trarily set to zero. If all the roots have not been found, then the unsolved deflated
polynomial is also returned.
Comments
Newton's method involves division by the derivative of the function. Should the
algorithm attempt to do any calculations at a point where the absolute values of the
derivative are less than TNNearlyZero, the routine stops and returns an error mes-
sage (Error = 2).
Convergence is determined with the Boolean function TestForRoot described at
the beginning of this chapter.
Sample Program
Input Files
It is possible to input the coefficients from a text file. The format for the text file is
as follows:
1. The degree of the polynomial
2. The coefficients in descending order, beginning with the leading coefficient
and decreasing to the constant term
Spaces or carriage returns can be used to separate the data. It does not matter
whether the file ends with a carriage return; for example, the polynomial
F(x) =x 3
- 2x
could be entered in a text file as
310 -20
Example
Problem. Determine the roots to the 7th degree polynomial:
Xfi + x5 - 49x·' + 69x3 + 120x2 + 98x - 240
Initial Polynomial:
Poly[6]: 1.00000000000000E+000
Poly[5]: 1.00000000000000E+000
Poly[4]: -4.90000000000000E+001
Poly[3]: 6.90000000000000E+001
Poly[2]: 1.20000000000000E+002
Poly[l]: 9.80000000000000E+001
Poly[O]: -2.40000000000000E+002
Initial approximation: O.OOOOOOOOOOOOOOE+OOO
Tolerance: 1.00000000000000E-008
Maximum number of iterations: 100
Number of calculated roots: 6
Root 1
Number of iterations: 7
Calculated root: 3.00000000000000E+000
Value of the function
at the calculated root: 4.83169060316868E-013
Value of the derivative
of the function at
the calculated root: -7.47999999999999E+002
Description
This example uses Muller's method (Burden and Faires 1985, 71-75) to find a
possibly complex root of a user-defined complex function. The algorithm finds a
root of a parabola defined by three distinct points of the given function. This
approximation to the root and its two predecessors are used to construct the next
parabola. This is repeated until the convergence criteria is satisfied. Muller's
method has the advantage of nearly always converging; however, it is slow because
it uses complex arithmetic. You must create a complex function, input an initial
guess (which need not be very accurate), the tolerance in the answer, and the
maximum number of iterations.
TNcomplex = record
Re, Im:Real;
end;
Input Parameters
Output Parameters
Comments
Muller's method involves constructing a parabola from three points. If they all lie
on a line whose slope in absolute value is less than TNNearlyZero, then a parabola
that intersects the x-axis cannot be constructed. Such a construction will halt the
algorithm and return Error = 2. Fortunately, this does not commonly occur.
Convergence is determined with the Boolean function TestForRoot described at
the beginning of this chapter. Complex arithmetic is used.
Example
Problem. Find a solution to the complex equation cos(x) = x.
1. First, code the following procedure TNTargetF into MULLERPAS:
(*------------- HERE IS THE FUNCTION ------------------*)
procedure TNTargetF(x : TNcomplex; var y : TNcomplex);
beg;n { this is the complex function y = Cos (x) - x }
y.Re := Cos(x.Re)*(Exp(-x.lm) + Exp(x.lm))/2 - x.Re;
y.lm := Sin(x.Re)*(Exp(-x.lm) - Exp(x.lm))/2 - x.lm;
end; { procedure TNTargetF }
(*-----------------------------------------------------*)
2. Run MULLERPAS:
Initial approximation to the root:
Re(Approximation)= -4
Im(Approximation)= 4
Tolerance (> 0, default = 1.000E-08): 1E-6
Maximum number of iterations (>= 0, default = 100): 100
Number of iterations: 18
Calculated root: -9.10998745393294EtOOO t 2.95017086170180EtOOOi
Value of the function
at the calculated root: -1.42534872793476E-011 t 3.75033337718378E-011i
Description
This example uses Laguerre's method (Ralston and Rabinowitz 1978, 380-383) and
linear deflation to find the possibly complex roots of a complex (or real) polynomial.
You must input the coefficients of the polynomial, an initial guess, the tolerance
with which to find the answer, and the maximum number of iterations.
User-Defined Types
TNcomplex = record
Re, Im:Real;
end;
TNIntVector = array[O •• TNArraySize] of Integer;
TNCompVector = array[O •• TNArraySize] of TNcomplex;
Input Parameters
4. degree ~ TNArraySize
Output Parameters
Comments
For some polynomials, certain starting values (Guess) will not yield convergence. If
the routine does not converge to a solution, try a different starting value. Note that
convergence is slower around multiple roots than around single roots.
Convergence is determined with the Boolean function TestForRaat described at
the beginning of this chapter.
Input Files
It is possible to input the coefficients from a text file. The format for the text file is
as follows:
1. The degree of the polynomial
2. The real and imaginary parts of the coefficients in descending order, begin-
ning with the leading coefficient and descending to the constant term
Spaces or carriage returns can be used to separate the data. It does not matter
whether the file ends with a carriage return; for example, the polynomial
F(x) = x-l - (2 + 2i)x3 + 4ix2 + (2 - 2i)x -1
where i represents the square root of - 1, could be entered in a text file like this:
410 -2 -2042 -2 -10
Example
Problem. Find all the roots to the complex polynomial
F(x) = l - (2 + 2i)x3 + 4ix2 + (2 - 2i)x - 1
where i is the square root of - 1.
Run LAGUERRE.PAS:
(K)eyboard or (F)ile input of data? K
Degree of the polynomial «= 30)? 4
Input the coefficients of the polynomial
where Poly[n] is the coefficients of x~n
Re(Poly[4]) = 1
Im(Poly[4]) = 0
Re(Poly[3]) = -2
Im(Poly[3]) = -2
Re(Poly[2]) 0
Im(Poly[2]) = 4
Initial approximation:
Re{Approximation) = 0
Im{Approximation) = 0
Tolerance (> 0, default = 1.000E-08): 1E-6
Maximum number of iterations (>= 0, default = 100): 100
Direct output to one of the following:
{S)creen
{P)rinter
(F) il e
Initial Polynomial:
Poly[4]: 1.00000000000000E+000 + O.OOOOOOOOOOOOOOE+OOOi
Poly[3]: -2.00000000000000E+000 + -2.00000000000000E+000i
Poly[2]: O.OOOOOOOOOOOOOOE+OOO + 4.00000000000000E+000i
Poly[l]: 2.00000000000000E+000 + -2.00000000000000E+000i
Poly[O]: -1.00000000000000E+000 + O.OOOOOOOOOOOOOOE+OOOi
Interpolation is useful when some values of a function are known but others are
required. For example, suppose values are known for a functionJ(x) at x = 2.3,2.4,
2.5, 2.6, 2.7, 2.8, and the value ofJ(x) is desired at x = 2.415. The routines in this
chapter provide the means to model to given values of J(x) with an appropriate
function, so that the function can be evaluated at other arbitrary points.
The goal of interpolation is to approximate the value of the function at a speci-
fied value of x, given N values of the function at N distinct points. This approxima-
tion will be a polynomial determined from the input data. The value of the
polynomial at x will be returned as the approximation to the value ofJ(x).
The Lagrange method (LAGRANGE.INC) accepts points in any order. The x-
values need not be equally spaced. An interpolating polynomial is explicitly calcu-
lated. Although an interpolating polynomial can be useful for computing deriva-
tives (and more), the Lagrange method is a lengthy process. Furthermore, high-
degree polynomials may cause significant round-off error in some interpolations.
Newton's general divided-difference algorithm (DIVDIF.INC) does not require
input to have equally spaced x-values, nor is it necessary that the x-values be
in either ascending or descending order. For large amounts of data, the divided-
difference routine (DIVDIF.INC) is more accurate than Lagrangian interpola-
tion (LAGRANGE.INC).
If there are many input points, the Lagrange (LAGRANGE.INC) and the
divided-difference (DIVDIF.INC) methods may result in high-degree polynomials
whose oscillatory nature can produce an inaccurate approximation. This is espe-
43
cially true if the interpolation occurs at a value near the midpoint between adjacent
input x-values. In such cases, the cubic spline methods (CUBEJ'RE.INC and
CUBE_CLA.INC) are preferable.
The cubic spline methods require that the x-values be entered in ascending
order. The clamped cubic spline method (CUBE_CLA.lNC) may yield more accu-
rate results than the free cubic spline method (CUBEJ'RE.lNC), but requires
knowledge of the first derivative of the function at the endpoints of the input data.
When this information is not available; the free cubic spline routine should be
used.
The values at which interpolation is to occur should lie in the closed interval
bounded by the extreme values of the input x-values. The preceding methods will
not give accurate approximations to values outside this interval (extrapolation).
Description
This example provides an interpolation algorithm (Burden and Faires 1985, 84 If;
Horowitz and Sahni 1984,429-430). Given a set of N data points (x,y), the routine
uses Lagrange polynomials to construct a polynomial to fit the data points. The
degree of the polynomial is at most N - l.
Note: The nature of high-degree polynomials may cause significant error if the
algorithm is used with large amounts of data (about N > 25). In such cases, DIV
DIF.INC, CUBEJ'RE.INC, or CUBE_CLA.INC should be used. You must sup-
ply the data points and the x-values at which interpolation will take place.
Input Parameters
Interpolation 45
TNArraySize fixes an upper bound on the number of elements in each vector. It
is used in the type definition of TNvector. TNArraySize is not a variable name and
is never referenced by the procedure; hence there is no test for condition 2. If
condition 2 is violated, the program will crash with an Index Out of Range error
(assuming the directive {$R + } is active).
Output Parameters
Sample Program
Input Files
Data may be entered from a text file. The x and y coordinates should be separated
by a space and followed by a carriage return. For example, data values of sqr(x)
could be entered in a text file as:
11
24
39
416
525
The data:
1.0000000 9.99847695156391E-001
2.0000000 9.99390827019096E-001
3.0000000 9.98629534754574E-001
4.0000000 9.97564050259824E-001
5.0000000 9.96194698091746E-001
6.0000000 9.94521895368273E-001
7.0000000 9.92546151641322E-001
8.0000000 9.90268068741570E-001
9.0000000 9.87688340595138E-001
10.0000000 9.84807753012208E-001
11.0000000 9.81627183447664E-001
12.0000000 9.78147600733806E-001
13.0000000 9.74370064785235E-001
14.0000000 9.70295726275996E-001
15.0000000 9.65925826289068E-001
16.0000000 9.61261695938319E-001
17.0000000 9.56304755963035E-001
18.0000000 9.51056516295154E-001
19.0000000 9.45518575599317E-001
20.0000000 9.39692620785908E-001
The polynomial:
Poly[19]= -1.72986376643586E-028
Poly[18]= 3.57504241395844E-026
Poly[17]= -3.43926153199199E-024
Poly[16]= 2.04507188280402E-022
Poly[15]= -8.41710490427928E-021
Poly[14]= 2.54454663251946E-019
Poly[13]= -5.85115478257286E-018
Poly[12]= 1.04567701944119E-016
Poly[ll]= -1.47131277721604E-015
Poly[10]= 1.64108936708876E-014
Poly[9]= -1.45382580196402E-013
Poly[8]= 1.02034999744286E-012
Poly[7]= -5.63354870368428E-012
Poly[6]= 2.41324946244329E-011
Poly[5]= -7.90989844502363E-011
Interpolation 47
Poly[4] = 4.05827439744465E-009
Poly[3]= -3.31023555675263E-010
Poly[2]= -1.52308331145220E-004
Poly[l]= -2.53687934078865E-010
Poly[O]= 1.00000000007368E+000
The data is taken from a function of which the derivative could be computed
exactly. Though the actual values in the right-hand column are not displayed on
screen, they are shown here to indicate the accuracy of the routine.
Description
This example provides an interpolation algorithm. Given a set of data points (x,y),
the routine uses Newton's interpolary divided-difference equation to interpolate
between the points (Burden and Faires 1985,100-102). The data points must have
unique x-values, but these values need not be evenly spaced nor set in any particu-
lar order. You must supply the data points and the x-values at which interpolation is
to take place.
User-Defined Types
Input Parameters
Interpolation 49
Output Parameters
Sample Program
The sample program DIVDIF.PAS provides I/O functions that demonstrate New-
ton's interpolary divided-difference algorithm.
Input Files
Data may be entered from a text file. The x and y coordinates should be separated
by a space and followed by a carriage return. For example, data values of sqr(x}
could be entered in a text file as
11
24
39
416
525
Example
Problem. Interpolate the cosine function between x = Ix and x = 20x.
Run DIVDIF.PAS:
(K)eyboard or (F)ile entry of the data points? F
File name? SAMPLE3C.DAT
(K)eyboard or (F)ile entry of the interpolated points? K
Point 1: 1.5
Point 2: 2.5
Point 3: 3.5
Point 4: 4.5
Point 5: 5.5
Point 6: 6.5
Point 7: 7.5
Point 8: 8.5
Point 9: 9.5
Point 10: 10.5
Point 11: 11.5
Point 12: 12.5
Point 13: 13.5
Point 14: 14.5
Poi nt 15: 15.5
X y
12.000 0.9781476
8.000 0.9902681
1.000 0.9998477
10.000 0.9848078
5.000 0.9961947
15.000 0.9659258
4.000 0.9975641
3.000 0.9986295
7.000 0.9925462
14.000 0.9702957
The data is taken from a function of which the derivative could be computed
exactly. Though the values in the right-hand column are not displayed on screen,
they are shown here to indicate the accuracy of the routine.
Interpolation 51
Free Cubic Spline Interpolation (CUBEJRE.INC)
Description
This example constructs a smooth curve through a given set of data points. The
curve is a cubic spline interpolant with the following properties:
1. It passes through every data point.
2. It is continuous.
3. Its first derivative is continuous.
4. Its second derivative is continuous.
The second derivative is assumed to be zero at both endpoints (thus the cubic
spline is "free") of the interval determined by the data (Burden and Faires 1985,
117 fI). Cubics that join adjacent data points are of the following form:
S[i](x) = CoefO[i] + Coefl[i](x - x[i]) + Coef2[i](x - x[i]?
+ Coef3[i](x- X[i])3
where i ranges between 1 and the number of data points minus 1, the x[i]' s are the
x-coordinates of the input data, and x[i] :5 x < x[i + 1]. The interpolated values of
f(x) are found by evaluating the ith cubic polynomial at x, where
x[i] :5 x :5 x[i + 1].
Input Parameters
Output Parameters
Sample Program
Interpolation 53
Input Files
Data may be entered from a text file. The x and y coordinates should be separated
by a space and followed by a carriage return. For example, data values of sqr(x)
could be entered in a text file as
11
24
39
416
525
Example
Problem. Construct an interpolating spline for the following figure:
2 3 4 5 6
Because a cusp occurs at x = 3.55, we will construct two splines, one for each
side of the cusp.
Run CUBE_FRE.PAS:
(K)eyboard or (F)ile entry of the data points? F
File name? SAMPLE3D.DAT
(K)eyboard or (F)ile entry of the interpolated points? F
File name? SAMPLE3E.DAT
Data: x y
1: 0.0000000000 2.8000000000
2: 0.1000000000 2.7000000000
3: 0.2000000000 2.6000000000
4: 0.6000000000 2.2000000000
5: 1.0000000000 1.8000000000
6: 1.4000000000 1.6000000000
7: 1.8000000000 1.4000000000
8: 2.0000000000 1.4200000000
9: 2.2000000000 1.4000000000
10: 2.6000000000 1.5000000000
11: 3.0000000000 1.8000000000
12: 3.4000000000 2.4000000000
13: 3.4500000000 2.6000000000
14: 3.5000000000 2.8000000000
15: 3.5500000000 2.9000000000
Interpolated Points: X y
1: 0.3000000000 2.5018157855
2: 0.5000000000 2.3042222482
3: 1.2000000000 1. 6916808945
4: 1.6000000000 1.4759529845
5: 2.1000000000 1.4132967676
6: 2.3000000000 1.3989477848
7: 2.5000000000 1.4480232575
8: 2.7000000000 1.5697457729
9: 2.9000000000 1.7293593063
10: 3.2000000000 1.9502390938
11: 3.3000000000 2.1142270171
Interpolation 55
Second half of the figure:
(K)eyboard or (F)ile entry of the data points? F
Data: X y
1: 3.5500000000 2.9000000000
2: 3.6000000000 2.8000000000
3: 3.6500000000 2.6500000000
4: 3.8000000000 2.5000000000
5: 4.0000000000 2.3500000000
6: 4.3000000000 2.2000000000
7: 4.8000000000 1.9500000000
8: 5.3000000000 1.6000000000
9: 5.6000000000 1.3000000000
10: 5.8000000000 1.2000000000
11: 6.0000000000 0.0000000000
Splines: CoefO Coefl Coef2 Coef3
1: 2.9000000000 -1. 6719664279 0.0000000000 -131.2134288293
2: 2.8000000000 -2.6560671441 -19.6820143244 256.0671441466
3: 2.6500000000 -2.7037649955 18.7280572976 -49.1308266290
4: 2.5000000000 -0.4016786037 -3.3808146854 8.1960385189
5: 2.3500000000 -0.7704798556 1.5368084259 -2.1173630243
6: 2.2000000000 -0.4200828166 -0.3688182960 0.4179678583
7: 1.9500000000 -0.4754252188 0.2581334916 -1.4145661079
8: 1.6000000000 -1.2782163082 -1.8637156703 9.3036778805
9: 1.3000000000 0.1155473174 6.5095944222 -47.9366550462
10: 1.2000000000 -3.0330135193 -22.2523986055 37.0873310092
Interpolated Points: X y
1: 3.7000000000 2.5554905401
2: 3.9000000000 2.4342200313
3: 4.1000000000 2.2862027357
4: 4.2000000000 2.2404374617
5: 4.5000000000 2.1045744477
6: 4.6000000000 2.0520666406
7: 5.0000000000 1.8539237670
8: 5.2000000000 1.7105990402
9: 5.5000000000 1.3442375346
10: 5.7000000000 1.3287140209
11: 5.9000000000 0.7112619930
Description
Tl1is example constructs a smooth curve through a given set of data points. The
curve is a cubic spline interpolant with the following properties:
1. It passes through every data point.
2. It is continuous.
3. Its first derivative is continuous.
4. Its second derivative is continuous.
The first derivative at the endpoints of the interval determined by the input data
is defined by the user (Burden and Faires 1985, l22 ff.). (This is what makes the
cubic spline "clamped.") The cubics that join adjacent data points are of the follow-
ing form:
S[i](x) = CoefO[i] + Coefl[i](x - x[i]) + Coef2[i](x - X[i])2
+ Coef3[i](x - x[i]?
where i ranges between 1 and the number of data points minus 1, the x[i]' s are the
x-coordinates of the input data, and x[i] S x < x[i + 1]. The interpolated values of
f(x) are found by evaluating the ith cubic polynomial at x, where x[i] S x S
x[i + 1].
User-Defined Types
Input Parameters
Interpolation 57
NumInter: Integer; Number of interpolations
XInter:TNvector; X-coordinates of points at which to interpolate
The preceding parameters must satisfy the following conditions:
1. X data points must be unique.
2. X data points must be in ascending order.
3. NumPoints, Numlnter S TNArraySize.
4. NumPoints > 1.
Output Parameters
Input Files
Data may be entered from a text file. The x- and y -coordinates should be separated
by a space and followed by a carriage return. The last two values in the file must be
the derivatives of the function at the endpoints. For example, data values of sqr(x)
could be entered in a text file as
11
24
39
416
525
210
Note that the last two values are the derivatives of sqr(x) at the endpoints x =1
and x = 5.
Example
Problem. Construct an interpolating spline for the following figure:
2 3 4 5 6
Interpolation 59
Because a cusp occurs at x = 3.55, we will construct two splines, one for each
side of the cusp.
Run CUBE_CLA.PAS:
(K}eyboard or (F}ile entry of the data points? F
Data: x y
1: 0.0000000000 2.8000000000
2: 0.1000000000 2.7000000000
3: 0.2000000000 2.6000000000
4: 0.6000000000 2.2000000000
5: 1.0000000000 1.8000000000
6: 1.4000000000 1.6000000000
7: 1.8000000000 1.4000000000
8: 2.0000000000 1.4200000000
9: 2.2000000000 1.4000000000
10: 2.6000000000 1.5000000000
11: 3.0000000000 1.8000000000
12: 3.4000000000 2.4000000000
13: 3.4500000000 2.6000000000
14: 3.5000000000 2.8000000000
15: 3.5500000000 2.9000000000
Data: X y
1: 3.5500000000 2.9000000000
2: 3.6000000000 2.8000000000
3: 3.6500000000 2.6500000000
4: 3.8000000000 2.5000000000
5: 4.0000000000 2.3500000000
6: 4.3000000000 2.2000000000
7: 4.8000000000 1.9500000000
8: 5.3000000000 1.6000000000
9: 5.6000000000 1.3000000000
10: 5.8000000000 1.2000000000
11: 6.0000000000 0.0000000000
Interpolation 61
Interpolated Points: X y
1: 3.7000000000 2.5490351248
2: 3.9000000000 2.4368630843
3: 4.1000000000 2.2844619846
4: 4.2000000000 2.2388141319
5: 4.5000000000 2.1097537107
6: 4.6000000000 2.0584802174
7: 5.0000000000 1.8354671712
8: 5.2000000000 1.6919117155
9: 5.5000000000 1.3872367766
10: 5.7000000000 1.2432752125
11: 5.9000000000 1.0138449575
x f(x)
1.0 45.0
1.1 49.2
1.2 54.5
1.3 59.8
1.4 65.1
1.5 70.4
The units are in hours and miles, and the data refers to a trip that started at
noon.f(l.O) = 45.0, so the distance traveled by one o'clock is 45.0 miles, andf(l.5)
= 70.4, so by half past one you will be 70.4 miles from where you were at noon.
The derivative of this distance function gives the velocity function. The car's
velocity at one o'clock is the value of the derivative at x = l.0. From the previous
data, it is impossible to compute the derivative exactly, but it is possible to approxi-
mate the derivative. The car traveled 49.2 - 45.0 = 4.2 miles in the six minutes
after one o'clock (1.1 - 1.0 = 0.1 hours = 6 minutes). Thus, the average velocity
of the car during thos~ six minutes is 4.2 / 0.1 = 42 miles per hour. This gives an
approximation to the velocity at one o'clock.
63
Each method described in this chapter approximates derivatives of a real func-
tion of one real variable.
The routines DERIV.lNC, DERIV2.1NC, and INTERDRV.lNC compute
derivatives of a function that is represented by tabular data. Consequently, their
accuracy depends heavily upon the precision and spacing of the data points.
The routines DERIVFN.lNC and DERIV2FN.lNC compute derivatives of a
user-defined function. Consequently, the accuracy of the values calculated with
these routines is limited by the precision of the computer.
Differentiation consists of subtracting two very close numbers and dividing by a
very small number; hence, it is extremely sensitive to round-off error. The accuracy
of the first derivative is approximately the square root of the precision with which
real numbers are represented; the accuracy of the second derivative is approxi-
mately equal to the fourth root. Thus, the precision of the first derivative will be
about IE - 8 when run with the 8087 math coprocessor, or about IE - 4 when run
without the 8087 math coprocessor. The precision of the second derivative will be
about IE - 4 with the 8087, or IE - 2 without it.
The first derivative of a function that is represented by a table of values can be
approximated in DERIV.lNC via a two-point formula, a three-point formula, or a
five-point formula. The accuracy of the formula increases with the number of
points used in the formula. In order to use the five-point formula, however, the
domain values of the data points (that is, the x-coordinates) must be equally spaced.
This is not required for the two-point and three-point formulas. Derivatives can
only be approximated at data points.
The second derivative of a function that is represented by a table of values can
be approximated in DERIV2.1NC via a three-point formula or a five-point for-
mula. The domain values of the data points must be equally spaced (regardless of
whether the three-point formula or five-point formula is used). Second derivatives
can only be approximated at data points.
The routine INTERDRV.lNC approximates a function by constructing a free
cubic spline to a set of data points. Cubic splines avoid the undesirable oscillatory
behavior of other interpolating polynomials. The derivative of the cubic spline at a
given domain value, which may be different from the input data values, will then
approximate the corresponding derivative of the function.
The first derivative of a user-supplied function is approximated in DERIV
FN.lNC via a three-point formula. The approximation is refined with Richardson
extrapolation. The derivative can be approximated at any point within the domain
of the function.
Numerical Differentiation 65
First Differentiation Using Two-Point, Three-Point, or
Five-Point Formulas (DERIV.INC)
Description
User-Defined Types
Input Parameters
Output Parameters
Comments
Numerical Differentiation 67
Since numerical differentiation is prone to round-off errors, TNNearlyZero is
different in this routine. The values of TNNearlyZero are TNNearlyZero = IE - 13
if using the 8087 math coprocessor and TNNearlyZero = IE - 6 if not using the
8087.
Sample Program
The sample program DERIV.PAS provides I/O functions that demonstrate differ-
entiation with two-point, three-point, and five-point formulas.
Input Files
Data points may be entered from a text file. The x- and y- coordinates should be
separated by a space and followed by a carriage return. For example, data values of
sqr(x) could be entered in a text file as
11
24
39
416
525
Derivative points may also be entered from a text file. Every derivative point
must be followed by a carriage return. For example, to determine the derivatives of
the preceding points, create the following file of derivative points:
1
2
3
4
5
Example
Problem. Approximate the first derivative off(x) = sqr(x) * cos (x) at several points
between one and two radians. The output from three runs is given. Actual values of
the derivatives to eight significant figures are also given.
Input Data: X y
1.0000000 5.40302305868140E-001
1.1000000 5.48851306924949E-001
1.2000000 5.21795166446410E-001
1.3000000 4.52073020375553E-001
1.4000000 3.33135600084472E-001
1.5000000 1.59158703752332E-001
1.6000000 -7.47507770912994E-002
1.7000000 -3.72360588514066E-001
1.8000000 -7.36134786805602E-001
1.9000000 -1.16707533637725E+000
2.0000000 -1.66458734618857E+000
Output using two-point differentiation:
Numerical Differentiation 69
Output using three-point differentiation:
Description
This example contains two algorithms that approximate the second derivative of a
functionf(x) when several data points (x,f(x)) are specified. You decide whether to
use a three-point or five-point formula (Gerald and Wheatley 1984, 236-237);
three points are used in the three-point formula, and five in the five-point formula.
You must supply the data points (x,j(x)) and the x-values of the data points at which
the second derivative is to be approximated. The second derivative may only be
approximated at x-values that were input as data points.
User-Defined Types
Input Parameters
Numerical Differentiation 71
5. XData points must be equally spaced.
6. XDeriv points must be a subset of the XData points.
7. NumPoints, NumDeriv :5 TNArraySize.
Output Parameters
Comments
If an x-value at which the second derivative is approximated is not among the data
points, the value - 9.9999999E35 is arbitrarily assigned to the derivative at that
point and Error = 1 is returned. When using five-point second differentiation with
only five data points, there is insufficient information for approximating the second
derivative at the second and fourth data points. Should an attempt be made to
approximate the second derivative at these points, the value 9.9999999E35 is arbi-
trarily assigned to the second derivative at that point and Error = 1 is returned.
Sample Program
Input Files
Data points may be entered from a text file. The x- and y-coordinates should be
separated by a space and followed by a carriage return. For example, data values of
sqr(x) could be entered in a text file as
11
24
39
416
525
Derivative points may also be entered from a text file. Every derivative point
must be followed by a carriage return. For example, to determine the second deriv-
atives of the preceding points, create the following file of derivative points:
1
2
3
4
5
Example
Problem. Approximate the second derivative of f(x) = sqr(x) * cos(x) at several
points between x = 1 and x = 2 radians. The output from two runs is given. Actual
values of the second derivatives to eight significant figures are also given.
Numerical Differentiation 73
Run DERIV2.PAS:
(K)eyboard or (F)ile entry of the data points? F
Point 1: 1.1
Point 2: 1.3
Point 3: 1.5
Point 4: 2.0
Point 5: 2.2
Input Data: X y
1.0000000 5.40302305868140E-001
1.1000000 5.48851306924949E-001
1.2000000 5.21795166446410E-001
1.3000000 4.52073020375553E-001
1.4000000 3.33135600084472E-001
1.5000000 1.59158703752332E-001
1.6000000 -7.47507770912994E-002
1.7000000 -3.72360588514066E-001
1.8000000 -7.36134786805602E-001
1.9000000 -1. 16707533637725E+000
2.0000000 -1.66458734618857E+000
Output using three-point second differentiation:
Numerical Differentiation 75
Differentiation with a Cubic Spline Interpolant
(INTERDRV.INC)
Description
This example contains an algorithm for approximating the first and second deriva-
tives of a function given several data points (x,f(x)). The algorithm assumes that a
free cubic spline interpolant (Burden and Faires 1985, 117-122) is an adequate ap-
proximation to the functionf(x), so that the slope of the interpolant at any value x.
is an adequate approximation tof'(x). See Chapter 3 (CUBE_FRE.lNC) for mor~
information on free cubic splines. The user must supply the data points (x,f(x)) and
the x-values at which to approximate the derivatives. Derivatives may be approxi-
mated at any x-value contained in the closed interval determined by the data
points. This routine will likely give significant errors if interpolation (Gerald and
Wheatley 1984, 227-231) is attempted outside the range of x-values (extrapo-
lation).
User-Defined Types
Input Parameters
Output Parameters
Sample Program
Numerical Differentiation 77
Input Files
Data points may be entered from a text file. The x- and y-coordinates should be
separated by a space and followed by a carriage return. For example, data values of
sqr(x) could be entered in a text file as
11
24
39
416
525
Derivative points may also be entered from a text file. Every derivative point
must be followed by a carriage return. For example, to determine the derivatives of
the preceding points, create the following file of derivative points:
1
2
3
4
5
Example
Problem. Determine the first and second derivative ofJ(x) = sqr(x) * cos (x) at
several points between one and two radians. Actual values of the derivatives to
eight significant figures are given here.
Run INTERDRY.PAS:
(K)eyboard or (F)ile entry of data points? F
File name? SAMPLE4B.DAT
(K)eyboard or (F)ile entry of derivative points? K
Number of derivative points (0-100)?5
Point 1: 1.1
Poi nt 2: 1. 3
Point 3: 1.55
Point 4: 1.95
Point 5: 2.20
The data is taken from a function of which the derivative could be computed
exactly. The actual values are shown here:
Numerical Differentiation 79
Differentiation of a User-Defined Function
(DERIVFN.INC)
Description
Given a user-defined function fix), this example will approximate the first deriva-
tive of the function at a set of x values. The formula
f' (x) = [f(x + ~X) - f(x - ~X)]/2*~
User-Defined Function
Input Parameters
Output Parameters
Comments
Sample Program
The sample program DERIVFN.PAS provides I/O functions that find the first
derivative of a function at a set of points.
Input Files
Derivative points may be entered from a text file. Every derivative point must be
followed by a carriage return. For example, to determine the derivatives at x-values
1 through 5, create the following file of derivative points:
1
2
3
4
5
Numerical Differentiation 81
Example
Problem. Determine the first derivative off(x) = sqr(x) * cos(x) at several points
between 1 and 2.2. Actual values of the derivatives to eight significant figures are
given here.
First, write the function into the DERIVFN.PAS program:
(* ----- here is the function to differentiate -------------------- *)
funct;on TNTargetF(X : Real) : Real;
beg;n
TNTargetF := Sqr(X)*Cos(X);
end; { function TNTargetF }
(* ---------------------------------------------------------------- *)
Run DERIVFN.PAS:
(K)eyboard or (F)ile entry of derivative points? K
Number of points (0-100)? 5
Poi nt 1: 1.1
Point 2: 1.3
Point 3: 1.55
Point 4: 1.95
Point 5: 2.2
Tolerance = 1.00000000000000E-004
Actual Values
X Derivative at X X Value at X 1st Deriv at X
1.100 -8.04494385380506E-002 1.1 0.5488513 -0.0804494
1.300 -9.32916380187812E-001 1.3 0.4520730 -0.9329164
1.550 -2.33751652942968E+000 1.55 0.0499596 -2.3375165
1.950 -4.97607456093019E+000 1.95 -1.4076126 -4.9760746
2.200 -6.50252751001358E+000 2.20 -2.8483454 -6.5025275
The data is taken from a function of which the derivative could be calculated
exactly. Though the values in the three right-hand columns (under "Actual Values")
are not displayed on screen, they are shown here to indicate the accuracy of the
routine.
Description
Given a user-defined function fix), this example will approximate the second deriv-
ative of the function at a set of x values. The three-point formula
f" (x) = [f(x + /lX) - 2f(x) + f(x - AX)]/AX2
gives a first approximation to the second derivative. Richardson extrapolation is
then used to refine the approximation (Burden and Faires 1985, 142-152).
User-Defined Types
Input Parameters
Numerical Differentiation 83
The preceding parameters must satisfy the following conditions:
1. NumDeriv $; TNArraySize
2. Tolerance ~ TNNearlyZero
Output Parameters
Comments
Sample Program
The sample program DERIV2FN.PAS provides I/O functions that find the second
derivative of a function at a set of points.
Example
Problem. Determine the second derivative of fix) = sqr(x) * cos (x) at several
points between 1 and 2.2. Actual values of the derivatives to eight significant fig-
ures are given here.
First, write the function into the DERIV2FN .PAS program:
(* ----- here is the function to differentiate -------------------- *)
(* ---------------------------------------------------------------- *)
Run DERIV2FN.PAS:
(K)eyboard or (F)ile entry of derivative pOints? K
Number of points (0-100)? 5
Point 1: 1.1
Point 2: 1.3
Point 3: 1.55
Point 4: 1.95
Point 5: 2.2
Tolerance (> 0, default = 1.000E-02)? 1E-4
Direct output to one of the following:
(S)creen
(P)rinter
(F)ile
Tolerance = 1.00000000000000E-004
Numerical Differentiation 85
Actual Values
X 2nd Derivative at X X Value at X 2nd Deriv at X
1.100 -3.56297143915941E+000 1.1 0.5488513 -3.5629715
1.300 -4.92757787674466E+000 1.3 0.4520730 -4.9275779
1.550 -6.20702925534123E+000 1.55 0.0499596 -6.2070293
1.950 -6.57863484485542E+000 1.95 -1.4076126 -6.5786348
2.200 -5.44342524529641E+000 2.20 -2.8483454 -5.4434252
The data is taken from a function of which the derivative could be calculated
exactly. Though the values in the three right-hand columns (under "Actual Values")
are not displayed on screen, they are shown here to indicate the accuracy of the
routine.
87
integrand. Simpson's method (ADAPSIMP.INC) and the Gaussian quadrature
method (ADAPGAUS.INC) are used with adaptive schemes. The Gaussian quad-
rature method permits, in some instances, the integrand to possess a singularity at
an endpoint of integration, since the function is evaluated at points that are not the
endpoints of the interval of integration.
The Romberg method (ROMBERG.INC) uses the trapezoid method and Rich-
ardson extrapolation to approximate the integral. It returns an approximation
within a user-specified accuracy. Except for extremely oscillatory functions or
functions that possess an endpoint singularity, this method is fastest and most
accurate. If the function oscillates substantially or possesses an endpoint singular-
ity, the adaptive Gaussian quadrature routine is preferred.
Description
This example uses Simpson's composite algorithm (Burden and Faires 1985, 156-
167) to approximate the definite integral of a functionJ\x) over an interval [a, b].
The interval is divided into N subintervals of equal length. The curve in each
subinterval is approximated by a second-degree Lagrange polynomial. The integral
of the resulting polynomial is then calculated. The sum of the integrals of the N
Lagrange polynomials approximates the integral of the function! over the interval
[a, bJ. You must supply the function, the limits of integration, and the number of
subintervals.
Input Parameters
Numerical Integration 89
Output Parameters
Sample Program
Example
Problem. Approximate the integral exp(3x) + sqr(x)/3 from 0 to 5 using Simpson's
composite algorithm.
1. Code function TNTargetF:
funct;on TNTargetF(x : Real) : Real;
(**************************************************************************)
(**** THIS IS THE FUNCTION TO INTEGRATE ****)
(**************************************************************************)
beg;n
TNTargetF .- Exp(3*X) + Sqr(X)/3;
end; { function TNTargetF }
Numerical Integration 91
Integration Using the Trapezoid Composite Rule
(TRAPZOID.INC)
Description
This example uses the trapezoid composite rule (Burden and Faires 1985,154-167)
to approximate the definite integral of a function fix) over an interval [a, b]. The
interval is divided into N subintervals of equal length. In each subinterval the
function is approximated by a straight line. The sum of the integrals of the result-
ing trapezoids approximates the integral of the function f over the interval [a, b].
You must supply the function, the limits of integration, and the number of subinter-
vals.
User-Defined Function
Input Parameters
Sample Program
The sample program T~APZOID.PAS provides I/O functions that demonstrate the
trapezoid composite rule.
Example
Problem. Approximate the integral exp(3x) + sqr(x)/3 from 0 to 5 using the trape-
zoid composite rule.
1. Code function TNTargetF:
function TNTargetF(x : Real) : Real;
(*****************************************************************************)
(**** THIS IS THE FUNCTION TO INTEGRATE ****)
(*****************************************************************************)
begin
TNTargetF .- Exp(3*X) + Sqr(X)/3;
end; { function TNTargetF }
Numerical Integration 93
2. Run TRAPZOID.PAS:
Lower limit of integration? 0
Upper limit of integration? 5
Number of intervals (> O)? 100
Direct output to one of the following:
(S)creen
(P)rinter
(F) il e
Description
Input Parameters
Numerical Integration 95
Output Parameters
Comments
Sample Program
The sample program ADAPSIMP.PAS provides I/O functions that demonstrate the
adaptive quadrature method with Simpson's rule.
Numerical Integration 97
Integration Using Adaptive Quadrature and Gaussian
Quadrature (ADAPGAUS.INC)
Description
Input Parameters
Comments
Numerical Integration 99
The following condition is satisfied by the numbers that follow it:
Integral from -1 to 1 ofJ(x) dx
equals
Sum from i = 1 to NumLegendreTerms of
Legendre[i].Weight * J(Legendre[i].Root)
for an arbitrary functionJ(x).
Legendre[I] ........................ ............... Root: 0.0950125098376370440185
Weight: 0.189450610455068496285
Legendre[2] ...................... ................. Root: 0.281603550778258913230
Weight: 0.182603415044923588867
Legendre[3] ... .................................... Root: 0.458016777657227386342
Weight: 0.169156519395002538189
Legendre[4] ................. ...................... Root: 0.617876244402643748447
Weight: 0.149595988816576732081
Legendre[5] ....................................... Root: 0.755404408355003Q33895
Weight: 0.124628971255533872052
Legendre[6] ....................................... Root: 0.865631202387831743880
Weight: 0.095158511682492784810
Legendre[7] ....................................... Root: 0.944575023073232576078
Weight: 0.062253523938647892863
Legendre[8] ... .................................... Root: 0.989400934991649932596
Weight: 0.027152459411754094852
Legendre[9] ....................................... Root: - 0.0950125098376370440185
Weight: 0.189450610455068496285
Legendre[10] ..................................... Root: - 0.281603550778258913230
Weight: 0.182603415044923588867
Legendre[ll] ..................................... Root: - 0.458016777657227386342
Weight: 0.169156519395002538189
Legendre[12] ..................................... Root: - 0.617876244402643748447
Weight: 0.149595988816576732081
Legendre[13] ..................................... Root: - 0.755404408355003033895
Weight: 0.124628971255533872052
Legendre[14] ..... ; ............................... Root: - 0.865631202387831743880
Weight: 0.095158511682492784810
Legendre[15] ..................................... Root: - 0.944575023073232576078
Weight: 0.062253523938647892863
Legendre[16] ..................................... Root: - 0.989400934991649932596
Weight: 0.027152459411754094852
Example
Problem. Approximate the integral exp(3x) + sqr(x)/3 from
quadrature with Gaussian quadrature algorithm.
°to 5 using adaptive
1. Code function TNTargetF:
funct;on TNTargetF(x : Real) : Real;
(*****************************************************************************)
(** THIS IS THE FUNCTION TO INTEGRATE ***)
(*****************************************************************************)
beg;n
TNTargetF := Exp(3*X) + Sqr(X)/3;
end; { function TNTargetF }
2. Run ADAPGAUS.PAS:
Lower limit of integration? 0
Upper limit of integration? 5
Tolerance (> 0, default = 1.000E-08): 1E-8
Maximum number of subintervals (> 0, default 1000): 1000
Direct output to one of the following:
(S)creen
(P)rinter
(F) il e
Integral: 1.08968601304609E+006
To eight significant figures, the correct answer is 1,089,686.0.
Description
This example contains an algorithm (Burden and Faires 1985,177-182) for approxi-
mating the integral of a functionf(x) over an interval [a, b] within a specified toler-
ance. The trapezoid rule is used to generate a preliminary approximation, and
Richardson extrapolation (Burden and Faires 1985,148-152) is subsequently used
to improve the approximation. Extrapolation continues until the fractional differ-
ence between successive approximations of the integral is less than the tolerance.
You must supply the function, the limits of integration, and the tolerance with
which to approximate the integral.
Input Parameters
Sample Program
The sample program ROMBERG.PAS provides I/O functions that demonstrate the
Romberg algorithm.
Example
Problem. Approximate the integral exp(3x} + sqr(x}/3 from 0 to 5 using the Rom-
berg algorithm.
1. Code function TNTargetF:
function TNT~rgetF(x : Real) : Real;
(*************************************************************************)
(**** THIS IS THE FUNCTION TO INTEGRATE ****)
(*************************************************************************)
begin
TNTargetF := Exp(3*X) + Sqr(X)/3;
end; { function TNTargetF }
This chapter provides routines for dealing with systems of linear equations. An
example of a system of linear equations is as follows:
2X+Y+Z=7
X-Y+Z=2
X+Y-Z=O
Matrix algebra is a collection of notations that constitutes a technique for han-
dling such systems. With matrix algebra, the preceding system would be written
Ax = b
where
[~
1
A= -1
1
105
cannot be solved. Even for different values of the right-hand side, the equations
can only be solved in certain exceptional cases. (If you change 4 and 5 to 3 and - 2,
then there are infinitely many solutions; but there are none if you change 4 and 5 to
3 and - 3.0001.)
Following is a description of several routines that operate on matrices and sys-
tems of linear equations.
The determinant of a square matrix is found via DET.INC.
The inverse of a nonsingular matrix is found via INVERSE.INC.
The direct techniques implemented to solve a system of N linear equations in N
unknowns are Gaussian elimination (GAUSELIM.INC), Gaussian elimination
with partial pivoting (PARTPIVT.INC), and direct factorization (DIRFACT.INC) ..
The Gauss-Seidel method (GAUSSIDL.INC), an iterative technique that con-
verges to the solution, is seldom used for solving small systems, since the time
required for sufficient accuracy exceeds that required for the preceding direct
techniques.
In general, Gaussian elimination with partial pivoting is the fastest, most accu-
rate algorithm (see Chapter 9, LEAST.INC, for an application of PARTPIVT.INC).
The following special cases may warrant the use of one of the other routines:
• If you are considering systems where round-off is minimal (that is, small systems
whose coefficients are all of nearly the same magnitude), Gaussian elimination
without pivoting may be used. It is somewhat faster than its pivoting counterpart
(PARTPIVT.INC).
• When considering sparse coefficient matrices, the Gaussian elimination routine
with partial pivoting is the most efficient and accurate routine. If the matrix is
small and the nonzero coefficients do not differ wildly from each other, regular
Gaussian elimination (GAUSELIM.INC) can usually be used safely.
• For large, dense matrices, the iterative technique (GAUSSIDL.INC) is the most
efficient; it creates less round-off error than the direct methods. However, the
Gauss-Seidel algorithm has' its own weaknesses (see the section, "Solving a Sys-
tem of Linear Equations with the Iterative Gauss-Seidel Method," for more
details).
• When it is necessary to solve several systems with the same coefficient matrix
but a different vector of· constant terms, the direct factorization method
(DIRFACT.INC) is the most efficient. This is because it does not require reduc-
tion of the coefficient matrix for each vector of constants. (See Chapter 7 for an
application of DIRFACT.INC.)
Description
Input Parameters
Sample Program
The sample program DET.PAS provides I/O functions that demonstrate how to
find the determinant of a matrix.
Input File
Data may be input from a text file. All entries in the text file should be separated
by a space or carriage return, and it does not matter if the text file ends with a
carriage return. The format of the text file should be like this:
1. The dimension of the matrix
2. The elements of the matrix in row order; that is,
[1, 1], [1, 2] ... [1, N], [2, 1] ... [2, N] ... [N, N],
where N is the dimension of the matrix
-i ! ~ =~:~ 1
[o 2 2
0
1
3
-3.0
-4.0
Run DET.PAS:
(K)eyboard or (F)ile input of data? F
The matrix:
1.00000000 2.00000000 0.00000000 -1.00000000
-1.00000000 4.00000000 3.00000000 -0.50000000
2.00000000 2.00000000 1.00000000 -3.00000000
0.00000000 0.00000000 3.00000000 -4.00000000
Determinant = -2.10000000000000E+001
Description
Input Parameters
Sample Program
The sample program INVERSE.PAS provides I/O functions that demonstrate how
to find the inverse of a matrix.
Input Files
Data may be input from a text file. All entries in the text file should be separated
by a space or carriage return, and it does not matter if the text file ends with a
carriage return. The format of the text file should be as follows:
1. The dimension of the matrix
2. The elements of the matrix in row order; that is,
[1, 1], [1, 2] ... [1, N], [2, 1] ... [2, N] ... [N, N],
where N is the dimension of the matrix
The matrix:
1.000000000 2.000000000 0.000000000 -1.000000000
-1.000000000 4.000000000 3.000000000 -0.500000000
2.000000000 2.000000000 1.000000000 -3.000000000
0.000000000 0.000000000 3.000000000 -4.000000000
Inverse:
-1.952380952 0.190476190 1.571428571 -0.714285714
0.761904762 0.047619048 -0.357142857 0.071428571
-1.904761905 0.380952381 1.142857143 -0.428571429
-1.428571429 0.285714286 0.857142857 -0.571428571
The matrix:
-1.952380952 0.190476190 1.571428571 -0.714285714
0.761904762 0.047619048 -0.357142857 0.071428571
-1.904761905 0.380952381 1.142857143 -0.428571429
-1.428571429 0.285714286 0.857142857 -0.571428571
Inverse:
1.000000000 2.000000000 0.000000000 -1.000000000
-1.000000000 4.000000000 3.000000000 -0.500000000
2.000000000 2.000000000 1.000000000 -3.000000000
-0.000000000 -0.000000000 3.000000000 -4.000000000
Description
Input Parameters
Output Parameters
Sample Program
Input File
Data may be input from a text file. All entries in the text file should be separated
by a space or carriage return, and it does not matter if the text file ends with a
carriage return. The format of the text file should be as follows:
1. The dimension of the coefficient matrix
2. The elements of the matrix in row order; that is,
[1, 1], [1, 2], ... , [1, N], [2, 1], "', [2, N], ... , [N, N],
where N is the dimension of the matrix
3. The elements of the constant vector, in the order [l], ... ,[N]
Example
Problem. Solve the following linear system:
w + 2x + Oy - z = 10.0
- w + 4x + 3y - 0.5z = 21.5
2w + 2x + y - 3z = 26.0
3y - 4z = 37.0
Run GAUSELIM.PAS:
(K)eyboard or (F)ile input of data? F
File name? 5AMPLE6A.DAT
Direct output to one of the following:
(5) creen
(P)rinter
(F) il e
The coefficients:
1.000000000 2.000000000 0.000000000 -1.000000000
-1.000000000 4.000000000 3.000000000 -0.500000000
2.000000000 2.000000000 1.000000000 -3.000000000
0.000000000 0.000000000 3.000000000 -4.000000000
The constants:
1.00000000000000E+001
2. 15000000000000E+001
2.60000000000000E+001
3.70000000000000E+001
The solution:
-1.00000000000000E+000
2.00000000000000E+000
3.00000000000000E+000
-7.00000000000000E+000
Description
User-Defined Types
Input Parameters
Output Parameters
Sample Program
The sample program PARTPIVT.PAS provides I/O functions that demonstrate how
to solve a system of linear equation with Gaussian elimination and partial pivoting.
Input File
Data may be input from a text file. All entries in the text file should be separated
by a space or carriage return, and it does not matter if the text file ends with a
carriage return. The format of the text file should be as follows:
1. The dimension of the matrix
2. The elements of the matrix in row order; that is,
[1, 1], [1, 2], ... , [1, N], [2, 1], ... , [2, N], ... , [N, N],
where N is the dimension of the matrix
3. The elements of the constant vector, in the order [l], ... ,[N]
Example
Problem. Solve the following linear system:
w + 2x + Oy - z = 10
-w + 4x + 3y - 0.5z = 21.5
2w + 2x + Y - 3z = 26
3y - 4z = 37
Run 'PARTPIVf.PAS:
{K)eyboard or {F)ile input of data? F
File name? SAMPLE6A.DAT
Direct output to one of the following:
{S)creen
{P)rinter
(F) i 1e
The coefficients:
1.000000000 2.000000000 0.000000000 -1.000000000
-1.000000000 4.000000000 3.000000000 -0.500000000
2.000000000 2.000000000 1.000000000 -3.000000000
0.000000000 0.000000000 3.000000000 -4.000000000
The constants:
1.00000000000000Et001
2. 15000000000000Et001
2.60000000000000Et001
3.70000000000000Et001
The solution:
-1.00000000000000EtOOO
2.00000000000000EtOOO
3.00000000000000EtOOO
-7.00000000000000EtOOO
Description
User-Defined Types
Sample Program
The sample program DIRFACT.PAS provides I/O functions that demonstrate how
to solve a system of linear equations with the method of direct factoring.
Input File
Data may be input from a text file. All entries in the text file should be separated
by a space or carriage return, and it does not matter if the text file ends with a
carriage return. The format of the text file should be as follows:
l. The dimension of the matrix
2. The elements of the matrix in row order; that is,
[1, 1], [1, 2], ... , [1, N], [2, 1], ... , [2, N], ... , [N, N],
where N is the dimension of the matrix
3. The elements of the first constant vector, in the order [1], ... ,[N], with each
element followed by a carriage return
2x + 3y = 10 2x + 3y = 1
-4x = 10 -4x = 2
a text file could be created to look like this:
2
2 3
-4 0
10
10
1
2
Example
Problem. Given the following set of coefficients:
2w + x + 5y - 8z
7w + 6x + 2y + 2z
- 1w - 3x - lOy + 4z
2w + 2w + 2y + z
compute solutions for each of the five constant vectors:
-15 14 -13
[ 17
0
-1~
50
-5 -12
1 84
-51
30
-15
5]
17 1 37 10
Run DIRFACT.PAS:
(K)eyboard or (F)ile input of data? F
File name? SAMPLE6C.DAT
Direct output to one of the following:
(S)creen.
(P)rinter
(F) il e
The constants:
O.OOOOOOOOOOOOOOE+OOO
1.70000000000000E+001
-1.00000000000000E+001
7.00000000000000E+000
The solution:
9.99999999999999E-001
1.00000000000000E+000
9.99999999999999E-001
9.99999999999999E-001
The constants:
-1.50000000000000E+001
5.00000000000000E+001
-5.00000000000000E+000
1.70000000000000E+001
The solution:
2.00000000000000E+000
4.99999999999999E+000
1.85358974546113E-015
3.00000000000000E+000
The constants:
1.40000000000000E+001
1.00000000000000E+000
-1.20000000000000E+001
1.00000000000000E+000
The solution:
1.00000000000000E+000
-1.00000000000000E+000
1.00000000000000E+000
-1.00000000000000E+000
The constants:
-1.30000000000000E+001
8.40000000000000E+001
-5.10000000000000E+001
3.70000000000000E+001
The solution:
3.99999999999999E+000
5.00000000000001E+000
6.00000000000000E+000
7.00000000000000E+000
The solution:
-1.01506105108586E-015
5.00000000000000E+OOO
O.OOOOOOOOOOOOOOE+OOO
O.OOOOOOOOOOOOOOE+OOO
Description
I
j=I
A[i,j] X",[j] - I
j=i+l
(A[i,j] Xm_Jj]) + B[i]
X,Ji] =
A[i,i]
The algorithm halts when the fractional difference for each element of the vector
X between two iterations is less than a specified tolerance.
If A is diagonally dominant (that is, each of the diagonal terms is greater than or
equal to the sum of the off-diagonal terms in the same row), then the sequence will
converge to the solution X. If the matrix A is not diagonally dominant, then the
sequence may converge to the solution, but more likely it will not. You must supply
the tolerance with which to approximate a solution.
Output Parameters
Sample Program
Input File
Data may be input from a text file. All entries in the text file should be separated
by a space or carriage return, and it does not matter if the text file ends with a
carriage return. The format of the text file should be as follows:
l. The dimension of the matrix
2. The elements of the matrix in row order; that is,
[1, 1], [1, 2], "', [1, N], [2, 1], ... , [2, N], ... , [N, N],
where N is the dimension of the matrix
3. The elements of the first constant vector, in the order [l], ... ,[N]
Example
Problem. Solve the following linear system to within a tolerance of IE - 12:
10v + w + 2x - 3y + 2z = 29
4v + SOw + x + z = 35
- 2v + 5w - 30x + Y + z = 25
6v + 4w + lOy + 3z = 46
- 3v - 2w - x + 6y + 25z = - 106
Run GAUSSIDL.PAS:
(K)eyboard or (F)ile input of data? F
The coefficients:
10.000000000 1.000000000 2.000000000 -3.000000000 2.000000000
4.000000000 50.000000000 1.000000000 0.000000000 1.000000000
-2.000000000 5.000000000 -30.000000000 1.000000000 1.000000000
6.000000000 4.000000000 0.000000000 10.000000000 3.000000000
-3.000000000 -2.000000000 -1.000000000 6.000000000 25.000000000
Number of iterations: 15
The result:
-2.99999999999997E+000
9.99999999999999E-001
9.99999999999998E-001
-1.99999999999999E+000
-4.00000000000000E+000
The routines in this chapter can find the eigenvalues and eigenvectors. A scalar c is
an eigenvalue (or characteristic value) of a square matrix A if there is a nonzero
vector v satisfying
Av = cv
The vector v is called the eigenvector corresponding to c.
The eigenvalues and eigenvectors of a matrix provide a lot of information about
the matrix. If a matrix is written in terms of a basis of eigenvectors, then it is
diagonal, meaning. that its only nonzero terms are on the main diagonal.
Each procedure in this chapter attempts to approximate at least' one real eigen-
value (and associated eigenvector) of a real square matrix. The eigenvector is nor-
malized so that the element with the largest magnitude is 1.
The power method (POWER.INC) approximates the eigenvalue that is largest
in magnitude (dominant eigenvalue). The iterative process will converge slowly or
not at all if the dominant eigenvalue is not simple or if it has nearly the same
magnitude as the next most-dominant eigenvalue.
The inverse power method (INVPOWER.INC) approximates the eigenvalue
nearest to a user-supplied real value. This process usually converges more rapidly
than the power method, and may be used to refine the approximate value of the
eigenvalue determined by the latter method (POWER.INC).
131
The Wielandt method (WIELANDT.lNC) attempts to approximate a user-
specified number of eigenvalues of a given matrix. The power method (POWER.
INC) is first used to approximate the dominant eigenvalue of the matrix. Deflation
is employed to form a deflated, square matrix (that is, a square matrix whose
dimension is one less than the original matrix). The eigenvalues of the deflated
matrix are identical to those of the original matrix except for the determined domi-
nant eigenvalue. Eigenvectors of the remaining eigenvalues from the original
matrix are also contained in the deflated matrix. The dominant eigenvalue of the
new deflated" matrix is then determined using the power method. Wielandt's
method is susceptible to round-off error, thus it may be desirable to use its results
as input to the inverse power method (INVPOWER.lNC).
The cyclic Jacobi method QACOBI.INC) approximates all the eigenvalues of a
symmetric matrix. The iterative process uses orthogonal plane rotations to reduce
the given matrix into a diagonal form. Although Jacobi's method is only applicable
to symmetric matrices, it is much more efficient and accurate than Wielandt's
method.
Description
The power method (Burden and Faires 1985, 452-456) approximates the dominant
real eigenvalue of a matrix and its associated eigenvector. The dominant eigen-
value is the eigenvalue of the largest absolute magnitude. Given a square matrix A
and a real nonzero vector v, a vector w is constructed by the matrix operation
Av = w. The vector w is normalized by dividing by its element of the largest
absolute magnitude q. If the absolute difference between each of the correspond-
ing elements in wand v is less than a specified tolerance, then the procedure halts.
Otherwise, v is set equal to w, and the operation repeats until a solution is found.
The magnitude q is the dominant eigenvalue, and w will be the associated eigen-
vector of the matrix A.
You must supply the matrix A, an initial approximation to the eigenvector v, and
the tolerance.
User-Defined Types
Input Parameters
Output Parameters
Comments
The power method will not converge if the initial approximation (Guess) to the
eigenvector is orthogonal to the dominant eigenvector. If the initial approximation
is orthogonal, then the power method will converge to a different eigenvector with-
out warning. If you suspect this has happened, run the routine with several differ-
ent initial approximations.
Sample Program
The sample program POWER.PAS provides I/O functions that demonstrate the
power method of approximating eigenvalues.
Input File
Data may be input from a text file. Entries in the text file should be separated by
spaces or carriage returns, and it does not matter if the text file ends with a carriage
return. The format of the text file should be as follows:
1. Dimension of the matrix
2. Elements of the matrix, in the order
[1, 1], [1, 2], ... , [1, N], ... , [N, 1], ... , [N, N],
where N is the dimension of the matrix
Example
Problem. Find the dominant eigenvalue of the matrix:
[ ~ 1~ ~]
024
using the initial guess (1, 2, 3).
Description
Where the power method converges to the dominant real eigenvalue of a matrix
(see POWER.INC), the inverse power method (Burden and Faires 1985, 459-462)
converges to the real eigenvalue nearest to a user-supplied real value. Given a
square matrix A, an initial approximation p to the eigenvalue, and an initial approx-
imation v to the eigenvector, the linear system (A - pI)w = v (where I is the
identity matrix) is solved via LV decomposition (see Chapter 6, "Solving a System
of Linear Equations with Direct Factoring"). The vector w is normalized by divid-
ing through by the element q with the largest absolute magnitude. If the absolute
difference between each of the corresponding elements in v and w is less than a
specified tolerance, then the procedure halts. Otherwise, v is set equal to w, and
the previous matrix equation is solved again. The process repeats until a solution is
reached. The eigenvalue of A closest to p will be (l/q + p), and w will be the
associated eigenvector.
You must supply the matrix A, the initial approximations p and v, and the toler-
ance.
User-Defined Types
Input Parameters
TNArray Size sets an upper bound on the number of elements in each vector. It
is used in the type definition of TNvector and TNmatrix. TNArraySize is not a
variable name and is never referenced by the procedure; hence there is no test for
condition 2. If condition 2 is violated, the program will crash with an Index Out of
Range error (assuming the directive {$R +} is active).
Output Parameters
The inverse power method approximates the solution of a system of linear equa-
tions. If the matrix (Mat - Eigenvalue * I) is singular, where I is the identity rnatrix,
the method will not converge to a solution and Error 5 will be returned. If this
occurs, run the routine again with a slightly different initial approximation,
ClosestVal.
The power method may not converge to repeated eigenvalues with linearly
dependent eigenvectors. Repeated eigenvalues with linearly independent eigen-
vectors do not pose a problem.
The inverse power method is sensitive to the initial approximation (ClosestVal).
If ClosestVal is not close to an eigenvalue or lies midway between two eigenvalues,
the algorithm will converge very slowly, if at all.
The eigenvectors are normalized such that the element of the largest absolute
magnitude in each vector is equal to one.
5mnple Program
Input File
Data may be input from a text file. Entries in the text file should be separated by
spaces or carriage returns, and it does not matter if the text file ends with a carriage
return. The format of the text file should be as follows:
1. Dimension of the matrix
2. Elements of the matrix, in the order
[1, 1], [1, 2], ... , [1, N], "', [N, 1], ... , [N, N],
where N is the dimension of the matrix
3. Elements of the initial guess, in the order [1], [2], ... , [N],
where N is the dimension of the matrix
with an initial guess of (11, 10), you could first create the following text file:
4
1
2
3
4
11
10
Example
Problem. Suppose you know that two of the eigenvalues of the matrix
2 10 0]
010
[
024
are approximately 1.999 and 0.7. Use the inverse power method with an initial
guess of (1, 2, 3) to refine these approximations.
Run INVPOWER.PAS with 1.999 as the approximate eigenvalue:
(K)eyboard or (F)i1e entry of data? K
Dimension of the matrix (1-30)? 3
Matrix[1. 1]: 2
Matrix[1. 2]: 10
Matrix[1. 3]: 0
Matrix[2. 1]: 0
Matrix[2. 2]: 1
Matrix[2. 3]: 0
Matrix[3. 1]: 0
Matrix[3. 2]: 2
Matrix[3. 3]: 4
Now input an initial guess for the eigenvector:
Vector[1]: 1
Vector[2]: 2
Vector[3]: 3
Approximate eigenvalue (default = 5.2857): 1.999
To1 erance (> O. default = 1.0.00E-06): 1E-8
Maximum number of iterations (> O. default = 200): 200
The matrix:
2.00000000000000E+000 1.00000000000000E+00I O.OOOOOOOOOOOOOOE+OOO
O.OOOOOOOOOOOOOOE+OOO 1.00000000000000E+000 O.OOOOOOOOOOOOOOE+OOO
O.OOOOOOOOOOOOOOE+OOO 2.00000000000000E+000 4.00000000000000E+000
Approximate eigenvalue: 1.99900000000000E+000
Tolerance: 1.00000000000000E-008
Maximum number of iterations: 200
Number of iterations: 4
The approximate eigenvector:
1.00000000000000E+000
9. 12736381850482E-014
-5. 13983108970145E-014
The associated eigenvalue: 2.00000000000091E+000
Run INVPOWERPAS with 0.7 as the approximate eigenvalue:
{K}eyboard or {F}ile entry of data? K
Dimension of the matrix {1-30}? 3
Matrix[l, 1]: 2
Matrix[l, 2]: 10
Matrix[l, 3]: 0
Matrix[2, 1]: 0
Matrix[2, 2]: 1
Matrix[2, 3]: 0
Matrix[3, 1]: 0
Matrix[3, 2]: 2
Matrix[3, 3]: 4
Now input an initial guess for the eigenvector:
Vector[l]: 1
Vector[2]: 2
Vector[3]: 3
Approximate eigenvalue {default = 5.2857}: 0.7
Tolerance {> 0, default = 1.000E-06}: 1E-8
Maximum number of iterations {> 0, default = 200}: 200
Direct output to one of the following:
{S}creen
{P}rinter
{F)ile
Number of iterations: 12
The approximate eigenvector:
1.00000000000000E+OOO
-1.00000002395103E-001
6.66666682633328E-002
Description
User-Defined Types
Input Parameters
Output Parameters
Comments
Input File
Data may be input from a text me. Entries in the text me should be separated by
spaces or carriage returns, and it does not matter if the text me ends with a carriage
return. The format of the text me should be as follows:
1. Dimension of the matrix
2. Elements of the matrix, in the order
[1, 1], [1, 2], ... , [1, N], ... , [N, 1], "', [N, N],
where N is the dimension of the matrix
Example
Problem. Find all real eigenvalues and eigenvectors of the matrix
2 100]
[
o 10
o 24
using an initial guess of (1, 2, 3).
The matrix:
2.00000000000000E+000 1.00000000000000E+001 O.OOOOOOOOOOOOOOE+OOO
O.OOOOOOOOOOOOOOE+OOO 1.00000000000000E+000 O.OOOOOOOOOOOOOOE+OOO
O.OOOOOOOOOOOOOOE+OOO 2.00000000000000E+000 4.00000000000000E+000
Tolerance: 1.00000000000000E-006
Maximum number of eigenvalues/eigenvectors to ,find: 3
Maximum number of iterations: 200
Number of iterations: 10
The approximate eigenvector:
-8.32731765653921E-007
4.60590249431080E-015
1.00000000000000E+000
The associated eigenvalue: 4.00000000000004E+000
Number of iterations: 0
The approximate eigenvector:
1.00000000000000E+000
O.OOOOOOOOOOOOOOE+OOO
O.OOOOOOOOOOOOOOE+OOO
Number of iterations: 0
The approximate eigenvector:
1.00000000000000E+000
-9.99999888969117E-002
6.66666592646070E-002
The associated eigenvalue: 9.99999999999991E-001
The exact solution is
Eigenvalue = 4; Eigenvector = (0, 0, 1)
Eigenvalue = 2; Eigenvector = (1, 0, 0)
Eigenvalue = 1; Eigenvector = (1, - 0.1, 2/30)
Description
The eigensystem of a symmetric matrix can be computed much more simply and
efficiently than the eigensystem of an asymmetric matrix. The cyclic Jacobi method
(Atkinson and Harley 1983, 154-160) is an iterative technique for approximating
the complete eigensystem of a symmetric matrix to within a given tolerance. It
consists of multiplying the matrix A by a series of rotation matrices RI . The rotation
matrices are chosen so that the elements of the upper triangular part of A (exclud-
ing the diagonal) are systematically annihilated; that is, Rl is chosen so that A[l, 2]
becomes zero, R2 is chosen so that A[l, 3] becomes zero, and so on. Since the matrix
is symmetric, this will also annihilate the lower triangular part of A. Because each
rotation wili probably change the value of elements annihilated in previous rota-
tions, the method is iterative. Eventually, the matrix will be diagonalized. The
eigenvalues will be the elements of the main diagonal of the diagonal matrix; the
eigenvectors will be the corresponding rows of the matrix created by the product of
the rotation matrices RI •
User-Defined Types
Input Parameters
Output Parameters
Sample Program
Input File
Data may be input from a text file. Entries in the text file should be separated by
spaces or carriage returns, and it does not matter if the text file ends with a carriage
return. The format of the text file should be as follows:
1. Dimension of the matrix
2. Elements of the matrix, in the order
[1, 1], [1, 2], "', [1, N], ... , [N, 1], ... , [N, N],
where N is the dimension of the matrix
[ ~ i ={ =~ 1
-3 -1
-1 -3
1
2
2
1
Run JACOBI.PAS:
{K)eyboard or {F)ile entry of data? F
The matrix:
1.000000000 2.000000000 -3.000000000 -1.000000000
2.000000000 1.000000000 -1.000000000 -3.000000000
-3.000000000 -1.000000000 1.000000000 2.000000000
-1.000000000 -3.000000000 2.000000000 1.000000000
Tolerance: 1.00000000000000E-008
Maximum number of iterations: 200
Number of iterations: 4
The approximate eigenvector:
1.00000000000000E+000
1.00000000000000E+000
-1.00000000000000E+000
-1.00000000000000E+000
155
boundary value problems for second-order ordinary differential equations.
Note that these routines work only with ordinary differentiai equations, not par-
tial differential equations. All of the routines in this chapter can solve problems
involving nonlinear equations (except the linear-shooting routine LINSHOT
2.1NC).
Two one-step techniques that solve initial value problems for first-order ordinary
differential equations are implemented. The first technique employs the
fourth-order Runge-Kutta method (RUNGE_I.lNC), also known as the classical
Runge-Kutta method. The second employs the Runge-Kutta-Fehlberg method
(RKF_l.INC).
Each one-step technique approximates the value of the dependent variable at a
mesh point, which is a value of the independent variable, by using only the infor-
mation obtained from the preceding mesh point. The Runge-Ku~ta method
employs equally spaced mesh points. On the other hand, the Runge-Kutta-
Fehlberg method varies the spacing of the mesh points in order to control the local
truncation error. This produces a corresponding bound on the global error.
The Adams-BashforthIAdams-Moulton predictor/corrector method (ADAMS_I.
INC) is a multistep method that uses information obtained at several preceding
mesh points to approximate the value of the dependent variable at the current
mesh point. The procedure employs the Adams-Bashforth four-step method to
obtain a predictor. It is subsequently used as input for the Adams-Moulton three-
step method to obtain a corrector. The corrector is the approximate value of the
solution. Mesh points are equally spaced, and the starting values for the process are
determined by the one step, fourth-order Runge-Kutta method.
The Runge-Kutta methods are the most reliable and should be used when you
are uncertain of the behavior of the differential equation (for example, if the solu-
tion to the differential equation is not very smooth). If you want the output to be
evenly spaced (in x) or do not require a high degree of accuracy, use the classical
Runge-Kutta method. Otherwise, the Runge-Kutta-Fehlberg method is the best
general purpose routine to use, since it provides control over the accuracy of the
solution.
The Adams-Bashforth/Adams-Moulton method achieves the same accuracy (for
equally spaced mesh points) as the fourth-order Runge-Kutta formula, but if is
significantly faster. Consequently, the Adams-Bashforth/Adams-Moulton method is
the most desirable method if you are reasonably certain that the differential equa-
tion is well-behaved.
Initial value problems for first-order ordinary differential equations are guaran-
teed to have a unique solution on the interval a, b if the function
x' = f(t, x)
The Runge-Kutta method can be generalized for any order ordinary differential
equation. The file RUNGE_N.INC contains an algorithm that can solve an initial
value problem for an nth-order differential equation with the fourth-order Runge-
Kutta formulas. The Lipshitz condition can be generalized for any order ordinary
differential equation. (For details, consult the reference book listed in the section,
"Solution to an Initial Value Problem for a First-Order Ordinary Differential Equa-
tion Using the Runge-Kutta Method.")
Although RUNGE-N.INC can be used to solve initial value problems for first-
order and second-order ordinary differential equations, we recommend that
RUNGE_l.INC and RUNGE--2.INC be used instead. The notation used by these
routines is somewhat simpler than the general case. There is no significant differ-
ence in computation time between the general program (RUNGE_N.INC) and the
specific programs (RUNGE_l.INC and RUNGE--2.INC).
Systems of coupled equations may also be solved with Runge-Kutta techniques.
A system of up to ten first-order ordinary differential equations can be solved with'
the file RUNGE_Sl.INC. A system of up to ten second-order ordinary differen-
tial equations can be solved with the file RUNGE_S2.INC. The algorithms in
both these files are based on the classical Runge-Kutta method with uniform spac-
ing between mesh points; hence, they do not allow for accuracy control (as
in the Runge-Kutta-Fehlberg method). (The Lipshitz condition for systems of
equations is given in the reference in the sections about RUNGE_Sl.INC and
RUNGE_S2.INC.)
Description
This example uses the Runge-Kutta method (Burden and Faires 1985,220-227) to
approximate the solution to a first-order ordinary differential equation with a speci-
fied initial condition.
Given a function of the form
dx/dt = TNTargetF(t, x)
which satisfies the conditions given at the beginning of this chapter, and an initial
condition
x[LowerLimit] = Xlnitial
and spacing
h = (UpperLimit - LowerLimit)/Numlntervals
the fourth-order Runge-Kutta method approximates x in the interval [LowerLimit,
UpperLimit].
The fourth-order Runge-Kutta formulas consist of the following:
F1 h * TNTargetF(t, x[t])
=
F2 h * TNTargetF(t + h/2, x[t] + Fl/2)
=
F3 h * TNTargetF(t + h/2, x[t] + F2/2)
=
F4 h * TNTargetF(t + h, x[t] + F3)
=
x[t + 1] = x[t] + (F1 + 2 * F2 + 2 * F3 + F4)/6
where t ranges from LowerLimit to UpperLimit in steps of h. These formulas give a
4
truncation error of order h •
You must supply LowerLimit, UpperLimit, Xlnitial, Numlntervals, and TNTar-
getF.
dx/dt = TNTargetF(t, x)
The function TNTargetF(t, x) is a user-defined function that calculates the deriv-
ative dx/dt.
Input Pararneters
Output Pararneters
Comments
This procedure will compute Numlntervals values in its calculations; however, you
will rarely need to use all the values. The vectors TValues and XValues will contain
only NumReturn values at roughly equally spaced t-values between the lower and
upper limits. (They will be equally spaced only when Numlntervals is a multiple of
NumReturn.) Thus, you can ensure a highly accurate solution (by making Numln-
tervals large) without generating an excessive amount of output (by making Num-
Return small).
The Runge-Kutta method uses the New/Dispose procedures to manipulate the
heap and should not be used in any program that uses Mark/Release to manipulate
the heap.
Warning: A stiff differential equation occurs when there are at least two very
different scales of the independent variable on which the dependent variable{s) is
changing; for example, y = x + e - lOOx. The Runge-Kutta method may generate a
numerical solution that bears no resemblance to the exact solution of the differen-
tial equation. This unstable numerical solution usually grows exponentially and
may be oscillatory. However, if the exact solution of the differential equation grows
as the independent variable increases, the instability may be difficult to detect. If a
suspected instability has been encountered, reduce the interval size (Numlnter-
vals).
Sample Program
The sample program RUNGE_l.PAS provides I/O functions that demonstrate the
Runge-Kutta method of solving initial value problems. Note that the file
RUNGE_l.INe is included after the function TNTargetF is defined.
Problem. Solve the following initial value problem with the Runge-Kutta method:
X' = xlt + t - 1 I:5t:52
x(l) = 1
1. Code the equation into the program RUNGE_l.PAS:
funct;on TNTargetF(t, X : Real) : Real;
(************************************************************************)
(**** THIS IS THE FIRST-ORDER DIFFERENTIAL EQUATION ****)
(************************************************************************)
beg;n
TNTargetF := x/t + t - 1
end; { function TNTargetF }
2. Run RUNGE_I.PAS:
Lower limit of interval? 1
X value at t = 1.00000000E+00:
t X
1.00000000 1.00000000000000E+000
1.10000000 1.10515880220649E+000
1.20000000 1.22121413182916E+000
1.30000000 1.34892645616477E+000
1.40000000 1.48893886869362E+000
1.50000000 1.64180233779216E+000
1.60000000 1.80799419315265E+000
1.70000000 1.98793197313186E+000
1.80000000 2. 18198400310574E+000
1.90000000 2.39047761619428E+000
2.00000000 2.61370563879444E+000
The exact solution is
x=l - t * In(t)
X(2) = 2.6137056
Description
This example uses the Runge-Kutta-Fehlberg method (Burden and Faires 1985,
230-235) to approximate a solution within a specified tolerance to a first-order
ordinary differential equation with a specified initial condition.
Where the Runge-Kutta method (see RUNGE_l.INe) uses a constant spacing
h, the Runge-Kutta-Fehlberg method varies the spacing so that the solution can be
approximated with accuracy.
Given a function of the form
dx/dt = TNTargetF(t, x)
which satisfies the conditions given at the beginning of this chapter, and an initial
condition
x[LowerLimit] = Xlnitial
both the fourth-order and fifth-order Runge-Kutta formulas are used to approxi-
mate x in the interval [LowerLimit, UpperLimit]. The number of subintervals is
continually increased until the fractional difference between the results of the
5
fourth-order and fifth-order formulas (which give a truncation error of h4 and h ,
respectively) in each subinterval is less than the specified tolerance.
You must supply LowerLimit, UpperLimit, Tolerance, and TNTargetF.
Output Parameters
This procedure will compute more values in its calculations than it will return in
the vectors TValues and XValues. The vectors TValues and XValues will contain
only NumReturn values at subintervals between the lower and upper limits. More
values will be returned in regions of large functional variation than in regions· of
small functional variation. Thus, you can ensure a highly accurate solution (by
making the Tolerance small) without generating an excessive amount of output (by
making NumReturn small).
The Runge-Kutta-Fehlberg method improves the accuracy in the solution by
reducing the spacing between successive values of t. However, if the Tolerance is
too small, the spacing required to reach Tolerance may be beyond the machine's
limit of precision. Consequently, the routine will not converge to a solution that
meets the required Tolerance and Error 5 will be returned.
The Runge-Kutta-Fehlberg method uses the New/Dispose procedures to manip-
ulate the heap and should not be used in any program that uses Mark/Release to
manipulate the heap.
Warning: A stiff differential equation occurs when there are at least two very
different scales of the independent variable on which the dependent variable(s) is
changing; for example, y = x + e- 1OOx • The Runge-Kutta-Fehlberg method may
generate a numerical solution that bears no resemblance to the exact solution of the
differential equation. This unstable numerical solution usually grows exponentially
and may be oscillatory. However, if the exact solution of the differential equation
grows as the independent variable increases, the instability may be difficult to
detect. If a suspected instability has been encountered, reduce the interval size
(NumI ntervals).
Sample Program
The sample program RKF_l.PAS provides I/O functions that demonstrate the
Runge-Kutta-Fehlberg method of solving initial value problems. Note that the file
RKF_l.INC is included after the function TNTargetF is defined.
begin
TNTargetF := x/t + t - 1;
end; { function TNTargetF }
2. Run RKF_l.PAS:
Lower limit of interval?
X value at t = 1.00000000E+00:
t X
1.00000000 1.00000000000000E+000
1.10000000 1.10515881708653E+000
1.20000000 1.22121416069278E+000
1.30000000 1.34892649817459E+000
1.40000000 1.48893892310351E+000
1.50000000 1.64180240395245E+000
1.60000000 1.80799427050390E+000
1.70000000 1.98793206119471E+000
1.80000000 2.18198410146987E+000
1.90000000 2.39047772450816E+000
2.00000000 2.61370575675625E+000
X value at t = 1.00000000E+00:
t X
1.00000000 1.00000000000000E+000
1.12208941 1. 12982837391732E+000
1.20585321 1.22836146826667E+000
1.29271260 1.33921121906568E+000
1.38286653 1.46405185209736E+000
1.47648998 1.60468229863568E+000
1.57374241 1.76304147973215E+000
1.67477301 1.94122165006705E+000
1.77972398 2. 14148082447423E+000
1.88873280 2.36625482837546E+000
2.00193373 2.61816928222327E+000
The exact solution is
x= t 2 - t In(t)
X(2) = 2.6137056
X(2.00193373) = 2.61B1693
In the first run, a solution could be approximated within tolerance with a spacing
of 0.1. In the second run, the algorithm had to vary the spacing in order to approxi-
mate a solution within the tolerance.
Description
User-Defined Types
dxldt = TNTargetF(t, x)
Input Parameters
Output Parameters
Comments
This procedure will compute Numlntervals values in its calculations; however, you
will rarely need to use the values. The vectors TValues and XValues will contain
only NumReturn values at roughly equally spaced t-values between the lower and
upper limits. (They will be equally spaced only when Numlntervals is a multiple of
NumReturn.) Thus, you can ensure a highly accurate solution (by making Numln-
tervals large) without generating an excessive amount of output (by making Num-
Return small).
The Adams-Bashforth/Adams-Moulton method uses the New/Dispose proce-
dures to manipulate the heap and should not be used in any program that uses
Mark/Release to manipulate the heap.
Warning: A stiff differential equation occurs when there are at least two very
different scales of the independent variable on which the dependent variable(s) is
changing; for example, y = x + e- • The Adams-Bashforth/Adams-Moulton
1OOx
method may generate a numerical solution that bears no resemblance to the exact
solution of the differential equation. This unstable numerical solution usually
grows exponentially and may be oscillatory. However, if the exact solution of the
differential equation grows as the independent variable increases, the instability
may be difficult to detect. If a suspected instability has been encountered, reduce
the interval size (Numlntervals).
Sample Program
The sample program ADAMS_l.PAS provides I/O functions that demonstrate the
Adams-Bashforth/Adams-Moulton predictor/corrector method of solving initial
value problems. Note that the file ADAMS_l.lNC is included after the function
TNTargetF is defined.
begin
TNTargetF := x/t + t - 1;
end; { function TNTargetF }
2. Run ADAMS_l.PAS:
Lower limit of interval? 1
X value at t = 1.00000000E+00:
t X
1.00000000 1.00000000000000E+000
1.10000000 1.10515880229293E+000
1.20000000 1.22121413201736E+000
1.30000000 1.34892645643801E+000
1.40000000 1.48893886904034E+000
1.50000000 1.64180233820416E+000
1.60000000 1.80799419362396E+000
1.70000000 1.98793197365806E+000
1.80000000 2. 18198400368348E+000
1.90000000 2.39047761682098E+000
2.00000000 2.61370563946811E+000
The exact solution is
X = t t In(t)
2
-
x(2) = 2.6137056
Initial Value and Boundary Value Methods 171
Solution to an Initial Value Problem/or a Second-Order
Ordinary Differential Equation
Using the Runge-Kutta Metlwd (RUNGE....2.INC)
Description
User-Defined Function
Input Parameters
Comments
This procedure will compute Numlntervals values in its calculations; however, you
will rarely need to use all these values. The vectors TValues, XValues, and XDeriv-
Values will contain only NumReturn values at roughly equally spaced t-values
between the lower and upper limits. (They will be equally spaced only when
Numlntervals is a multiple of NumReturn.) Thus, you can ensure a highly accurate
solution (by making Numlntervals large) without generating an excessive amount
of output (by making NumReturn small).
Warning: A differential equation occurs when there are at least two very differ-
ent scales of the independent variable on which the dependent variable(s) is chang-
ing; for example, y = x + e -IOOx. The Runge-Kutta method may generate a
numerical solution that bears no resemblance to the exact solution of the differen-
tial equation. This unstable numerical solution usually grows exponentially and
may be oscillatory. However, if the exact solution of the differential equation grows
as the independent variable increases, the instability may be difficult to detect. If a
suspected instability has been encountered, reduce the interval size (Numlnter-
vals).
The sample program RUNGE-2.PAS provides I/O functions that demonstrate the
Runge-Kutta method of solving initial value problems for second-order ordinary
differential equations. Note that the file RUNGE-2.INC is included after the func-
tion TNTargetF is defined.
Example
Problem. A weight with mass m lies on a frictionless table and is connected to a
spring with spring constant k:
'4-Wall
F( w)
k
m
. . It ~
FnctlOn ess surlace
If the weight is subject to a driving force F sin(oo t) (00 represents the frequency
of the driving force and t is time), the equation of motion of the mass is as follows:
m d2x/dt2 + k x = F sin(oo t)
Given
m = 2kg
F=9N
k = 32 N/m
00 = 5 cycles/sec
x(O) = 0 m
dx(O)/dt = -2.5 m/sec
find the position and velocity of the block from t = 0 second to t = 2 seconds.
1. Rewrite the preceding second-order differential equation:
d 2x/dt2 = F/m sin(oo t) - kim x
begin
TNTargetF := 9/2 * Sin (5 * t) - 32/2 * x;
end; { function TNTargetF }
3. Run RUNGE--2.PAS:
Lower limit of interval? 0
Upper limit of interval? 2
t Value of X Derivative of X
0.00000000 O.OOOOOOOOOOOOOOE+OOO -2.50000000000000E+000
0.20000000 -4.20735284275848E-001 -1.35075642830665E+000
0.40000000 -4.54648724216734E-001 1.04036531118478E+000
0.60000000 -7.05605786993375E-002 2.47497991717220E+000
0.80000000 3.78400378699554E-001 1.63411037473655E+000
1.00000000 4.79461767300631E-001 -7.09151289407567E-001
1.20000000 1.39708469016311E-001 -2.40042152228323E+000
1.40000000 -3.28491796183335E-001 -1.88475529635974E+000
1.60000000 -4.94677974769030E-001 3.63745224811839E-001
1.80000000 -2.06059519715175E-001 2.27781864414105E+000
2.00000000 2. 72008842396951E-001 2.09767516082021E+000
F., cos{oo t)
dx/dt = - - - - - - -
2 2
m (oo - (0 ) 0
where 00
0
is the natural frequency of the system
2
000 = kim
The period of oscillation is given by
t = 2 'Tr/oo = 1.257 sec
The data is taken from a function of which the derivative could be computed
exactly. Following are the actual values:
t Values of X Derivative of X
0.0 O.OOOOOOOOOOOOE + 000 - 2.500000000000E + 000
0.2 - 4.207354924039E - 001 ...:. 1.350755764670E + 000
0.4 - 4.546487134128E - 001 1.040367091367E + 000
0.6 -7.056000402993E - 002 2.474981241501E + 000
0.8 3.784012476539E - 001 1.634109052159E + 000
1.0 4.794621373315E - 001 -7.091554636580E - 001
1.2 1.397077490994E - 001 - 2.400425716625E + 000
1.4 - 3.284932993593E - 001 - 1.884755635858E + 000
1.6 -4.946791233116E - 001 3.637500845215E - 001
1.8 - 2.060592426208E - 001 2.277825654711E + 000
2.0 2.720105554446E - 001 2.097678822691E + 000
Description
x[LowerLimit] = a 1
x(l)[LowerLimit] = a2
x(n-l)[LowerLimit] = an
and spacing
h = (UpperLimit - LowerLimit)/Numlnteroals
rewrite the nth-order differential equation as n first-order differential equations:
X(l) = Y
I
X(2) = y(I)I = Y2
X(3) = y(I)2 = Y3
x(n-I) = Y
(1)
n-2 = Yn - I
FIx = h * yJt]
FIY l = h * Yz[t]
FIy z
lI
_ = h * y _Jt] lI
F2Y n_l = h * TNTargetF(t + h/2, x[t] + Flx/2, yl[t] + Fly/2, ... , yn_Jt]
+ FIy n_/2)
F3x = h * (yJt] + F2y/2)
F3Yl = h * (Yz[t] + F2y/2)
F3Y l
II
_ = h * TNTargetF(t + h/2, x[t] + F2x/2, yJt] + F2y/2, ... , yn_Jt]
+ F2Yn_/2)
F4x = h * (yJt] + F3y)
F4Yl = h * (Yz[t] + F3yz)
User-Defined Types
TNRowSize is an upper bound for the number of values returned for a particular
variable (NumReturn). TNColumnSize is an upper bound for the order of the differ-
ential equation (Order).
User-Defined Function
Input Parameters
Output Parameters
Comments
The first row of SolutionValues will be the values of t between the limits, the
second row of SolutionValues will be the values of x between the limits, the third
row of SolutionValues will be the values of X(l) between the limits, and so on.
This procedure will compute Numlntervals values in its calculations; however,
you will rarely need to use all those values. The rows of SolutionValues will contain
only NumReturn values at roughly equally spaced t-values between the lower and
upper limits. (They will be equally spaced only when Numlntervals is a multiple of
NumReturn.) Thus, you can ensure a highly accurate solution (by making Numln-
tervals large) without generating an excessive amount of output (by making Num-
Return small).
There are no bounds on the order of the differential equation.
Sample Program
The sample program RUNGE_N.PAS provides I/O functions that demonstrate the
Runge-Kutta method of solving initial value problems for high-order ordinary dif-
ferential equations. Note that the file RUNGE-N.lNC is included after the func-
tion TNTargetF is defined.
Example
Problem. Find the solution to the following fourth-order ordinary differential equa-
tion from t = 0 to t = 1:
4
d x(t)/dt = - 4 x(t) d x(t)/dt
3 3
x(O) = 1
dx(O)/dt = - 1
d x(O)/dt = 2
2 2
3 3
d x(O)/dt = - 6
begin
TNTargetF := -4 * V[l] * V[4];
end; { function TNTargetF }
2. Run RUNGE-N.PAS:
Order of the equation (1-10)? 4
Lower limit of interval? 0
Upper limit of interval?
Enter X value at t O.OOOOOOOE+OOO: 1
Derivative 1 of X at t O.OOOOOOOE+OOO: -1
Derivative 2 of X at t O.OOOOOOOE+OOO: 2
Derivative 3 of X at t O.OOOOOOOE+OOO: -6
Number of values to return (1-100)? 10
Number of intervals (>= 10, default = 10)? 100
Direct output to one of the following:
(S)creen
(P)rinter
(F) il e
t Val ue X[2]
0.00000000 -1.00000000000000E+000
0.10000000 -8.26446283273189E-001
0.20000000 -6.94444446826215E-001
0.30000000 -5.91715977923112E-001
0.40000000 -5.10204082090465E-001
0.50000000 -4.44444443661452E-001
0.60000000 -3.90624997971428E-001
0.70000000 -3.46020758007957E-001
0.80000000 -3.08641970911504E-001
0.90000000 -2.77008304743045E-001
1.00000000 -2.49999993429933E-001
t Value X[3]
0.00000000 2.00000000000000E+000
0.10000000 1.50262961438149E+000
0.20000000 1. 15740742373768E+000
0.30000000 9.10332288053840E-001
0.40000000 7.28862989793594E-001
0.50000000 5.92592607536865E-001
0.60000000 4.88281263842229E-001
0.70000000 4.07083261374878E-001
0.80000000 3.42935540127152E-001
0.90000000 2.91587706310718E-001
1.00000000 2.50000010753535E-001
t Value X[4]
0.00000000 -6.00000000000000E+000
0.10000000 -4.09808076056272E+000
0.20000000 -2.89351855059016E+000
0.30000000 -2.10076680857258E+000
0.40000000 -1.56184925333600E+000
0.50000000 -1. 18518520443061E+000
0.60000000 -9.15527359078898E-001
0.70000000 -7. 18382215400418E-001
0.80000000 -5.71559223064178E-001
0.90000000 -4.60401631119694E-001
1.00000000 -3.75000005740567E-001
x(I) = 0.5
dx(I}jdt = - 0.25
d2x(I}jdt2 = 0.25
d3x(I}jdt3 = - 0.375
Description
which satisfies the Lipshitz condition (the Lipshitz condition for first-order and
second-order ordinary differential equations is given at the beginning of this chap-
ter; consult the previous book reference for details of the Lipshitz condition for
systems), and initial conditions
xJLowerLimit] = a l
xJLowerLimit] = a 2
x.JLowerLimit] = am
and spacing
h = (UpperLimit - LowerLimit)/Numlnteroals
the fourth-order general Runge-Kutta method can be used to approximate simulta-
neously the xj's.
User-Defined Types
User-Defined Functions
V[10] = XlO
Input Parameters
Output Parameters
SolutionValues : TNmatrix; Values oft, Xl' X2 ' ... xm between the limits
Error: Byte; 0: No errors
1: NumReturn < 1
2: Numlnteroals < NumReturn
3: NumEquations < 1
4: LowerLimit = UpperLimit
Comments
The first row of SolutionValues will be the values of t between the limits, the
second row of SolutionValues will be the values of Xl between the limits, the third
row of SolutionValues will be the values of x2 between the limits, and so on.
All ten user-defined functions are called from the procedure. If your system has
less than ten equations, you must still define all ten functions or the program will
not compile. The superfluous functions should be defined as follows (TNTargetFIO
is used as an example):
funct;on TNTargetF10(V : TNvector} : Real;
beg;n
end: { function TNTargetF10 }
If you need to solve a system with more than ten equations, then edit the include
file RUNGE_Sl.lNC. The following line should be added to the end of procedure
Step:
F[ll] := Spacing * TNTargetFll(CurrentValues};
More statements (for F[12], and so on) may be added as necessary. All new
functions (for example, TNTargetFll) must be defined in your top-level program.
Note: Before making any changes to the include file, make sure you have a backup
copy.
This procedure will compute Numlntervals values in its calculations; however,
you will rarely need to use these values. The rows of SolutionValues will contain
only NumReturn values at roughly equally spaced t-values between the lower and
upper limits. (They will be equally spaced only when Numlntervals is a multiple of
NumReturn.) Thus, you can ensure a highly accurate solution (by making Numln-
tervals large) without generating an excessive amount of output (by making Num-
Return small).
This routine stores much information on the heap. If you try to accurately solve a
large system (that is, if both NumEquations and Numlntervals are large), you may
get run-time error $FF, Heap/Stack collision. If this happens, the dimension of
TNvector and TNmatrix should be reduced as much as possible. If this is not possi-
ble, then remove any RAM-resident software (for example, SideKick, SuperKey, or
a print buffer).
Sample Program
The sample program RUNGE_S1.PAS provides I/O functions that demonstrate the
Runge-Kutta method of solving initial value problems for systems of first-order
ordinary differential equations. Note that the file RUNGE_S1.INC is included
after the TNTargetF functions are defined.
Example
Problem. A weight with mass m lays on a frictionless table and is connected to a
spring with spring constant k:
~Wall
F( w)
k
m
Frictionlts surface
beg;n
TNTargetFl := V[2];
end; { function TNTargetF1 }
funct;on TNTargetF2{V TNvector): Real;
begin
TNTargetF2 := 9/2 * Sin(5 * V[O]) - 32/2 * V[l];
end; { function TNTargetF2 }
function TNTargetF3(V : TNvector) : Real;
(**************************************************)
(* THIS IS THE THIRD DIFFERENTIAL EQUATION *)
(**************************************************)
(* *)
(* dx[3] *)
(* TNTargetF3(t, x[l], x[2], ..• x[m]) *)
(* dt *)
(* *)
(* The vector V is defined: *)
(* V[0] t *)
(* V[l] X[l] *)
(* V[2] X[2] *)
(* *)
(* *)
(* *)
(* V[m] X[m] *)
(* *)
(* where m is the number of coupled equations. *)
(**************************************************)
begin
end; { function TNTargetF3 }
Functions TNTarget4 to TNTargetlO should be defined like the function
TNTargetF3.
F", cos{w t)
dxldt = -------
2 2
m (wo - ( )
Description
2 2
d x",Idt = TNTargetFm(t, Xl' X'I' X2' X' 2' .. ,' X rn ' X' J
where x') indicates dx/dt, which satisfies the Lipshitz condition (the Lipshitz con-
dition for first-order and second-order ordinary differential equations is given at
the beginning of this chapter; consult the previous book reference for details of the
Lipshitz condition for systems), and initial condition
xJLowerLimit] =a l
x' JLowerLimit] = hI
X' 2[LowerLimit] = h2 .
xJLowerLimit] = am x'JLowerLimit] = hm
and spacing
h = (UpperLimit - LowerLimit)/Numlntervals
dx)dt = Ym
dx)dt = TNTargetFm(t, Xl' YI, X2' Y2' "', Xm' yJ
Then the fourth-order general Runge-Kutta method can be used to approximate
the xj's and the yj's simultaneously,
The general Runge-Kutta formulas for these equations are as follows:
Flxl = h * Yl
FIY I = h * TNTargetFI(t, xl[t], yJt], x2[t], y2[t], "', xJt], yJt])
Flx2 = h * Y2
FIY2 = h * TNTargetF2(t, xJt], yJt], x2[t], y2[t], "', xJt], yJt])
Flxm = h * Ym
FIY m = h * TNTargetFm(t, xJt], yJt], x2[t], y2[t], .. ,' xJt], yJt])
User-Defined Types
TNData = record
x : Real;
xDeriv : Real;
end; { TNData record}
TNvector = array[O •• TNRowSize] of TNData;
TNmatrix = array[O •• TNColumnSize] of TNvector;
TNRowSize is an upper bound for the number of values returned for a particular
variable (NumReturn). TNColumnSize is an upper bound for the number of second-
order differential equations (NumEquations).
V[O].X =t
V[l].x = x[l]
V[l].xDeriv = x'[l]
V[2].x = x[2]
V[2].xDeriv = x ' [2]
V[10].x = x[10]
V[10].xDeriv = x ' [10]
The procedure defined in RUNGE_S2.INC solves this system of coupled differ-
ential equations (a maximum of ten equations). All ten functions must be defined,
even if your system contains less than ten equations.
Input Parameters
Output Parameters
SolutionValues : TNmatrix; Values oft, Xi' and x'} between the limits
Error: Byte; 0: No errors
1: NumReturn < 1
2: Numlntervals < NumReturn
3: NumEquations < 1
4: LowerLimit = UpperLimit
Comments
The first row of SolutionValues will be the values of t between the . limits, the
second row of Solution Values will be the values of Xl and XiI between the limits, the
third row of Solution Values will be the values of x2 and x' 2 between the limits, and
so on.
All ten user-defined functions are called from the procedure. If your system has
less than ten equations, you must still define all ten functions or the program will
not compile. The superfluous functions should be defined as follows (TNTargetF10
is used as an example):
function TNTargetFIO(V : TNvector) : Real;
begin
end; { function TNTargetFIO }
numerical solution that bears no resemblance to the exact solution of the differen-
tial equation. This unstable numerical solution usually grows exponentially and
may be oscillatory. However, if the exact solution of the differential equation grows
as the independent variable increases, the instability may be difficult to detect. If a
suspected instability has been encountered, reduce the interval size (Numlnter-
vals).
Example
Problem.Two weights of mass m each hang from a pendulum of length l and are
connected by a spring with spring constant k:
Ceiling
k
®-~
y UUUUUUUUUUUU X
where g is the acceleration due to gravity, t is time, and x and yare the displace-
ments of the two weights from their rest positions. Given
m = 2kg
k = 32 N/m
2
g = 9.8 m/sec
l = 0.6l25 m
x(O) = 1 m
y(O) = -1 m
dx(O)/dt = 0 m/sec
dy(O)/dt = 0 m/sec
var
t Real;
begin
t := v[O] .x;
TNTargetFl .- -9.8 * V[1].x/O.6125 - 32/2 * (V[l].x - V[2].x);
end; { function TNTargetFl }
var
t Real;
begin
t := v[O] .x;
TNTargetF2 := -9.8 * V[2].x/O.6125 + 32/2 * (V[l].x - V[2].x);
end; { function TNTargetF2 }
var
t Real;
begin
end; { function TNTargetF3 }
Functions TNTargetF 4 to TNTargetF10 should be defined like function TNTargetF3.
3. Run RUNGE_S2.PAS:
Number of second order equations: (1-20)? 2
Lower limit of interval? 0
Upper limit of interval?
Enter X[l] value at t O.OOOOOOOOE+OO: 0.01
Enter X'[l] value at t O.OOOOOOOOE+OO: 0.00
Enter X[2] value at t = O.OOOOOOOOE+OO: -0.01
Enter X'[2] value at t = O.OOOOOOOOE+OO: 0.00
Number of values to return (1-100)? 10
Number of intervals (>= 10, default=10)? 100
Description
This example uses the shooting method to approximate the solution to a second-
order ordinary differential equation with specified boundary conditions (Burden
and Faires 1985,526-531).
Given a second-order differential equation (Burden and Faires 1985, 261-269) of
the form
2 2
d y/dx = TNTargetF(x, y, y')
where y' represents dy/dx, which satisfies the conditions given at the beginning of
this chapter, boundary conditions
y[LowerLimit] = Lowerlnitial
y[UpperLimit] = Upperlnitial
and spacing
h = (UpperLimit - LowerLimit)/Numlntervals
and an initial approximation (guess) to the slope at LowerLimit
y'[LowerLimit] = InitialSlope
the shooting method first solves the second-order initial value problem (using the
method described in RUNGE-2.1NC). Based on a comparison of the solution at
UpperLimit with the boundary condition Upperlnitial, a new approximation to the
slope at LowerLimit is made. In this way, a new "shot" at the solution is made by
observing the result of the previous "shot." Subsequent iterations use information
from two previous shots and the secant method (see Chapter 2, "Roots of a Func-
tion Using the Secant Method") to approximate the slope at LowerLimit. This pro-
cess is repeated until the fractional difference between successive approximations
to the boundary condition at UpperLimit is less than a specified tolerance.
You must supply the LowerLimit, UpperLimit, LowerInitial, Upperlnitial,
InitialSlope, Numlntervals, Tolerance, and TNTargetF.
User-Defined Function
Input Parameters
Comments
The parameter Tolerance can be misleading. The shooting method converges to the
initial slope, which satisfies the upper boundary condition. Convergence is
achieved when the fractional difference between Upperlnitial and the upper
boundary approximation is less than Tolerance. This does not mean that every
value between the boundaries has been approximated with the same degree of
accuracy. To improve the accuracy of the entire approximation, increase the num-
ber of intervals. The example demonstrates the different effects of Tolerance and
Numlntervals.
The shooting algorithm will compute Numlntervals values in its calculations.
However, you will rarely need to use all those values. The vectors XValues,
YValues, and YDerivValues will contain only NumReturn values at roughly equally
spaced t-values between the lower and upper limits. (They will be equally spaced
only when Numlntervals is a multiple of NumReturn.) Thus, you can ensure a
highly accurate solution (by making Numlntervals large) without generating an
excessive amount of output (by making NumReturn small).
cal solution that bears no resemblance to the exact solution of the differential equa-
tion. This unstable numerical solution usually grows exponentially and may be
oscillatory. However, if the exact solution of the differential equation grows as the
independent variable increases, the instability may be difficult to detect. If a sus-
pected instability has been encountered, reduce the interval size (Numlntervals).
Sample Program
The sample program SHOOT2.PAS provides I/O functions that demonstrate the
shooting method of solving boundary value problems. Note that the file SHOOT
2.1NC is included after the function TNTargetF is defined.
Example
Problem. Use the nonlinear shooting method to solve the following boundary value
problem:
Y 1/ = 192 sqr(y/y') 0:5x:51
y(l) = 1
y(2) = 16
begin
TNTargetF := 192 * Sqr(y/yPrime);
end; {function TNTargetF}
2. Run SHOOT2.PAS:
Lower limit of interval? 0
Upper limit of interval? 1
Enter Y value at X = O.OOOOOOOOE+OO: 1
Enter Y value at X = 1.00000000E+00: 16
Enter a guess for the slope at X = O.OOOOOOOOE+OO (default=1.50E+01): 15
Number of points returned (1-500)? 10
Number of intervals (>= 10, default=10)? 10
Tolerance (> 0, default = 1.000E-06)? 1E-12
Maximum number of iterations (> 0, default 100)? 100
Direct output to one of the following:
(S)creen
(P) ri nter
(F) il e
X Y Value Derivative of Y
O.OOOOOOOOOOOOOOE+OOO 1.00000000000000E+000 4.00053795390884E+000
1.00000000000000E-001 1.46417721408153E+000 5.32386904044879E+000
2.00000000000000E-001 2.07370562259973E+000 6.91162114244397E+000
3.00000000000000E-001 2.85621262766442E+000 8.78752756627335E+000
4.00000000000000E-001 3.84170902091389E+000 1.09754927855527E+001
5.00000000000000E-001 5.06259931530967E+000 1.34994802016423E+001
6.00000000000000E-001 6.55368547624580E+000 1.63834750611955E+001
7.00000000000000E-001 8.35216836918581E+000 1.96514712240017E+001
8.00000000000000E-001 1.04976483580762E+001 2.33274661179548E+001
9.00000000000000E-001 1.30321255669365E+001 2.74354587043771E+001
1.00000000000000E+000 1.60000000000094E+001 3. 19994486182108E+001
Now solve the same problem using a smaller spacing, but with a greater tolerance:
Lower limit of interval? 0
Upper limit of interval? 1
Description
This example uses the linear shooting method to approximate the solution to a
second-order linear ordinary differential equation with specified boundary condi-
tions (Burden and Faires 1985, 519-524).
Given a second-order differential equation (Burden and Faires 1985, 261-264) of
the form
d2y/dx2 = TNTargetF(x, y, y')
which is linear in both y and y', where y' represents dy/dx, and which satisfies the
conditions given at the beginning of this chapter, boundary conditions
y[LowerLimit] = Lowerlnitial
y[UpperLimit] = Upperlnitial
and spacing
h = (UpperLimit - LowerLimit)/Numlntervals
the shooting method solves the two initial value problems (see RUNGE-2.INC):
y'[LowerLimit] =0 y[LowerLimit] = Lowerlnitial
y'[LowerLimit] =1 y[LowerLimit] = Lowerlnitial
(These values are particular to this implementation; any other nonidentical set of
initial conditions will suffice.) Since neither of these initial values of y' is likely to
be correct, the solutions generated are not likely to satisfy the boundary condition
at Upperlnitial. However, because of the linearity of the equation, an appropriate
linear combination of these two solutions will be a solution to the boundary value
problem. The linear shooting method requires that only two initial value problems
be solved, where the ordinary shooting method (SHOOT2.INC) requires that an
unknown number of initial value problems be solved before the method converges
to a solution.
You must supply the LowerLimit, UpperLimit, Lowerlnitial, Upperlnitial,
Numlntervals, and TNTargetF.
User-Defined Function
Input Parameters
Output Parameters
Comments
The sample program LINSHOT2.PAS provides I/O functions that demonstrate the
linear shooting method of solving boundary value problems. Note that the file
LINSHOT2.1NC is included after the function TNTargetF is defined.
Example
Problem.Use the linear shooting method to solve the following boundary value
problem:
y" = y' /x - y/sqr(x) +1 1 =x S 10
y(1) 1
y(10) 76.974149
beg;n
TNTargetF := yPrime/x - y/Sqr(x) + 1;
end; { function TNTargetF }
2. Run LINSHOT2.PAS:
Lower limit of interval? 1
Upper limit of interval? 10
Enter Y value at X = 1.00000000E+00: 1
Enter Y value at X = 1.00000000E+01: 76.974149
Number of points returned (1-500)? 9
Number of intervals (>= 9, default = 9)? 9
Direct output to one of the following:
(S)creen
(P)rinter
(F)ile
Lower limit: 1.00000000000000E+00
Upper limit: 1.00000000000000E+01
Value of Y at 1.0000: 1.00000000000000E+00
Value of Y at 10.0000: 7.69741490000000E+01
NumIntervals: 9.00000000000000E+00
X Y Value Derivative of Y
1.00000000000000E+000 1.00000000000000E+000 1.00000001942594E+000
2.00000000000000E+000 2.61370547174514E+000 2.30685275847028E+000
3.00000000000000E+000 5.70416298088411E+000 3.90138768358927E+000
4.00000000000000E+000 1.04548224122436E+001 5.61370562650429E+000
5.00000000000000E+000 1.69528103026793E+001 7.39056208402440E+000
6.00000000000000E+000 2.52494430584438E+001 9.20824053324639E+000
6.99999999999999E+000 3.53786288412165E+001 1.10540898579729E+001
7.99999999999999E+000 4.73644675641047E+001 1.29205584690303E+001
8.99999999999999E+000 6.12249787166508E+001 1.48027754364805E+001
9.99999999999998E+000 7.69741490000000E+001 1.66974149235206E+001
The exact solution is
y = x * x - x * In(x)
y(l) = 1 y'(l) = 1
y(10) = 7.6974149 y'(10) = 16.6974149
The second approximation is more accurate.
Given a set of data points, this chapter provides routines to model the data with a
function of a given type. The most common application of this concept is linear
regression.
With linear regression, there is some control variable, say X, and some observed
variable, say Y. X and Yare known or suspected to have some linear relationship,
say
Y=a*X+b
but the parameters a and b are unknown. Usually there is some experimental error
or some other nonlinear influence on Y, so that there are no values of a and b for
which the preceding equation holds exactly. The method of regression provides a
formula for a and b in terms of the values of X and Y such that the error is mini-
mized. The error is the sum of squares of the errors (a * X + b - Y) on each data
point. Except in certain unusual cases, there is exactly one value for a and one
value for b that makes this sum of squares the least possible. This is called the
least-squares solution.
This concept of least squares also applies' when more variables are present-
then it is often called multiple regression. Using this method, you can find the best
model for a given set of data that is linear in a given set of other data sets or
functions. Models that are nonlinear variables can also be treated as long as the
unknown parameters appear linearly.
221
Least-Squares Approximation (LEAST.INC)
Description
where ~(X) is the jth basis vector evaluated at the data value Xli]. A vector Y is
constructed that contains the y-values of the data points. The coefficients of the
basis vectors that form the least-squares approximation will be the n vector C, such
that the Euclidean norm of(AC - Y) (represented by II AC - Y liz) is a minimum.
This requirement is converted to the requirement that
Input Parameters
So 1ut ion: TNRowVector; Coefficients of the basis vectors that form the least-
squares approximation
YFit : TNColumnVector; Values of the least-squares fit evaluated at the XData
values
Res i dua 1s : TN Co 1umnVector; Difference between YData and YFit values
StandardDevi ation : Real; Square root of the variance-indicates the goodness of fit
Error: Byte; 0: No error
1: NumPoints < 2
2: NumTerms < 1
3: NumTerms > NumPoints
4: Least-squares solution does not exist (see "Com-
ments")
Comments
The least-squares routine is kept in two modules (include files). One is called
LEAST.lNC and must always be included in your top-level program. The choice of
the second module will depend upon the functional form (basis vectors) to which
you fit the data. Following are the five basis modules included in this package:
POLY.LSQ
This module uses Chebyshev polynomials to fit a polynomial to the data points.
NumTerms must be one greater than the degree of the polynomial (for example, to
fit a fourth-degree polynomial, input NumTerms = 5). To get a straight-line least-
squares fit, use this module and fit a curve with only two coefficients. The elements
of the Solution vector will be as follows:
Solution[j] = aJ 1 s j S NumTerms
where a. is the coefficient of Xl-I.
}
POWER.LSQ
This module will fit the function
y = al
where a and b are real numbers to the data points. A linear equation is obtained by
taking the log of both sides, like so:
In{y) = In{a) + b * In{x)
and expanding on basis vectors 1 and In{x). The x-values of the data points must all
be positive, and the y-values of the data points must all have the same sign. The
number of coefficients in the approximation will be two regardless of the value of
NumTerms (unless NumTerms > NumPoints, in which case Error 3 will occur).
The elements of the Solution vector will be as follows:
Solution[l] =a
Solution[2] =b
EXP.LSQ
This module will fit the function
y = aebx
where a and b are real numbers to the data points. A linear equation is obtained by
taking the log of both sides, like so:
In{y) = In{a) + bx
WG.LSQ
This module will fit the function
y = a In(bx)
where a and b are real numbers to the data points. A linear equation is obtained by
rewriting the equation:
y = a In(b) + a In(x)
and expanding on basis vectors 1 and In(x). The x-values of the data points must all
have the same sign. The number of coefficients in the approximation will be two
regardless of the value of NumTerms (unless NumTerms > NumPoints, in which
case Error 3 will occur). The elements of the Solution vector will be as follows:
Solution[l] =a
Solution[2] =b
USERLSQ
This module is included if you need a least-squares approximation on a set of basis
vectors different from the ones listed earlier. This module allows you to create your
own set of basis vectors. The source code contains detailed instructions of how to
flesh out the skeleton contained in USER.LSQ.
A least-squares solution may not exist for some input data and choice of basis
vectors (Error 4). The reasons for this will depend on the module you are using. For
example, it is impossible to fit an exponential function (module EXP.LSQ) to data
with y-values of differing signs; Error 4 will occur if you try. The same data could
be fit with a polynomial and no error would result. Error 4 will also occur if all the
x-values of the data are identical.
The demonstration program LEAST.PAS contains I/O routines that allow you to
run the least-squares approximation routine. Note that there are two include com-
mands:
{$I POLY.LSQ} (* load the basis vectors *)
{$I LEAST. INC} (* load procedure LeastSquares *)
The LEAST.lNC file must always be included after the basis module. To change
the basis vectors of the approximation, simply load a different basis module with
the first INCLUDE command.
Input Files
Data may be entered from a text file. The x- and y-coordinates should be separated
by a space and followed by a carriage return. For example, data values of sqr(x)
could be entered in a text file as
11
24
39
416
525
Example
Problem. Given the following data (contained in the file SAMPLE9A.DAT), fit a
fourth-degree polynomial and a logarithmic function to the data:
O.oooooooOOOOOOOE+OO 1.33830225764886E-03
0.10000000000000E+00 4.43184841193803E-02
0.20000000000000E+00 5.39909665131879E-01
0.30000000000000E+00 2.41970724519143E+00
0.40000000000000E+00 3.98942280401433E+00
0.02000000000000E+00 2.91946925791461E-03
0.04000000000000E+00 6.11901930113775E-03
0.06000000000000E+00 1.23221916847303E-02
0.08000000000000E+00 2. 38408820146486E-02
0.12000000000000E+00 7.91545158298001E-02
0.14000000000000E+00 1.35829692336855E-01
0.16000000000000E+00 2.23945302948430E-01
0.18000000000000E+00 3.54745928462313E-01
0.22000000000000E+00 7.89501583008939E-01
0.24000000000000E+00 1.10920834679455E+00
0.26000000000000E+00 1.49727465635745E+00
0.28000000000000E+00 1.94186054983213E+00
0.32000000000000E+00 2.89691552761483E+00
0.34000000000000E+00 3.33224602891800E+00
0.36000000000000E+00 3.68270140303323E+00
0.38000000000000E+00 3.91042693975456E+00
*----------------------------------------*
Polynomial Least Squares Fit
*----------------------------------------*
*----------------------------------------*
Logarithmic Least Squares Fit
*----------------------------------------*
Fourier transforms are used to analyze periodic phenomena such as waves. A con-
tinuous functionfthat has period 2'TT (= 2 * 3.14159265...); that is, satisfies
f(x + 2'TT) = f(x)
for all x, can be decomposed into sines and cosines:
f(x) = a[O] + a[l] * cos(x) + b[l] * sin(x) + a[2] * cos(2x)
+ b[2] * sin(2x) + ...
This is an infinite series where the coefficients get closer and closer to zero. The
routines in this chapter can be used to calculate the coefficients.
The Fast Fourier Transform (FFT) is a particular algorithm for computing Fou-
rier transforms efficiently.
This chapter includes two kinds of units. One group consists of four variations of
the FFT method of calculating discrete Fourier transforms (FFTB2.INC, FFT
B4.INC, FFT87B2.INC, FFT87B4.INC), each optimized for certain conditions.
All are variations of the original Cooley-Tukey method. The second group consists
of six applications (COMPFFT.INC, REALFFT.INC, COMPCNVL.lNC, REAL
CNVL.INC, COMPCORR.INC, REALCORR.INC) that can each be used with
any of the FFT methods. You can select the FFT method most appropriate to the
circumstances and combine it with the appropriate application or integrate it into
another program (Brigham 1974; Nussbaumer 1982).
In each FFT unit the procedure calls have exactly the same form (although there
are different restrictions on the data) so that anyone FFT unit can be combined
233
with any of the application units without rewriting code. Each of these algorithms
will compute either a forward or an inverse transform.
Each unit contains two procedures needed to prepare for the FFf calculation:
procedure Testlnput and procedure MakeSinCosTable. TestInput examines the
input data to ensure that it satisfies certain conditions (for example, that there is
more than 1 data point). MakeSinCosTable precalculates a table of the nth roots of
unity for look up in the FFf calculation.
FFfB2.1NC contains a procedure that implements the Cooley-Tukey powers-of
two (radix 2 or base 2) Fast Fourier Transforms. It is optimized to reduce the
number of real multiplications by taking advantage of the symmetries of certain
roots of unity and by using a complex multiplication that requires three real multi-
plications and three real additions. This algorithm is appropriate when the time for
real multiplications is significantly greater than the time for real additions; for
example, when running on an 8088 machine with no numeric coprocessor.
FFf87B2.1NC implements the same algorithm as FFTB2.1NC. The difference
is that the complex multiplications are done with four real multiplications and two
real additions. By using this standard form of complex multiplication, storage over-
head and assignment statements are reduced. This algorithm is appropriate when
the time for a real multiplication is close to the time for a real addition; for example,
when running an 8088 machine with an 8087 numeric coprocessor, or an 80286
machine with an 80287 numeric coprocessor.
Standard Turbo Pascal uses 6-byte reals; Turbo-87 Pascal (which utilizes the
8087 coprocessor) uses 8-byte reals. Consequently, given the same amount of
memory, more data points can be manipulated in Turbo Pascal than in Turbo-87
Pascal. Both FFfB2.INC and FFf87B2.1NC require the number of data points to
be a power of two up to a maximum of 4,096 points when in Turbo-87, or 8,192
points when in standard Turbo Pascal.
FFfB4.1NC and FFf87B4.1NC implement powers-offour (radix 4 or base 4)
Fast Fourier Transforms. The powers-of-four method is the same as the Cooley-
Tukey algorithm except at each stage of reduction a given transform is converted
into four transforms each with one fourth the data points of its predecessor (Nus-
sbaumer 1982). When this algorithm is optimized, there are about 25 percent
fewer multiplications and slightly fewer additions than in a radix-2 algorithm. The
algorithm has the disadvantage of only being applicable to data sets where the
number of points is a power of four up to a maximum of 4,096 points whether in
Turbo-87 or standard Turbo. A reduction in execution time of about 20 percent is
accomplished when FFfB4.1NC or FFf87B4.1NC is used over its B2 counterpart.
There are six application programs that use the basic FFT routines
contained in the previously mentioned include files (COMPFFT.INC,
REALFFT.lNC, COMPCNVL.lNC, REALCNVL.INC, COMPCORR.INC, and
REALCORR.INC).
Fast Fourier Transforms are particularly useful for analyzing periodic signals.
Such a signal is represented by a function J satisfying
J(t + T) = J(t)
where t is time and T is the period. Under mild hypotheses,J can be expanded into
a Fourier series such as the following:
oc
J(t) = N- 112 I
11= - oc
F(n) exp (2'Ti' i n tiT)
where i is the square root of - 1. The term exp (2'Ti' i n tiT) is a sinusoid of period
Tin and frequency niT, and its coefficient F(n) gives the strength of that frequency
component in the original signal.
(The factor of N- I12 is a matter of convention. Some books do not include it in this
formula, resulting in a factor of liN in the formula for X that follows.)
The formula for the coefficients X(k) is as follows:
I
N-I
1/2
X(k) = N- x(n) exp (-i 2'TT n kiN)
n'" 0
This formula for X makes sense for any integer k. X is then periodic, satisfying
X(k + N) = X(k)
for all k. In formulas and programs, it is more convenient to let k run from t9 °
N - 1, but for analyzing signals it makes more sense to think of k as running from
(- N/2) to (N/2 - 1). This is because values of k near zero represent the low
frequency information, and values of k near or greater than N12 represent frequen-
cies that are so high that the discretization is too coarse to realize them accurately
anyway. Therefore, if k is between N12 and N, X(k) should be thought of as the
coefficient of
exp (2'TT i (k - N) tiT)
rather than
exp (2'TT i k tiT)
In other words, negative frequencies are represented on the right half of the
transform.
m = 0,1, ... N - 1
,,= 0
where subscripts are taken modulo N (circular convolution). The basic theorem that
allows us to calculate these effectively using FFTs is shown in the following:
y In
= X HIn In
m = 0,1, ... N - 1
where capital letters indicate the transforms of the functions represented by lower-
case letters. Thus the procedure 'for convolution works like this:
1. Transform both given data sets using FFTs.
2. Multiply the resulting transforms.point by point.
3. Find the inverse transform of this product using FFTs.
Cm =L
" xn hn+m m = 0,1, ... N - 1
n = 0
where subscripts are taken modulo N (circular convolution). This can be computed
using FFTs with a method analogous to that used in COMPCNVL.INC:
c m
= Xm H N-m m = 0,1, ... N - 1
Commonly x and h are real functions; in which case the preceding formula
reduces to C = X Ii , where the bar stands for complex conjugation. Thus the
procedure fo; correGtio~ works like this:
1. Transform both given data sets using FFTs.
2. Multiply each element of the transform of the first data set by the appropri-
ate element of the transform of the second data.
3. Find the inverse transform of this product using FFTs.
Anyone of the FFT include files can be used with any of the applications.
Data Sampling
While sampling theory is beyond the scope of this Toolbox, we would like to men-
tion several common problems associated with data sampling (Brigham 1974; Press
et al. 1986, Ch.12). The most common frustration is aliasing. A Fourier transform
only represents frequencies up to a certain limit (called the Nyquist limit, or
Nyquist frequency), which is given by half the sampling rate. (For example, if a
signal is sampled sixty times a second, the Nyquist frequency will be 30 Hz.) 'A
sample containing frequencies greater than this limit will not be properly trans-
formed. The high frequencies will falsely contribute to the transform. This contri-
bution will be indistinguishable from a contribution of a frequency below the
Nyquist frequency.
There are several ways to combat aliasing. Increasing the sampling rate will
increase the Nyquist frequency and thus reduce aliasing effects. It is also possible
to pass the signal through a low pass filter, thus eliminating the high frequencies
before sampling. If the Fourier transform of a signal does not converge to zero at
the Nyquist frequency, the transform has very likely been aliased.
The Fourier transform assumes that the sample represents a periodic function
and that the sample is an integer multiple of one period. If the latter condition is
not true, spurious frequencies will show up in the transform. For example, if a sine
wave is sampled from 0 to 270 degrees (instead of the full period), a sharp bound-
ary is created because the sine of 0 does not equal the sine of 270. High frequen-
cies will be introduced into the transform to account for that sharp boundary.
User-Defined Types
These user-defined types are different from others in this Toolbox, because they
involve pointers. Pointers are used to transcend the limitations imposed by the 64K
data segment size of Turbo Pascal. One array of 8,000 elements uses the entire data
segment (in Turbo-87). However, it is possible to store these arrays on the heap,
and to point to them with pointers that only require 4 bytes. The size of the heap
(and hence the maximum size and number of TNvectors) is determined by the
amount of memory in the machine.
Procedure TestInput
Description
This example determines the number of data points in terms of a power of two
[four]. If the number of data points is not a power of two [four], then an error is
returned.
Input Parameters
NumPoi nts : Integer; Number of data points
Output Parameters
NumberOfBits : Byte; Number of data points as a power of two [four]
Error: Byte; 0: No errors
1: NumPoints < 2
2: NumPoints not a power of two [four]
Description
This example creates a look-up table ofNumPoints/2 [3/4 NumPoints] roots of unity.
The roots of unity are defined as follows:
Rootn = exp ( - i 2'TT n/NumPoints n = O..NumPoints/2 [3/4 NumPoints]
where i is the square root of -1. These values are stored in two tables: SinTable,
containing the imaginary parts of the roots of unity, and CosTable, containing the
real parts of the roots of unity. It is faster to look up these values in a table than to
calculate them in the FFf procedure.
Input Pararneters
NumPoints : Integer; Number of data points
Output Pararneters
Si nTabl e : TNvectorPtr; Table of sine values
CosTabl e : TNvectorPtr; Table of cosine values
Procedure FFT
Description
This example implements the particular variation of the Cooley-Tukey algorithm.
The Fast Fourier Transform of the data XReal, Xlmag is made in place and is thus
returned in the vectors XReal, Xlmag. The inverse transform of the data can also
be calculated with this procedure.
Input Parameters
NumberOfBi ts : Byte; Number of data points as a power of two [four]
NumPoi nts : Integer; Number of data points
Inverse: Boolean; FALSE equals forward transform; TRUE equals inverse
transform
XRea 1 ': TNvectorPt r; Pointer to real values of the data points
XImag : TNvectorPtr; Pointer to imaginary values of the data points
SinTable: TNvectorPtr; Table of sine values
CosTable : TNvectorptr: Table of cosine values
Output Parameters
XReal : TNvectorPtr; Pointer to real values of the discrete Fourier transform of the
input data
XImag : TNvectorPtr; Pointer to imaginary values of the discrete Fourier transform
of the input data
Each of the six application programs calls the three procedures contained within
FFT algorithm files.
COMPFFT.INC
Description
This example is the most basic application, performing a complex Fast Fourier
Transform. It simply calls Testlnput, MakeSinCosTable, and FFT sequentially; thus
accomplishing an in-place transformation of the complex data XReal, Xlmag.
Input Parameters
NumPoints : Integer; Number of data points
Inverse: Boolean; FALSE equals forward transform; TRUE equals inverse trans-
form
XReal : TNvectorPtr; Pointer to real values of the data points
XImag : TNvectorPtr; Pointer to imaginary values of the data points
Output Parameters
XReal : TNvectorptr; Pointer to real values of the discrete Fourier transform of the
input data
XImag : TNvectorPtr; Pointer to imaginary values of the discrete Fourier transform
of the input data
Error: Byte; 0: No errors
1: NumPoints < 2
2: NumPoints not a power of two [four]
Comments
The complex Fast Fourier method uses the New/Dispose procedures to manipulate
the heap and should not be used in any program that uses Mark/Release to manip-
ulate the heap.
REALFFT.INC
Description
This example performs a complex Fast Fourier Transform of real data. The Num-
Points real data points are first mapped onto NumPoints/2 complex data points. A
complex Fast Fourier Transform of these complex points is performed by calling
TestInput, Make Sin Cos Table , and FFT. The NumPoints/2 transform is then mapped
onto NumPoints complex points. The real part of the transformation will be even,
and the imaginary part of the transformation will be odd. If you are implementing
this application with a radix-4 algorithm, be sure that the number of real data
points (NumPoints) is twice the power of four.
Input Parameters
NumPoints : Integer; Number of data points
Inverse: Boo1 ean; FALSE equals forward transform; TRUE equals inverse trans-
form
XReal : TNvectorPtr; Pointer to real values of the data points
At least four data points are required, because this algorithm transforms the real
vector to a complex vector half the size. If only two real data points were entered,
the routine would have to take the transform of a single complex point.
Comments
This method uses the New/Dispose procedures to manipulate the heap and should
not be used in any program that uses Mark/Release to manipulate the heap.
COMPCNVL.INC
Description
The calculation of the convolution of two complex vectors is facilitated with a Fast
Fourier Transform routine. The discrete convolution of two functions x and h is
defined by
N-I
Ym = IX
1\ = 0
n hm _ n m = 0,1, ... N - 1
where subscripts are taken modulo N (circular convolution). The basic theorem
that allows us to calculate these effectively using FFTs is as follows:
m = 0,1, ... N - 1
Input Parameters
NumPoints : Integer; Number of data points
XReal : TNvectorPtr; Pointer to real values of the first set of data points
XImag : TNvectorptr; Pointer to imaginary values of the first set of data points
HReal : TNvectorptr; Pointer to real values of the second set of data points
HImag : TNvectorptr; Pointer to imaginary values of the second set of data points
Output Parameters
XReal : TNvectorptr; Pointer to real values of the convolution of XReal, Xlmag and
HReal, Hlmag
XImag : TNvectorptr; Pointer to imaginary values of the convolution of XReal, Xlmag
and HReal, Hlmag
Error: Byte; 0: No errors
1: NumPoints < 2
2: NumPoints not a power of two [four]
Comments
This method uses the New/Dispose procedures to manipulate the heap and should
not be used in any program that uses Mark/Release to manipulate the heap.
Description
The calculation of the convolution of two real vectors is facilitated with a Fast
Fourier Transform routine. This procedure is exactly the same as the previous
procedure for complex convolution (COMPCNVL.lNC) except that only one
Fourier transform need be performed. The procedure is as follows:
1. Given two real vectors XReal and HReal, combine them into a complex vec-
tor XReal + iHReal, where i is the square root of -1.
2. Transform this complex vector.
3. Extract the transforms of the two real functions from the transform of the
complex function (using the symmetry Xf = X_ f , where the bar stands for
complex conjugation).
4. Multiply the resulting transforms point by point.
5. Find the inverse transform of this product using FFTs. REALCNVL.lNC is
about 25 percent faster than its complex counterpart for the same set of real
data.
Input Parameters
NumPoints : Integer; Number of data points
XReal : TNvectorPtr; Pointer to real values of the first set of data points
HReal : TNvectorptr; Pointer to real values of the second set of data points
Output Parameters
XReal : TNvectorPtr; Pointer to real values of the convolution of XReal and HReal
XImag : TNvectorPtr; Pointer to imaginary values of the convolution of XReal and
HReal
Error: Byte; 0: No errors
1: NumPoints < 2
2: NumPoints not a power of two [four]
Comments
This method uses the New/Dispose procedures to manipulate the heap and should
not be used in any program that uses Mark/Release to manipulate the heap.
COMPCORR.INC
Description
The calculation of the correlation of two complex vectors is facilitated with a Fast
Fourier Transform routine. The discrete correlation of two complex functions x and
h is defined by
N-l
Ym = "'xh m+n
L n
m = 0,1, ... N - 1
11 = 0
where subscripts are taken modulo N (circular correlation). The basic theorem that
allows us to calculate these effectively using FFTs is as follows:
ym =X m
H N-m m = 0,1, ... N - 1
where capital letters indicate the transforms of the functions represented by lower-
case letters and - indicates the complex conjugate. (Commonly x and h are real
functions, in which case the preceding formula reduces to the more familiar
expression Grn = Xm H,1l' where the bar stands for complex conjugation. (See REAL
CORR.INC for a real version of correlation.) Thus the procedure for correlation
works like this:
1. Transform both given data sets using FFTs.
2. Multiply each element of the transform of the first data set by the appropri-
ate element of the transform of the second data.
3. Find the inverse transform of this product using FFTs.
Output Parameters
XReal : TNvectorptr; Pointer to real values of the correlation of XReal, Xlmag and
HReal, Hlmag (or the autocorrelation of XReal, Xlmag if
Auto = TRUE)
Xlmag : TNvectorPtr; Pointer to imaginary values of the correlation of XReal, Xlmag
and HReal, Hlmag (or the autocorrelation of XReal, Xlmag if
Auto = TRUE)
Error: Byte; 0: No errors
1: NumPoints < 2
2: NumPoints not a power of two [four]
Comments
If you are performing an autocorrelation of the vector XReal, Xlmag, then set
Auto = TRUE. In this case, the vector HReal, Hlmag will not contain any informa-
tion but must still be passed into the procedure. Autocorrelations are faster to
compute, since only one forward transformation must be made, as opposed to two
for crosscorrelation.
REALCORR.INC
Description
The calculation of the convolution of two real vectors is facilitated with a Fast
Fourier Transform routine. This procedure is exactly the same as the previous
procedure for complex correlation (COMPCORR.INC) except that only one for-
ward Fourier transform need be performed. The procedure is as follows:
1. Given two real vectors XReal and HReal, combine them into a complex
vector XReal + iHReal, where i is the square root of· - 1.
2. Transform this complex vector.
3. Extract the transforms of the two real vectors from the transform of the
complex vector (using the symmetry Xf = X_ f , where the bar stands for
complex conjugation).
4. Multiply each element of the transform of the first data set by the appro-
priate element of the transform of the second data.
5. Find the inverse transform of this product using FFTs.
Input Parameters
NumPoi nts : Integer; Number of data points
Auto: Bool ean; FALSE equals crosscorrelation; TRUE equals autocorrelation
XRea 1 : TNvectorPtr; Pointer to real values of the first set of data points
HRea 1 : TNvectorPtr; Pointer to real values of the second set of data points (for cross-
correlation)
The preceding parameters must satisfy the following conditions:
1. NumPoints 2! 2.
2. NumPoints must be a power of two [four].
Comments
If you are performing an autocorrelation of the vector XReal, then set Auto equal
to TRUE. In this case, the vector HReal will not contain any information but must
still be passed into the procedure. Autocorrelations are faster to compute, since
only one forward transformation must be made, as opposed to two for crosscorrela-
tion. This method uses the New/Dispose procedures to manipulate the heap and
should not be used in any program that uses Mark/Release to manipulate the heap.
Sample Program
Example
Problem. Perform a Fourier transform and an autocorrelation of a 32-point square
wave. Also, convolve and crosscorrelate this square wave with a saw-tooth wave
(assume you are working in Turbo-87).
1. First, make sure that the FFT file FFT87B2.1NC is included in
FFTPROGS.PAS:
(*{$I FFTB2. INC} (* radix 2, 8088 version *)
(*{$I FFTB4.INC} (* radix 4, 8088 version *)
{$I FFT87B2.INC} (* radix 2, 8087 version *)
(*{$I FFT87B4.INC} (* radix 4, 8087 version *)
The input data file SAMPI0A.DAT is as follows (note that this is in real
format):
o
o
o
o
o
o
o
o
Results of Autocorrelation:
1.94454364826301E+OOO -6.89287897017933E-016
1.76776695296637E+OOO -5.49532360539321E-016
1.59099025766973E+OOO -6.55792559089139E-016
1.41421356237309E+OOO -4. 17055809337878E-016
1.23743686707646E+OOO -4.44089209850063E-016
1.06066017177982E+OOO -2.74766180269661E-016
8.83883476483184E-001 -3.84393694788862E-016
7.07106781186547E-001 -1. 17756934401283E-016
5.30330085889910E-001 5.37799391805346E-017
3.53553390593274E-001 -1.96261557335472E-016
1.76776695296637E-001 -1.88132137453390E-016
O.OOOOOOOOOOOOOOE+OOO -9.32242397343491E-017
-3.92523114670944E-016 4.44089209850063E-016
2.35513868802566E-016 O.OOOOOOOOOOOOOOE+OOO
4.71027737605132E-016 -1.06260198549818E-016
4.71027737605132E-016 -1.66822323735151E-016
1.57009245868377E-016 5.81728018656864E-016
-3.14018491736755E-016 3.92523114670944E-016
O.OOOOOOOOOOOOOOE+OOO 6.00281407857881E-016
-1.57009245868377E-016 2.89485797069821E-016
-3.92523114670944E-016 4.44089209850063E-016
-3. 14018491736755E-016 1. 17756934401283E-016
1.76776695296636E-001 4.39904846020120E-016
3.53553390593273E-001 -1.96261557335472E-017
5.30330085889910E-001 5.37799391805346E-017
7.07106781186547E-001 3.53270803203849E-016
8.83883476483184E-001 2.43643288684648E-016
1.06066017177982E+OOO 2.20794252002406E-016
1.23743686707646E+OOO -4. 44089209850063E-016
1.41421356237309E+OOO 1.57009245868377E-016
1.59099025766973E+OOO 5.07490473185598E-017
1.76776695296637E+OOO 3.04205413869981E-016
Keeping in mind that this is a periodic function (see "Data Sampling"), you
can see that this is a triangle wave.
Let's now convolve the square wave with a saw-tooth wave. The input file for
the saw-tooth wave (SAMPIOC.DAT) is as follows:
o
o
o
o
o
o
The programs in this chapter can only be run by those users with PC-DOS.
There are some programs that graphically demonstrate the usefulness of the
least-squares routines in Chapter 9 and the Fourier transforms in Chapter 10. A
graphics monitor is required. Each program reads a data set from an input file, and
uses this Toolbox to display the results. You will see curves being fitted to data
using the least-squares routines and also see a signal being transformed into its
Fourier spectrum.
The programs LSQIBM.COM and FFTIBM.COM graphically illustrate the
power and utility of the Turbo Pascal Numerical Methods Toolbox. To run them,
you'll need an IBM Color Graphics Adapter (CGA) or a suitable clone. (The pro-
grams LSQHERC.COM AND FFTHERC.COM can be used to graphically illus-
trate the Toolbox on a Hercules Monochrome Graphics Adapter or compatible.)
And to print, you'll need an Epson or IBM compatible dot-matrix printer. As
explained in this chapter, these programs can be recompiled to run on other hard-
ware, including the IBM Enhanced Graphics Adapter (EGA) in its high resolution
mode. The programs can also be recompiled to take advantage of the 8087 (or
80287) math coprocessor.
261
Function of the Least-Squares Graphics Demonstration
Program
The input data from LSQIN.DAT along with the least-squares fits and coeffi-
cients are output to a text file called LSQOUT.DAT.
A default output data file can easily be arranged by changing a constant
WriteToFile in LSQDEMO.PAS and recompiling it. (See the section, "Rebuilding
LSQIBM.COM; to recompile it, and the comment in LSQDEMO.PAS next to the
constant WriteToFile.)
Printing
'.11
Rebuilding LSQIBM.COM
The sources to LSQIBM.COM are provided. To recompile, you need Turbo Pascal
(version 3.0) and Turbo Graphix Toolbox (version 1.06A or later), as well as this
Toolbox.
To rebuild LSQIBM.COM, the following files are needed:
From CHAP9:
EXP.LSQ (for the exponential model)
FOURIER.LSQ (Fourier series)
LOG.LSQ (logarithmic)
POLY.LSQ (polynomial)
POWER.LSQ (power law)
From CHAPll:
GENERIC.LSQ
IOCHECK.lNC (a modification of COMMON.lNC)
LEAST.MOD (a modification of LEAST.lNC from Chapter 9)
LSQDEMO.PAS
From Turbo Pascal:
(These files are not included in this Toolbox.)
TURBO.COM or TURBO-87.COM (the Turbo Pascal compiler)
TURBO.MSG (compiler error messages)
Rebuilding FFTIBM.COM
The sources to FFTIBM.COM are provided. To recompile, you need Turbo Pascal
(version 3.0) and Turbo Graphix Toolbox (version 1.06A or later), as well as this
Toolbox.
To rebuild FFTIBM.COM, the following files are needed:
From CHAPIO:
COMPFFT.lNC
FFTB2.1NC
REALFFT.lNC
From CHAPll:
4X6.FON (originally from the Graphix Toolbox)
8X8.FON (originally from the Graphix Toolbox)
ERROR.MSG (originally from the Graphix Toolbox)
FFTDEMO.PAS
IOCHECK.lNC (a modification of COMMON.INC)
From Turbo Pascal:
(These files are not included in this Toolbox.)
TURBO.COM or TURBO-87.COM (the Turbo Pascal compiler)
TURBO.MSG (compiler error messages)
To recompile for the IBM Enhanced Graphics Adapter (EGA), you must copy
GRAPHIX.EGA to GRAPHIX.SYS. A copy of GRAPHIX.EGA is included with
this Toolbox since many users purchased copies of the Turbo Graphix Toolbox
before support for the EGA was added.
LSQIBM.COM and FFTIBM.COM will run on machines with an EGA card.
However, to take advantage of the higher resolution offerred by this card, you must
copy GRAPHIX.EGA and recompile.
To recompile to take advantage of the math coprocessor, you must use TURBO-
87.COM instead of TURBO.COM. Some increased performance in the FFr demo
program will be obtained ifFFfB2.1NC is replaced by FFT87B2.INC from Chap-
ter 10 of the Toolbox.
271
Gerald, Curtis F., and Patrick O. Wheatley. Applied Numerical Analysis, 3rd ed.
Reading: Addison-Wesley Publishing Co., 1984. This is another excellent source for
learning numerical analysis.
Press, William H., Brian P. Flannery, Saul A. Teukolsky, and William T. Vetter-
ling. Numerical Recipes: The Art of Scientific Computing. New York: Cambridge
University Press, 1986. This book is user-oriented, discusses many of the subtle
problems encountered when implementing numerical methods, and has program
listings in Turbo Pascal.
Ralston, Anthony, and Philip Rabinowitz. A First Course in Numerical AnalysiS.
New York: McGraw-Hill Book Co., 1978. A well-written mathematics text that is a
step more sophisticated than the preceding ones.
273
$FF error, 145, 190,202 Cubic spline methods
clamped, 44, 57-62
A free, 44,52-56, 64,76-79
ADAMS_l.INC, 156, 168-170 Cyclic Jacobi method, 132, 149-153
ADAMS_l.PAS, 170-171
Adams-Bashforth/Adams-Moulton D
predictor-corrector method, 156, Data sampling, 239-240
168-171 Data types, 12
Adams-Bashforth formula, 168 Defined constants, 12
Adams-Moulton formula, 168 Deflation, 16
ADAPGAUS.INC, 88, 98-100 and Laguerre, 37-41
ADAPGAUS.PAS, 101 and Newton-Horner, 28-32
ADAPSIMP.INC, 88, 95-96 of a matrix, 132, 143-148
ADAPSIMP.PAS, 96-97 DERIV2FN.lNC, 64-65, 83-84
Adaptive schemes, 87 DERIV2FN.PAS, 84-86
quadrature, 95-101 DERIV2.1NC, 64, 71-73
Aliasing, 239 DERIV2.PAS, 73-75
Derivative, 63
B approximation of, 64-86
Backward substitution, 114, 120 DERIVFN .INC, 64, 80-81
Basis vectors, 222 DERIVFN .PAS, 81-82
BISECT.INC, 18-19 DERIV.INC, 64, 66-68
Bisection method, 15 DERIV.PAS, 68-70
root of a function using, 18-20 Determinant of a matrix, 105-109
BISECT.PAS, 19-20 DET.INC, 106-108, 128
Boundary value problems, 155-156, 158 DET.PAS, 108-109
using Linear Shooting/Runge-Kutta Diagonally dominant, 126
methods, 215-219 Diagonal matrix, 131
using Shooting/Runge-Kutta methods, Differential equations
208-214 first-order, 159-171
coupled,186-195
C linear, 155
nth order, 157, 178-185
Chebyshev polynomials, 224 ordinary, 155
Color Graphics Adapter, 3, 12 second-order, 158, 172-177, 192,
for graphics demos, 261, 267 208-219
COMMON.INC,7, 11, 13 coupled, 196-207
COMPCNVL.INC, 233, 235, 237-238 stiff, 161
application, 246-247 . systems of, 157
COMPCORR.INC, 233, 235, 238 Direct factorization of matrices, 106,
application, 249-250 120-125
COMPFFT.INC, 233, 235, 237 DIRFACT.lNC. 106, 120-122
application, 244-245 DIRFACT.PAS, 122-125
Compiler directives, 13 Distribution disks, 7-11
Convergence, rate of, 15-16 DIVDIF.lNC, 49-50
Cooley-Tukey method, 233-234 DIVDIF.PAS, 50-51
Coprocessor, 3, 270
CUBE_CLA.lNC, 57-58 E
CUBE_CLAPAS, 59-62
CUBEJRE.lNC, 52-53 Eigensystem, 149
CUBEJRE.PAS, 53-56 Eigenvalue, 131-132
Eigenvector, 131
Index 275
Laguerre's method, 16 Mesh points, 156
finding roots of complex polynomial, MULLERINC, 33-35
37-41 MULLERPAS, 35-36
LEAST.IN C, 106, 222-226 Muller's method, 16
LEAST. MOD, 267 finding roots of complex function,
LEAST.PAS, 227-231 33-36
Least-squares approximation, 222-231
Least-squares solution N
graphics demo, 262-263 New/Dispose, 96, 99, 145
linear regression, 221 in Adams-Bashforth/Adams-Moulton,
multiple regression, 221 170
Linear equations, 105-106 in complex fast Fourier, 245-252
differential, 155 in Linear Shooting/Runge-Kutta, 217
with direct factoring, l20-125 in Runge-Kutta, 161, 182, 191,202,211
with Gaussian elimination, 114-119 in Runge-Kutta-Fehlberg, 165
LINSHOT2.1NC, 156, 158, 215-218 NEWTDEFL.INC, 28-30
LINSHOT2.PAS, 218-219 NEWTDEFL.PAS, 30-32
Lipshitz condition, 157 Newton-Horner method, 15-16
LOG.LSQ, 226 with deflation, 28-32
LSQDEMO.COM, 268 Newton-Raphson method, 15-16
LSQDEMO.PAS, 263 root of a function using, 21-24
rebuilding for EGA, 269 Newton's general divided-difference
rebuilding for printing, 265-267 algorithm, 43, 49-51
LSQHERC.COM, 261 Nonlinear shooting method, 158,
rebuilding for Hercules, 269 208-214
LSQIBM.COM, 261-263 Numerical differentiation, 63-65
graphics demo, 265-267 five-point formulas, 64, 66-75
rebuilding, 267 three-point formulas, 64, 66-75
rebuilding for EGA, 269 two-point formulas, 64, 66-70
LUJ)ecompose, l20-l21 Nyquist frequency, 239, 264
LU_Solve, l20-l22
p
M
Partial pivoting, 106, 117
MakeSinCosTable, 242-243 and direct factoring, l20
Mark/Release, 96, 99, 145 PARTPIVT.INC, 106, 117-118
in Adams-Bashforth/Adams-Moulton, PARTPIVT. PAS , 118-119
170 POLY.LSQ, 224, 267
in complex fast Fourier, 245-252 Polynomials
in Linear Shooting/Runge-Kutta, 217 Lagrange, 89
in Runge-Kutta, 161, 182, 191, 202, 211 Legendre, 98-100
in Runge-Kutta-Fehlberg, 165 methods to approximate roots of, 16,
Matrix 25-41
algebra, 105 POWERINC, 131-135
diagonal, 131 POWERLSQ, 225, 267
direct factorization, 106, l20 Power method, 131-136
identity, 139 and Wielandt's deflation, 143-148
nonsingular, 106, l20-l25 POWERPAS, 135-136
orthogonal, 134 Powers-of-four, 234
permutation, l20 Powers-of-two, 234
rotation, 149
square, 106, 131-133, 137
symmetric, 132, 149-153
Index 277
278 Turbo Numerical Methods Toolbox
Borland
Software
=i 3:-0
¥J'
BORLAND
INTERNATIONAL 4585 Scotts Valley Drive, Scotts Valley, CA 95066
All the SideKick windows stacked up over Lotus 1-2-3.- Here's SideKick running over Lotus 1-2-3. In the
From bottom to top: SideKick's "Menu Window," ASCII SideKick Notepad you'll notice data that's been imported
Table, Notepad, Calculator, Appointment Calendar, Monthly directly from the Lotus screen. In the upper right you can
Calendar, and Phone Dialer. see the Calculator.
Minimum system configuration: IBM PC. Xl. AT. PCjr. and true compatibles: PC-DOS (MS-DOS)
2.0 or greater. 128K RAM. One disk drive.
SuperKey is a registered trademark 01 Borland Internationat. Inc. IBM. XT, AT. and PCjr are
registered trademarks 01 International Business Machines Corp. MS-DOS is a registered
trademark 01 Microsoft Corp. BOR 0062C
If you use an IBM® PC, you need
T U R B 0
Lightning®
Turbo Lightning teams up II you ever write a You can teach Turbo
with the Random House word, think a word, or Lightning new words
Concise Word List to say a word, you need You can teach your new Turbo
check your spelling as Turbo Lightning Lightning your name, business
you type! associates' names, street
Turbo Lightning, using the names, addresses, correct
BO,OOO-word Random House capitalizations, and any
Dictionary, checks your spelling specialized words you use
as you type. If you misspell a frequently. Teach Turbo
word, it alerts you with a Lightning once, and it
"beep." At the touch of a key, knows forever.
Turbo Lightning opens a
window on top of your Turbo Lightning is the
application program and engine that powers
suggests the correct spelling.
Just press one key and the
"
. Borland's Turbo Lightning
Libraryf'J
misspelled word is instantly Turbo Lightning brings
replaced with the correct word. electronic power to the
The Turbo Lightning Proofreader
Random House Concise Word
Turbo Lightning works List and Random House
hand-in-hand with the Thesaurus. They're at your
Random House Thesaurus fingertips-even while you're
to give you instant access running other programs. Turbo
to synonyms Lightning will also "drive"
soon-to-be-released
Turbo Lightning lets you
encyclopedias, extended
choose just the right word from
thesauruses, specialized
a list of alternates, so you
dictionaries, and many other
don't say the same thing the
popular reference works. You
same way every time. Once
get a head start with this
Turbo Lightning opens the
first volume in the Turbo
Thesaurus window, you see a Lightning Library.
list of alternate words; select
the word you want, press
ENTER and your new word will The Turbo Lightning Thesaurus
instantly replace the original
word. Pure magic!
BORLAND
INTERNATIONAL
Turbo Lightning and Turbo Lightning Library are registered trademarks of Bortand tnternational, Inc.
IBM. XT, AT, and PCjr are registered trademarks of International Business Machines Corp, Random
House is a registered trademark of Random House, Inc, Copyright 1987 Borland International
BOR 0070B
Your Development Toolbox and Technical Reference Manual for Thrbo lightning@)
l I G H T N I N G
"The real future of Lightning clearly lies not with the spelling checker and thesaurus currently
included, but with other uses of its powerful look-up engine." Ted Silveira, Profiles
"This newest product from Borland has it all." Don Roy, Computing Now!
Minimum system configuration: IBM PC. XT. AT. PCjr. Portable. and true compatibles. 256K RAM minimum. PC·DOS (MS·DOS) 2.0
or greater. Turbo Lightning software required. Optional-Turbo Pascal 3.0 or greater to edit and compile Turbo Pascal source code.
Turbo Pascal, Turbo Lightning and Turbo Lightning Library are registered trademarks and Lightning Word Wizard is a trademark of Borland International, Inc. Random
House is a registered trademark of Random House, Inc. IBM, XT, AT, and PCjr are registered trademarks of International Business Machines Corp. MS-DOS is a
registered trademark of Microsoft Corp. CqJyright 1987 Borland International BOR00878
'EEI Ell
I TIE DATABASE
~r'~II: .ANABEI
The high-performance database manager
®
Minimum system configuration: IBM PC, Xl, AT, and true compatibles. 384K RAM minimum. IBM Color Graphics Adapter, Hercules
Monochrome Graphics CArd, or equivalent. PC·DOS (MS· DOS) 2.0 or greater. Hard disk and mouse optional. Lotus 1·2·3, dBASE,
or PFS: File optional.
Reflex is a trademark of Bortand/Analytica tnc. Lotus 1-2-3 is a registered trademark of Lotus
Development Corporation. dBASE is a registered trademark of Ashton-Tate. PFS: Fite is a
registered trademark of Software Publishing Corporation. IBM. XT. AT, and IBM Color Graphics
Adapter are registered trademarks of International BUSiness Machines Corporation. Hercules
Graphics Card is a trademark of Hercules Computer Technology. MS-DOS is a registered
trademark of Microsoft Corp. Copyright 1987 Borland International BOA 0066C
REILEX: TIE WIIISII"·
Includes 22 "instant templates" covering a broad range of
business applications (listed below). Also shows you how to
customize databases, graphs, crosstabs, and reports. It's an invaluable
analytical tool and an important addition to another one of
our best sellers, Reflex: The Database Manager.
Fast-start tutorial examples:
Learn Reflex@) as you work with practical business applications. The Reflex Workshop Disk supplies
databases and reports large enough to illustrate the power and variety of Reflex features. Instructions in each
Reflex Workshop chapter take you through a step-by-step analysis of sample data. You then follow simple
steps to adapt the files to your own needs.
22 practical business applications:
Workshop's 22 "instant templates" give you a wide range of analytical tools:
Administration • Tracking Manufacturing Quality Assurance
• Scheduling Appointments • Analyzing Product Costs
• Planning Conference Facilities Accounting and Financial Planning
• Managing a Project • Tracking Petty Cash
• Creating a Mailing System • Entering Purchase Orders
• Managing Employment Applications • Organizing Outgoing Purchase Orders
Sales and Marketing • Analyzing Accounts Receivable
• Researching Store Check Inventory • Maintaining Letters of Credit
• Tracking Sales Leads • Reporting Business Expenses
• Summarizing Sales Trends • Managing Debits and Credits
• Analyzing Trends • Examining Leased Inventory Trends
• Tracking Fixed Assets
Production and Operations • Planning Commercial Real Estate Investment
• Summarizing Repair Turnaround
Whether you're a newcomer learning Reflex basics or an experienced "power user" looking for tips, Reflex:
The Workshop will help you quickly become an expert database analyst.
Minimum system configuration: IBM PC, AT, and XT, and true compatibles. PC-DOS (MS-DOS) 2.0 or greater. 384K RAM minimum. Requires Rellex:
The Database Manager, and IBM Color Graphics Adapter, Hercules Monochrome Graphics Card or equivalent.
Reflex is a registered trademark and Reflex: The Workshop is a trademark of Borland/Analytica, Inc. IBM, AT, and Xl are registered trademarks 01 International Business
Machines Corp. Hercules is a trademark 01 Hercules Computer Technology. MS-DOS is a registered trademark 01 Microsoft Corp. Copyright 1987 Borland International
BOA 0088B
Version 3.0 with 8087 support and BCD reals
MicroCalc: A sample spreadsheet on your disk "What I think the computer industry is headed
with ready-to-compile source code. for: well-documented, standard, plenty of
good features, and a reasonable price."
IBM~ PC Version: Supports Turtle Graphics, -Jerry Pournelle, BYTE
color, sound, full tree directories, window
routines, input/output redirection, and
much more.
LOOK AT TURBO NOW!
@' More than 500,000 users worldwide. @' Turbo Pascal named "Most
Significant Product of the Year" by
@' Turbo Pascal is the de facto industry PC WEEK.
standard.
@' Turbo Pascal 3.0-the fastest Pascal
@' Turbo Pascal wins PC MAGAZINE'S development environment on the
award for technical excellence. planet, period.
Suggested Retail Price: $99.95; CP/M~-80 version without 8081 and BCD: $69.95
Features lor 16-bit Systems: 8087 math co-processor support for intensive calculations.
Binary Coded Decimals (BCD): eliminates round-off error! Amust for any serious business application.
Minimum system configuration: 128K RAM minimum. Includes 8087 &BCD features for 16-bit MS-DOS 2.0 or later and
CP/M-86 1.1 or later. CP/M-80 version 2.2 or later 48K RAM minimum (8087 and BCD features not available). 8087
version requires 8087 or 80287 co-processor.
Turbo Pascal is a registered trademark of Bor1and International. Inc. CP/M is a registered trademark
of Digital Research Inc. IBM is a registered tradermrk of International Business Machines Corp.
MS-DOS is a registered tradermrk of Micr~ft Corp. VlbrdStar is a registered trademark of
MicroPro International. Copyright 1987 Borland International BOA 00619
VERSION 2.0
And if you're already proficient, Turbo Tutor can sharpen up the fine paints.
The manual and program disk focus on the whole spectrum of Turbo
Pascal programming techniques.
• For the Novice: It gives you a concise history of Pascal, tells you how to write a
simple program, and defines the basic programming terms you need to know.
• Programmer's Guide: The heart of Turbo Pascal. The manual covers the fine points
of every aspect of Turbo Pascal programming: program structure, data types, control
structures, procedures and functions, scalar types, arrays, strings, pointers, sets, files,
and records.
• Advanced Concepts: If you're an expert, you'll love the sections detailing such topics as
linked lists, trees, and graphs. You'll also find sample program examples for PC-DOS and
MS-DOS.@)
10,000 lines of commented source code, demonstrations of 20 Turbo Pascal features, multiple-
choice quizzes, an interactive on-line tutor, and more!
Turbo Tutor may be the only reference work about Pascal and programming you'll ever need!
Minimum system configuration: Turbo Pascal 3.0. PC-DOS (MS-DOS) 2.0 or later. 192K RAM minimum (CP/M-80
version 2.2 or later: 64K RAM minimum).
Turbo Pascal and Turbo Tutor are regislered tademarks of Borland International Inc. CP1M is a
registered trademark of Digital Research Inc. MS-DOS is a registered trademark of Microsoft Corp.
Copyright 1987 Borland International BOR 0064C
Is The Perfect Complement To Turbo Pascal®
It contains a complete library of Pascal procedures that
allows you to sort and search your data and build powerful database
applications. It's another set of tools from Borland that will give
even the beginning programmer the expert's edge.
Minimum system conliguration: 128K RAM and one disk drive (CP/M·80: 48K). 16·bit systems: Turbo Pascal 2.0 or greater lor
MS·DOS or PC·DOS 2.0 or greater. Turbo Pascal 2.1 or greater lor CP/M·86 1.0 or greater. 8·bit systems: Turbo Pascal 2.0 or
greater lor CP/M·80 2.2 or greater.
BORLAND
INTERNATIONAL
Turbo Pascal and Turbo Database Toolbox are registered trademarks of Borland International
Inc. CP1M is a registered trademark of Dignal Research. Inc. MS-DOS is a registered
trademark of Microsoft Corp. Copyright 1987 Borland International BOR 00630
7VRBO MS'CAl.
GRAPHlx200tsox®
A Library of Graphics Routines for Use with Turbo Pasca/®
High-resolution graphics for your IBM" PC, AT," XT," PCjr", true PC compatibles, and the Heath
Zenith Z-100:" Comes complete with graphics window management.
Even if you're new to Turbo Pascal programming, the Turbo Pascal Graphix Toolbox will get you started
right away. It's a collection of tools that will get you right into the fascinating world of high-resolution
business graphics, including graphics window management. You get immediate, satisfying results. And
we keep Royalty out of American business because you don't pay any-even if you distribute your own
compiled programs that include all or part of the Turbo Pascal Graphix Toolbox procedures.
"While most people only talk about low-cost personal computer software, Borland has been doing
something about it. And Borland provides good technical support as part of the price."
John Markov & Paul Freiberger, syndicated columnists.
II you ever plan to create Turbo Pascal programs that make use of business graphics or scientific
graphics, you need the Turbo Pascal Graphix Toolbox.
Suggested Retail Price: $69.95 (not copy protected)
Minimum system configuration: IBM PC, Xl, AT, PCjr, true compatibles and the Heath Zenith Z-100. Turbo Pascal 3.0 or later. 192K
RAM minimum. Two disk drives and an IBM Color Graphics Adapter (CGA), IBM Enhanced Graphics Adapter (EGA), Hercules Graphics
Card or compatible.
BORLAND
INTERNATIONAL
Turbo Pascal and Turbo Graphix Toolbox are registered trademarks 01 Borland International.
Inc. IBM. XT. AT. and PCjr are registered trademarks 01 International Business Machines
Corporation. Hercules Graphics Card is a trademark 01 Hercules Computer Technology. Heath
Zenith Z-1 00 is a trademark 01 Zenith Data Systems Epson is a registered trademark 01
Epson Corp. Copyright 1987 Borland Intemational BOR 0068C
2VRBO PASCAl.
.EDJ2OR 2totsox®
It's All You Need To Build Your Own Text Editor
Or Word Processor
Build your own lightning-fast editor and incor- Create your own word processor. We provide all
porate it into your Turbo Pascale programs. the editing routines. You plug in the features you want.
Turbo Editor Toolbox gives you easy-to-install You could build aWordStar e -like editor with pull-down
modules. Now you can integrate a fast and powerful menus like Microsoft'se Word, and make it work as fast
editor into your own programs. You get the source as WordPerfect.e
code, the manual, and the know-how.
To demonstrate the tremendous power of Turbo Editor Toolbox, we give you the source code for
two sample editors:
Simple Editor A complete editor ready to include in your programs. With windows, block commands, and
memory-mapped screen routines.
MicroStar A full-blown text editor with a complete pull-down menu user interface, plus a lot more.
Modify MicroStar's pull-down menu system and include it in your Turbo Pascal programs.
And Turbo Editor Toolbox has features that word processors selling for several hundred dollars can't begin to match.
Just to name a few:
[3" RAM-based editor. You can edit very large [3" Multiple windows. See and edit up to eight
files and yet editing is lightning fast. documents-or up to eight parts of the same
[3" Memory-mapped screen routines. In- document-all at the same time.
stant paging, scrolling, and text display. [3" Multitasking. Automatically save your
[3" Keyboard installation. Change control text. Plug in a digital clock, an appointment
keys from WordStar-like commands to any that alarm-see how it's done with MicroStar's
you prefer. "background" printing.
Best of all, source code is included for everything in the Editor Toolbox.
Suggested Retail Price: $69.95 (not copy protected)
Minimum system configuration: IBM PC, XT, AT, 3270, PCjr, and true compatibles. PC-DOS (MS-DOS) 2.0 or greater. 192K RAM.
You must be using Turbo Pascal 3.0 for IBM and compatibles.
Turbo Pascal and Turbo Editor Toolbox are registered trademarks of Borland International. Inc.
WordStar is a registered trademark of MicroPro International Corp. Word and MS-DOS are
registered trademarks of Microsoft Corp. WordPerfect is a trademark of Satellite Software
International. IBM, XT, AT, and PCjr are registered trademarks of International Business Machines
Corp. BOR 0067B
®
TURBO BRIDGE
Now play the world's most popular card game-bridge. Play one-on-one with your computer or against up to
three other opponents. With Turbo Pascal source code, you can even program your own bidding or scoring
conventions.
"There has never been a bridge program written which plays at the expert level, and the ambitious user will
enjoy tackling that challenge, with the format already structured in the program. And for the inexperienced player,
the bridge program provides an easy-to-follow format that allows the user to start right out playing. The user can
'play bridge' against real competition without having to gather three other people."
-Kit Woolsey, writer of several articles and books on bridge,
and twice champion of the Blue Ribbon Pairs.
TURBO GO-MOKU
Prepare for battle when you challenge your computer to a game of Go-Moku-the exciting strategy game also
known as Pente." In this battle of wits, you and the computer take turns placing X's and D's on a grid of 19X19
squares until five pieces are lined up in a row. Vary the game if you like, using the source code available on your
disk.
Suggested Retail Price: $69.95 (not copy protected)
Minimum system configuration: IBM PC, XT, AT, Portable, 3270, PClr, and true compatibles. PC·DOS (MS·DOS) 2.0 or later. 192K
RAM minimum. To edit and compile the Turbo Pascal source code, you must be using Turbo Pascal 3.0 for IBM PCs and
compatibles.
~ " BORLAND Turbo Pascal and Turbo GameWorks are registered trademarks of Borland International. Inc.
Pente is a registered trademark of Parker Brothers. IBM. XT. AT. and PCjr are registered
'# I N T ERN A T ION A L trademarks of International Business Machiles Corporation. MS-DOS is a registered trademark
of Microsoft Corporation. Copyright 1987 Borland International BOR0065C
TURBO
The Turbo Prolog Toolbox enhances Turbo Prolog-our 5th-generation computer programming
language that brings supercomputer power to your IBM PC and compatibles-with its more than 80
tools and over 8,000 lines of source code that can be incorporated into your programs, quite easily.
Minimum system configuration: iBM PC, XT, AT or true compatibles. PC-DOS (MS-DOS) 2.0 or later. Requires Turbo Prolog 1.10
or higher. Dual-floppy disk drive or hard disk. 512K.
Turbo Prolog Toolbox and Turbo Prolog are trademarks of Borland International, Inc. Reflex
is a registered trademark of Borland/Analytica, Inc. dBASE III is a registered trademark of
Ashton-Tate. Lotus 1-2-3 and Symphony are registered trademarks of Lotus Development
Corp. IBM, XT. and AT are registered trademarks of International Business Machines Corp.
MS-DOS is a registered trademark of Microsoft Corp. BOR 0240
1UIIII IABIC®
The high-speed BASIC you've been waiting for!
You probably know us for our Turbo Pascale and Turbo Pr%g.'" Well, we've done
it again! We've created Turbo Basic, because BASIC doesn't have to be slow.
If BASIC taught you how to walk, Turbo Basic will teach you how to run!
With Turbo Basic, your only speed is "Full Speed Ahead"! Turbo Basic is a complete development
environment with a lightning fast compiler, an interactive editor and a trace debugging system. And
because Turbo Basic is also compatible with BASICA, chances are that you already know how to use
Turbo Basic. .
Turbo Basic ends the basic confusion
There's now one standard: Turbo Basic. And because Turbo Basic is a Borland product, the price is
right, the quality is there, and the power is at your fingertips. Turbo Basic is part of the fast-growing
Borland family of programming languages we call the "Turbo Family." And hundreds of thousands of
users are already using Borland's languages. So, welcome to a whole new generation of smart PC
users!
Free spreadsheet included with source code!
Yes, we've included MicroCalc, our sample spreadsheet, complete with source code. So you can get
started right away with a "real program." You can compile and run it "as is," or modify it.
Minimum system configuration: IBM PC. AT. XT or true compatibles. 256K. One floppy drive. PC-DOS (MS-DOS) 2.0 or lalel.
Turbo Basic and Turbo Pascal are registered trademarks and Turbo Prolog is a trademark of
Borland International. Inc. IBM. AT. and XT are regislered trademarks of International Business
Machines Corp. MS-DOS is a registered trademark of Microsolt Corp.
Copyright 1987 Borland International BOR 0265A
Includes tree
MicroCalc spreadsheet
with source code
Suggested Retail Price: $99.95* (not copy protected) ·Introductory offer good through July 1. 1987.
Minimum system configuration: IBM PC, XT, AT and true compatibles. PC-DOS (MS-DOS) 2.0 or later. One floppy drive. 320K.
Turbo C and Turbo Pascal are registered traoomarks and Turbo Prolog is a trademark 01 Borland
International. Inc. Microsoft C and MS-DOS ere registered trademarks 01 Microsoft Corp. lattice C
is a registered trademark 01 lattice. Inc. IBM, Xl, and AT are registered trademarks 01 International
Business Machines Corp. BOA 0243
EIIIEIA: "E BllfJlE'u
The solution to your most complex
equations-in seconds!
If you're a scientist, engineer, financial analyst, student, teacher, or any other professional working with
equations, Eureka: The Solver can do your Algebra, Trigonometry and Calculus problems in a snap.
Eureka also handles maximization and minimization problems, plots functions, generates reports, and
saves an incredible amount of time. Even if you're not a computer specialist, Eureka can help you
solve your real-world mathematical problems fast, without having to learn numerical approximation
techniques. Using Borland's famous pull-down menu design and context-sensitive help screens, Eureka
is easy to learn and easy to use-as simple as a hand-held calculator.
X + exp(X) = 10 solved instantly instead 01 eventually!
Imagine you have to "solve for X," where X + exp(X) = 10, and you don't have Eureka: The Solver.
What you do have is a problem, because it's going to take a lot of time guessing at "X." With Eureka,
there's no guessing, no dancing in the dark-you get the right answer, right now. (Ps: X = 2.0705799,
and Eureka solved that one in .4 of a second!)
How to use Eureka: The Solver
It's easy. You can then tell Eureka to
1. Enter your equation into the • Evaluate your solution
full-screen editor • Plot a graph
2. Select the "Solve" command • Generate a report, then send the output
3. Look at the answer to your printer, disk file or screen
4. You're done , • Or all of the above
Minimum system configuration: IBM PC, AT, XT, Portable, Suggested Retail Price: $99.95*
3270 and true compatibles. PC-DOS (MS-DOS) 2.0 and (not copy protected)
later. 384K.
Eureka: The Solver is a trademark of Borla1d International, Inc. IBM, AT, and XT are registered
trademarks of International Business Machnes Corp. MS-DOS is a registered trademark of
Microsoft Corp. Copyright 1987 Borland International BOR 0221A
'Introductory price expires July 1, 1987
"II'EII'PI®
IJJ Ift'lIJ : IIRSAIIIIER
THE BEBKTllfI
Release 2.0
Macintosh™
The most complete and comprehensive collection of
desk accessories available for your Macintosh!
Thousands of users already know that SideKick is the best collection of desk accessories available
for the Macintosh. With our new Release 2.0, the best just got better.
We've just added two powerful high-performance tools to SideKick-Outlook": The Outliner
and MacPlan'": The Spreadsheet. They work in perfect harmony with each other and while you
run other programs!
Outlook: The Outliner
• II's the desk accessory with more power than a stand-alone outliner
• A great desktop publishing tool, Outlook lets you incorporate both text and graphics
into your outlines
• Works hand-in-hand with MacPlan
• Allows you to work on several outlines at the same time
MacPlan: The Spreadsheet
• Integrates spreadsheets and graphs
• Does both formulas and straight numbers
• Graph types include bar charts, stacked bar charts, pie charts and line graphs
• Includes 12 example templates free!
• Pastes graphics and data right into Outlook creating professional memos and reports, complete
with headers and footers.
~ Calendar III
•
'"
em ExplIf'Is.c
~
[] 0319iLJ,bor
Phone Log []I
~
<4.E.6S1 H.-.... iII's
6.,,1a Ovtrtlt"
~ Analog clock 0" 18'£ Tot.,ExpfOSK
~ Alarm system IJ
•
'"
1845'1: ".tPro(11
~ Calculator
~ Report generator MacPlan does both spreadsheets and business
~ Telecommunications (new version now graphs. Paste them into your Outlook files and
supports XModem file transfer protocol) generate professional reports.
SideKick is a registered trademark and Outlook and MacPlan are trademarks of Borland
International, Inc. Macintosh is a trademark of Mcintosh Laboratory, Inc. licensed to Apple
Computer, Inc. Copyright 1987 Borland International BOR 00690
'I £IIEI;• IAIIASEI
~.]
111 IIA1AIA"
The easy-to-use relational database that thinks like a spreadsheet.
Rellex lor the Mac lets you crunch numbers by entering lormulas and link
databases by drawing on-screen lines.
5 free ready-to-use templates are included on the examples disk:
• A sample 1040 tax application with Sched- • A client billing application set up for a law
ule A, Schedule 8, and Schedule 0, each office, but easily customized by any
contained in a separate report document. professional who bills time.
• A portfolio analysis application with linked • A parts explosion application that breaks
databases of stock purchases, sales, and down an object into its component parts
dividend payments. for cost analysis.
• A checkbook application.
Reflex for the Mac accomplishes all of these tasks without programming-using
spreadsheet-like formulas. Some other Reflex for the Mac features are:
• Visual database design. • Data types which include variable length text, number,
• "What you see is what you get" report and form layout integer, automatically incremented sequence number,
with pictures. date, time, and logical.
• Automatic restructuring of database files when data • Up to 255 fields per record.
types are changed, or fields are added and deleted. • Up to 16 files simultaneously open.
• Display formats which include General, Decimal, • Up to 16 Mac fonts and styles are selectable for
Scientific, Dollars, Percent. individual fields and labels.
• file Edit form.t OllCrlb. D,"'''' S•• n .. Milt Window
OltlB . . IOUflr"llw
<om""'''I.''''
Compon.ntOfl.llg IO_ltnembl
~:;: ___
.1
MII_ .. p,rt
ISP.~t I.to.
"Peru.S
After opening the "Overview" window, you The link lines you draw establish both You can have multiple windows open
draw link lines between databases directly visual and electronic relationships between simultaneously to view all members 01 a
onto your Macintosh screen. your databases. linked set-which are interactive and truly
relational.
Critic's Choice
"... a powerful relational database ... uses a visual approach to information management." InfoWorld
"... gives you a lot of freedom in report design; you can even import graphics." A+ Magazine
". '.' bridges the gap between the pretty programs and the power programs." Stewart Alsop, PC Letter
Minimum system configuration: Macintosh 512K or Macintosh Plus with one disk drive.
Turbo Pascal and SideKick are registered trademarks of Borland International, Inc. and Reflex is a
registered trademark of Borlandl Analytica, Inc. Macintosh is a trademark 01 Mcintosh Laboratories, Inc. licensed
to Apple Computer with its express permission. Lisa is a registered trademark of Apple Computer, Inc. Inside
Macintosh is a copyright of Apple Computer, Inc.
Copyright 1987 Borland International BOR 0167A
-
Borland
Software
OllD.Bll fODAY I
_________ ~ I
j
j
==:115,~O
-: L
BORLAND INTERNATIONAL
IIIIi
J 4585 Scotts Valley Drive Scotts Valley, California 95066 I'
J In Ii
I To Orde~ ,.- ..., California
By Credit call
r I'
I'
--------11
I
I
BOA 0234 I
I
NOTES
NOTES
NOTES
ew from Borland's Scientific & Engineering Division, Turbo
, Pascal Numerical Methods Toolbox implements the latest
high-level mathematical methods to solve common scientific
and engineering problems. Fast.
So every time you need to calcu- • Least squares approximations
late an integral, work with Fourier • Fourier transforms
Transforms or incorporate any of the
classical numerical analysis tools into 5 free ways to look at
your programs, you don't have to "Least Squares Fit"!
reinvent the wheel. Because the As well as a free demo "Fast
Numerical Methods Toolbox is a Fourier Transforms," you also get
complete collection of Turbo Pascal "Least Squares Fit" in 5 different
routines and programs that gives you forms-which gives you 5
applied state-of-the-art math tools. different methods of fitting curves
It also includes two graphics demo to a collection of data pOints.
programs, Least Squares Fit and Fast You instantly get the picture 1 The
Fourier Transforms, to give you the
5 different forms are:
picture along with the numbers.
The Numerical Methods Toolbox is 1. Power 4. 5-term Fourier
a must for you if you're involved with 2. Exponential 5. 5-term
any type of scientific or engineering 3. Logarithm Polynomial
computing. Because it comes with They're all ready to compile
complete source code, you have total and run "as is." To modify or
control of your application. add graphics to your own
programs, you simply add Turbo
What Numerical Methods
Toolbox will do for you now: Graphix Toolbox" to your soft-
• Find solutions to equations ware library. Our Numerical
• Interpolations I
Methods Toolbox is designed to
• Calculus: numerical derivatives work hand-in-hand with our
and integrals Turbo Graphix Toolbox to make
• Matrix operations: inversions, professional graphics in your
determinants and eigenvalues own programs an instant part of
• Differential equations the picture!
Minimum s,stem conllgul.don: IBM PC. XT, AT and true compatibles. PC-DOS (MS-OOSI 2.0 or later 256K. Turbo
Pascal 2.0 or later. The graphics modules reQuire a graphics monilor with an IBM CGA. IBM EGA, or Hercules compatible
adapter card, and reQuire the Turbo Graph~ Toolbox. MS-DOS generic version will not support Turbo Graphix Toolbox
routines. An 8087 or 80287 numeric co-processor is not reQuired, bul recommended lor optimal perlormance.
Turbo Pascal Numerical Methods Toolbox ~ a trademark and Turbo Pascal and Turbo Graphix Toolbox are registered
trademarks 01 Borland International, Inc. IBM, XT, and AT are registered trademarks of International Business Machines
Corp. MS-DOS is a registered trademark 01 Microso" Corp. Hercules is a trademark 01 Hercules Compuler TechnOlogy.
Apple is a registered trademark 01 Apple Computer, Inc. Macintosh is a trademark 01 Mcintosh Laboralory. Inc. licensed
10 Apple Computer. Copyright 1986 Borland Inlernational BOA 0224