C PCPRB8.FOR 27 February 1991 C PROGRAM PCPRB8 C C The Freudenstein-Roth function. C C Reference C C F Freudenstein, B Roth, C Numerical Solutions of Nonlinear Equations, C Journal of the Association for Computing Machinery, C Volume 10, 1963, Pages 550-556. C C The function F(X) is of the form C C FX(1) = X1 - X2**3 + 5*X2**2 - 2*X2 - 13 + 34*(X3-1) C FX(2) = X1 + X2**3 + X2**2 - 14*X2 - 29 + 10*(X3-1) C C Starting from the point (15,-2,0), the program is required to produce C solution points along the curve until it reaches a solution point C (*,*,1). It also may be requested to look for limit points in the C first or third components. C C The correct value of the solution at X3=1 is (5,4,1). C C Limit points in the first variable occur at: C C (14.28309, -1.741377, 0.2585779) C (61.66936, 1.983801, -0.6638797) C C Limit points for the third variable occur at: C C (20.48586, -0.8968053, 0.5875873) C (61.02031, 2.230139, -0.6863528) C INTEGER LIW INTEGER LRW INTEGER NVAR C PARAMETER (NVAR=3) PARAMETER (LIW=NVAR+29) PARAMETER (LRW=29+(6+NVAR)*NVAR) C EXTERNAL DENSLV EXTERNAL DFROTH EXTERNAL FXROTH EXTERNAL PITCON C REAL FPAR(1) INTEGER I INTEGER IERROR INTEGER IPAR(1) INTEGER ITRY INTEGER IWORK(LIW) INTEGER J INTEGER LOUNIT CHARACTER NAME*12 REAL RWORK(LRW) REAL XR(NVAR) C LOUNIT=6 C WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'PCPRB8:' WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'PITCON sample program.' WRITE(LOUNIT,*)'Freudenstein-Roth function' WRITE(LOUNIT,*)'2 equations, 3 variables.' C C Carry out continuation for each method of approximating the jacobian: C ITRY=0 99 CONTINUE C C Set work arrays to zero: C DO 10 I=1,LIW IWORK(I)=0 10 CONTINUE DO 20 I=1,LRW RWORK(I)=0.0 20 CONTINUE C C Set some entries of work arrays. C C IWORK(1)=0 ; This is a startup C IWORK(2)=2 ; Use X(2) for initial parameter C IWORK(3)=0 ; Program may choose parameter index C IWORK(4)=0 ; Update jacobian every newton step C IWORK(5)=3 ; Seek target values for X(3) C IWORK(6)=1 ; Seek limit points in X(1) C IWORK(7)=1 ; Control amount of output. C IWORK(8)=6 ; Output unit C IWORK(9)=* ; Jacobian choice. Will vary from 0, 1, 2. C 0=Use user's jacobian routine C 1=Forward difference approximation C 2=Central difference approximation C IWORK(1)=0 IWORK(2)=2 IWORK(3)=0 IWORK(4)=0 IWORK(5)=3 IWORK(6)=1 IWORK(7)=1 IWORK(8)=LOUNIT IWORK(9)=ITRY C C RWORK(1)=0.00001; Absolute error tolerance C RWORK(2)=0.00001; Relative error tolerance C RWORK(3)=0.01 ; Minimum stepsize C RWORK(4)=20.0 ; Maximum stepsize C RWORK(5)=0.3 ; Starting stepsize C RWORK(6)=1.0 ; Starting direction C RWORK(7)=1.0 ; Target value (Seek solution with X(3)=1) C RWORK(1)=0.00001 RWORK(2)=0.00001 RWORK(3)=0.01 RWORK(4)=20.0 RWORK(5)=0.3 RWORK(6)=1.0 RWORK(7)=1.0 C C Set starting point. C XR(1)=15.0 XR(2)=-2.0 XR(3)=0.0 C WRITE(LOUNIT,*)' ' IF(ITRY.EQ.0)THEN WRITE(LOUNIT,*)'Use user jacobian routine.' ELSEIF(ITRY.EQ.1)THEN WRITE(LOUNIT,*)'Approximate jacobian with forward difference.' ELSEIF(ITRY.EQ.2)THEN WRITE(LOUNIT,*)'Approximate jacobian with central difference.' ENDIF WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'Step Type of Point '// *'X(1) X(2) X(3)' WRITE(LOUNIT,*)' ' I=0 NAME='Start point ' WRITE(LOUNIT,'(1X,I3,2X,A12,2X,3G14.6)')I,NAME,(XR(J),J=1,NVAR) C C Take another step. C DO 30 I=1,30 C CALL PITCON(DFROTH,FPAR,FXROTH,IERROR,IPAR,IWORK,LIW, * NVAR,RWORK,LRW,XR,DENSLV) C IF(IWORK(1).EQ.1)THEN NAME='Corrected ' ELSEIF(IWORK(1).EQ.2)THEN NAME='Continuation ' ELSEIF(IWORK(1).EQ.3)THEN NAME='Target point ' ELSEIF(IWORK(1).EQ.4)THEN NAME='Limit point ' ENDIF C WRITE(LOUNIT,'(1X,I3,2X,A12,2X,3G14.6)')I,NAME,(XR(J),J=1,NVAR) C IF(IWORK(1).EQ.3)THEN WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'We have reached the point we wanted.' GO TO 40 ENDIF C IF(IERROR.NE.0)THEN WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'An error occurred.' WRITE(LOUNIT,*)'We will terminate the computation now.' GO TO 40 ENDIF C 30 CONTINUE C 40 CONTINUE WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'Jacobians ',IWORK(19) WRITE(LOUNIT,*)'Factorizations ',IWORK(20) WRITE(LOUNIT,*)'Solves ',IWORK(21) WRITE(LOUNIT,*)'Functions ',IWORK(22) ITRY=ITRY+1 IF(ITRY.LE.2)GO TO 99 STOP END SUBROUTINE FXROTH(NVAR,FPAR,IPAR,X,F,IERROR) C C Evaluate the function at X. C C ( X1 - ((X2-5.0)*X2+2.0)*X2 - 13.0 + 34.0*(X3-1.0) ) C ( X1 + ((X2+1.0)*X2-14.0)*X2 - 29.0 + 10.0*(X3-1.0) ) C INTEGER NVAR C REAL F(*) REAL FPAR(*) INTEGER IERROR INTEGER IPAR(*) REAL X(NVAR) C F(1)=X(1) * -((X(2)-5.0)*X(2)+2.0)*X(2) * -13.0 * +34.0*(X(3)-1.0) C F(2)=X(1) * +((X(2)+1.0)*X(2)-14.0)*X(2) * -29.0 * +10.0*(X(3)-1.0) C RETURN END SUBROUTINE DFROTH(NVAR,FPAR,IPAR,X,FJAC,IERROR) C C Evaluate the Jacobian: C C ( 1.0 (-3.0*X(2)+10.0)*X(2)-2.0 34.0 ) C ( 1.0 (3.0*X(2)+2.0)*X(2)-14.0 10.0 ) C INTEGER NVAR C REAL FJAC(NVAR,NVAR) REAL FPAR(*) INTEGER IERROR INTEGER IPAR(*) REAL X(NVAR) C FJAC(1,1)=1.0 FJAC(1,2)=(-3.0*X(2)+10.0)*X(2)-2.0 FJAC(1,3)=34.0 C FJAC(2,1)=1.0 FJAC(2,2)=(3.0*X(2)+2.0)*X(2)-14.0 FJAC(2,3)=10.0 C RETURN END