home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Phoenix CD 2.0
/
Phoenix_CD.cdr
/
01e
/
libry31a.zip
/
LIBRY4.DOC
< prev
next >
Wrap
Text File
|
1987-01-20
|
42KB
|
1,126 lines
.pa
MATHEMATICAL
These mathematical procedures provide a wide variety of functions that
are useful for engineering applications. I began developing this
library after my first experience with SLIM (name changed to protect the
guilty) back in 1976. So I guess I have them to thank. If you've ever
used SLIM perhaps you've had a similar experience...
Let's say you want to solve what one would think to be a simple problem
of simultaneous linear equations, and due to some unfortunate turn of
events, you are given SLIM with which to do this. You are confronted
with no less than 60 subroutines from which to choose. As it turns out,
none will do what you want. Anyway, you narrow it down to the
Smith-Jones-Simpson-Wilkerson-Feinstein-Gidley method for partially
symmetric matrices with strange sparceness stored in convoluted-
compacted form (not to be confused with space-saving form) and the
Kleidowitz-Yonkers-Lakeville-Massey-Perkins modification of the Orwell-
Gauss-Newton-Schmeltz algorithm for ill-conditioned Hilbert matrices
stored in space-saving form (not to be confused with convoluted-
compacted form). An obvious choice indeed (any fool with 6 Ph.D.s in
math could see that the latter is the preferred method). After numerous
attempts to link your program with all 347 of the necessary SLIM
libraries you are finally ready to run your program. Of course, it
bombs out because you were reading tuesday afternoon's manual and it's
now wednesday morning. Tough luck...
Actually, SLIM is so bad that even I can't achieve hyperbole. Language
fails me. The above could easily be a random selection from one of the
six volumes that are supposed to be a "user's guide."
Anyway, I have a completely different philosophy of programming. I hope
that it works well for you: if there is a clearly preferred method of
doing anything, then use it and discard the others.
An excellent example of this is Gauss quadrature. There are only two
methods of numerical integration that will (at least in theory) converge
given an infinite order: trapezoidal rule and Gauss quadrature. Gauss
quadrature is so far superior in accuracy to any other method (e.g.
Newton-Cotes) it's practically incredible. Therefore, why have 50 other
methods from which to choose? An interesting note about this is that,
contrary to some rumors 10 subdivisions of 20-point Gauss quadrature
(200 point evaluation) is not as accurate as 96-point Gauss quadrature
with less than half the function evaluations.
Another example is Gauss-Jordan elimination. SLIM has subroutines that
provide a "high accuracy solution" (If Andy Rooney understood this, he'd
probably ask, "Who would want a low accuracy solution anyway?"). And of
course, there's the iterative improvement type solutions. What they
don't tell you is that the residual must be computed in greater
precision than the matrix is reduced in (say single to double precision)
which is OK; but why not just use double precision all the way through?
If you are already using double precision and that's not good enough for
you... well, you're just out of luck. If your matrix is ill-
conditioned (as is typically the case with Vandermondes - that's what
you get when you're curve-fitting) about the best thing you can use is
Gauss-Jordan elimination with full column pivoting. It is fast enough
and I have used vector emulation throughout. I have tried all sorts of
things like QR decompositions, Givens rotations, and a whole bunch of
others; but it all boils down to the precision of the math coprocessor.
I only use one algorithm for solving full matrices. I have tested it
using a series of standard Hilbert matrices against SLIM's "high
accuracy solution" and found mine to be both faster and more accurate.
One last example is the solution of ordinary differential equations. If
you read numerical mathematics texts, they always give examples of the
error in such-and-such a method as compared to the exact solution. If
you look closely at the microscopic print in the margin of the legend
under the fly wing that got stuck on the copier, you will probably find
the pithy little statement "exact solution determined using 4-th order
Runge-Kutta with 0.0000000001 step size." Runge-Kutta is self-starting
(which is a bigger pain in the neck than you might think if you use a
method that isn't), convenient to use, and works very well. So why not
use it, I ask? Of course there are certain (unknown before the fact)
cases where other methods are much faster; but I would rather spend
less time experimenting with zillions of algorithms and let the machine
sweat a little more. I've been through it once; and once is enough.
Runge-Kutta is the one for me. A word about automatic stepsize control.
It has been my experience that this takes so long at each step that you
are better off to use a smaller step from the start. I had two such
subroutines that I spent a lot of time on; but I purged them one day in
a fit of frustration.
One final note: WATCH YOUR VARIABLE TYPES! ESPECIALLY DOUBLE
PRECISION! I've seen many a programmer bewildered by a program that
crashes because they have passed "0." instead of "0.D0" to a subroutine
expecting a double precision value or forgotten to put the "REAL*8" out
in front of a function.
REAL*8 FUNCTION DOFX(D)
It is really best to use IMPLICIT statements like
IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z)
than to individually type each variable. You're bound to miss one.
Another biggy is "O" instead of "0" or vise versa.
.pa
QUICK LIST OF MATHEMATICAL FUNCTIONS
AERF: arc-error function
APOLY: area of a polygon
BETA: complete beta function
BETAR: incomplete beta function
BINML: binomial coefficient
CHEBY: Chebyshev polynomial of N-th degree
DWSN: Dawson's integral
ERF: error function
ERFC: complementary error function
FACT: factorial
FGQ0I: numerical integration from zero to infinity
FGQ0ID: double precision form of FGI0I
FGQ1: numerical integration over definite interval
FGQ1D: double precision form of FGQ1
FGQ2: numerical integration over definite interval in 2 dimensions
FGQ3: numerical integration over definite interval in 3 dimensions
FPG: normal probability distribution in percent
FRF: cubic error function
FTC95: Student's T-distribution for 95% confidence
GAMAL: natural log of the Gamma function for large values
GAMMA: Gamma function
GXNYM: numerical integration of X**N*Y**M*dXdY over polygonal region
PLG0: Legendre polynomial of N-th degree
PLG1: Legendre polynomial of N-th degree (first derivitive)
PLG2: Legendre polynomial of N-th degree (second derivitive)
PLG3: Legendre polynomial of N-th degree (third derivitive)
PPLG0: Legendre polynomial of N-th degree in 2 dimensions
RPQ: rational polynomial evaluation
RPQD: double precision form of RPQ
SEVAL: cubic spline evaluation
QUICK LIST OF MATHEMATICAL SUBROUTINES
BANDG: banded matrix solver
BRYDN: modified Broyden's method for nonlinear simultaneous equations
CUBIC: solve a cubic equation
CUBICD: double precision form of CUBIC
FMIN: find the minimum of a function within an interval
GAUSP: Gauss-Jordan elimination of simultaneous equations
GAUSPD: double precision form of GAUSP
JLDAY: convert month/day/year to Julian (annual) day
LLSLC: linear least-squares with linear constraints
LNREG: linear regression
MDAY: convert Julian (annual) day to month/day/year
MINV: matrix inversion
MINVD: double precision form of MINV
PROOT: polynomial root finder
QUAD: solve a quadratic
QUADD: double precision form of QUAD
RK4: 4th order Runge-Kutta (solution of ordinary differential eqns.)
RK6D: 6th order Runge-Kutta (solution of ordinary differential eqns.)
SPLNE: determine coefficients for a cubic spline
SVD: singular value decomposition using Householder transformations
SWMLR: stepwise multiple linear regression
TDIAG: tridiagonal matrix solver
TREND: compute trend of periodic function
ZERO: find the zero (root) of a function within an interval
.pa
NAME: AERF
PURPOSE: arc-error function
TYPE: REAL*4 function (far external)
SYNTAX: X=AERF(E)
INPUT: E (REAL*4) the erf of something
OUTPUT: X (REAL*4) the arc-erf
NAME: APOLY
PURPOSE: area of a polygon
TYPE: REAL*4 function (far external)
SYNTAX: A=APOLY(X,Y,N)
INPUT: X,Y (REAL*4) an array of points (X,Y)
OUTPUT: A (REAL*4) the area enclosed
NOTE: points must be in the order you draw a "connect-the-dots"
picture, the connection between the last point and the first
point is assumed (e.g. for a triangular region N=3)
NAME: BETA
PURPOSE: complete beta function
TYPE: REAL*4 function (far external)
SYNTAX: B=BETA(X,Y)
INPUT: X,Y (REAL*4)
OUTPUT: B (REAL*4)
NAME: BETAR
PURPOSE: incomplete beta function
TYPE: REAL*4 function (far external)
SYNTAX: B=BETAR(X,Y,R)
INPUT: X,Y,R (REAL*4)
OUTPUT: B (REAL*4)
NAME: BINML
PURPOSE: binomial coefficient
TYPE: REAL*4 function (far external)
SYNTAX: B=BINML(N,I)
INPUT: N,I (INTEGER*2)
OUTPUT: B (REAL*4)
NAME: CHEBY
PURPOSE: Chebyshev polynomial of N-th degree
TYPE: REAL*4 function (far external)
SYNTAX: C=CHEBY(N,X)
INPUT: X (REAL*4), N (INTEGER*2)
OUTPUT: C (REAL*4)
NAME: DWSN
PURPOSE: Dawson's integral
TYPE: REAL*4 function (far external)
SYNTAX: D=DWSN(X)
INPUT: X (REAL*4)
OUTPUT: D (REAL*4)
NAME: ERF
PURPOSE: error function
TYPE: REAL*4 function (far external)
SYNTAX: E=ERF(X)
INPUT: X (REAL*4)
OUTPUT: E (REAL*4)
NOTE: this is an unusually fast full machine precision algorithm
NAME: ERFC
PURPOSE: complementary error function
TYPE: REAL*4 function (far external)
SYNTAX: E=ERFC(X)
INPUT: X (REAL*4)
OUTPUT: E (REAL*4)
NAME: FACT
PURPOSE: factorial
TYPE: REAL*4 function (far external)
SYNTAX: F=FACT(N)
INPUT: N (INTEGER*2)
OUTPUT: F (REAL*4)
NAME: FGQ0I (note the "0" (zero))
PURPOSE: numerical integration from zero to infinity
TYPE: REAL*4 function (far external)
SYNTAX: F=FGQ0I(SPARE)
INPUT: SPARE (REAL*4) whatever you like (put 0. if no spare)
FOFX (REAL*4 function) that YOU MUST SUPPLY and must be
called by the sequence F=FOFX(X,SPARE)
OUTPUT: F (REAL*4)
NOTE: for double precision use FGQ0ID and name your function DOFX
20-point Gauss quadrature is used (96-point for double
precision)
NAME: FGQ1
PURPOSE: numerical integration over definite interval
TYPE: REAL*4 function (far external)
SYNTAX: F=FGQ1(A,B,SPARE)
INPUT: A,B (REAL*4) interval over which to integrate FOFX
SPARE (REAL*4) whatever you like (put 0. if no spare)
FOFX (REAL*4 function) that YOU MUST SUPPLY and must be
called by the sequence F=FOFX(X,SPARE)
OUTPUT: F (REAL*4)
NOTE: your function will be integrated from A to B
for double precision use FGQ1D and name your function DOFX
20-point Gauss quadrature is used (96-point for double
precision)
NAME: FGQ2
PURPOSE: numerical integration over definite interval in 2-D
TYPE: REAL*4 function (far external)
SYNTAX: F=FGQ3(A,B,SPARE)
INPUT: A,B (REAL*4) interval over which to integrate FOFX
SPARE (REAL*4) whatever you like (put 0. if no spare)
FOFX (REAL*4 function) that YOU MUST SUPPLY and must be
called by the sequence F=FOFX(X,SPARE)
FY1OFX (REAL*4 function) that YOU MUST SUPPLY and must be
called by the sequence Y1=FY1OFX(X)
FY2OFX (REAL*4 function) that YOU MUST SUPPLY and must be
called by the sequence Y2=FY2OFX(X)
OUTPUT: F (REAL*4)
NOTE: your function will be integrated in X from A to B and Y
from FY1OFX(X) to FY2OFX(X), 20*20=400 point Gauss quadrature
is used
NAME: FGQ3
PURPOSE: numerical integration over definite interval in 3-D
TYPE: REAL*4 function (far external)
SYNTAX: F=FGQ3(A,B,SPARE)
INPUT: A,B (REAL*4) interval over which to integrate FOFX
SPARE (REAL*4) whatever you like (put 0. if no spare)
FOFX (REAL*4 function) that YOU MUST SUPPLY and must be
called by the sequence F=FOFX(X,SPARE)
FY1OFX (REAL*4 function) that YOU MUST SUPPLY and must be
called by the sequence Y1=FY1OFX(X)
FY2OFX (REAL*4 function) that YOU MUST SUPPLY and must be
called by the sequence Y2=FY2OFX(X)
FZ1OFXY (REAL*4 function) that YOU MUST SUPPLY and must be
called by the sequence Z1=FZ1OFXY(X,Y)
FZ2OFXY (REAL*4 function) that YOU MUST SUPPLY and must be
called by the sequence Z2=FZ2OFXY(X,Y)
OUTPUT: F (REAL*4)
NOTE: your function will be integrated in X from A to B and Y
from FY1OFX(X) to FY2OFX(X) and Z from FZ1OFXY(X,Y) to
FZ2OFXY(X,Y), 20*20*20=8000 point Gauss quadrature
is used
NAME: FPG
PURPOSE: normal probability distribution in percent
TYPE: REAL*4 function (far external)
SYNTAX: F=FPG(XBAR,SIGMA,X,DX)
INPUT: XBAR (REAL*4) mean
SIGMA (REAL*4) standard deviation
X (REAL*4) independent variable
DX (REAL*4) increment in X (if you want to know the
probability of 0,5,10,15,20,25,...,95,100%, then X=0,5,10,...
and DX=5)
OUTPUT: F (REAL*4) probability in %
NAME: FRF
PURPOSE: cubic error function
TYPE: REAL*4 function (far external)
SYNTAX: F=FRF(X)
INPUT: X (REAL*4)
OUTPUT: F (REAL*4)
NOTE: FRF is similar to ERF and also varies from -1 to 1 as X varies
from -infinity to +infinity, it pops up in some problems like
the error function (e.g. statistical thermodynamics)
NAME: FTC95
PURPOSE: Student's T-distribution for 95% confidence
TYPE: REAL*4 function (far external)
SYNTAX: F=FTC95(N)
INPUT: N (INTEGER*2) number of degrees of freedom
OUTPUT: F (REAL*4)
NOTE: "Student" is just the guy's hokey pseudonym - it doesn't mean
anything
NAME: GAMAL
PURPOSE: natural log of the Gamma function for large values
TYPE: REAL*4 function (far external)
SYNTAX: G=GAMAL(X)
INPUT: X (REAL*4)
OUTPUT: G (REAL*4)
NAME: GAMMA
PURPOSE: Gamma function
TYPE: REAL*4 function (far external)
SYNTAX: G=GAMMA(X)
INPUT: X (REAL*4)
OUTPUT: G (REAL*4)
NAME: GXNYM
PURPOSE: numerical integration of X**N*Y**M*dXdY over polygonal region
TYPE: REAL*4 function (far external)
SYNTAX: G=GXNYM(X,Y,N,M,L)
INPUT: X,Y (REAL*4) the points defining the boundary of a polygon
N (INTEGER*2) the exponent on X (may be -+0)
M (INTEGER*2) the exponent on Y (may be -+0)
L (INTEGER*2) the number of points
OUTPUT: G (REAL*4)
NOTE: 8-point Gauss quadrature and Green's Lemma are used,
integration will be exact (machine precision) for N+M<16.
If N=M=0 then you get the area and you should use APOLY.
To get the location of the centroid and the moment of inertia
use the following
XC=GXNYM(X,Y,1,0,L)/GXNYM(X,Y,0,0,L)
YC=GXNYM(X,Y,0,1,L)/GXNYM(X,Y,0,0,L)
AI=GXNYM(X,Y,1,1,L)
You can get radii of gyration and other useful things like FEM
coefficients from GXNYM.
NAME: PLG0
PURPOSE: Legendre polynomial of N-th degree
TYPE: REAL*4 function (far external)
SYNTAX: P=PLG0(N,X)
INPUT: N (INTEGER*2) degree
X (REAL*4)
OUTPUT: P (REAL*4)
NAME: PLG1
PURPOSE: Legendre polynomial of N-th degree (first derivitive)
TYPE: REAL*4 function (far external)
SYNTAX: P=PLG1(N,X)
INPUT: N (INTEGER*2) degree
X (REAL*4)
OUTPUT: P (REAL*4)
NAME: PLG2
PURPOSE: Legendre polynomial of N-th degree (second derivitive)
TYPE: REAL*4 function (far external)
SYNTAX: P=PLG2(N,X)
INPUT: N (INTEGER*2) degree
X (REAL*4)
OUTPUT: P (REAL*4)
NAME: PLG3
PURPOSE: Legendre polynomial of N-th degree (third derivitive)
TYPE: REAL*4 function (far external)
SYNTAX: P=PLG3(N,X)
INPUT: N (INTEGER*2) degree
X (REAL*4)
OUTPUT: P (REAL*4)
NAME: PPLG0
PURPOSE: Legendre polynomial of N-th degree (2 dimensions)
TYPE: REAL*4 function (far external)
SYNTAX: P=PPLG0(N,X,Y)
INPUT: N (INTEGER*2) degree
X,Y (REAL*4)
OUTPUT: P (REAL*4)
NAME: RPQ
PURPOSE: rational polynomial evaluation
TYPE: REAL*4 function (far external)
SYNTAX: R=RPQ(P,N,Q,M,X)
INPUT: P (REAL*4) coefficients in numerator polynomial
N (INTEGER*2) number of terms in P
Q (REAL*4) coefficients in denominator polynomial
M (INTEGER*2) number of terms in Q
X (REAL*4)
OUTPUT: R (REAL*4)
NOTE: RPQ evaluates a rational polynomial as shown below using
Horner's method
R=(P1+P2*X+P3*X**2+P4*X**3+...)/(Q1+Q2*X+Q3*X**2+Q4*X**3...)
for double precision use RPQD and make sure that you declare
R,P,Q,RPQD all to be REAL*8
NAME: SEVAL
PURPOSE: cubic spline evaluation
TYPE: REAL*4 function (far external)
SYNTAX: Y=SEVAL(XI,YI,N,C,X)
INPUT: XI,YI (REAL*4) specified points
N (INTEGER*2) number of points
C (REAL*4) coefficients determined by first calling SPLNE
X (REAL*4) point where you want to interpolate
OUTPUT: Y (REAL*4) interpolated value
NOTE: you need to first call SPLNE, you only need to do this once
NAME: BANDG
PURPOSE: banded matrix solver
TYPE: subroutine (far external)
SYNTAX: CALL BANDG(A,B,X,N,M,IER)
INPUT: A,B (REAL*4) matrices
N (INTEGER*2) number of equations
M (INTEGER*2) bandwidth
OUTPUT: X (REAL*4) solution vector
IER (INTEGER*2) error indicator
IER=0 indicates normal return
IER=1 indicates M too small (M<3)
IER=2 indicates m too big (M>N)
IER=3 indicates M not odd
IER=4 indicates [A] singular
NOTE: The Gauss-Seidel method is used.
The system of equations must be of the form
[A][X]=[B]
[A] should look something like the following (note zeroes)...
if N=3 if N=5 if N=7
dimension A(3,N) A(5,N) A(7,N)
0 2 -2 0 0 3 -2 -1 0 0 0 6 -3 -2 -1
-1 2 -1 0 -2 5 -2 -1 0 0 -3 9 -3 -2 -1
-1 2 -1 -1 -2 6 -2 -1 0 -2 -3 11 -3 -2 -1
-1 2 -1 -1 -2 6 -2 -1 -1 -2 -3 12 -3 -2 -1
: : : : : : : : : : : : : : :
-1 2 1 -1 -2 6 -2 -1 -1 -2 -3 12 -3 -2 -1
-1 2 1 -1 -2 6 -2 -1 -1 -2 -3 11 -3 -2 0
-1 2 1 -1 -2 5 -2 0 -1 -2 -3 9 -3 0 0
-2 2 0 -1 -2 3 0 0 -1 -2 -3 6 0 0 0
NAME: BRYDN
PURPOSE: modified Broyden's method for nonlinear simultaneous equations
TYPE: subroutine (far external)
SYNTAX: CALL BRYDN(XMIN,XMAX,X,F,N,M,W,MW,PRINT,IER)
INPUT: XMIN,XMAX (REAL*4) arrays, upper and lower limits on X
N (INTEGER*2) number of elements in X
M (INTEGER*2) number of residuals
W (REAL*4) working space of dimension MW
MW (INTEGER*2) dimension of working space MW>=4N+3M+N*N+N*M
PRINT (LOGICAL*2) set to .true. for printer output (also you
need to open file 6 on the PC)
RESID (a far external subroutine) which YOU MUST SUPPLY that
may be called by the following
CALL RESID(X,F)
and returns "M" different Fs given "N" Xs
OUTPUT: X (REAL*4) solution vector
F (REAL*4) final residual
IER (INTEGER*2) error indicator
IER=0 indicates normal return
IER=1 indicates N<2
IER=2 indicates M<N
IER=3 indicates XMIN>XMAX conflict
IER=4 indicates insufficient work space (increase MW)
NOTE: This is a FAR more difficult problem than most people even
imagine! You can compare it to finding the lowest spot on
the face of the Earth given the rules of the game "Battle
Ship" where all you can do is call out coordinates and your
opponent answers with the elevation. You may find a rut or a
ditch somewhere; but it's highly unlikely that you will find
death valley, let alone some trench in the Pacific. Now extend
this analogy to many-dimensional space... get the point? So
don't get too upset if BRYDN can't find the solution you want.
I've tried all sorts of algorithms and this seems to work the
best for general problems.
If you have 6 equations and 6 unknowns then M=N=6. If you
have 12 equations and 6 unknowns then M=12 and N=6.
If all of your equations are linear, by all means don't use
BRYDN, you want LLSLC.
If you have only one unknown then use FMIN.
NAME: CUBIC
PURPOSE: solve a cubic equation
TYPE: subroutine (far external)
SYNTAX: CALL SUBROUTINE CUBIC(A1,A2,A3,A4,R1,C1,R2,C2,R3,C3)
INPUT: A1,A2,A3,A4 (REAL*4) coefficients in cubic equation
OUTPUT: (R1,C1),(R2,C2),(R3,C3) (REAL*4) roots
NOTE: the equation has the form A1*X**3+A2*X**2+A3*X+A4=0
for double precision use CUBICD and don't forget to declare
A1,A2,A3,A4,R1,C1,R2,C2,R3,C3 all to be REAL*8 and if you
use constants, don't forget the 1.D0 or whatever.
NAME: FMIN
PURPOSE: find the minimum of a function within an interval
TYPE: subroutine (far external)
SYNTAX: CALL FMIN(A,B,X,F,SPARE)
INPUT: A,B (REAL*4) the search interval
SPARE (REAL*4) whatever you want to pass to your function
FOFX (a far external REAL*4 function) that you must supply
that may be called by the following
F=FOFX(X,SPARE)
OUTPUT: X (REAL*4) location of minimum
F (REAL*4) function at minimum
NAME: GAUSP
PURPOSE: Gauss-Jordan elimination of simultaneous equations
TYPE: subroutine (far external)
SYNTAX: CALL GAUSP(A,B,X,JPIVOT,N,D,C,IER)
INPUT: A,B (REAL*4) arrays containing the equations to be solved
having dimensions (N,N) and (N) respectively
JPIVOT (INTEGER*2) working space of dimension "N"
N (INTEGER*2) the number of equations
OUTPUT: X (REAL*4) solution vector of dimension N
D (REAL*4) ALOG10(ABS(DETERMINANT)) log of the determinant
C (REAL*4) ALOG10(AMAX1(ABS(PIVOT))/AMIN1(ABS(PIVOT))) a
measure of the gradedness or condition of the matrix
IER (INTEGER*2) error indicator (IER=0 is normal)
NOTE: For double precision use GAUSPD and don't forget to declare
A,B,X,D,C all to be REAL*8.
Note that A and B will be destroyed upon return.
Note that A is in row-order as indicated below
A(1,1)*X(1)+A(1,2)*X(2)+...A(1,N)*X(N)=B(1)
A(2,1)*X(1)+A(2,2)*X(2)+...A(2,N)*X(N)=B(2)
A(3,1)*X(1)+A(3,2)*X(2)+...A(3,N)*X(N)=B(3)
..........................................
A(N,1)*X(1)+A(N,2)*X(2)+...A(N,N)*X(N)=B(N)
When you really have to be careful about this is with DATA
statements as FORTRAN fills arrays in column-order.
Full column pivoting is used.
Vector emulation is used throughout to give maximum speed.
NAME: JLDAY
PURPOSE: convert month/day/year to Julian (annual) day
TYPE: subroutine (far external)
SYNTAX: CALL JLDAY(MONTH,IDAY,IYEAR,JDAY)
INPUT: MONTH,IDAY,IYEAR (INTEGER*2)
OUTPUT: JDAY (INTEGER*2)
NOTE: I know that "Julian" day is a misnomer; but that's what
everyone else calls it.
NAME: LLSLC
PURPOSE: linear least-squares with linear constraints
TYPE: subroutine (far external)
SYNTAX: CALL LLSLC(X,Y,A,NE,NC,W,NW,IOPT,IER)
INPUT: X,Y,A,W (REAL*8) see example below
NE (INTEGER*2) number of variables (X or A)
NC (INTEGER*2) number of constraints (may be zero)
NW (INTEGER*2) working space, NW>=(NE+NC)*(NE+NC+3))
IOPT (INTEGER*2) option, see example below
OUTPUT: A (REAL*8)
IER (INTEGER*2) error indicator
IER=0 normal
IER=1 dimension error
IER=2 singular matrix - no inverse
IER=3 too few constants+constraints (minimum=2)
IER=4 too few equations (minimum=1)
IER=5 too few constraints (minimum=0)
IER=6 too many constraints (maximum=ne-1)
IER=7 IOPT not defined
IER=8 matrices not initialized
IER=9 too few points for least-squares
IER=10 constraints added before least-squares points
IER=11 too many constraints were added
IER=12 too few constraints were specified
IER=13 working space too small (increase NW)
NOTE: This is a little tricky, so I have included an example.
Sorry, there is no REAL*4 equivalent.
eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
PROGRAM EXMPL
C
C THIS IS AN EXAMPLE OF HOW TO USE LLSLC
C
C IN THIS EXAMPLE A POWER SERIES (1,X,X**2,X**3,...) IS USED TO
C APPROXIMATE THE FUNCTION F(T)=ERFC(T), THE APPROXIMATION
C IS OVER T=0,1 AND THE CONSTRAINTS ARE THAT THE APPROXIMATION
C FIT EXACTLY AT THE END POINTS
C
IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z)
REAL*4 ERFC
PARAMETER (NE=5,NC=2,NS=NE+NC,NF=NE-NC,NW=NS*(NS+3),M=20)
DIMENSION X(NE),A(NE),WORK(NW),T(M),F(M)
C
OPEN(6,FILE='PRN')
C
C DEFINE THE FUNCTION F(T)=ERFC(T)
C
DO 100 I=1,M
T(I)=DBLE(I-1)/DBLE(M-1)
100 F(I)=DBLE(ERFC(SNGL(T(I))))
C
WRITE(6,1000) NE,NC,NF,NS
1000 FORMAT('1LINEAR LEAST-SQUARES WITH LINEAR CONSTRAINTS',//,
&' NUMBER OF TERMS.................... ',I3,/,
&' NUMBER OF CONSTRAINTS.............. ',I3,/,
&' NUMBER OF FREE CONSTANTS........... ',I3,/,
&' NUMBER OF SIMULTANEOUS EQUATIONS... ',I3)
C
C STEP#1 INITIALIZE (NOTE: IOPT=1)
C
CALL LLSLC(X,Y,A,NE,NC,WORK,NW,1,IER)
IF(IER.NE.0) GO TO 900
C
C STEP#2 FEED THE DATA TO LLSLC (NOTE: IOPT=2)
C
DO 210 I=1,M
TI=T(I)
FI=F(I)
X(1)=1.D0
DO 200 J=2,NE
200 X(J)=X(J-1)*TI
Y=FI
C
CALL LLSLC(X,Y,A,NE,NC,WORK,NW,2,IER)
210 IF(IER.NE.0) GO TO 900
C
C STEP#3.1 ADD CONSTRAINT#1 (NOTE: IOPT=3)
C NOTE: IF YOU HAVE NO CONSTRAINTS SET NC=0 AND SKIP SECTION 3
C
I=1
X(I)=1.D0
TI=T(I)
FI=F(I)
DO 310 J=2,NE
310 X(J)=X(J-1)*TI
Y=FI
C
CALL LLSLC(X,Y,A,NE,NC,WORK,NW,3,IER)
IF(IER.NE.0) GO TO 900
C
C STEP#3.2 ADD CONSTRAINT#2 (NOTE: IOPT=3)
C
I=M
TI=T(I)
FI=F(I)
X(I)=1.D0
DO 320 J=2,NE
320 X(J)=X(J-1)*TI
Y=FI
C
CALL LLSLC(X,Y,A,NE,NC,WORK,NW,3,IER)
IF(IER.NE.0) GO TO 900
C
C STEP#4 SOLVE FOR COEFFICIENTS (NOTE: IOPT=4)
C
CALL LLSLC(X,Y,A,NE,NC,WORK,NW,4,IER)
IF(IER.NE.0) GO TO 900
C
C LIST Q-FACTOR, COEFFICIENTS, AND SIGNIFICANCE LEVELS
C
WRITE(6,5000) Y
5000 FORMAT(/,' Q-FACTOR=',1PG12.5,' (Q<1 INDICATES STABLE, ',
&'Q>1 INDICATES UNSTABLE)',//,' I A S')
C
DO 500 I=1,NE
500 WRITE(6,5100) I,A(I),X(I)
5100 FORMAT(1X,I3,2(1X,1PG12.5))
C
C LIST AGREEMENT
C
WRITE(6,5200)
5200 FORMAT(/,' AGREEMENT',/,
&' I X F Y E')
C
EMAX=0.D0
EAVE=0.D0
C
DO 520 I=1,M
X(1)=1.D0
Y=A(1)*X(1)
DO 510 J=2,NE
X(J)=X(J-1)*T(I)
510 Y=Y+A(J)*X(J)
C
E=Y-F(I)
EAVE=EAVE+DABS(E)
IF(DABS(E).GT.DABS(EMAX)) EMAX=E
520 WRITE(6,5300) I,T(I),F(I),Y,E
5300 FORMAT(1X,I3,4(1X,1PG12.5))
C
C LIST AVERAGE AND MAXIMUM ERROR
C
EAVE=EAVE/DBLE(M)
WRITE(6,5400) EAVE,EMAX
5400 FORMAT(30X,'AVERAGE ERROR=',1PG12.5,/,
&30X,'MAXIMUM ERROR=',1PG12.5)
GO TO 999
C
900 WRITE(6,9000) IER
9000 FORMAT(' ERROR CODE:',I5)
C
999 CLOSE(6)
STOP
END
eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
NAME: LNREG
PURPOSE: linear regression
TYPE: subroutine (far external)
SYNTAX: CALL LNREG(X,Y,N,A,B,R,IER)
INPUT: X,Y (REAL*4) data points
N (INTEGER*2) number of points
OUTPUT: A,B (REAL*4) slope, intercept
R (REAL*4) regression coefficient in %
IER (INTEGER*2) error indicator (IER=0 is normal)
NOTE: LNREG only fits a straight line of the form Y=A*X+B to a set
of points. If you want something more elaborate use LLSLC.
NAME: MDAY
PURPOSE: convert Julian (annual) day to month/day/year
TYPE: subroutine (far external)
SYNTAX: CALL MDAY(JDAY,IYEAR,MONTH,IDAY)
INPUT: JDAY,IYEAR (INTEGER*2)
OUTPUT: MONTH,IDAY (INTEGER*2)
NOTE: I know that "Julian" day is a misnomer; but that's what
everyone else calls it.
NAME: MINV
PURPOSE: matrix inversion
TYPE: subroutine (far external)
SYNTAX: CALL MINV(A,IPIVOT,JPIVOT,N,D,C,IER)
INPUT: A (REAL*4) array containing the matrix to be inverted
IPIVOT,JPIVOT (INTEGER*2) working spaces of dimension "N"
N (INTEGER*2) rank (the number of equations)
OUTPUT: A (REAL*4) inverted matrix
D (REAL*4) ALOG10(ABS(DETERMINANT)) log of the determinant
C (REAL*4) ALOG10(AMAX1(ABS(PIVOT))/AMIN1(ABS(PIVOT))) a
measure of the gradedness or condition of the matrix
IER (INTEGER*2) error indicator (IER=0 is normal)
NOTE: For double precision use MINVD and don't forget to declare
A,D,C all to be REAL*8.
Note that A is inverted "in place."
Note that A is in row-order (see GAUSP).
NAME: PROOT
PURPOSE: polynomial root finder
TYPE: subroutine (far external)
SYNTAX: CALL PROOT(P,R,C,N,IER)
INPUT: P (REAL*8) array containing the polynomial coefficients
N (INTEGER*2) number of terms in P
OUTPUT: R (REAL*8) array containing the real part of the roots
C (REAL*8) array containing the complex part of the roots
IER (INTEGER*2) error indicator
IER=0 indicates no error
IER=1 indicates input error
IER=2 indicates failure to converge
NOTE: Bairstow's method is used.
The polynomial must be of the form
P1+P2*X+P3*X**2+P4*X**3+....=0
Also note that this is a VERY difficult problem for large N,
so don't be too upset if you get complex roots when you
thought you should be getting real roots. That's life in the
real world.
NAME: QUAD
PURPOSE: solve a quadratic
TYPE: subroutine (far external)
SYNTAX: CALL SUBROUTINE QUAD(A1,A2,A3,R1,C1,R2,C2)
INPUT: A1,A2,A3 (REAL*4) coefficients in quadratic equation
OUTPUT: (R1,C1),(R2,C2) (REAL*4) roots
NOTE: the equation has the form A1*X**2+A2*X+A3=0
for double precision use QUADD and don't forget to declare
A1,A2,A3,R1,C1,R2,C2 all to be REAL*8 and if you use
constants, don't forget the 1.D0 or whatever.
NAME: RK4
PURPOSE: 4th order Runge-Kutta (solve ordinary differential equations)
TYPE: subroutine (far external)
SYNTAX: CALL RK4(X,DX,Y,DY,N,W,V,SPARE)
INPUT: X (REAL*4) independent variable
DX (REAL*4) increment in X (step size)
Y (REAL*4) array, dependent variables at X
DY (REAL*4) array, change in Y
N (INTEGER*2) number of dependent variables (Y)
W (REAL*4) working space of dimension (N,4)
V (REAL*4) working space of dimension (N)
SPARE (REAL*4) whatever you want to pass to DYDX
DYDX (far external) a subroutine that YOU MUST PROVIDE and
must be called by the sequence
CALL DYDX(X,Y,DY,N,SPARE)
OUTPUT: Y (REAL*4) new values of dependent variable
NOTE: this isn't too tricky; but I have provided an example that may
help
eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
PROGRAM EXMPL
C
C THIS IS AN EXAMPLE OF HOW TO USE RK4
C IN THIS EXAMPLE A 3-rd ORDER DIFFERENTIAL EQUATION WILL BE SOLVED
C
IMPLICIT INTEGER*2 (I-N),REAL*4 (A-H,O-Z)
PARAMETER (N=3)
DIMENSION Y(N),DY(N),W(N,4),V(N)
C
C DEFINE THE INITIAL CONDITIONS, THE STEP SIZE, AND THE NUMBER OF STEPS
C
DATA Y/0.,0.,0./
DATA DX,NSTEP/.1,100/
C
OPEN(6,FILE='PRN')
C
WRITE(6,1000)
1000 FORMAT('1SOLUTION OF 3-RD ORDER EQUATION',//,
&' X DDY/DXX DY/DX Y')
C
DO 100 ISTEP=1,NSTEP
CALL RK4(X,DX,Y,DY,N,W,V,SPARE)
100 WRITE(6,1100) X,Y
1100 FORMAT(4(1X,1PG12.5))
C
CLOSE(6)
STOP
END
SUBROUTINE DYDX(X,Y,DY,N,SPARE)
C
C YOU MUST PROVIDE THIS SUBROUTINE!
C
IMPLICIT INTEGER*2 (I-N),REAL*4 (A-H,O-Z)
DIMENSION Y(N),DY(N)
C
DY(1)=X*COS(X)-SIN(X)
DY(2)=Y(1)
DY(3)=Y(2)
C
RETURN
END
eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
NAME: RK6D
PURPOSE: 6th order Runge-Kutta (solve ordinary differential equations)
TYPE: subroutine (far external)
SYNTAX: CALL RK6D(X,DX,Y,DY,N,W,V,SPARE)
INPUT: X (REAL*8) independent variable
DX (REAL*8) increment in X (step size)
Y (REAL*8) array, dependent variables at X
DY (REAL*8) array, change in Y
N (INTEGER*2) number of dependent variables (Y)
W (REAL*8) working space of dimension (N,8) <-- note that the
"8" here is not a misprint - it takes 8 steps to get 6th order
even though it only takes 4 to get 4-order Runge-Kutta
V (REAL*4) working space of dimension (N)
SPARE (REAL*8) whatever you want to pass to DYDX
DYDX (far external) a subroutine that YOU MUST PROVIDE and
must be called by the sequence
CALL DYDX(X,Y,DY,N,SPARE)
OUTPUT: Y (REAL*8) new values of dependent variable
NOTE: This is just like RK4 except for the "8" above and that you
need to declare everything REAL*8 instead of REAL*4.
NAME: SPLNE
PURPOSE: determine coefficients for a cubic spline
TYPE: subroutine (far external)
SYNTAX: CALL SPLNE(XI,YI,N,C)
INPUT: XI,YI (REAL*4) points to be fit
N (INTEGER*2) number of points
OUTPUT: C (REAL*4) coefficient array of dimension (3,N)
NOTE: you need to call this once, then use SEVAL to evaluate spline
NAME: SVD
PURPOSE: singular value decomposition using Householder transformations
TYPE: subroutine (far external)
SYNTAX: CALL SVD(A,N,M,S,U,V,W,IER)
INPUT: A (REAL*8) the matrix (N,M) to be decomposed - row order
N (INTEGER*2) number of rows in A
M (INTEGER*2) number of columns in A
W (REAL*8) working space of dimension MAX0(N,M)
OUTPUT: S (REAL*8) singular values
U (REAL*8) left eigenvectors (M,N)
V (REAL*8) right eigenvectors (N,M)
IER (INTEGER*2) error indicator (IER=0 is normal)
NAME: SWMLR
PURPOSE: stepwise multiple linear regression
TYPE: subroutine (far external)
SYNTAX: CALL SWMLR(MA,MF,IER,ND,NF,R,X,Y,A,ATA,ATB,JPIVOT,
&F,FTF,FTR,LF,SPARE)
INPUT: MA (INTEGER*2) maximum number of terms in regression
MF (INTEGER*2) maximum number of functions
IER (INTEGER*2) error indicator (IER=0 is normal)
ND (INTEGER*2) number of data points (X,Y)
X,Y (REAL*8) data points to fit
ATA (REAL*8) working space of dimension (MA,MA)
ATB (REAL*8) working space of dimension (MA)
ATB (INTEGER*2) working space of dimension (MA)
F (REAL*8) working space of dimension (MA)
FTF (REAL*8) working space of dimension (MA)
FTR (REAL*8) working space of dimension (MA)
OUTPUT: NF (INTEGER*2) number of functions used
R (REAL*8) residual error after approximation
A (REAL*8) coefficients
LF (INTEGER*2) order in which the functions are used
NOTE: This is a little tricky, so an example is provided below.
Note that SWMLR will find the set of functions providing the
most compact fit. This is different from the fit that has
the least error or has the largest F-value and a LOT HARDER
TO DETERMINE.
eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
PROGRAM EXAMPL
C
C THIS IS AN EXAMPLE OF HOW TO USE SWMLR
C
IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z)
PARAMETER (MA=8,MF=8,ND=20)
REAL*4 ERF
CHARACTER NAME*10
DIMENSION X(ND),Y(ND),A(MA),ATA(MA,MA),ATB(MA),JPIVOT(MA),F(MF),
&FTF(MF),FTR(MF),LF(MF),R(ND),NAME(MF)
DATA NAME/
&'1. ',
&'X ',
&'X*X ',
&'X*X*X ',
&'EXP(X) ',
&'EXP(-X) ',
&'X*EXP(X) ',
&'X*EXP(-X) '/
C
OPEN(6,FILE='PRN')
WRITE(6,1000)
1000 FORMAT('1EXAMPLE OF STEPWISE LINEAR REGRESSION')
C
C DEFINE THE DATA TO BE FIT (THE ERROR FUNCTION)
C
DO 100 I=1,ND
X(I)=DBLE(I-1)/DBLE(ND-1)
100 Y(I)=DBLE(ERF(SNGL(X(I))))
C
CALL SWMLR(MA,MF,IER,ND,NF,R,X,Y,A,ATA,ATB,JPIVOT,
&F,FTF,FTR,LF,SPARE)
C
WRITE(6,2000)
2000 FORMAT(/,' I J A NAME')
C
DO 200 I=1,NF
200 WRITE(6,2010) I,LF(I),A(I),NAME(LF(I))
2010 FORMAT(2(2X,I1),2X,1PG12.5,2X,A10)
C
WRITE(6,2100)
2100 FORMAT(/,' X Y APPROX ERROR')
C
EAVE=0.
EMAX=0.
C
DO 211 I=1,ND
CALL SOFX(X(I),F,SPARE)
C
APPROX=0.D0
DO 210 J=1,NF
210 APPROX=APPROX+A(J)*F(LF(J))
C
E=APPROX-Y(I)
EAVE=EAVE+DABS(E)
EMAX=DMAX1(EMAX,DABS(E))
211 WRITE(6,2110) X(I),Y(I),APPROX,E
2110 FORMAT(4(1X,1PG12.5))
C
EAVE=EAVE/DBLE(ND)
WRITE(6,2200) EAVE,EMAX
2200 FORMAT(/,
&' AVERAGE ERROR=',1PG12.5,/,
&' MAXIMUM ERROR=',1PG12.5)
C
CLOSE(6)
STOP
END
SUBROUTINE SOFX(X,F,SPARE)
C
C YOU MUST PROVIDE THIS SUBROUTINE!
C
IMPLICIT INTEGER*2 (I-N),REAL*8 (A-H,O-Z)
PARAMETER (MF=8)
DIMENSION F(MF)
C
F(1)=1.
F(2)=X
F(3)=X*X
F(4)=X*X*X
F(5)=EXP(X)
F(6)=EXP(-X)
F(7)=X*EXP(X)
F(8)=X*EXP(-X)
C
RETURN
END
eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
NAME: TDIAG
PURPOSE: tridiagonal matrix solver
TYPE: subroutine (far external)
SYNTAX: CALL TDIAG(D,B,X,W,N,IER)
INPUT: D (REAL*4) matrix of dimension (N,3)
B (REAL*4) right-hand-side of dimension (N)
W (REAL*4) working space of dimension (N,3)
N (INTEGER*2) number of equations
OUTPUT: X (REAL*4) solution vector of dimension (N)
IER (INTEGER*2) error indicator (IER=0 is normal)
NOTE: the tridiagonal system of equations is defined by [D][X]=[B]
[D] should look something like the following
0 2 -2
-1 2 -1
-1 2 -1
. . .
-2 2 0
For double precision use TDIAGD and don't forget to declare
D,B,X,W all to be of type REAL*8.
NAME: TREND
PURPOSE: compute trend of periodic function
TYPE: subroutine (far external)
SYNTAX: CALL TREND(X,Y,W,NX,N,C,NC)
INPUT: X,Y (REAL*4) points to be fit
W (REAL*4) working space of dimension (NX)
NX (INTEGER*2) number of points
N (INTEGER*2) number of known points (N<NX)
NC (INTEGER*2) number of terms in expansion (try 6)
OUTPUT: Y (REAL*4) upon return, the N+1,N+2,N+3,...NX values of Y will
be trended
C (REAL*4) coefficients in trend series
NOTE: Chebyshev approximation is used
NAME: ZERO
PURPOSE: find the zero (root) of a function within an interval
TYPE: subroutine (far external)
SYNTAX: CALL ZERO(A,B,X,IER,SPARE)
INPUT: A,B (REAL*4) interval to look for root in
SPARE (REAL*4) whatever you want to pass to your function
FOFX (REAL*4 function) that YOU MUST SUPPLY and must be
called by the sequence F=FOFX(X,SPARE)
OUTPUT: X (REAL*4) the location of the root within (A<=X<=B)
IER (INTEGER*2) error indicator (IER=0 is normal)