混沌理论chaos + chaos = order（Can two chaotic systems give rise to order论文实验复现）

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt 
import  random as r

p=0.495
p1=0.095
p2=0.745

a=3
b=2
s=0
a1=[0]*100
a2=[0]*100
a32=[0]*100
a22=[0]*100
a44=[0]*100
cob=[a32,a22,a44]
vv=[0]*100

ite=10000 #迭代次数


for time in range(1,ite+1):
    #A1
    s=0
    for i in range(100):
        a1[i]+=s
        x=r.random()
        if x<p:
            s+=1
        else:
            s-=1
        
    
    #A2 
    s=0
    for i in range(100):
        a2[i]+=s 
        x=r.random()
        if abs(s)%3==0:
            if x<p1:
                s+=1
            else:
                s-=1
        else:
            if x<p2:
                s+=1
            else:
                s-=1
        

    #(a,b)
    para=[[3,2],[2,2],[4,4]]
    for k in range(3):
        pa=para[k]
        s=0
        for zuhe in range(int(100*1.0/(pa[0]+pa[1])+0.5)): 
            
            for aa1 in range(pa[0]):
                if zuhe*(pa[0]+pa[1])+aa1>=100:
                    break                
                cob[k][zuhe*(pa[0]+pa[1])+aa1]+=s
                x=r.random()
                if x<p:
                    s+=1
                else:
                    s-=1 

                    
            for aa2 in range(pa[1]):
                if zuhe*(pa[0]+pa[1])+pa[0]+aa2>=100:
                    break
                cob[k][zuhe*(pa[0]+pa[1])+pa[0]+aa2]+=s
                x=r.random()
                if abs(s)%3==0:
                    if x<p1:
                        s+=1
                    else:
                        s-=1
                else:
                    if x<p2:
                        s+=1
                    else:
                        s-=1
    #v
    oo=0.4
    s=0
    for i in range(100):
        choose=r.random()
        if choose<0.4:
            vv[i]+=s
            x=r.random()
            if x<p:
                s+=1
            else:
                s-=1
        else:
            vv[i]+=s
            x=r.random()
            if abs(s)%3==0:
                if x<p1:
                    s+=1
                else:
                    s-=1
            else:
                if x<p2:
                    s+=1
                else:
                    s-=1
                 

for i in range(100):
    a1[i]=a1[i]*1.0/ite
    a2[i]=a2[i]*1.0/ite
    vv[i]=vv[i]*1.0/ite
    for j in range(3):
        cob[j][i]=cob[j][i]*1.0/ite
plt.plot(a1)
plt.plot(a2)
plt.plot(vv)
for j in range(3):
    plt.plot(cob[j])
plt.show()