Sprott 式ストレンジアトラクタ発見機
米国の物理学者である J. C. Sprott は、今から 30 年ほど前の 1993 年に "Automatic Generation of Strange Attractors" なる題名の論文を発表した。
https://sprott.physics.wisc.edu/pubs/paper203.pdf
題名どおり、この論文ではストレンジアトラクタを自動的に生成する具体的な方法が提示されている。ストレンジアトラクタというのは数学や物理における用語であり詳しい説明は専門書に譲るが、簡単な仕組みであるにも関わらず非常に美しい図形を描く（場合もある）ことで有名だ。
具体的な例としては本稿冒頭に表示されている画像がそれであり、また以下のような画像も Sprott の方法により発見したストレンジアトラクタの図である：
これらの図形は以下の式で表される、x と y の 2 次式により生成されたものである。
$$
\left\{ \begin{array}{ccl} x'&=&a_1+a_2x+a_3x^2+a_4xy+a_5y+a_6y^2 \\\ y'&=&a_7+a_8x+a_9x^2+a_{10}xy+a_{11}y+a_{12}y^2 \end{array} \right.
$$
ちょっと項数が多いので、これが単純な式？と思われるかもしれないが、ここには三角関数も平方根も出てきてはいない。中学生でも計算できるし、変数の意味さえ分かれば小学生でも計算は可能である。
Sprott の論文では BASIC というプログラミング言語で発見用のプログラムが記述されていた。30 年前は BASIC もまだまだ現役だったのだ。もっとも、今でも BASIC そのものは現役ではあるが、行番号が不要になっていたり、構造化言語になっていたりと、論文に記載されている BASIC とは別の言語になっている感がある。
BASIC は手軽な言語ではあるが、秒進分歩と言われる計算機の世界においては、やはり古いという印象は否めない。すべて大域変数となっている Sprott のプログラムもアルゴリズムの構造が分かりづらい。
そのため、私自身がプログラムの構造を理解するため、また、より多くの人々に気軽にストレンジアトラクタの発見過程を楽しんでもらうため、今回、私の方でこの Sprott のプログラムを Processing (Python mode) に移植した。
このプログラムは全部で 110 行弱であるが、移植のついでに、少しばかり魅力的な図になるように描画方法を加色混合としたり、高解像度画像を画像ファイルをして保存できる機能を追加している。
まずは、そのプログラムを以下に示す：
# strange attractor generater with Sprott's algorithm
# program by KojiSaito ( Twitter @KojiSaito )
InitX=InitY=.05
def init():
global XE,YE,Lsum,Xmin,Xmax,Ymin,Ymax
YE=.05;XE=InitX+.000001
Lsum=0
Xmin=Ymin=1000000;Xmax=Ymax=-Xmin
id=[int(random(25)) for i in range(12)]
coef=[(t-12)*.1 for t in id]
idStr=""
for c in id:idStr+=chr(ord('A')+c)
return coef,idStr
def next(coef,x,y):
xNew=coef[0]+x*(coef[1]+coef[2]*x+coef[3]*y)+y*(coef[ 4]+coef[ 5]*y)
yNew=coef[6]+x*(coef[7]+coef[8]*x+coef[9]*y)+y*(coef[10]+coef[11]*y)
return xNew,yNew
def updateBound(x,y):
global Xmin,Xmax,Ymin,Ymax
Xmin=min(x,Xmin);Xmax=max(x,Xmax)
Ymin=min(y,Ymin);Ymax=max(y,Ymax)
def LyapunovExponent(coef,x,y,n):
global XE,YE,Lsum
tx,ty=next(coef,XE,YE)
dLx=tx-x;dLy=ty-y
dL2=dLx*dLx+dLy*dLy
df=1e12*dL2;rs=1/sqrt(df)
XE=x+rs*(tx-x);YE=y+rs*(ty-y)
Lsum+=log(df)
return .721347*Lsum/n # L=Lyapunov exponent
def search():
while True:
coef,idStr=init()
x,y=InitX,InitY
found=True
for i in range(1,11000):
x,y=next(coef,x,y);updateBound(x,y)
l=LyapunovExponent(coef,x,y,i)
if abs(x)+abs(y)>1e6 or (i>1000 and l<.002): #Sprott 's original is .005
found=False; break
if found:break
print('ID="'+idStr+'"'," L=",l)
return coef,idStr
def calcScale(coef,w,h):
init()
x,y=InitX,InitY
for i in range(1000):x,y=next(coef,x,y);updateBound(x,y)
dx=.1*(Xmax-Xmin);dy=.1*(Ymax-Ymin)
xl=Xmin-dx;xh=Xmax+dx
yl=Ymin-dy;yh=Ymax+dy
sx=w/(xh-xl);sy=h/(yh-yl)
return sx,sy,xl,yl
def setColor(idStr):
colorH=538
for c in idStr:colorH=ord(c)+colorH*33
colorH%=251
s=[0]*3;t=colorH%3
s[t]=3;s[(t+1)%3]=1.5;s[(t+2)%3]=.8
colorMode(HSB,250);c=color(colorH,190,200) # make color here
colorMode(RGB,255);blendMode(ADD);
plotColor=color(s[0]*(red(c)+40),s[1]*(green(c)+40),s[2]*(blue(c)+40))
stroke(plotColor)
return plotColor
IdStt=""
Coef=[]
PlotColor=0
N=90000
def drawSA(w,h):
global Coef,IdStr,PlotColor
clear()
Coef,IdStr=search()
PlotColor=setColor(IdStr)
sx,sy,tx,ty=calcScale(Coef,w,h)
x,y=InitX,InitY
for i in range(N):x,y=next(Coef,x,y);point((x-tx)*sx,(y-ty)*sy)
def render(fileName):
w=h=2000
hiResImage=createImage(w,h,RGB)
sx,sy,tx,ty=calcScale(Coef,w,h)
x,y=InitX,InitY
for i in range(N*w/width*h/height):
x,y=next(Coef,x,y)
u,v=int(sx*(x-tx)),int(sy*(y-ty))
if u<0 or w<=u or v<0 or h<=v:continue
c=hiResImage.pixels[w*v+u]
hiResImage.pixels[w*v+u]=color(red(c)+red(PlotColor),
green(c)+green(PlotColor),
blue(c)+blue(PlotColor))
hiResImage.save(fileName)
def setup():size(500,500);drawSA(width,height)
def keyPressed():
if key==' ':drawSA(width,height)
if key=='s':
print("Generating Hi-Res Image: "+IdStr+".png ...")
render(IdStr+'.png')
# save(IdStr+'.jpg') # if you need...
print("Done.")
def draw():random(1)使い方は至って単純である。まず、Processing を起動し、「ファイル」から「新規」を選び、Python モードになっていることを確認して、上記のプログラムをコピー＆ペーストする。あとは適当な名前をつけて、このプログラムを保存し、その後、実行するだけである（後述する高解像度画像の保存機能を使わなければ、プログラムの保存は不要かもしれない）。
実行すると、Processing のウィンドウが現れ、運が良ければそこに魅力的なストレンジアトラクタが表示される。
もしあまり魅力的でない図形が表示されても慌てることはない。スペースバーを押下すれば、次の「魅力的であろう」ストレンジアトラクタを検索して表示する。
この「魅力的であろう」という部分であるが、現状は Sprott の論文に記載されていたプログラムと同様にリアプノフ指数による選別しか行っていない。Sprott 自信は論文の中でフラクタル次元も援用する方法を記述しているが、その具体的な計算プログラムは掲載していない。
フラクタル次元の具体的な計算方法については、Sprott の "Strange Attractors: Creating Patterns in Chaos" という本の中で言及されているが、その処理は上に掲載したプログラムには入っていない。
ちなみに Sprott の本は PDF 版が以下の URL で公開されている。大変ありがたいことである。
https://sprott.physics.wisc.edu/fractals/booktext/SABOOK.PDF
この本を頑張って読解し、フラクタル次元についても、できればこの発見プログラムにその機能を導入したいと思っている。しかし現時点では「今後の課題」という状況にある。
Sprott のプログラムではリアプノフ指数の閾値を 0.005 としていたが、本稿のプログラムでは 0.002 としている。これは、こちらの方が - つまり、リアプノフ指数が 0.004 や 0.003 のアトラクタでも十分魅力的なものも存在するため、このような変更を行っている。
もしかすると、厳密には本稿で得られるアトラクタは「ストレンジ」なアトラクタではないかもしれない。が、より面白い図形が表示されるのを期待して、あえてこのような処理にしている。気になる人は元の値に戻して実行してみるのも面白いと思う。
最後に、魅力的な図形が得られた場合の保存方法を紹介する。s キーを押下すれば、Processing のコンソールエリアに Generating Hi-Res Image: KBSGOUWREMES.png ... などと生成される PNG 画像のファイル名を含んだ文字列が表示される。.png の前の文字列は Sprott のストレンジアトラクタ ID である。
高解像度画像の生成にはしばらくの時間がかかる。M1 の mac book air で 10 秒弱程度かかるようだ。処理が終わるまで少し待って欲しい。無事、画像ファイルが生成・保存できると、コンソールエリアに Done と表示される。
高解像度画像生成時に、まれにエラーとなることがある。これは計算により、ストレンジアトラクタの軌道が非常に大きな値となってしまうことが原因のようだ。この場合、本ソフトウェアでは高解像度画像画像が生成できない。諦めて再度ストレンジアトラクタ発見の旅に出かけてもらいたい。
高解像度画像ファイルのファイル名に使用される暗号のような記号列を私は Sprott のストレンジアトラクタ ID と呼んでいる。というのも、この ID さえあれば、Sprott 式の 2 次式からなるストレンジアトラクタが再現できるからである。この ID は、発見したストレンジアトラクタの情報共有において極めて重要なものなのである。
この ID については、国外でも利用している事例が散見される。
もともと英語圏で発表された識別の方法なので、海外で使用されているのは当たり前といえばそのとおりである。とにかく、国外でも通用する ID であることをここでは強調しておきたい。私のプログラムで適当に作っている ID では決してないので、もし面白いストレンジアトラクタを見つけたら積極的に情報発信していただければと思う。
なお、この ID がどのように決められているのかについては、以下のページに簡単な解説を書いておいた。もし興味を持たれたかたは読んでいただきたい：
話を生成される高解像度画像ファイルに戻そう。画像ファイルの保存場所は、本稿で紹介した Processing のスケッチをあなたが保存したディレクトリとなる。s キーを押下した後、当該ディレクトリを見てもらえれば、そこに画像ファイルが保存されているはずである。
もし本プログラムによるストレンジアトラクタ探索が面白いと思ったら、見つけた面白いストレンジアトラクタの画像と ID をハッシュタグ #StrAttrCat と共に tweet して頂けたら幸いである（ちなみに #StrAttrCat は Strange Attractor Catalog の略である）。
また、今回のプログラムに関して、質問や機能追加の要望、などなど、何か聞きたいことがあれば Twitter の @KojiSaito まで連絡していただければと思う。それでは楽しいストレンジアトラクタライフを。
この記事が気に入ったらサポートをしてみませんか？