pitcon.f90

求解非线性方程组的一组源代码,FORTRAN90.用于解决N个未知数,N-1个方程.

subroutine pitcon ( df, fpar, fx, ierror, ipar, iwork, liw, nvar, rwork, &
  lrw, xr, slname )
!
!*******************************************************************************
!
!! PITCON is the user-interface routine for the continuation code.
!
!
!  A) Introduction:
!  ---------------
!
!  PITCON solves nonlinear systems with one degree of freedom.
!
!  PITCON is given an N dimensional starting point X, and N-1 nonlinear
!  functions F, with F(X) = 0.  Generally, there will be a connected
!  curve of points Y emanating from X and satisfying F(Y) = 0.  PITCON
!  produces successive points along this curve.
!
!  The program can be used to study many sorts of parameterized problems,
!  including structures under a varying load, or the equilibrium
!  behavior of a physical system as some quantity is varied.
!
!  PITCON is a revised version of ACM TOMS algorithm 596.
!
!  Both versions are available via NETLIB, the electronic software
!  distribution service.  NETLIB has the original version in its TOMS
!  directory, and the current version in its CONTIN directory.
!  For more information, send the message "send index from contin"
!  to "[e-mail]".
!
!
!  B) Acknowledgements:
!  -------------------
!
!  PITCON was written by
!
!    Professor Werner C Rheinboldt and John Burkardt,
!    Department of Mathematics and Statistics
!    University of Pittsburgh,
!    Pittsburgh, Pennsylvania, 15260, USA.
!
!    E-Mail: [e-mail]
!            [e-mail]
!
!  The original work on this package was partially supported by the National
!  Science Foundation under grants MCS-78-05299 and MCS-83-09926.
!
!
!  C) Overview:
!  -----------
!
!  PITCON computes a sequence of solution points along a one dimensional
!  manifold of a system of nonlinear equations F(X) = 0 involving NVAR-1
!  equations and an NVAR dimensional unknown vector X.
!
!  The operation of PITCON is somewhat analogous to that of an initial value
!  ODE solver.  In particular, the user must begin the computation by
!  specifying an approximate initial solution, and subsequent points returned
!  by PITCON lie on the curve which passes through this initial point and is
!  implicitly defined by F(X) = 0.  The extra degree of freedom in the system is
!  analogous to the role of the independent variable in a differential
!  equations.
!
!  However, PITCON does not try to solve the algebraic problem by turning it
!  into a differential equation system.  Unlike differential equations, the
!  solution curve may bend and switch back in any direction, and there may be
!  many solutions for a fixed value of one of the variables.  Accordingly,
!  PITCON is not required to parametrize the implicitly defined curve with a
!  fixed parameter.  Instead, at each step, PITCON selects a suitable variable
!  as the current parameter and then determines the other variables as
!  functions of it.  This allows PITCON to go around relatively sharp bends.
!  Moreover, if the equations were actually differentiated - that is, replaced
!  by some linearization - this would introduce an inevitable "drift" away from
!  the true solution curve.  Errors at previous steps would be compounded in a
!  way that would make later solution points much less reliable than earlier
!  ones.  Instead, PITCON solves the algebraic equations explicitly and each
!  solution has to pass an acceptance test in an iterative solution process
!  with tolerances provided by the user.
!
!  PITCON is only designed for systems with one degree of freedom.  However,
!  it may be used on systems with more degrees of freedom if the user reduces
!  the available degrees of freedom by the introduction of suitable constraints
!  that are added to the set of nonlinear equations.  In this sense, PITCON may
!  be used to investigate the equilibrium behavior of physical systems with
!  several degrees of freedom.
!
!  Program options include the ability to search for solutions for which a
!  given component has a specified value.  Another option is a search for a
!  limit or turning point with respect to a given component; that is, of a
!  point where this particular solution component has a local extremum.
!
!  Another feature of the program is the use of two work arrays, IWORK and
!  RWORK.  All information required for continuing any interrupted computation
!  is saved in these two arrays.
!
!
!  D) PITCON Calling Sequence:
!  --------------------------
!
!  subroutine PITCON(DF,FPAR,FX,IERROR,IPAR,IWORK,LIW,NVAR,RWORK,LRW,XR,SLVNAM)
!
!  On the first call, PITCON expects a point XR and a routine FX defining a
!  nonlinear function F.  Together, XR and FX specify a curve of points Y
!  with the property that F(Y) = 0.
!
!  On the first call, PITCON simply verifies that F(XR) = 0.  If this is not
!  the case, the program attempts to correct XR to a new value satisfying
!  the equation.
!
!  On subsequent calls, PITCON assumes that the input vector XR contains the
!  point which had been computed on the previous call.  It also assumes that
!  the work arrays IWORK and RWORK contain the results of the prior
!  calculations.  PITCON estimates an appropriate stepsize, computes the
!  tangent direction to the curve at the given input point, and calculates a
!  predicted new point on the curve.  A form of Newton's method is used to
!  correct this point so that it lies on the curve.  If the iteration is
!  successful, the code returns with a new point XR.  Otherwise, the stepsize
!  may be reduced, and the calculation retried.
!
!  Aside from its ability to produce successive points on the solution curve,
!  PITCON may be asked to search for "target points" or "limit points".
!  Target points are solution vectors for which a certain component has a
!  specified value.  One might ask for all solutions for which XR(17) = 4.0, for
!  instance.  Limit points occur when the curve turns back in a given
!  direction, and have the property that the corresponding component of the
!  tangent vector vanishes there.
!
!  If the user has asked for the computation of target or limit points, then
!  PITCON will usually proceed as described earlier, producing a new
!  continuation point on each return.  But if a target or limit point is
!  found, such a point is returned as the value of XR, temporarily interrupting
!  the usual form of the computation.
!
!
!  E) Overview of PITCON parameters:
!  --------------------------------
!
!  Names of routines:
!
!    DF     Input,        external DF, evaluates the Jacobian of F.
!    FX     Input,        external FX, evaluates the function F.
!    SLVNAM Input,        external SLVNAM, solves the linear systems.
!
!  Information about the solution point:
!
!    NVAR   Input,        integer NVAR, number of variables, the dimension of XR.
!    XR     Input/output, real XR(NVAR), the current solution point.
!
!  Workspace:
!
!    LIW    Input,        integer LIW, the dimension of IWORK.
!    IWORK  Input/output, integer IWORK(LIW), work array.
!
!    LRW    Input,        integer LRW, the dimension of RWORK.
!    RWORK  Input/output, real RWORK(LRW), work array.
!
!    FPAR   "Throughput", real FPAR(*), user defined parameter array.
!    IPAR   "Throughput", integer IPAR(*), user defined parameter array.
!
!  Error indicator:
!
!    IERROR Output,       integer IERROR, error return flag.
!
!
!  F) Details of PITCON parameters:
!  -------------------------------
!
!  DF     Input, external DF, the name of the Jacobian evaluation routine.
!         This name must be declared external in the calling program.
!
!         DF is not needed if the finite difference option is used
!         (IWORK(9) = 1 or 2). In that case, only a dummy name is needed for DF.
!
!         Otherwise, the user must write a routine which evaluates the
!         Jacobian matrix of the function FX at a given point X and stores it
!         in the FJAC array in accordance with the format used by the solver
!         specified in SLVNAM.
!
!         In the simplest case, when the full matrix solverDENSLV solver
!         provided with the package is used, DF must store  D F(I)/D X(J) into
!         FJAC(I,J).
!
!         The array to contain the Jacobian will be zeroed out before DF is
!         called, so that only nonzero elements need to be stored.  DF must
!         have the form:
!
!       subroutine DF(NVAR,FPAR,IPAR,X,FJAC,IERROR)
!
!           NVAR   Input, integer NVAR, number of variables.
!
!           FPAR   Input, real FPAR(*), vector for passing parameters.
!                  This vector is not used by the program, and is only provided
!                  for the user's convenience.
!
!           IPAR   Input, integer IPAR(*), vector for passing integer
!                  parameters.  This vector is not used by the program, and is
!                  only provided for the user's convenience.
!
!           X      Input, real X(NVAR), the point at which the
!                  Jacobian is desired.
!
!           FJAC   Output, real FJAC(*), array containing Jacobian.
!
!                  If DENSLV is the solver:  FJAC must be dimensioned
!                  FJAC(NVAR,NVAR) as shown above, and DF sets
!                  FJAC(I,J) = D F(I)/DX(J).
!
!                  If BANSLV is the solver:  the main portion of the Jacobian,
!                  rows and columns 1 through NVAR-1, is assumed to be a banded
!                  matrix in the standard LINPACK form with lower bandwidth ML
!                  and upper bandwidth MU.  However, the final column of the
!                  Jacobian is allowed to be full.
!
!                  BANSLV will pass to DF the beginning of the storage for
!                  FJAC, but it is probably best not to doubly dimension FJAC
!                  inside of DF, since it is a "hybrid" object.  The first
!                  portion of it is a (2*ML+MU+1, NEQN) array, followed by a
!                  single column of length NEQN (the last column of the
!                  Jacobian).  Thus the simplest approach is to declare FJAC to
!                  be a vector, and then then to store values as follows:
!
!                    If J is less than NVAR, then
!                      if I-J .LE. ML and J-I .LE. MU,
!                        set K = (2*ML+MU)*J + I - ML
!                        set FJAC(K) = D F(I)/DX(J).
!                      else
!                        do nothing, index is outside the band
!                    endif
!                    Else if J equals NVAR, then
!                      set K = (2*ML+MU+1)*(NVAR-1)+I,
!                      set FJAC(K) = D F(I)/DX(J).
!                  endif.
!
!           IERROR Output, integer IERROR, error return flag.  DF should set
!                  this to 0 for normal return, nonzero if trouble.
!
!  FPAR   Input/output, real FPAR(*), a user defined parameter array.
!
!         FPAR is not used in any way by PITCON.  It is provided for the user's
!         convenience.  It is passed to DF and FX, and hence may be used to
!         transmit information between the user calling program and these user
!         subprograms. The dimension of FPAR and its contents are up to the
!         user.  Internally, the program declares DIMENSION FPAR(*) but never
!         references its contents.
!
!  FX     Input, external FX, the name of the routine to evaluate the function.
!         FX computes the value of the nonlinear function.  This name must be
!         declared external in the calling program.  FX should evaluate the
!         NVAR-1 function components at the input point X, and store the result
!         in the vector FVEC.  An augmenting equation will be stored in entry
!         NVAR of FVEC by the PITCON program.  FX should have the form:
!
!           subroutine FX ( NVAR, FPAR, IPAR, X, FVEC, IERROR )
!
!           Input, integer NVAR, number of variables.
!           Input/output, real FPAR(*), user parameters.
!           Input/output, integer IPAR(*), array of user parameters.
!           Input, real X(NVAR), the point of evaluation.
!           Output, real FVEC(NVAR-1), the value of the function at X.  
!           Output, integer IERROR, 0 for no errors, nonzero for an error.
!
!  IERROR Output, integer IERROR, error return flag.
!
!         On return from PITCON, a nonzero value of IERROR is a warning of some
!         problem which may be minor, serious, or fatal.
!
!         0, No errors occurred.
!
!         1, Insufficient storage provided in RWORK and IWORK, or NVAR is less
!            than 2.  This is a fatal error, which occurs on the first call to
!            PITCON.
!
!         2, A user defined error condition occurred in the FX or DF
!          subroutines.  PITCON treats this as a fatal error.
!
!         3, A numerically singular matrix was encountered.  Continuation
!            cannot proceed without some redefinition of the problem.  This is
!            a fatal error.
!
!         4, Unsuccessful corrector iteration.  Loosening the tolerances
!            RWORK(1) and RWORK(2), or decreasing the maximum stepsize RWORK(4)
!            might help.  This is a fatal error.
!
!         5, Too many corrector steps.  The corrector iteration was proceeding
!            properly, but too slowly.  Increase number of Newton steps
!            IWORK(17), increase the error tolerances RWORK(1) or RWORK(2), or
!            decrease RWORK(4).  This is a fatal error.
!
!         6, Null tangent vector.  A serious error which indicates a data
!            problem or singularity in the nonlinear system.  This is a fatal
!            error.
!
!         7, Root finder failed while searching for a limit point.
!            This is a warning.  It means that the limit point
!            computation has failed, but the main computation (computing the
!            continuation curve itself) may continue.
!
!         8, Limit point iteration took too many steps.  This is a warning
!            error.  It means that the limit point computation has failed, but
!            the main computation (computing the continuation curve itself) may
!            continue.
!
!         9, Target point calculation failed.  This generally means that
!            the program detected the existence of a target point, set up
!            an initial estimate of its value, and applied the corrector
!            iteration, but that this corrector iteration failed.
!            This is a only a warning message.  PITCON can proceed to compute
!            new points on the curve.  However, if the target point was
!            really desired, PITCON has no automatic recovery method to
!            retry the calculation.  The best prescription in that case
!            is to try to guarantee that PITCON is taking smaller steps
!            when it detects the target point, which you may do by reducing
!            HMAX, stored as RWORK(4).
!
!         10, Undiagnosed error condition.  This is a fatal error.
!
!  IPAR   Input/output, integer IPAR(*), user defined parameter array.
!
!         IPAR is not used in any way by PITCON.  It is provided for the user's
!         convenience in transmitting parameters between the calling program