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

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

Step0. 使用するパッケージ

import numpy as np
import matplotlib.pyplot as plt
import random

Step1. 写像を定義

写像を定義

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. 時系列の生成

時系列の生成

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. クラスター数の判定

クラスター数の判定

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. プロットのために整理する

プロット用に整列させる

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. すべてのプログラムを実行しプロット

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)

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()

実行結果の例