# parametri
alpha = 4.15
mu = 0.001
sigma = -1
|
|
# početni uslovi
x_0 = 0
y_0 = -2.8
# broj iteracija (+1, zbog početnog uslova)
N = 4999
|
|
# Rulkovljeva mapa
def x_new(x_old, y_old, alpha):
return alpha/(1+x_old**2) + y_old
def y_new(y_old, x_old, mu, sigma):
return y_old - mu*(x_old - sigma)
|
|
# kako naći fiksne tačke? how to find all the roots? (http://ask.sagemath.org/question/263/numerically-find-all-roots-in-an-interval)
def find_all_roots(f, a, b, eps=0.0000000001):
roots = []
intervals_to_check = [(a,b)]
while intervals_to_check:
start, end = intervals_to_check.pop()
try:
root = find_root(f, start, end)
except RuntimeError:
continue
if root in roots:
continue
if abs(f(root)) < 1:
roots.append(root)
intervals_to_check.extend([(start, root-eps), (root+eps, end)])
roots.sort()
return roots
|
|
# iteracije + naći fiksne tačke
var('x_solution')
x = [x_0]
y = [y_0]
fixed_points = [] # (x_solution, y)
max_iterates = []
min_iterates = []
for i in range(N):
x_temp = x_new(x[i], y[i], alpha)
y_temp = y_new(y[i], x[i], mu, sigma)
x_solutions = find_all_roots(x_new(x_solution, y[i], alpha) - x_solution, -4, 4)
fixed_points += zip([y[i]]*len(x_solutions), x_solutions)
max_iterates.append(alpha + y[i])
min_iterates.append(alpha/(1+(alpha+y[i])**2)+y[i])
x.append(x_temp)
y.append(y_temp)
|
|
# plot
plotovi = list_plot(fixed_points, color='black') + list_plot(zip(y, max_iterates), color='red') + list_plot(zip(y, min_iterates), color='blue')
plotovi += list_plot(zip(y,x), plotjoined=True, color='green')
plotovi += list_plot(zip(y,[sigma]*N), color='pink')
show(plotovi)
|
|
# first return map
first_map = []
for i in range(N):
first_map.append((x[i],x[i+1]))
#plot
list_plot(first_map, plotjoined=False)
|
|
# a derivative? šta je to nullcline?
prvi_izvod = []
for i in range(N):
prvi_izvod.append((x[i], x[i+1]-x[i]))
# plot
list_plot(prvi_izvod, plotjoined=False)
|
|
|
|