I am studying bifurcations of limit cycles in symmetric systems that exhibit branching at nonsimple singular points (e.g. where two real Floquet multipliers simultaneously pass through +1).
Neither the BP function nor the SPB function changes signs at such points, so I have been trying to define an overloaded parameter value to locate them. In particular, I use GETP('EIG',IND,U) in my PVALS functions to get the Floquet multipliers, and I set PAR(I) to the third smallest real eigenvalue sorted by distance to 1 in absolute value (since the first is trivially 1 and the second is already well detected by BP and LP functions). I then set UZR={I:1} to try to identify these special points where two real multipliers change stability.
The issue is that when the value of PAR(I) changes sign, AUTO enters the routine LCSPBV to locate the special point, but the Floquet multipliers are not recomputed between the steps in the Mueller iterations. Thus, the value of PAR(I) does not change, and the point is poorly detected.
I implemented a few changes to src/bvp.f90 (see the diff below) to resolve this for my case, and it works quite well as I can branch switch from those equivariant bifurcations that were previously undetected. However, I'm not sure if this resolves the issue for all problem types or breaks anything else. Also, I'm not sure if such technical details are still being maintained. If there is interest, I'd be happy to provide more information.
diff --git a/src/bvp.f90 b/src/bvp.f90
index 3062fab..6c59e24 100644
--- a/src/bvp.f90
+++ b/src/bvp.f90
@@ -280,7 +280,31 @@ CONTAINS
CALL PVLI(AP,ICP,UPS,NDIM,PAR,FNCI)
IF(IID.GE.2)WRITE(9,*)
ENDIF
- DO ITEST=1,AP%NTEST
+ ! ZGN: Check for bifurcations first, before new call to FNCI on 301
+ DO ITEST=NUZR+1,AP%NTEST
+ ! Check for special points
+ CALL LCSPBV(AP,DSOLD,DSTEST,PAR,ICP,ITEST,FUNI,BCNI,ICNI,FNCI, &
+ TEST(ITEST),RLCUR,RLOLD,RLDOT,NDIM,UPS,UOLDPS,UDOTPS, &
+ UPOLDP,TM,DTM,P0,P1,EV,THL,THU,IUZ,VUZ,NITPS,ATYPE,STEPPED, &
+ NTSTNA)
+ IF(STEPPED)ISTEPPED=ITEST
+ IF(LEN_TRIM(ATYPE)>0)THEN
+ IFOUND=ITEST
+ AP%ITP=LBITP(ATYPE,.TRUE.)
+ AP%ITP=AP%ITP+SIGN(10,AP%ITP)*ITPST
+ ENDIF
+ ENDDO
+
+ ! ZGN: Before checking UZRS, recalculate Floquet
+ IF(IAM==0)THEN
+ AP%IID=0 ! Disable repeated print
+ SP1 = FNCI(AP,ICP,UPS,NDIM,PAR,3,ATYPEDUM) ! New FNCI call
+ AP%IID=IID ! Restore IID
+ CALL PVLI(AP,ICP,UPS,NDIM,PAR,FNCI)
+ IF(IID.GE.2)WRITE(9,*)
+ ENDIF
+ ! ZGN: Now check UZR
+ DO ITEST=1,NUZR
! Check for special points
CALL LCSPBV(AP,DSOLD,DSTEST,PAR,ICP,ITEST,FUNI,BCNI,ICNI,FNCI, &
TEST(ITEST),RLCUR,RLOLD,RLDOT,NDIM,UPS,UOLDPS,UDOTPS, &
@@ -1264,7 +1288,7 @@ CONTAINS
INTEGER I,IAM,IID,ITMX,IBR,NTOT,NTOP,NITSP1,ISTOP,NCOL,NFPR,NITPSS
LOGICAL SEARCH
- DOUBLE PRECISION DS,DSMAX,EPSS,Q0,Q1,DQ,RDS,RRDS,S0,S1
+ DOUBLE PRECISION DS,DSMAX,EPSS,Q0,Q1,DQ,RDS,RRDS,S0,S1,SP1
DOUBLE PRECISION DETS,FLDFS,DSOLDS,DSTESTS
CHARACTER(4) :: ATYPEDUM
@@ -1399,7 +1423,8 @@ CONTAINS
ENDIF
! Check for zero.
-
+ !ZGN: Recalculate FLOWKM before PVLI
+ SP1 = FNCI(AP,ICP,UPS,NDIM,PAR,3,ATYPEDUM) ! New FNCI call
CALL PVLI(AP,ICP,UPS,NDIM,PAR,FNCI)
IF(IID.GE.2)WRITE(9,*)
Q=FNCS(AP,ICP,UPS,PAR,ATYPE,IUZ,VUZ,ITEST,FNCI)