# -*- coding: utf8 -*- from numpy import * print(" nr=1,5 r=2,0 r=2,5 ") M=zeros([201,3]) a=[1.5,2.0,2.5] for j in arange(0,3,1): M[0,j]=0.5 for j in arange(0,3,1): for i in arange(1,201,1): M[i,j]=a[j]*M[i-1,j]*(1-M[i-1,j]) for i in arange(0,201,1): if 0<=i<=2: print(" {0: 7.0f} {1: 7.4f} {2: 7.4f} {3: 7.4f} " . format(i,M[i,0],M[i,1],M[i,2])) elif 2<i<=5: print(".") elif 197<=i<=200: print(" {0: 7.0f} {1: 7.4f} {2: 7.4f} {3: 7.4f} " . format(i,M[i,0],M[i,1],M[i,2]))
nr=1,5 r=2,0 r=2,5 0 0.5000 0.5000 0.5000 1 0.3750 0.5000 0.6250 2 0.3516 0.5000 0.5859 . . . 197 0.3333 0.5000 0.6000 198 0.3333 0.5000 0.6000 199 0.3333 0.5000 0.6000 200 0.3333 0.5000 0.6000
# -*- coding: utf8 -*- from numpy import * print(" nr=3,1 r=3,25 r=3,5 ") M=zeros([1008,3]) a= [3.1,3.25,3.5] for j in arange(0,3,1): M[0,j]=0.5 for j in arange(0,3,1): for i in arange(1,1008,1): M[i,j]=a[j]*M[i-1,j]*(1-M[i-1,j]) for i in arange(0,1008,1): if 0<=i<=3: print(" {0: 7.0f} {1: 7.4f} {2: 7.4f} {3: 7.4f} " . format(i,M[i,0],M[i,1],M[i,2])) elif 4<i<=7: print(".") elif 1000<=i<=1007: print(" {0: 7.0f} {1: 7.4f} {2: 7.4f} {3: 7.4f} " . format(i,M[i,0],M[i,1],M[i,2]))
nr=3,1 r=3,25 r=3,5 0 0.5000 0.5000 0.5000 1 0.7750 0.8125 0.8750 2 0.5406 0.4951 0.3828 3 0.7699 0.8124 0.8269 . . . 1000 0.5580 0.4953 0.5009 1001 0.7646 0.8124 0.8750 1002 0.5580 0.4953 0.3828 1003 0.7646 0.8124 0.8269 1004 0.5580 0.4953 0.5009 1005 0.7646 0.8124 0.8750 1006 0.5580 0.4953 0.3828 1007 0.7646 0.8124 0.8269
# -*- coding: utf8 -*- import matplotlib.pyplot as plt from numpy import * M=zeros([1008,3]) a= [3.1,3.25,3.5] for j in arange(0,3,1): M[0,j]=0.5 for j in arange(0,3,1): for i in arange(1,1008,1): M[i,j]=a[j]*M[i-1,j]*(1-M[i-1,j]) x=arange(987,999,1) y=M[987:999,0] y1=M[987:999,1] y2=M[987:999,2] plt.title('Удвоение цикла роста популяции для r=3,1;3,25;3,5') plt.plot(x,y, label="T=1,ymax=%s,ymin=%s"%(round(max(y),3),round(min(y),3))) plt.plot(x,y1, label="T=2,ymax=%s,ymin=%s"%(round(max(y1),3),round(min(y1),3))) plt.plot(x,y2, label="T=4,ymax=%s,ymin=%s"%(round(max(y2),3),round(min(y2),3))) plt.grid() plt.legend(loc="best") plt.ylabel("x(n)") plt.xlabel("n") plt.show()
# -*- coding: utf8 -*- import matplotlib.pyplot as plt from numpy import * print(" nr=3,57 ") M=zeros([1041,1]) a= [3.57] for j in arange(0,1,1): M[0,j]=0.5 for j in arange(0,1,1): for i in arange(1,1041,1): M[i,j]=a[j]*M[i-1,j]*(1-M[i-1,j]) for i in arange(0,1041,1): if 1000<=i<=1015: print(" {0: 7.0f} {1: 7.4f}" . format(i,M[i,0])) x=arange(999,1040,1) y=M[999:1040,0] plt.title(' Не детерминированный хаос для r=3,57') plt.plot(x,y) plt.grid() plt.ylabel("x(n)") plt.xlabel("n") plt.show()
nr=3,57 1000 0.4751 1001 0.8903 1002 0.3487 1003 0.8108 1004 0.5477 1005 0.8844 1006 0.3650 1007 0.8275 1008 0.5096 1009 0.8922 1010 0.3434 1011 0.8050 1012 0.5604 1013 0.8795 1014 0.3784 1015 0.8397
# -*- coding: utf8 -*- import matplotlib.pyplot as plt from numpy import* N=1000 y=[] y.append(0.5) for r in arange(2.8,4.0,0.0001): for n in arange(1,N,1): y.append(round(r*y[n-1]*(1-y[n-1]),4)) y=y[N-250:N] x=[r ]*250 plt.plot( x,y, color='black', linestyle=' ', marker='.', markersize=1) plt.title("Диаграмма ветвления при 2,8<= r<=4,0,0<=x<=1") plt.xlabel("r") plt.ylabel("y") plt.axvline(x=3.01,color='black',linestyle='--') plt.axvline(x=3.45, color='black',linestyle='--') plt.axvline(x=3.6,color='black',linestyle='--') plt.axvline(x=3.7,color='black',linestyle='--') plt.axvline(x=3.8,color='black',linestyle='--') plt.axvline(x=3.9,color='black',linestyle='--') plt.show()
# -*- coding: utf8 -*- import matplotlib.pyplot as plt from numpy import* N=1000 y=[] y.append(0.5) for r in arange(3.8,3.9,0.0001): for n in arange(1,N,1): y.append(round(r*y[n-1]*(1-y[n-1]),4)) y=y[N-250:N] x=[r ]*250 plt.plot( x,y, color='black', linestyle=' ', marker='.', markersize=1) plt.title("Диаграмма ветвления при 3,8<= r<=3,9,0<=x<=1") plt.xlabel("r") plt.ylabel("y") plt.axvline(x=3.83,color='black',linestyle='--') plt.show()
# -*- coding: utf8 -*- import matplotlib.pyplot as plt from numpy import * a=2.7 x1=0.62 def ff(x): return a*x*(1-x) b=a*x1*(1-x1)/x1 def fl(x): return b*x x=0.1 y=0 Y=[] X=[] for i in arange(1,30,1): X.append(x) Y.append(y) y=ff(x) X.append(x) Y.append(y) x=y/b plt.title('Паутинная диаграмма логистической \n функции λx(1-x) при λ = 2.7') plt.plot(X,Y,'r') x1=arange(0,1,0.001) y1=[ff(x) for x in x1] y2=[fl(x) for x in x1] plt.plot(x1,y1,'b') plt.plot(x1,y2,'g') plt.grid(True) plt.show()
# -*- coding: utf8 -*- from numpy import * from scipy. integrate import odeint import matplotlib.pyplot as plt for F in [0.6,0.7,0.75,0.8]: def f(y,t): y1,y2=y return [y2,-y2-y1**3+y1+F*cos(t)] t=arange(100,200,0.001) y0=[1.0,0.0] [y1,y2]=odeint(f, y0, t, full_output=False,rtol=1e-12).T if F==0.6: plt.subplot(221) plt.title('Фазовая плоскость F=0.6,T=2'r'$\pi$') plt.plot(y1,y2, color='black', linestyle=' ', marker='.', markersize=0.1) plt.grid(True) plt.subplot(222) plt.title('Решение x(t): F=0.6,T=2'r'$\pi$') plt.plot(t,y1, color='black', linestyle=' ', marker='.', markersize=0.1) plt.grid(True) elif F==0.7: plt.subplot(223) plt.plot(y1,y2, color='black', linestyle=' ', marker='.', markersize=0.1, label='Фазовая плоскость \n F=0.7,T=4'r'$\pi$') plt.legend(loc='upper left') plt.grid(True) plt.subplot(224) plt.plot(t,y1, color='black', linestyle=' ', marker='.', markersize=0.1, label='Решение x(t): F=0.7,T=4'r'$\pi$') plt.legend(loc='upper left') plt.grid(True) plt.show() if F==0.75: plt.subplot(221) plt.title('Фазовая плоскость F=0.75,T=8'r'$\pi$') plt.plot(y1,y2, color='black', linestyle=' ', marker='.', markersize=0.1) plt.grid(True) plt.subplot(222) plt.title('Решение x(t): F=0.75,T=8'r'$\pi$') plt.plot(t,y1, color='black', linestyle=' ', marker='.', markersize=0.1) plt.grid(True) elif F==0.8: plt.subplot(223) plt.plot(y1,y2, color='black', linestyle=' ', marker='.', markersize=0.1, label='Фазовая плоскость\n F=0.8,Хаос') plt.legend(loc='upper left') plt.grid(True) plt.subplot(224) plt.plot(t,y1, color='black', linestyle=' ', marker='.', markersize=0.1, label='Решение x(t): F=0.8,Хаос') plt.legend(loc='upper left') plt.grid(True) plt.show()
# -*- coding: utf8 -*- import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt s,r,b=10,28,8/3 def f(y, t): y1, y2, y3 = y return [s*(y2-y1), -y2+(r-y3)*y1, -b*y3+y1*y2] t = np.linspace(0,20,2001) y0 = [-8, 8, 27] [y1,y2,y3]=odeint(f, y0, t, full_output=False).T plt.plot(y1,y3, color='black', linestyle=' ', marker='.', markersize=2) plt.xlabel('x') plt.ylabel('z') plt.grid(True) plt.title("Проекция траектории Лоренца на плоскость xz") plt.show()
# -*- coding: utf8 -*- import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D #Создаем функцию правой части системы уравнений. s,r,b=10,25,3 def f(y, t): y1, y2, y3 = y return [s*(y2-y1), -y2+(r-y3)*y1, -b*y3+y1*y2] #Решаем систему ОДУ и строим ее фазовую траекторию t = np.linspace(0,20,2001) y0 = [1, -1, 10] [y1,y2,y3]=odeint(f, y0, t, full_output=False).T fig = plt.figure(facecolor='white') ax=Axes3D(fig) ax.plot(y1,y2,y3,linewidth=2) plt.xlabel('y1') plt.ylabel('y2') plt.title("Начальные условия: y0 = [1, -1, 10]") y0 = [1.0001, -1, 10] [y1,y2,y3]=odeint(f, y0, t, full_output=False).T fig = plt.figure(facecolor='white') ax=Axes3D(fig) ax.plot(y1,y2,y3,linewidth=2) plt.xlabel('y1') plt.ylabel('y2') plt.title("Начальные условия: y0 = [1.0001, -1, 10]") plt.show()
# -*- coding: utf8 -*- import numpy as np from scipy.integrate import odeint import matplotlib.pyplot as plt a,b,c=0.398,2.0,4.0 def f(y, t): y1, y2, y3 = y return [(-y2-y3), y1+a*y2, b+y3*(y1-c)] t = np.linspace(0,50,5001) y0 = [0,0, 0] [y1,y2,y3]=odeint(f, y0, t, full_output=False).T plt.plot(y1,y2, color='black', linestyle=' ', marker='.', markersize=2) plt.xlabel('x') plt.ylabel('y') plt.grid(True) plt.title("Проекция ленты Ройсслера на плоскость xy") plt.show()
Source: https://habr.com/ru/post/436014/