PythonによるGlobally Coupled Mapの相図を描画する

Globally coupled map に関する特徴や数学的背景に関しては以下を参照

https://qiita.com/kojiwada/items/c3f423ca9c18448400c6

### Step0. 使用するパッケージ
~~~python:
import numpy as np
import matplotlib.pyplot as plt
import random
~~~

### Step1. 写像を定義
~~~python:写像を定義
def f(x, alpha):
  y = 1 - alpha*x*x
  return y

def coupled_f(x_list, alpha, e):
  total = 0
  y_list = []
  num = len(x_list)

  for n in range(num):
    total += f(x_list[n], alpha)
  for n in range(num):
    y_list.append((1 - e) * f(x_list[n], alpha) + (e / num) * total)
  return y_list
~~~

### Step2. 時系列の生成
~~~python:時系列の生成
def series(alpha, e, initial, length, trans_period):
  num_coupled = len(initial)  # 結合数
  x_list = [n + trans_period for n in range(length)]
  y_list = []

  for n in range(num_coupled):
    y_list.append([initial[n]])
  length -= 1
  for n1 in range(length + trans_period):
    pre_y_list = []
    for n2 in range(len(initial)):
      pre_y_list.append(y_list[n2][n1])

    after_y_list = coupled_f(pre_y_list, alpha, e)
    for n2 in range(len(initial)):
      y_list[n2].append(after_y_list[n2])

  for n1 in range(num_coupled):
    y_list[n1] = y_list[n1][trans_period:]
  return x_list, y_list
~~~

### Step3. クラスター数の判定
~~~python:クラスター数の判定
def count_num_cluster(y_list):
  couple = np.size(np.array(y_list).T[-1,:])
  num_cluster = 0
  cluster_list = []
  for n in range(couple):
    check = 0
    for index, cluster in enumerate(cluster_list):
      if abs(cluster - np.array(y_list).T[-1,n]) < 10**(-5):
        check += 1
    if check == 0:
      cluster_list.append(np.array(y_list).T[-1,n])
      num_cluster += 1

  num_element = np.zeros(num_cluster, dtype=int)
  for n in range(couple):
    for index, cluster in enumerate(cluster_list):
      if abs(cluster - np.array(y_list).T[-1,n]) < 10**(-5):
        num_element[index] += 1

  return num_cluster
~~~

### Step4. プロットのために整理する
~~~python:プロット用に整列させる
def arrange_list(alpha_range, e_range):
  alpha_range_ar = []
  e_range_ar     = []

  for e in e_range:
    for alpha in alpha_range:
      alpha_range_ar.append(alpha)
      e_range_ar.append(e)
  return alpha_range_ar, e_range_ar

# プロットのためのデータを取得
def cal_plot(e_range, alpha_range, num_mean, couple, length, trans_period, sample_range):
  plot_list = []
  for n in range(num_mean):
    print("n :", n)
    plot_list.append([])
    for e in e_range:
      for alpha in alpha_range:
        initial = random.choices([random.random() for _ in range(sample_range)], k = couple)
        x_list, y_list = series(alpha, e, initial, length, trans_period)
        num_cluster = count_num_cluster(y_list)
        if num_cluster == 1:
          plot_list[n].append(0)
        elif num_cluster < 10:
          plot_list[n].append(10)
        elif num_cluster < couple:
          plot_list[n].append(20)
        elif num_cluster == couple:
          plot_list[n].append(30)
  plot_list = np.sum(np.array(plot_list), axis=0)
  plot_list = [n / num_mean for n in plot_list]
  return plot_list
~~~

### Step5. すべてのプログラムを実行しプロット
~~~python:
e_range = np.linspace(0.0, 0.6, 61)
alpha_range = np.linspace(1.4, 2.0, 61)

couple = 50
length = 10
trans_period = 500
sample_range = 10**4

num_mean = 10
plot_list = cal_plot(e_range, alpha_range, num_mean, couple, length, trans_period, sample_range)
alpha_range, e_range = arrange_list(alpha_range, e_range)
~~~

~~~python:
marker_size = 13
plt.figure(figsize=(5, 4))
plt.scatter(alpha_range, e_range, c=plot_list, cmap='viridis', s=marker_size, marker="s")
plt.xlabel(r'$\alpha$')
plt.ylabel(r'$\varepsilon$')
plt.show()
~~~

#### 実行結果の例
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3127285/70b8f797-81ed-db89-2210-380cf7b96545.png)