Solves a system of nonlinear equations using a modified Powell hybrid algorithm with a user-supplied Jacobian.
FCN
User-supplied SUBROUTINE to evaluate
the system of equations to be solved. The usage is
CALL FCN (X, F, N), where
X The point at which
the functions are evaluated. (Input)
X should not be
changed by FCN.
F The computed function values at the point X. (Output)
N Length of X, F. (Input)
FCN must be declared EXTERNAL in the calling program.
LSJAC
User-supplied SUBROUTINE to evaluate
the Jacobian at a point X. The usage
is
CALL LSJAC
(N, X,
FJAC), where
N Length of X. (Input)
X The point at which
the function is evaluated. (Input)
X should not be
changed by LSJAC.
FJAC The computed N by N Jacobian at the point X. (Output)
LSJAC must be declared EXTERNAL in the calling program.
X A vector of
length N.
(Output)
X
contains the best estimate of the root found by NEQNJ.
ERRREL Stopping
criterion. (Input)
The root is accepted if the relative error
between two successive approximations to this root is less than ERRREL.
Default:
ERRREL = 1.e-4
for single precision and 1.d-8 for double precision.
N The number of
equations to be solved and the number of unknowns.
(Input)
Default: N = size (X,1).
ITMAX The
maximum allowable number of iterations. (Input)
Suggested value
= 200.
Default: ITMAX = 200.
XGUESS A vector
of length N. (Input)
XGUESS
contains the initial estimate of the root.
Default: XGUESS =
0.0.
FNORM A scalar that has the value F(1)2 + + F(N)2 at the point X. (Output)
Generic: CALL NEQNJ (FCN, LSJAC, X [, ])
Specific: The specific interface names are S_NEQNJ and D_NEQNJ.
Single: CALL NEQNJ (FCN, LSJAC, ERRREL, N, ITMAX, XGUESS, X, FNORM)
Double: The double precision name is DNEQNJ.
Routine NEQNJ is based on the MINPACK subroutine HYBRDJ, which uses a modification of M.J.D. Powells hybrid algorithm. This algorithm is a variation of Newtons method, which takes precautions to avoid large step sizes or increasing residuals. For further description, see More et al. (1980).
1. Workspace may be explicitly provided, if desired, by use of N2QNJ/DN2QNJ. The reference is:
CALL N2QNJ (FCN, LSJAC, ERRREL, N, ITMAX, XGUESS, X, FNORM, FVEC, FJAC, R, QTF, WK)
The additional arguments are as follows:
FVEC A vector of length N. FVEC contains the functions evaluated at the point X.
FJAC An N by N matrix. FJAC contains the orthogonal matrix Q produced by the QR factorization of the final approximate Jacobian.
R A vector of length N * (N + 1)/2. R contains the upper triangular matrix produced by the QR factorization of the final approximate Jacobian. R is stored row-wise.
QTF A vector of length N. QTF contains the vector TRANS(Q) * FVEC.
WK A work vector of length 5 * N.
2. Informational errors
Type Code
4 1 The number of calls to FCN has exceeded ITMAX. A new initial guess may be tried.
4 2 ERRREL is too small. No further improvement in the approximate solution is possible.
4 3 The iteration has not made good progress. A new initial guess may be tried.
The following 3 Χ 3 system of nonlinear equations
is solved with the initial guess (4.0, 4.0, 4.0).
USE NEQNJ_INT
USE UMACH_INT
IMPLICIT NONE
! Declare variables
INTEGER N
PARAMETER (N=3)
!
INTEGER K, NOUT
REAL FNORM, X(N), XGUESS(N)
EXTERNAL FCN, LSJAC
! Set values of initial guess
! XGUESS = ( 4.0 4.0 4.0 )
!
DATA XGUESS/4.0, 4.0, 4.0/
!
!
CALL UMACH (2, NOUT)
! Find the solution
CALL NEQNJ (FCN, LSJAC, X, XGUESS=XGUESS, FNORM=FNORM)
! Output
WRITE (NOUT,99999) (X(K),K=1,N), FNORM
99999 FORMAT (' The roots found are', /, ' X = (', 3F5.1, &
')', /, ' with FNORM = ',F5.4, //)
!
END
! User-supplied subroutine
SUBROUTINE FCN (X, F, N)
INTEGER N
REAL X(N), F(N)
!
REAL EXP, SIN
INTRINSIC EXP, SIN
!
F(1) = X(1) + EXP(X(1)-1.0) + (X(2)+X(3))*(X(2)+X(3)) - 27.0
F(2) = EXP(X(2)-2.0)/X(1) + X(3)*X(3) - 10.0
F(3) = X(3) + SIN(X(2)-2.0) + X(2)*X(2) - 7.0
RETURN
END
! User-supplied subroutine to
! compute Jacobian
SUBROUTINE LSJAC (N, X, FJAC)
INTEGER N
REAL X(N), FJAC(N,N)
!
REAL COS, EXP
INTRINSIC COS, EXP
!
FJAC(1,1) = 1.0 + EXP(X(1)-1.0)
FJAC(1,2) = 2.0*(X(2)+X(3))
FJAC(1,3) = 2.0*(X(2)+X(3))
FJAC(2,1) = -EXP(X(2)-2.0)*(1.0/X(1)**2)
FJAC(2,2) = EXP(X(2)-2.0)*(1.0/X(1))
FJAC(2,3) = 2.0*X(3)
FJAC(3,1) = 0.0
FJAC(3,2) = COS(X(2)-2.0) + 2.0*X(2)
FJAC(3,3) = 1.0
RETURN
END
The roots found are
X = ( 1.0 2.0
3.0)
with FNORM =.0000
PHONE: 713.784.3131 FAX:713.781.9260 |