C PCPRB4.FOR 27 February 1991 C PROGRAM PCPRB4 C C Ermentrout-Aronson-Kopell problem. C C The function F(X) is the right hand side of an ODE, and is of interest C in studying the equilibrium behavior of the system. C C F(X) is of the form: C C FX(1)=R*(1-GAMMA-R*R)+S*(-MU*SIN(PHI)+GAMMA*COS(PHI)) C FX(2)=S*(1-GAMMA-S*S)+R*( MU*SIN(PHI)+GAMMA*COS(PHI)) C FX(3)=Q*(R*R-S*S)-GAMMA*SIN(PHI)*(R*R+S*S)/(R*S) C +MU*COS(PHI)*(R*R-S*S)/(R*S) C C GAMMA is the parameter of interest. C Q may be fixed at -2.0, MU at 1.0, for this problem. C R, S and PHI are the fundamental variables. C C One family of solutions is: C (1, 1, 0, GAMMA) which becomes unstable as GAMMA increases from zero. C Two branches of asymmetric solutions diverge from this solution at one C point. C C Another family is: C (1-2*GAMMA, 1-2*GAMMA, PI, GAMMA) C INTEGER LIW INTEGER LRW INTEGER NVAR C PARAMETER (NVAR=4) PARAMETER (LIW=NVAR+29) PARAMETER (LRW=29+(6+NVAR)*NVAR) C EXTERNAL DENSLV EXTERNAL FXBARD EXTERNAL FPBARD EXTERNAL PITCON C REAL FPAR(2) INTEGER I INTEGER IERROR INTEGER IPAR(1) INTEGER IWORK(LIW) INTEGER J INTEGER LOUNIT CHARACTER NAME*12 REAL PARMU REAL PARQ REAL RWORK(LRW) REAL XR(NVAR) C LOUNIT=6 C C Set parameters C DO 4 I=1,LIW IWORK(I)=0 4 CONTINUE DO 5 I=1,LRW RWORK(I)=0.0 5 CONTINUE PARQ=-2.0 PARMU=1.0 FPAR(1)=PARQ FPAR(2)=PARMU C C Set some entries of work arrays. C C IWORK(1)=0 ; This is a startup C IWORK(2)=4 ; Use X(4) for initial parameter C IWORK(3)=0 ; Program may choose parameter index C IWORK(4)=0 ; Update jacobian every newton step C IWORK(5)=1 ; Seek target values for X(1) C IWORK(6)=0 ; Seek no limit points. 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)=4 IWORK(3)=0 IWORK(4)=0 IWORK(5)=1 IWORK(6)=0 IWORK(7)=1 IWORK(8)=LOUNIT IWORK(9)=0 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.10 ; Maximum stepsize C RWORK(5)=0.02 ; 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.0001 RWORK(2)=0.0001 RWORK(3)=0.01 RWORK(4)=0.10 RWORK(5)=0.02 RWORK(6)=1.0 RWORK(7)=1.0 C C Set starting point: C XR(1)=0.71055 XR(2)=1.05053 XR(3)=0.46070 XR(4)=0.5 C WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'PCPRB4' WRITE(LOUNIT,*)'PITCON sample program.' WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'Ermentrout-Aronson-Kopell function' WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'Step Type of Point '// *'X(1) X(2) X(3) X(4)' WRITE(LOUNIT,*)' ' I=0 NAME='Start point ' WRITE(LOUNIT,'(1X,I3,2X,A12,2X,4G14.6)')I,NAME,(XR(J),J=1,NVAR) C DO 10 I=1,30 CALL PITCON(FPBARD,FPAR,FXBARD,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,4G14.6)')I,NAME,(XR(J),J=1,NVAR) C IF(IWORK(1).EQ.3)THEN WRITE(LOUNIT,*)'We have reached the desired point.' 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) STOP END SUBROUTINE FXBARD(NVAR,FPAR,IPAR,X,FX,IERROR) C C*********************************************************************** C C Evaluate the nonlinear function C INTRINSIC COS INTRINSIC SIN C INTEGER NVAR C REAL COS REAL FPAR(2) REAL FX(NVAR) REAL GAMMA INTEGER IERROR INTEGER IPAR(*) REAL PARMU REAL PARQ REAL PHI REAL R REAL S REAL SIN REAL X(NVAR) C PARQ=FPAR(1) PARMU=FPAR(2) R=X(1) S=X(2) PHI=X(3) GAMMA=X(4) C FX(1)=R*(1.0-GAMMA-R*R) + S*(-PARMU*SIN(PHI)+GAMMA*COS(PHI)) FX(2)=S*(1.0-GAMMA-S*S) + R*( PARMU*SIN(PHI)+GAMMA*COS(PHI)) FX(3)=PARQ*(R*R-S*S)-GAMMA*SIN(PHI)*(R*R+S*S)/(R*S) * +PARMU*COS(PHI)*(R*R-S*S)/(R*S) RETURN END SUBROUTINE FPBARD(NVAR,FPAR,IPAR,X,FPRIME,IERROR) C C*********************************************************************** C C Evaluate the jacobian of the nonlinear function. C INTRINSIC COS INTRINSIC SIN C INTEGER NVAR C REAL COS REAL FPAR(2) REAL FPRIME(4,4) REAL GAMMA INTEGER IERROR INTEGER IPAR(1) REAL PARMU REAL PARQ REAL PHI REAL R REAL S REAL SIN REAL X(NVAR) C PARQ=FPAR(1) PARMU=FPAR(2) R=X(1) S=X(2) PHI=X(3) GAMMA=X(4) C FPRIME(1,1)=1.0-GAMMA-3.0*R*R FPRIME(1,2)=-PARMU*SIN(PHI)+GAMMA*COS(PHI) FPRIME(1,3)=S*(-PARMU*COS(PHI)-GAMMA*SIN(PHI)) FPRIME(1,4)=-R+S*COS(PHI) FPRIME(2,1)=(PARMU*SIN(PHI)+GAMMA*COS(PHI)) FPRIME(2,2)=1.0-GAMMA-3.0*S*S FPRIME(2,3)=R*( PARMU*COS(PHI)-GAMMA*SIN(PHI)) FPRIME(2,4)=-S+R*COS(PHI) FPRIME(3,1)=2.0*PARQ*R - GAMMA*SIN(PHI)*(R*R-S*S)/(R*R*S) * +PARMU*COS(PHI)*(R*R+S*S)/(R*R*S) FPRIME(3,2)=-2.0*PARQ*S - GAMMA*SIN(PHI)*(S*S-R*R)/(S*S*R) * -PARMU*COS(PHI)*(R*R+S*S)/(S*S*R) FPRIME(3,3)=-GAMMA*COS(PHI)*(R*R+S*S)/(R*S) * -PARMU*SIN(PHI)*(R*R-S*S)/(R*S) FPRIME(3,4)=-SIN(PHI)*(R*R+S*S)/(R*S) RETURN END