「自然界の黄金比」をプロットしてみる

話題になった黄金率の話にうすーく乗っかってみたくなりました。

モチベーション

www.watto.nagoya

正直、連分数に関してちゃんと理解できていないのですが、とにかく、

(略)これは、黄金比の有理数による近似は常に精度が最も悪くなることを意味するのである。 これまでと同様の手順で、100枚の葉をつけた状態が、下図である。有理数による近似の精度が悪い＝腕が見えづらい＝葉の重なりが少ないってことで、葉の重なりが最も少なくなる配置が得られた。

というところに素直に感動してしまいました。自然界、まじやべぇ。*1

というわけで、これは matplotlib を使って動画化してみてみるしかないと奮い立ちました。先に言ってしまうと、自分の理解が薄かったせいで最初に考えていたよりも作業は大変でした。

課題のおさらい

幹から1枚目の葉を出したら、1枚目の葉からθの角度で2枚目の葉を出す。3枚目の葉は2枚目の葉からやはりθの角度で出す。以下n枚目の葉まで、同様に繰り返す。 幹からn枚目の葉の中心までの長さは、√nとする。

こうして幹の周りから同じ角度でずらしながらで直径１の葉っぱを配置していきます。中心からのきょりも順番の数の平方根で遠ざかっていって、全体に効率よく太陽の光をあびることのできる配置を考えるんですね。そうすると黄金比を使ったやつが一番となるらしいです。今回は元エントリーの例に倣って 無理数の代表　θ = π×π も作って比較してみます。

プロット結果

とっとと実行結果からお見せします。葉っぱの枚数は300枚としました。 まずは　π×π

そして　黄金比を使った場合。

きちんと前の円の群れで欠けたるところにすっぽりはまっていくように見えますね。不思議。

コード

jupyterで書いていますが、plt.show()使えば、jupyter以外でもできると思います。未確認。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib import animation

from IPython.display import HTML

T = 360 * 1.0 / (1 + (1+math.sqrt(5))/2)
#T = 360 * 1.0 / (math.pi * math.pi) 

fig = plt.figure()
ax = plt.gca()
ax.set_xlim(-20, 20)
ax.set_ylim(-20, 20)

def ccl_plt(n, d):
    r = n * d * math.pi / 180 
    x = math.sqrt(n) * math.cos(r)
    y = math.sqrt(n) * math.sin(r)
    return x, y

def init():
    return []

def animate(j):
    plt.cla()
    patches = []
    for i in range(j):
        if i == 0:
            continue
        patches.append(ax.add_patch(plt.Circle(ccl_plt(i, T), 1.0, alpha=0.1, color='g')))

    return patches

    
anim = animation.FuncAnimation(fig, animate,
                               init_func=init,
                               frames=300,
                               interval=80,
                               blit=True)


#plt.show()
HTML(anim.to_html5_video())

もう一つ、簡単な静止画バージョン。

import math
import matplotlib.pyplot as plt
%matplotlib inline

def ccl_plt(n, d):
    r = n * d * math.pi / 180 
    x = math.sqrt(n) * math.cos(r)
    y = math.sqrt(n) * math.sin(r)
    return x, y

fig = plt.figure()
ax = fig.add_subplot(111)
plt.xlim(-20, 20)
plt.ylim(-20, 20)
for i in range(300):

    if i == 0:
        continue
    t = 360 * 1.0 / (1 + (1+math.sqrt(5))/2)
    #t = 360 * 1.0 / (math.pi * math.pi)
    x, y = ccl_plt(i, t)
    ax.add_patch(plt.Circle((x, y), 1.0, alpha=0.1, color='g'))
plt.show()

蛇足

今回の課題は簡単かなとたかをくくっていたのですが、いざやってみると黄金比の角度をどう設定すればいいのかわからず苦労しました。検索して見つけたのですが、黄金角っちゅうのがあるのですね。入力が度数である場合は以下の計算で黄金角を導き出されるようです。

t = 360 * 1.0 / (1 + (1+math.sqrt(5))/2)

参考：http://gakuen.gifu-net.ed.jp/~contents/museum/golden/page62.html

分母の部分が黄金比の式になってますね。

いやはや、時間がかかりすぎてしまった。眠い。

*1:というかこうしたネタを華麗に持ち出してくるwattoさん、しゅごい。