NEQBJ

Solves a system of nonlinear equations using factored secant update with a user-supplied Jacobian.

Required Arguments

FCN — User-supplied SUBROUTINE to evaluate the system of equations to be solved. The usage is CALL FCN (N, X, F), where

N – Length of X and F.   (Input)

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)

FCN must be declared EXTERNAL in the calling program.

JAC — User-supplied SUBROUTINE to evaluate the Jacobian at a point X. The usage is CALL JAC (N, X, FJAC, LDFJAC), where

N – Length of X.   (Input)

X – Vector of length N at which point the Jacobian is evaluated.   (Input)
X should not be changed by JAC.

FJAC – The computed N by N Jacobian at the point X.   (Output)

LDFJAC – Leading dimension of FJAC.   (Input)

JAC must be declared EXTERNAL in the calling program.

X — Vector of length N containing the approximate solution.   (Output)

Optional Arguments

N — Dimension of the problem.   (Input)
          Default: N = size (X,1).

XGUESS — Vector of length N containing initial guess of the root.   (Input)
Default: XGUESS = 0.0.

XSCALE — Vector of length N containing the diagonal scaling matrix for the variables.   (Input)
XSCALE is used mainly in scaling the distance between two points. In the absence of other information, set all entries to 1.0. If internal scaling is desired for XSCALE, set IPARAM(6) to 1.
Default: XSCALE = 1.0.

FSCALE — Vector of length N containing the diagonal scaling matrix for the functions.   (Input)
FSCALE is used mainly in scaling the function residuals. In the absence of other information, set all entries to 1.0.
Default: FSCALE = 1.0.

IPARAM — Parameter vector of length 6.   (Input/Output)
Set IPARAM (1) to zero for default values of IPARAM and RPARAM.
See Comment 4.
Default: IPARAM = 0.

RPARAM — Parameter vector of length 5.   (Input/Output)
See Comment 4.

FVEC — Vector of length N containing the values of the functions at the approximate solution.   (Output)

FORTRAN 90 Interface

Generic:          CALL NEQBJ (FCN, JAC, X [,…])

Specific:         The specific interface names are S_NEQBJ and D_NEQBJ.

FORTRAN 77 Interface

Single:            CALL NEQBJ (FCN, JAC, N, XGUESS, XSCALE, FSCALE, IPARAM, RPARAM, X, FVEC)

Double:          The double precision name is DNEQBJ.

Description

Routine NEQBJ uses a secant algorithm to solve a system of nonlinear equations, i. e.,

F (x) = 0

where F : Rn® Rn, and x Ξ Rn.

From a current point, the algorithm uses a double dogleg method to solve the following subproblem approximately:

subject to ||s||2 ≤ dc

to get a direction sc, where F(xc) and J(xc) are the function values and the approximate Jacobian respectively evaluated at the current point xc. Then, the function values at the point xn = xc + sc are evaluated and used to decide whether the new point xn should be accepted.

When the point xn is rejected, this routine reduces the trust region dc and goes back to solve the subproblem again. This procedure is repeated until a better point is found.

The algorithm terminates if the new point satisfies the stopping criterion. Otherwise, dc is adjusted, and the approximate Jacobian is updated by Broyden's formula,

where Jn = J(xn), Jc = J(xc), and y = F (xn) - F (xc). The algorithm then continues using the new point as the current point, i.e. xc ¬ xn.

For more details, see Dennis and Schnabel (1983, Chapter 8).

Comments

1.         Workspace may be explicitly provided, if desired, by use of N2QBJ/DN2QBJ. The reference is:

CALL N2QBJ (FCN, JAC, N, XGUESS, XSCALE, FSCALE, IPARAM, RPARAM, X, FVEC, WK, LWK)

The additional arguments are as follows:

WK — A work vector of length LWK. On output WK contains the following information: The third N locations contain the last step taken. The fourth N locations contain the last Newton step. The final N2 locations contain an estimate of the Jacobian at the solution.

LWK — Length of WK, which must be at least 2 * N2 + 11 * N.   (Input)

2.         Informational errors

Type           Code

   3                  1       The last global step failed to decrease the 2-norm of F(X) sufficiently; either the current point is close to a root of F(X) and no more accuracy is possible, or the secant approximation to the Jacobian is inaccurate, or the step tolerance is too large.

   3                  3       The scaled distance between the last two steps is less than the step tolerance; the current point is probably an approximate root of F(X) (unless STEPTL is too large).

   3                  4       Maximum number of iterations exceeded.

   3                  5       Maximum number of function evaluations exceeded.

   3                  7       Five consecutive steps of length STEPMX have been taken; either the 2-norm of F(X) asymptotes from above to a finite value in some direction or the maximum allowable stepsize STEPMX is too small.

3.         The stopping criterion for NEQBJ occurs when the scaled norm of the functions is less than the scaled function tolerance (RPARAM(1)).

4.         If the default parameters are desired for NEQBJ, then set IPARAM(1) to zero and call routine NEQBJ. Otherwise, if any nondefault parameters are desired for IPARAM or RPARAM, then the following steps should be taken before calling NEQBJ:

CALL N4QBJ (IPARAM, RPARAM)
Set nondefault values for desired IPARAM, RPARAM elements.

Note that the call to N4QBJ will set IPARAM and RPARAM to their default values, so only nondefault values need to be set above.

The following is a list of the parameters and the default values:

IPARAM — Integer vector of length 6.

IPARAM(1) = Initialization flag.

IPARAM(2) = Number of good digits in the function.
Default: Machine dependent.

IPARAM(3) = Maximum number of iterations.
Default: 100.

IPARAM(4) = Maximum number of function evaluations.
Default: 400.

IPARAM(5) = Maximum number of Jacobian evaluations.
Default: not used in NEQBJ.

IPARAM(6) = Internal variable scaling flag.
If IPARAM(6) = 1, then the values of XSCALE are set internally.
Default: 0.

RPARAM — Real vector of length 5.

RPARAM(1) = Scaled function tolerance.
The scaled norm of the functions is computed as

where fi is the i-th component of the function vector F, and fsi is the i-th component of FSCALE.
Default:

where ɛ is the machine precision.

RPARAM(2) = Scaled step tolerance. (STEPTL)
The scaled norm of the step between two points x and y is computed as

where si is the i-th component of XSCALE.

Default: ɛ2/3, where ɛ is the machine precision.

RPARAM(3) = False convergence tolerance.
Default: not used in NEQBJ.

RPARAM(4) = Maximum allowable step size. (STEPMX)

Default:                  1000 * max(ɛ1, ɛ2), where

ɛ2 = ||s||2, s = XSCALE, and t = XGUESS.

RPARAM(5) = Size of initial trust region.
Default: based on the initial scaled Cauchy step.

If double precision is desired, then DN4QBJ is called and RPARAM is declared double  precision.

5.         Users wishing to override the default print/stop attributes associated with error messages issued by this routine are referred to “Error Handling” in the Introduction.

Example

The following 3 ΄ 3 system of nonlinear equations

is solved with the initial guess (4.0, 4.0, 4.0).

 

      USE NEQBJ_INT

      USE UMACH_INT

 

      IMPLICIT   NONE

!                                 Declare variables

      INTEGER    N

      PARAMETER  (N=3)

!

      INTEGER    K, NOUT

      REAL       X(N), XGUESS(N)

      EXTERNAL   FCN, JAC

!                                 Set values of initial guess

!                                 XGUESS = (  4.0  4.0  4.0 )

!

      DATA XGUESS/3*4.0/

!                                 Find the solution

      CALL NEQBJ (FCN, JAC, X, XGUESS=XGUESS)

!                                 Output

      CALL UMACH (2, NOUT)

      WRITE (NOUT,99999) (X(K),K=1,N)

99999 FORMAT ('  The solution to the system is', /, '  X = (', 3F8.3, &

            ')')

!

      END

!                                 User-defined subroutine

      SUBROUTINE FCN (N, X, F)

      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 JAC (N, X, FJAC, LDFJAC)

      INTEGER    N, LDFJAC

      REAL       X(N), FJAC(LDFJAC,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

Output

 

The solution to the system is
X = (   1.000   2.000   3.000)

Visual Numerics, Inc.
Visual Numerics - Developers of IMSL and PV-WAVE
http://www.vni.com/
PHONE: 713.784.3131
FAX:713.781.9260