SECTION 11
SOLVING EQUATIONS
See Acton, F.S., Numerical Methods that Work, New York, Harper & Row, 1970.
This method is designed to find the root of the equation F(X) = 0, where X is a scalar variable (not a vector of variables). It is essential that the equation possess a single root between A and B. In general, one should employ this algorithm only if it can be verified by mathematical analysis that this condition holds. If it does, convergence is guaranteed. If it does not, unpredictable result will follow. The MAIN program must contain the following:
IMPLICIT REAL*8 (A-H,O-Z)
EXTERNAL FUNC
A= lower limit of region containing root
B= upper limit of region containing root
ITERL= iteration limit
ACC= relative accuracy of root (e.g. 1.D-6)
CALL REGFAL(A,B,SOL,FUNC,ACC,ITERL,IER)
where FUNC must be a function subprogram that evaluates the left hand side of the equation. It must have the following structure
FUNCTION FUNC(X)
IMPLICIT REAL*8 (A-H,O-Z)
FUNC= some expression involving X
RETURN
END
Upon returning from REGFAL, SOL will contain the solution. Error conditions are indicated by
IER = 0 Convergence attained
IER =-1 Iteration limit exceeded
IER =-3 A is not less than B
IER =-42 FUNC(A) and FUNC(B) do not have opposite signs; hence
there is either no root or there are several roots in the
interval.
See Acton, F.S., Numerical Methods that Work, New York, Harper & Row, 1970.
NEWRAP is a Newton-Raphson type algorithm with various features to en- hance the probability of convergence. It is designed to solve the system of equations G(X) = 0, where G is a vector valued function and X is the vector of variables. The MAIN program must contain the following
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION X(np)
CHARACTER*8 ALABEL(np)
C where np is the actual number of variables (.le.100)
COMMON/BSTACK/AINT(nq)
C where nq is a number greater than or equal to the scratch storage needed
COMMON/BSTAK/NQ,NTOP
EXTERNAL FQQQQ
CALL DFLT
NP=
NQ=
ITERL=
ACC=
X(1)=
. . .
X(NP)=
CALL LABEL(ALABEL,NP)
CALL NEWRAP(X,NP,ACC,ITERL,ALABEL,FQQQQ,IER)
where
FQQQQ is the name of the subroutine that evaluates G(X) (see
Section 11.2 for the form of this; unlike 11.2, the name here is
arbitrary and up to the user)
NP = np; set variable NP equal to actual number
NQ = nq; set variable NQ equal to number nq
ITERL = iterl; set variable ITERL to actual iteration limit
ACC = acc; set variable ACC to actual accuracy desired
X(1 )= set starting values of X-vector
. . .
X(NP) =
Use CALL LABEL(ALABEL,NP) if automatic labeling of the variables
is desired (see Section 1.5)
STORAGE REQUIREMENT NQ = nq = 3*NP*NP+7*NP
TERMINATION FLAGS
IER = 1 Convergence reached, stepsize meets accuracy
IER = 2 Convergence reached, G'G meets accuracy
IER =-1 Iteration limit exceeded
IER =-3 Input error
IER =-4 Not enough scratch storage
IER =-66 Error in numerical evaluation of Jacobian
IER =-67 Jacobian is singular
IER =-68 Error in evaluating equations
IER =-69 The norm G(X)'G(X) is diverging
Other options.
A. In addition to the above, the user may include
COMMON/BPRINT/IPT,NFILE,NDIF,NPUNCH,JPT,MFILE.
If IPT or JPT are greater than 0, detailed printout by iteration will be produced: IPT,JPT=1 the X-vector is printed out at each iteration IPT,JPT=2 in addition the G-vector is printed out. In addition, the the type of iteration is indicated as follows: NORM normal iteration DIVR dd where dd is an integer, shows that the norm was diverging and required dd shrinkages of the step size to result in diminution of the norm; SING dd (dd is an integer) shows that the Jacobian was singular; required dd increases of the main diagonal to become non singular.
B. The following common blocks may be important to tailor the algorithm to a particular problem:
COMMON/BNWRP1/IVFLIM,INVLIM
COMMON/BNWRP2/EPCOR,EPSVRT
where
IVFLIM = number of step shrinkages permitted in the presence
of errors in evaluating equations (Default=7) INVLIM = number of
consecutive times that Jacobian may be singular. When the
Jacobian is singular, its diagonal is augmented in absolute value
by EPCOR at most INVLIM times (Default=7)
EPCOR = used to augment diagonal of a singular Jacobian
(Default=0.1)
EPSVRT = accuracy required in matrix inversion (Default=1.D-15)
C. The derivatives needed for the Jacobian are evaluated numerically. The intervals over which the finite approximation to the derivative is taken can be altered by including COMMON/BFIDIF/FDFRAC,FDMIN (see Section 3.1)
D. The algorithm terminates if either the step size for each variable becomes less than ACC or if the norm G'G becomes less than ACC. One but not both termination criteria may be suppressed by including in in the MAIN program
COMMON/BSTOP/NVAR1,ISTOP(3)
and setting ISTOP(1)=1 if the stepsize criterion is not to be checked and setting ISTOP(2)=1 if the G'G criterion is not to be checked (see also Section 1.1). ISTOP(3) plays no role in NEWRAP.
E. If NP > 100 is needed, the user should reset the dimensions of the arrays IP and JO in SUBROUTINE MATVRT and should change the integer 100 in statement 5 in SUBROUTINE NEWRAP.
F. EXAMPLE
IMPLICIT REAL*8 (A-H,O-Z)
EXTERNAL FQQQQ
COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
COMMON/BNWRP1/IVFLIM,INVLIM
COMMON/BNWRP2/EPCOR,EPSVRT
COMMON/BFIDIF/FDFRAC,FDMIN
COMMON/BSTACK/AINT(1000)
COMMON/BSTAK/NQ,NTOP
DIMENSION X(25),FPD(25),SPD(25,25)
CHARACTER*8 ALABEL(25)
open (unit=8,file='newtst.out',status='unknown')
call dflt
NQ=1000
X(1)=1.5
X(2)=1.5
X(3)=1.5
X(4)=1.5
X(5)=1.5
ACC=1.D-6
FDFRAC=1.D-4
FDMIN=1.D-6
NP=5 ITERL=50
CALL LABEL(ALABEL,NP)
CALL NEWRAP(X,NP,ACC,ITERL,ALABEL,FQQQQ,IER)
WRITE (6,1234) IER
1234 FORMAT(' IER = ',I3)
STOP
END
SUBROUTINE FQQQQ(X,NP,FPD,*)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION X(NP),FPD(NP)
IF(X(1).LE.0..OR.X(2).LE.0..OR.X(3).LE.0..OR.X(4).LE.0.)
RETURN 1
TEMP=DEXP(X(2))*DLOG(X(3))+X(4)
IF(TEMP.LE.0.)
RETURN 1
FPD(1)=X(1)*TEMP**X(5)-1.
FPD(2)=X(1)+X(2)*X(3)+X(4)*X(5)-3.
FPD(3)=X(1)-X(2)+2.*X(3)**2-X(4)-X(5)
FPD(4)=X(1)**X(2)+X(2)**X(1)-X(4)**X(5)-X(3)**X(4)
FPD(5)=X(1)+X(2)+X(3)+X(4)+X(5)-5.
RETURN
END
G. AN ALTERNATIVE ROUTINE
A second version of NEWRAP is also provided; this version is called NEWRAP1 and has exactly the same calling sequence. This version uses Fortran-90 features, which make it unnecessary for the user to include the BSTACK and BSTAK common blocks and to set the value of NQ.
H. AN ALERT
When linking a program that calls NEWRAP or NEWRAP1, some linkers may complain that FQQQQ has not been defined. This should cause no problems and an executable program will nevertheless be generated.
See Acton, F.S., Numerical Methods that Work, New York, Harper & Row, 1970.
Consider a vector X and a vector-valued function F(X) and suppose that we wish to find a solution such that F(X)=0. SLVMIN solves this problem by converting it into a minimization problem. It creates an objective function that is the sum of the squares of the elements of F(X) and uses OPT to minimize that function with respect to the elements of the X-vector. The artificial function may have local minima at points other than solutions of F(X)=0 and any final results should be examined closely to see that the value of the objective function is adequately close to zero. Also, there may be several solutions to F(X)=0 and SLVMIN will find only one from a particular given starting point. The MAIN program must contain the following
IMPLICIT REAL*8 (A-H,O-Z)
C plus all the statements required for using SUBROUTINE OPT;
CALL OPT(X,NP,F,METHOD,ITERL,MAX,IER,ACC,SLVEQ1,ALABEL)
where X is the vector of variables, SLVEQ1 is a required subroutine name which the user may not change, and where all the other parameters have their usual interpretation; in particular NP is the number of variables and equations in the system. The user must also provide a subroutine as follows:
SUBROUTINE FQQQQ(X,NP,FPD,*)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION X(NP),FPD(NP)
FPD(1)=
FPD(2)=
. . .
FPD(NP)=
RETURN
END
where each FPD( )= is to be set equal to one of the NP functions of X that must equal zero. Note that the user has no choice over the naming of this subroutine: it must be called FQQQQ. Inside the subroutine the user may check for valid X-values; if an invalid X- value is encountered, the user should issue the command RETURN 1.
EXAMPLE: Let the system of equations to be solved be given by
LOG(X(1))+LOG(X(2)) = 0
X(1)+X(2) = 2
which has the solution X(1)=X(2)=1. The subroutine FQQQQ should be as follows:
SUBROUTINE FQQQQ(X,NP,FPD,*)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION X(NP),FPD(NP)
IF(X(1).LE.0..OR.X(2).LE.0.) RETURN 1
FPD(1)=DLOG(X(1))+DLOG(X(2))
FPD(2)=X(1)+X(2)-2.
RETURN
END
STORAGE REQUIREMENTS: Same as for the relevant optimization algorithm used.
TERMINATION FLAGS: Same as for the relevant optimization algorithm used.
Note especially that if GRADX is employed and the termination message is BACKTRACK CYCLE EXCEEDED, a solution has almost certainly been found if the value of the objective function is sufficiently small (e.g. less than 1.D-6).