Sprott 式ストレンジアトラクタ発見機

2022年3月17日 14:02

　米国の物理学者である 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 まで連絡していただければと思う。それでは楽しいストレンジアトラクタライフを。

この記事が参加している募集

#私の作品紹介

107,196件

この記事が気に入ったらサポートをしてみませんか？