2016-12-02 3 views
-1

В моем задании домашней работы я написал код, который моделирует фигуры в классической механике по учебнику Тейлора. Рисунки 12.40 и 12.41, если кому-то интересно узнать.Bifurcation Diagram not plotting/Plot просто не отображается

Я был в состоянии воспроизвести одну из них, которая является ниже код (который может быть хорошим справочником для кода я ФАКТИЧЕСКИ возникли проблемы):

import nympy as np 
import matplotlib.pyplot as plt 

# We need to calculate the first fixed point 
r1=np.array(np.arange(0,4,0.09)) 
x1 = np.zeros((len(r1),1)) 

# Now calculating the second fixed point 
r2=np.array(np.arange(1,4,0.1)) 
x2 = (r2 -1)/r2 

# Now finding when the fixed points split up again 
r3=np.array(np.arange(3,4,0.1)) 
y1 = (((r3**2 - 2*r3 - 3)**0.5) + 1 + r3)/(2*r3) 
y2 = ((-(r3**2 - 2*r3 - 3)**0.5) + 1 + r3)/(2*r3) 

# Now finding the experimental values for 1/2 of a split 
x3 = [] 
for r in np.arange(0,4,0.09): 
    x = 0.666 
    for i in range(100): 
     x = (r**2) * x * (1.0 -x) - (r**3) * (x**2)*((1-x)**2) 
    x3.append(x) 

# Doing the same as above second 1/2 
x4 = [] 
for r in np.arange(0,4,0.09): 
    x = 0.8 
    for i in range(100): 
     x = (r**2) * x * (1.0 -x) - (r**3) * (x**2)*((1-x)**2) 
    x4.append(x) 

plt.plot(r1,x3,'bo', label='Experimental') 
plt.plot(r1,x4,'bo') 
plt.plot(r3,y2,'k-') 
plt.plot(r3,y1,'k-') 
plt.plot(r1,x1,'k-', label='Theoretical') 
plt.plot(r2,x2,'k-') 
plt.legend(loc=2) 
plt.show() 

А вот код для второго изображения что, похоже, не работает. И я не знаю, почему. Любая помощь будет оценена по достоинству. Фигура просто не замышляет, и я не знаю, почему.

import numpy as np 
import matplotlib.pyplot as plt 

for r in n.arange(2.8,4,0.01): 
    x = 0.5 
    for i in range(150): 
     x = r*x*(1-x) 
     if i >= 125: 
      plt.plot(r,x,'k') 
plt.xlim (2.8,4) 
plt.show() 
+0

Добро пожаловать так стека переполнения. Пожалуйста, прочитайте [как задать хороший вопрос] (http://stackoverflow.com/help/how-to-ask). Кроме того, вставьте свой код вместо вставки изображений. В сообщении есть кнопка '{}', которая отбрасывает все на четыре пробела и будет отображаться как код –

+0

Используете ли вы IPython ноутбук? – Xevaquor

ответ

0

Существует хорошая причина, по которой ваш код не работает.

Вы звоните plt.plot неоднократно со скалярными значениями, когда массив, как ожидается, для й и у: это довольно трудно построить линию только один точки. Попробуйте это, вы увидите, вы получите пустой участок:

plt.plot(0,1) 
plt.show(); 

Я понятия не имею, если это то, что вы ищете в конечном счете, но вот рабочая версия, которую вы можете по крайней мере использовать для отправная точка. Обратите внимание, что я изменил цикл на I: вместо того, чтобы делать i in range(150) с последующим if i >=125, лучше просто i in range(125,150), вы будете избегать ненужных повторений:

fig, ax = plt.subplots(nrows=1, ncols=1) 
r1 = [] 
x1 = [] 
for r in np.arange(2.8,4,0.01): 
    x=0.5 
    # Instead of range(150) then checking if i >= 125... 
    for i in range(125,150): 
     x = r*x*(1-x) 
     r1.append(r) 
     x1.append(x) 

ax.plot(r1,x1,'k') 
plt.show() 

Result graph

0

Я наткнулся на этот один год старый пост и понял, что есть простой ответ на проблему, о которой сообщает ОП.

Недостающий бит был marker аргумент.

Замените plot(r,x,'k') на plot(r,x,'.',color='k') и неожиданно вы увидите участок бифуркации.

Bifurcation for logistic regression

import numpy as np 
import matplotlib.pyplot as plt 

plt.figure() 
for r in n.arange(2.8,4,0.01): 
    x = 0.5 
    for i in range(150): 
     x = r*x*(1-x) 
     if i >= 125: 
      plt.plot(r,x,'.',color='black',markersize=2) 
plt.xlim (2.8,4) 
plt.show() 

Этот код прекрасно, хотя и довольно неэффективно. Гораздо быстрее позвонить plot() только один раз.

plt.figure() 
def mark(r,x,diagram,v): 
    N,M = np.shape(diagram) 
    rmin,rmax = r[0],r[-1] 
    for (i,j) in zip(r,x): 
     nx = int(j*(M-1)) 
     nr = int((i-rmin)/(rmax-rmin)*(N-1)) 
     diagram[nx,nr] = v 

r = np.arange(2.8,4,0.01) 
diagram = np.zeros((200,200)) 
x0 = 0.5 
x = np.ones_like(r)*x0 
for i in range(150): 
    x = np.multiply(np.multiply(r,x),(1-x)) 
    if i >= 125: 
     mark(r,x,diagram,1) 

plt.imshow(np.flipud(diagram), extent=[r[0],r[-1],0,1], 
      aspect=(r[-1]-r[0]),cmap=plt.cm.Greys) 
plt.show() 

Another bifurcation diagram

Последний код использует

CPU times: user 224 ms, sys: 29 µs, total: 224 ms 
Wall time: 261 ms 

Хотя бывший менее эффективный код занимает более 30 раз больше

CPU times: user 8.59 s, sys: 91.7 ms, total: 8.68 s 
Wall time: 8.68 s