コラッツ予想：統計的振る舞いを確かめるためのプログラム（pythonコード）

コラッツの操作の統計的振る舞いの数値シミュレーションと、それを予測する算出式をpythonでプログラムしたものです。

グラフのプロット部分がまだ未完成なので、後日修正します。

# -*- coding: utf-8 -*-
"""
Created on Thu Nov 9 21:29:36 2017

@author: HobbyMath
"""
import numpy as np
import matplotlib.pylab as plt
from mpl_toolkits.mplot3d import Axes3D

count = 0
all_num = 10**5
maxval = 1
stopping_time = np.ones(all_num , dtype=int)
collatz_max = np.ones(all_num , dtype=int)
for num in range(1,all_num + 1):
tmp = num
collatz_max[num]=maxval
print(num)
step = 0
for time in range(10**3):

if np.mod(tmp,2) == 0:
tmp = int(tmp / 2)
step += 1
if tmp == 1:
stopping_time[num] = step
count+=1
break
if np.mod(tmp,2) == 1:
tmp = 3 * tmp + 1
step += 1
if(tmp > maxval):
maxval = tmp
collatz_max[num] = tmp

size = 500
ave_stopping_time = np.zeros(all_num, dtype=float)
var_stopping_time = np.zeros(all_num, dtype=float)
for num in range(10**3, all_num - size):
tmp_stopping_time = stopping_time[num - size : num + size]
ave_stopping_time[num] = np.average(tmp_stopping_time)
var_stopping_time[num] = np.var(tmp_stopping_time)

x = np.arange(1, all_num + 1)
y = np.log(x.astype(np.float64))

D = 2.95/3
#D = 1
v = np.log(4 / 3.0) / (np.log(3) * D)
scale = 2.0 / np.log(3.0)
alpha=1.5 * scale / v
beta = (1.5 + 0.5 * v)**2 * scale / v**3
predict_ave = alpha * y
predict_var = beta * y

#Plot average of stopping times
plt.hold(True)
plt_range = np.arange(10**3, all_num - size)
plt.plot(y[plt_range], ave_stopping_time[plt_range])
plt.plot(y[plt_range], predict_ave[plt_range])
plt.ylabel("Average of stopping times")
plt.xlabel("log(n)")
plt.hold(True)

#Plot variance of stopping times
plt.plot(y[plt_range], var_stopping_time[plt_range])
plt.plot(y[plt_range], predict_var[plt_range])
plt.ylabel("Variance of stopping times")
plt.xlabel("log(n)")

#Plot maximum value
plt.plot(y, np.log(collatz_max))
plt.ylabel("Maximum Value")
plt.xlabel("log(n)")

fig = plt.figure()
ax = fig.add_subplot(1,1,1)

ax.hist(stopping_time, bins=60)
plt.ylabel("Variance of stopping times")
plt.xlabel("Stopping times")

predict_hist = np.zeros(400, dtype=float)
count = 0
delta = 10**(-3)
maxnum = int(scale * np.log(all_num) / delta)
for Ts in range(1, 400):
print(Ts)
predict_hist[count] = 0.0

for Xind in range(1,maxnum):
# X=scale * np.log(num)
X = Xind * delta
predict_hist[count] += 2 * delta * X * (3 + v) **(0.5) / (scale * (2 * np.pi * (2 * Ts + X)**3) ** (0.5)) * np.exp(X / scale) * np.exp(-(3 * X - 2 * v* Ts)**2 / (2 *(3 + v) * (2 * Ts + X)))
count+=1

plt.plot(predict_hist)）