# Contributions from Stephanie Perez, Jessica Phung, and Imani Chatman

#1a
f(x)=2.5*x*(1-x) #This line gives the function with r=2.5.
count=srange(0,20) #This line gives the range of numbers and thus the amount of iterations we want.
initial=[0.9] #This gives the initial value of 0.9 we chose.
for i in count: #This creates the for loop.
    new=f(initial[i]) #This changes every new output using the previous output as an input. (whatever 'initial' currently is, it will be input into                        the function f(x) and the [i] indicates that it will be the 1st, 2nd, 3rd, etc. iteration of that). The output will be                              called 'new'.
    initial.append(new) #This adds the new values to the list, making whatever 'new' is, the new current value of 'initial'.
print(initial) #This prints the list.
list_plot(initial, plotjoined=True, axes_labels=("time", "f(x)")) # this plots the lis
#This shows equilibrium behavior after time 5.

#1b
f(x)=3.1*x*(1-x) #Provides the function with r=3.1
count2=srange(0,20) #This line states the amount of iterations by setting a range
initial2=[0.9] #States the initial condition
for i in count:# Establishes the for loop
    new2=f(initial2[i])# Using old output as the new input
    initial2.append(new2) #Places all the values into a list
print(initial2)# Prints the list
list_plot(initial2, plotjoined=True, axes_labels=["time", "f(x)"])
#This time series demonstrates unstable oscillatory behavior.

#1c
f(x)=3.5*x*(1-x) # given function with r=3.5
count3=srange(0,20) # establishes the amount of iterations that we want
initial3=[0.9] # creates a list called initial3 where the initial condition is 0.9
for i in count3: # for every integer in the 'count3' list...
    new3=f(initial3[i]) # we will run initial3 condition=0.9 through the f(x) function and call the result 'new3'
    initial3.append(new3) # every 'new3' result will be appended to the 'initial3' list and become the new initial condition value
print(initial3) # prints out the list
list_plot(initial3, plotjoined=True, axes_labels=["time", "f(x)"]) # plots everything on 'initial3'
# We are observing chaotic behavior in this time series.

#1d
f(x)=3.9*x*(1-x) # given function with r=3.9
count4=srange(0,20) # establishes the amount of iterations that we want
initial4=[0.9] # creates a list called initial4 where the initial condition is 0.9
for i in count4: # for 20 iterations, we will...
    new4=f(initial4[i]) # we will run whatever the current 'initial4' is through the f(x) function and call the result 'new4'
    initial4.append(new4) # every 'new4' result will be appended to the 'initial4' list and become the new initial condition
print(initial4) # prints out the list
list_plot(initial4, plotjoined=True, axes_labels=["time", "f(x)"]) # plots everything on 'initial4'
# We are observing more apparent chaotic behavior in this time series.

#2
var("r") #Declares "r" as a symbolic symbol
r_values=[2.5, 3.1, 3.5, 3.9] #Creates a list for "r" values so we can run them all at the same time
count= srange(0,20) #Gives a range for the amount of interations
for r in r_values:# Establishes a for loop where the actions listed below will be done for every r value in the list above
    newg=r*(1-x)*(x) # States the function
    init=[0.9] #Makes the list "init" with an initial condition of 0.9 for x
    for i in count: #Creates a for loop where for 20 iterations...
        newg2=newg(init[i]) #Places all new values into a list... we will plug whatever current value 'init' is into the function 'newg' 20x
        init.append(newg2) #Append all new values in order to create a list, becomes new 'init'
    show(list_plot(init, plotjoined=True, axes_labels=("time", "different r values"))) # indent one time only or else it will fall under both for                                                                                            loops and run 80x
# Here, we observe equilibrium, unstable (increasing) oscillations, chaos, and chaos, as before.

#3
h(x)=3.9*x*(1-x) # discrete logistic function w/ r=3.9 and ic=0.4
count=srange(0,20) # we want to run it 20x
newinitial=[0.4] # initial condition so that new current conditions can be added
for i in count: # for 20 iterations, we will...
    n=h(newinitial[i]) #...plug whatever 'newinitial' currently is into h(x) and call it 'n'
    newinitial.append(n) # append 'n' to the 'newinitial' list and it becomes the new current condition
print(newinitial) # prints 'newinitial' list

# Basically the same exact function, except with a different ics and we're making it green on the time series so we can differentiate
g(x)=3.9*x*(1-x) # discrete logistic function w/ r=3.9 and ic=0.401
count=srange(0,20)
newinitial1=[0.401] 
for i in count: 
    n=g(newinitial1[i]) 
    newinitial1.append(n)
print(newinitial1)

list_plot(newinitial, plotjoined=True, axes_labels=("time", "h(x)+g(x)"), legend_label=("h(x)"))+list_plot(newinitial1, plotjoined=True, color="green", legend_label=("g(x)")) # plots both of them together overlaying

# The relationship between the time series starts off similarly to each other and following the same pattern, but over time becomes chaotic and has different bifurcations. This is an example of the butterfly effect.

#4a
a(x)=3.7*x*(1-x) # discrete logistic function w/ r=3.7 and ic=0.7
count=srange(0,20)
newinit1=[0.7]
for i in count: # for 20 iterations...
    egg=a(newinit1[i]) # whatever 'newinit1' currently is will be plugged into a(x) and be called 'egg'
    newinit1.append(egg) # whatever 'egg' is will be appended to 'newinit1' and become the new current value3
print(newinit1)
b(x)=3.7*x*(1-x) # discrete logistic function w/ r=3.7 and ic=0.701
count=srange(0,20)
newinit2=[0.701]
for i in count:
    nice=b(newinit2[i])
    newinit2.append(nice)
print(newinit2)
list_plot(newinit1, plotjoined=True, axes_labels=("time", "a(x)+b(x)"), legend_label=("a(x)"))+ list_plot(newinit2, plotjoined=True, color="green", legend_label=("b(x)")) # plots both of them together overlaying

#4b
e(x)=3.8*x*(1-x) # discrete logistic function w/ r=3.8 and ic=0.8
count=srange(0,20)
einitial=[0.8] 
for i in count: 
    newe=e(einitial[i]) 
    einitial.append(newe)
print(einitial)
k(x)=3.8*x*(1-x) # discrete logistic funciton w/ r=3.8 and ic=0.801
count=srange(0,20)
kinitial=[0.801] 
for i in count: 
    newk=k(kinitial[i]) 
    kinitial.append(newk)
print(kinitial)
list_plot(einitial, plotjoined=True, axes_labels=("time", "e(x)+k(x)"), legend_label=("e(x)"))+list_plot(kinitial, plotjoined=True, color="green", legend_label=("k(x)"))

#4c
c(x)=3.9*x*(1-x) # discrete logistic function w/ r=3.9 and ic=0.9
count=srange(0,20)
c_initial=[0.9] 
for i in count: 
    n=c(c_initial[i]) 
    c_initial.append(n)
print(c_initial)
d(x)=3.9*x*(1-x) # discrete logistic function w/ r=3.9 and ic=0.901
count=srange(0,20)
d_initial=[0.901] 
for i in count: 
    n=d(d_initial[i]) 
    d_initial.append(n)
print(d_initial)
list_plot(c_initial, plotjoined=True, axes_labels=("time", "c(x)+d(x)"), legend_label=("c(x)"))+list_plot(d_initial, plotjoined=True, color="green", legend_label=("d(x)")) # plots both of them together overlaying

# Overall, the relationship between the pairs of time series shows that a slight change in the initial condition causes bifurcations in the long term states despite similar initial transient states.

#5
count=srange(0,8) # establishes count for first 8 pairs
rvalue=[] # empty list for the values of the difference in r-values from iteration to iteration
for i in count:
    rvalue.append(abs((newinitial[i])-(newinitial1[i]))) # calculates the difference between each functions from each other (which is g(x)-h(x)                                                                from Exercise 3), gives difference in absolute value for each iteration and appends that                                                            value to the 'rvalue' list
print(rvalue)
list_plot(rvalue, plotjoined=True, axes_labels=("time", "difference in r values"))

#6
p=list_plot(list(zip(newinitial[:-1], newinitial[1:])), axes_labels=("h(x_N)", "h(x_N+1)")) # we are using data from Exercise 3, where the ic=0.4,                                                                                               where the values were placed in a list called                                                                                                       'newinitial'.
# 'newinitial[:-1]' means all elements up until the last one, which is h(x_N), and 'newinitial[1:]' means the first element up until the last one, which is h(x_N+1). We plot h(x_N) and h(x_N+1) against each other in order to make the Poincare plot, and label the axes accordingly. We will name all of this p.

show(p) # plots 'p'

#7
func=3.9*x*(1-x) # the function we were told to plot
q=plot(func, (x, 0,1)) # make sure plot it for 0<=x<=1 and name it 'q'
show(p+q) # plots p and q together on the same space

#The function lines up exactly with the points in the previous exercise. This makes sense because both exercises use the same r value, but #7 is a continuous function (line) and #6 is a discrete time function (points).

#8a
r11=3.5   # our r-value will be 3.5
r(x)=r11*x*(1-x) # logistic equation with r-value
xvals1=[0.9] # our chosen initial condition
for i in srange(0,250,1): # we want to iterate this 250 times by a timestep of 1. for every iteration...
    xnew=r(xvals1[-1]) #... we will plug in the current 'xvals1' into r(x) and call it 'xnew'
    xvals1.append(xnew) # we will then append whatever 'xnew' is to the 'xvals' list and have it be the new current value
show(xvals1) # shows the list

#8b
xvaluesnew=xvals1[100:] #takes the 'xvalues' list from the previous part of the question and indexes them to only count for iterations 100-250. we                          rename this new list 'xvaluesnew'
xvaluesnew # shows the list (we don't have to put show() since it's just numbers either way)

#8c
formlistrx=150*[3.5] # makes a list of just 3.5 iterated 150x (because out of 200, we threw away the first 100)
coord=list(zip(formlistrx, xvaluesnew)) # this zips together the values of r and our function X (in this case, r(x)) as points. we name it 'coord'
coord # shows us the points we've zipped together as a list

#There are 150 different values of X.

#9 We threw away the first 100 values in the previous exercise because we are interested in observing the long-term behavior. This way, we are able to observe the behavior past the transient states.

#10
masterlist=[]
r10=srange(2.0, 4.0, 0.01) # our r-values will range from 2 to 4 by time steps of 0.01 (aka, 200 iterations)
for r in r10: # for every iteration of r...
    func1(x)=r*x*(1-x) # logistic function
    xval=[0.9] # initial condition in a list called 'xval'
    for o in srange(0,250,1): # for every iteration (in this case, 250x)...
        new1=func1(xval[-1]) #... whatever the current 'xval' is will be plugged into func1(x) and be called 'new1'
        xval.append(new1) # watever 'new1' is will be plugged into the 'xval' list and become the new current value
    coord1=list(zip(150*[r], xval[100:])) # put them together, and our x coordinates will be 150 iterations of r-values ranging from 2.0-4.0.
    masterlist=masterlist+coord1
list_plot(masterlist, axes_labels=["r", "Xn"], size=1)

#11
# The approximate value that produces chaotic behavior are about r=3.6

#12
# As r increases, the behavior of the discrete logistic model goes through a Hopf bifurcation at about r=3.0, then a period doubling bifurcation at about r=3.45, then chaotic behavior at about r=3.6.
# This plot aligns with the results from Exercise 1 because it reflects the changing nature of the oscillations. At r=3, there is 2-point oscillatory behavior, at r=3.5, there's 4-point oscillatory behavior, and as said before, around 3.6, the behavior becomes chaotic.

diff(0.15*y - 0.6*z*x - 0.01*z*x, x)