source('iterate_logistic.R')

# Part (a): Generate solutions with the requested values for rho:
#
N = 400
u0 = 0.2
rhos = c( 2.9, 2.95, 2.99 )
for( ri in 1:length(rhos) ){
    rho = rhos[ri]
    u = iterate_logistic(rho, u0, N)
    if( ri==1 ){
        sols = u
    }else{
        sols = rbind( sols, u )
    }
}
rownames(sols) = c('sol1', 'sol2', 'sol3')

#postscript("../../WriteUp/Graphics/Chapter2/chap_2_sect_9_prob_17_plot_part_a.eps", onefile=FALSE, horizontal=FALSE)
matplot(0:(N-1), t(sols), type='l', xlab='iterate', ylab='u_n')
grid()
#dev.off()

# Compute the limiting value for each value of rho:
#
lvs = ( rhos - 1 ) / rhos
print(lvs)

# Find the iterate at which each solution is first 1e-3 close to the limiting value:
#
first_iterate = c()
for( ri in 1:length(rhos) ){
    lv = lvs[ri]
    close = abs( sols[ri,] - lv ) < 1.e-3
    n = which(close)[1] - 1
    first_iterate = c( first_iterate, n )
}
print('first iterate that is close to the limiting value')
print(first_iterate)


# Part (b):
#
N = 100
u0 = 0.2
rhos = c( 3.01, 3.05, 3.1 )
for( ri in 1:length(rhos) ){
    rho = rhos[ri]
    u = iterate_logistic(rho, u0, N)
    if( ri==1 ){
        sols = u
    }else{
        sols = rbind( sols, u )
    }
}
rownames(sols) = c('sol1', 'sol2', 'sol3')

#postscript("../../WriteUp/Graphics/Chapter2/chap_2_sect_9_prob_17_plot_part_b.eps", onefile=FALSE, horizontal=FALSE)
matplot(0:(N-1), t(sols), type='l', xlab='iterate', ylab='u_n')
grid()
#dev.off()


# Compute the limiting value for each value of rho:
#
# Could solve the system:
#
# u_2 = rho u_1 (1-u_1)
# u_1 = rho u_2 (1-u_2)
#
# but that's what the above fixed points are finding thus we can just grab the numbers above.
#
fixed_points = rbind( sols[ , dim(sols)[2] ], sols[ , dim(sols)[2]-1 ] )
print('limit cycle fixed points for each solution:')
print(fixed_points)


# Find the iterate at which each solution is first 1e-3 close to these limiting value:
#
first_iterate = c()
for( ri in 1:length(rhos) ){
    u1 = fixed_points[1, ri]
    u2 = fixed_points[2, ri]
    for( ti in 1:(N-1) ){
        if( ( abs( sols[ri,ti]-u1 ) < 1.e-3 ) & ( abs( sols[ri,ti+1]-u2 ) < 1.e-3 ) ){
            first_iterate = c( first_iterate, ti )
            break
        }
    }
}
print('first iterate that is close to the limiting value')
print(first_iterate)