Empirical Feigenbaum Diagram in Sine Map: Coding with R

1 Objectives:

• Modeling the Pitchfork bifurcation of Sine Map with R
• constructing Feigenbaum Diagram in Sine Map
• numerical verification of Feigenbaum constant \[\delta = \frac{r_{n}-r_{n-1}}{r_{n+1}-r_{n}} = 4.669.... \]

2 Sine Map Model:

• r: bifurcation parameter
• x0: initial value
• n: number of iteration
• M: number of latest iteration points used for plot

\[x_{n+1} = f(x_{n}) = r sin(\pi x_{n}) \],\[r \in \left[0, 1 \right]\]and\[0\leq x_{0}\leq 1 \]

3 Period: n=60

nperiod=60 
r.change = seq(0, 1, by=0.0001)
Orbit = sapply(r.change, sine.map,  x0=0.01, n= nperiod, M= 55)
Orbit_dim = (Orbit/Orbit)
r_dim = t(t(Orbit_dim)*r.change)
mycolor = wes_palette("Darjeeling1", dim(Orbit_dim)[1], type = c( "continuous"))
plot(r_dim, Orbit, xlim = c(0.7, 1.0), ylim = c(0, 1), pch = ".", ylab = "x", 
     xlab = "r",  cex=1.5, col= mycolor , lwd=1)
title(main= paste("Sine Map: nth period = ", nperiod, sep=" "), col.main="red")
grid()

4 Period: n=1000, r.change = seq(0, 1, by=0.0001)

• Feigenbaum constant is the limiting ratio of each bifurcation interval to the next between every period doubling: Numerical verification from using my R code programming:

nperiod=1000 
r.change = seq(0, 1, by=0.0001)
Orbit = sapply(r.change, sine.map,  x0=0.01, n= nperiod, M= 300)
Orbit_dim = (Orbit/Orbit)
r_dim = t(t(Orbit_dim)*r.change)
mycolor = wes_palette("Darjeeling1" , dim(Orbit_dim)[1], type = "continuous")
plot(r_dim, Orbit, xlim = c(0.7, 1.0),ylim = c(0, 1), pch = ".", 
     ylab = "x", xlab = "r",  cex=1.5, col= mycolor , lwd=1)
title(main= paste("Sine Map: nth period = ", nperiod ,  sep=" "), col.main="red")
grid()

Period_1_r = 0.7179
unique((round(Orbit[, which(r.change ==  Period_1_r)],2)))

## [1] 0.64

Period_2_r = 0.8324
unique((round(Orbit[, which(r.change ==  Period_2_r)],2)))

## [1] 0.82 0.45

Period_4_r = 0.8568
unique((round(Orbit[, which(r.change ==  Period_4_r)],2)))

## [1] 0.79 0.52 0.86 0.38

Period_8_r = 0.8620
unique((round(Orbit[, which(r.change ==  Period_8_r)],2)))

## [1] 0.80 0.50 0.86 0.36 0.78 0.54 0.85 0.38

• Theoretical Feigenbaum constant in Sine Map based on nperiod = 1000 \[\delta = \frac{r_{n}-r_{n-1}}{r_{n+1}-r_{n}} = 4.669.... \]

deltaFC_1 = (Period_2_r - Period_1_r)/ (Period_4_r - Period_2_r )
deltaFC_1

## [1] 4.692623

deltaFC_2 = (Period_4_r - Period_2_r)/ (Period_8_r - Period_4_r )
deltaFC_2 

## [1] 4.692308

5 Conclusion :

• The empirical Feigenbaum constant estimated to be \(\delta\) = 4.692… which is very closed to the theoretical one: 4.669…. . The small computational discrepancy is mainly due to 1) the how small the digitization of r interval grid 2) increasing the nth period to 2000 . It is very computational expensive to re-fine the r interval grid to be order of \(10^{-5}\) or less and to increase the nth period iterations simultaneously.

6 Appendix: Sine Map nperiod = 50 : Online Interactive (drag me with your mouse to see)

nperiod=50
r.change = seq(0.1, 1, by=0.0001)
Orbit = sapply(r.change, sine.map,  x0=0.01, n= nperiod, M= 10)
Orbit_dim = (Orbit/Orbit)
r_dim = t(t(Orbit_dim)*r.change)
mycolor = wes_palette("Darjeeling1" , dim(Orbit_dim)[1], type = "continuous")
mydata = data.frame(melt(r_dim), melt(Orbit) )
mydata070 = mydata[which(mydata[,"value"] >= 0.68),]
plot_ly(data = mydata070,  x= mydata070[,"value"], y = mydata070[,"value.1"], 
        type = 'scatter', marker = list(size = 2, color = 'rgba(255, 152, 143, 1)'))%>%
layout(title = 'Sine Map nperiod = 50 :
       Interactive (drag me with your mouse to see)' )