C PCPRB5.FOR 27 February 1991 C PROGRAM PCPRB5 C C Discretized two point boundary value problem with parameter LAMBDA: C C Y''+LAMBDA*EXP(Y)=0 C Y(0)=0 Y'(1)=0 C C Expecting limit point in LAMBDA at roughly LAMBDA=0.878. C C This test illustrates the use of the full Newton, modified Newton, C and "cheap" Newton methods, which update the Jacobian on every iteration, C or only once per point, or only when convergence fails. C INTEGER LIW INTEGER LRW INTEGER NVAR C PARAMETER (NVAR=22) PARAMETER (LRW=29+(NVAR+6)*NVAR) PARAMETER (LIW=NVAR+29) C EXTERNAL DENSLV EXTERNAL FPFULL EXTERNAL FXBEND EXTERNAL PITCON C REAL FPAR(1) INTEGER I INTEGER IERROR INTEGER IPAR(2) INTEGER ITRY INTEGER IWORK(LIW) INTEGER J INTEGER LOUNIT INTEGER MODCON CHARACTER NAME*12 INTEGER NIT REAL RWORK(LRW) REAL XR(NVAR) C LOUNIT=6 ITRY=0 C WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'PCPRB5' WRITE(LOUNIT,*)'PITCON sample program.' WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'Two point BVP, limit point in Lambda' WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'This program will be run THREE times.' C 19 CONTINUE C ITRY=ITRY+1 WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'This is run number ',ITRY NIT=0 C WRITE(LOUNIT,*)' ' DO 4 I=1,LIW IWORK(I)=0 4 CONTINUE DO 5 I=1,LRW RWORK(I)=0.0 5 CONTINUE IERROR=0 C C Set input quantities C IF(ITRY.EQ.1)MODCON=0 IF(ITRY.EQ.2)MODCON=1 IF(ITRY.EQ.3)MODCON=2 C C IWORK(1)=0 ; This is a startup C IWORK(2)=NVAR ; Use X(NVAR) for initial parameter C IWORK(3)=0 ; Program may choose parameter index C IWORK(4)=MODCON ; Control frequency of jacobian update in newton iteration C IWORK(5)=NVAR ; Seek target values for X(NVAR) C IWORK(6)=NVAR ; Seek limit points in X(NVAR) C IWORK(7)=1 ; Control amount of output. C IWORK(8)=6 ; Output unit C IWORK(9)=0 ; Jacobian choice. C IWORK(1)=0 IWORK(2)=NVAR IWORK(3)=0 IWORK(4)=MODCON IWORK(5)=NVAR IWORK(6)=NVAR IWORK(7)=1 IWORK(8)=LOUNIT IWORK(9)=0 C C Jacobian upper/lower bandwidths: C IPAR(1)=1 IPAR(2)=1 C C RWORK(1)=0.0001; Absolute error tolerance C RWORK(2)=0.0001; Relative error tolerance C RWORK(3)=0.01 ; Minimum stepsize C RWORK(4)=0.25 ; Maximum stepsize C RWORK(5)=0.05 ; Starting stepsize C RWORK(6)=1.0 ; Starting direction C RWORK(7)=0.80 ; Target value (Seek solution with X(NVAR)=0.80 C RWORK(1)=0.0001 RWORK(2)=0.0001 RWORK(3)=0.01 RWORK(4)=0.25 RWORK(5)=0.05 RWORK(6)=1.0 RWORK(7)=0.80 C C Set starting point C DO 55 I=1,NVAR XR(I)=0.0 55 CONTINUE C IF(ITRY.EQ.1)THEN WRITE(LOUNIT,*)'Using full Newton method.' WRITE(LOUNIT,*)'Jacobian updated on every iteration.' ELSEIF(ITRY.EQ.2)THEN WRITE(LOUNIT,*)'Using modified Newton method.' WRITE(LOUNIT,*)'Jacobian held fixed while correcting a point.' ELSE WRITE(LOUNIT,*)'Using "cheap" Newton method.' WRITE(LOUNIT,*)'Jacobian held fixed as long as possible,' WRITE(LOUNIT,*)'perhaps over multiple points, until' WRITE(LOUNIT,*)'convergence fails.' ENDIF WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'Step Type of Point Lambda' WRITE(LOUNIT,*)' ' I=0 NAME='Start point ' WRITE(LOUNIT,'(1X,I3,2X,A12,2X,G14.6)')I,NAME,XR(NVAR) C DO 10 I=1,50 CALL PITCON(FPFULL,FPAR,FXBEND,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(NVAR) C IF(IWORK(1).EQ.3)THEN NIT=NIT+1 WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'Complete value for target point ',NIT WRITE(LOUNIT,*)' ' WRITE(LOUNIT,'(1X,5G14.6)')(XR(J),J=1,NVAR) IF(NIT.GE.2)GO TO 20 ENDIF C IF(IERROR.NE.0)THEN WRITE(LOUNIT,*)'Iteration halted, IERROR=',IERROR GO TO 20 ENDIF C 10 CONTINUE C 20 CONTINUE WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'Jacobians ',IWORK(19) WRITE(LOUNIT,*)'Factorizations ',IWORK(20) WRITE(LOUNIT,*)'Solves ',IWORK(21) WRITE(LOUNIT,*)'Functions ',IWORK(22) IF(ITRY.LT.3)GO TO 19 STOP END SUBROUTINE FXBEND(NVAR,FPAR,IPAR,X,FX,IERROR) C C*********************************************************************** C C Function evaluator for two point boundary value problem. C INTRINSIC EXP C INTEGER NVAR C REAL EXP REAL FPAR(1) REAL FX(NVAR) REAL HSQ INTEGER I INTEGER IERROR INTEGER IPAR(2) REAL X(NVAR) C HSQ=NVAR-2 HSQ=1.0/HSQ HSQ=HSQ*HSQ FX(1)=X(1) DO 10 I=2,NVAR-2 FX(I)=X(I-1)-2.0*X(I)+X(I+1)+X(NVAR)*EXP(X(I))*HSQ 10 CONTINUE FX(NVAR-1)=X(NVAR-1)-X(NVAR-2) RETURN END SUBROUTINE FPFULL(NVAR,FPAR,IPAR,X,FPRIME,IERROR) C C*********************************************************************** C C Full storage jacobian of function C INTRINSIC EXP C INTEGER NVAR C REAL EXP REAL FPAR(1) REAL FPRIME(NVAR,NVAR) REAL HSQ INTEGER I INTEGER IERROR INTEGER IPAR(2) REAL X(NVAR) C HSQ=NVAR-2 HSQ=1.0/HSQ HSQ=HSQ*HSQ FPRIME(1,1)=1.0 DO 10 I=2,NVAR-2 FPRIME(I,I-1)=1.0 FPRIME(I,I)=-2.0+X(NVAR)*EXP(X(I))*HSQ FPRIME(I,I+1)=1.0 FPRIME(I,NVAR)=EXP(X(I))*HSQ 10 CONTINUE FPRIME(NVAR-1,NVAR-2)=-1.0 FPRIME(NVAR-1,NVAR-1)=1.0 RETURN END