pitcon.f90

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

!                  a point.  RWORK(7) does not have a default value.  The
!                  program does not set it, and it is not referenced unless
!                  IWORK(5) has been set.
!
!  RWORK(8)        EPMACH, the value of the machine precision.  The computer
!                  can distinguish 1.0+EPMACH from 1.0, but it cannot
!                  distinguish 1.0+(EPMACH/2) from 1.0. This number is used
!                  when estimating a reasonable accuracy request on a given
!                  computer.  PITCON computes a value for EPMACH internally.
!
!  RWORK(9)        STEPX, the size, using the maximum-norm, of the last
!                  step of Newton correction used on the most recently
!                  computed point, whether a starting point, continuation
!                  point, limit point or target point.
!
!  RWORK(10)       A minimum angle used in the steplength computation,
!                  equal to 2.0*ARCCOS(1-EPMACH).
!
!  RWORK(11)       Estimate of the angle between the tangent vectors at the
!                  last two continuation points.
!
!  RWORK(12)       The pseudo-arclength coordinate of the previous continuation
!                  pointl; that is, the sum of the Euclidean distances between
!                  all computed continuation points beginning with the start
!                  point.  Thus each new point should have a larger coordinate,
!                  except for target and limit points which lie between the two
!                  most recent continuation points.
!
!  RWORK(13)       Estimate of the pseudo-arclength coordinate of the current
!                  continuation point.
!
!  RWORK(14)       Estimate of the pseudoarclength coordinate of the current
!                  limit or target point, if any.
!
!  RWORK(15)       Size of the correction of the most recent continuation
!                  point; that is, the maximum norm of the distance between the
!                  predicted point and the accepted corrected point.
!
!  RWORK(16)       Estimate of the curvature between the last two
!                  continuation points.
!
!  RWORK(17)       Sign of the determinant of the augmented matrix at the
!                  last continuation point whose tangent vector has been
!                  calculated.
!
!  RWORK(18)       This quantity is only used if the jacobian matrix is to
!                  be estimated using finite differences.  In that case,
!                  this value determines the size of the increments and
!                  decrements made to the current solution values, according
!                  to the following formula:
!
!                    Delta X(J) = RWORK(18) * (1.0 + ABS(X(J))).
!
!                  The value of every entry of the approximate jacobian could
!                  be extremely sensitive to RWORK(18).  Obviously, no value
!                  is perfect.  Values too small will surely cause singular
!                  jacobians, and values too large will surely cause inaccuracy.
!                  Little more is certain.  However, for many purposes, it
!                  is suitable to set RWORK(18) to the square root of the
!                  machine epsilon, or perhaps to the third or fourth root,
!                  if singularity seems to be occuring.
!
!                  RWORK(18) defaults to SQRT(SQRT(RWORK(8))).
!
!                  The user may set RWORK(18).  If it has a nonzero value on
!                  input, that value will be used.  Otherwise, the default
!                  value will be used.
!
!  RWORK(19)       Not currently used.
!
!  RWORK(20)       Maximum growth factor for the predictor stepsize based
!                  on the previous secant stepsize.  The stepsize algorithm
!                  will produce a suggested step that is no less that the
!                  previous secant step divided by this factor, and no greater
!                  than the previous secant step multiplied by that factor.
!                  RWORK(20) defaults to 3.
!
!  RWORK(21)       The (Euclidean) secant distance between the last two
!                  computed continuation points.
!
!  RWORK(22)       The previous value of RWORK(21).
!
!  RWORK(23)       A number judging the quality of the Newton corrector
!                  convergence at the last continuation point.
!
!  RWORK(24)       Value of the component of the current tangent vector
!                  corresponding to the current continuation index.
!
!  RWORK(25)       Value of the component of the previous tangent vector
!                  corresponding to the current continuation index.
!
!  RWORK(26)       Value of the component of the current tangent vector
!                  corresponding to the limit index in IWORK(6).
!
!  RWORK(27)       Value of the component of the previous tangent vector
!                  corresponding to the limit index in IWORK(6).
!
!  RWORK(28)       Value of RWORK(7) when the last target point was
!                  computed.
!
!  RWORK(29)       Sign of the determinant of the augmented matrix at the
!                  previous continuation point whose tangent vector has been
!                  calculated.
!
!  RWORK(30)       through RWORK(30+4*NVAR-1) are used by the program to hold
!                  an old and new continuation point, a tangent vector and a
!                  work vector.  Subsequent entries of RWORK are used by the
!                  linear solver.
!
!
!  I) Programming Notes:
!  --------------------
!
!  The minimal input and user routines required to apply the program are:
!
!    Write a function routine FX of the form described above.
!    Use DENSLV as the linear equation solver by setting SLVNAM to DENSLV.
!    Skip writing a Jacobian routine by using the finite difference option.
!    Pass the name of FX as the Jacobian name as well.
!    Declare the name of the function FX as external.
!    Set NVAR in accordance with your problem.
!
!  Then:
!
!    Dimension the vector IWORK to the size LIW = 29+NVAR.
!    Dimension the vector RWORK to the size LRW = 29+NVAR*(NVAR+6).
!    Dimension IPAR(1) and FPAR(1) as required in the function routine.
!    Dimension XR(NVAR) and set it to an approximate solution of F(XR) = 0.
!
!  Set the work arrays as follows:
!
!    Initialize IWORK to 0 and RWORK to 0.0.
!
!    Set IWORK(1) = 0 (Problem startup)
!    Set IWORK(7) = 3 (Maximum internally generated output)
!    Set IWORK(9) = 1 (Forward difference Jacobian)
!
!  Now call the program repeatedly, and never change any of its arguments.
!  Check IERROR to decide whether the code is working satisfactorily.
!  Print out the vector XR to see the current solution point.
!
!
!  The most obvious input to try to set appropriately after some practice
!  would be the error tolerances RWORK(1) and RWORK(2), the minimum, maximum
!  and initial stepsizes RWORK(3), RWORK(4) and RWORK(5), and the initial
!  continuation index IWORK(2).
!
!  For speed and efficiency, a Jacobian routine should be written. It can be
!  checked by comparing its results with the finite difference Jacobian.
!
!  For a particular problem, the target and limit point input can be very
!  useful.  For instance, in the case of a discretized boundary value problem
!  with a real parameter it may be desirable to compare the computed solutions
!  for different discretization-dimensions and equal values of the parameter.
!  For this the target option can be used with the prescribed values of the
!  parameter. Limit points usually are of importance in connection with
!  stability considerations.
!
!
!  The routine REPS attempts to compute the machine precision, which in practice
!  is simply the smallest power of two that can be added to 1 to produce a
!  result greater than 1.  If the REPS routine misbehaves, you can replace
!  it by code that assigns a constant precomputed value, or by a call to
!  the PORT/SLATEC routine R1MACH.  REPS is called just once, in the
!  PITCON routine, and the value returned is stored into RWORK(8).
!
!  In subroutines DENSLV and BANSLV, the parameter "EPS" is set to RWORK(18)
!  and used in estimating jacobian matrices via finite differences.  If EPS is
!  too large, the jacobian will be very inaccurate.  Unfortunately, if EPS is
!  too small, the "finite" differences may become "infinitesmal" differences.
!  That is, the difference of two function values at close points may be
!  zero.  This is a very common problem, and occurs even with a function
!  like F(X) = X*X.  A singular jacobian is much worse than an inaccurate one,
!  so we have tried setting the default value of RWORK(18) to
!  SQRT(SQRT(RWORK(8)).  Such a value attempts to ward off singularity at the
!  expense of accuracy.  You may find for a particular problem or machine that
!  this value is too large and should be adjusted.  It is an utterly arbitrary
!  value.
!
!
!  J) References:
!  -------------
!
!  1.
!  Werner Rheinboldt,
!  Solution Field of Nonlinear Equations and Continuation Methods,
!  SIAM Journal of Numerical Analysis,
!  Volume 17, 1980, pages 221-237.
!
!  2.
!  Cor den Heijer and Werner Rheinboldt,
!  On Steplength Algorithms for a Class of Continuation Methods,
!  SIAM Journal of Numerical Analysis,
!  Volume 18, 1981, pages 925-947.
!
!  3.
!  Werner Rheinboldt,
!  Numerical Analysis of Parametrized Nonlinear Equations
!  John Wiley and Sons, New York, 1986
!
!  4.
!  Werner Rheinboldt and John Burkardt,
!  A Locally Parameterized Continuation Process,
!  ACM Transactions on Mathematical Software,
!  Volume 9, Number 2, June 1983, pages 215-235.
!
!  5.
!  Werner Rheinboldt and John Burkardt,
!  Algorithm 596, A Program for a Locally Parameterized Continuation Process,
!  ACM Transactions on Mathematical Software,
!  Volume 9, Number 2, June 1983, Pages 236-241.
!
!  6.
!  J J Dongarra, J R Bunch, C B Moler and G W Stewart,
!  LINPACK User's Guide,
!  Society for Industrial and Applied Mathematics,
!  Philadelphia, 1979.
!
!  7.
!  Richard Brent,
!  Algorithms for Minimization without Derivatives,
!  Prentice Hall, 1973.
!
!  8.
!  Tony Chan,
!  Deflated Decomposition of Solutions of Nearly Singular Systems,
!  Technical Report 225,
!  Computer Science Department,
!  Yale University,
!  New Haven, Connecticut, 06520,
!  1982.
!
!
!  L)  Sample programs:
!  -------------------
!
!
!  A number of sample problems are included with PITCON.  There are several
!  examples of the Freudenstein-Roth function, which is a nice example
!  with only a few variables, and nice whole number starting and stopping
!  points.  Other sample problems demonstrate or test various options in
!  PITCON.
!
!
!  1:
!  pcprb1.f
!  The Freudenstein-Roth function.  3 variables.
!
!  This is a simple problem, whose solution curve has some severe bends.
!  This file solves the problem with a minimum of fuss.
!
!
!  2:
!  pcprb2.f
!  The Aircraft Stability problem.  8 variables.
!
!  This is a mildly nonlinear problem, whose solution curve has some
!  limit points that are difficult to track.
!
!
!  3:
!  pcprb3.f
!  A boundary value problem.  22 variables.
!
!  This problem has a limit point in the LAMBDA parameter, which we seek.  We
!  solve this problem 6 times, illustrating the use of full and banded
!  jacobians, and of user-generated, or forward or central difference
!  approximated jacobian matrices.
!
!  The first 21 variables represent the values of a function along a grid
!  of 21 points between 0 and 1.  By increasing the number of grid points,
!  the problem can be set up with arbitrarily many variables.  This change
!  requires changing a single parameter in the main program.
!
!
!  4:
!  pcprb4.f
!  The Freudenstein Roth function.  3 variables.
!
!  This version of the problem tests the parameterization fixing option.
!  The user may demand that PITCON always use a given index for continuation,
!  rather than varying the index from step to step.  The Freudenstein-Roth
!  curve may be parameterized by the variable of index 2, although this
!  may increase the number of steps required to traverse the curve.
!
!  This file carries out 6 runs, with and without the parameterization fixed,
!  and with no limit point search, or with limit points sought for index 1
!  or for index 3.
!
!
!  5:
!  pcprb5.f
!  A boundary value problem.  21 variables.
!
!  This problem has a limit point in the LAMBDA parameter, which we do not
!  seek.  We do seek target points where LAMBDA = 0.8, which occurs twice,
!  before and after LAMBDA "goes around the bend".
!
!  We run this test 3 times, comparing the behavior of the full Newton
!  iteration, the Newton iteration that fixes the Jacobian during the
!  correction of each point, and the Newton iteration that fixes the
!  Jacobian as long as possible.
!
!
!  6:
!  pcprb6.f
!  Freudenstein-Roth function.  3 variables.
!
!  This version of the problem demonstrates the jacobian checking option.
!  Two runs are made.  Each is allowed only five steps.  The first
!  run is with the correct jacobian.  The second run uses a defective
!  jacobian, and demonstrates not only the jacobian checker, but also
!  shows that "slightly" bad jacobians can cause the Newton convergence
!  to become linear rather than quadratic.
!
!
!  7:
!  pcprb7.f
!  The materially nonlinear problem, up to 71 variables.
!