daruma3940の日記

理解や文章に間違い等あればどんなことでもご指摘お願いします

f:id:daruma3940:20160520223745p:plain
SDT5に向けて将棋盤アプリ作ったのじぇ
play.google.com
f:id:daruma3940:20160521003616p:plain
SDT5関係ある???
f:id:daruma3940:20160520223745p:plain
二日目暇な開発者さんにダウンロードしてもらってお金を稼ぐ算段なのじぇ
f:id:daruma3940:20160521003616p:plain
あんたねぇ...
f:id:daruma3940:20160520223745p:plain
どうぶつしょうぎも作ったのじぇ(公開はしてない もうちょっと機能追加したい)
f:id:daruma3940:20171027232525p:plain
f:id:daruma3940:20160709192554j:plain
でメインである僕は強くなってるの????
f:id:daruma3940:20160520223745p:plain
もうダメダメなのじぇ.....


f:id:daruma3940:20160709192554j:plain
僕の開発もがんばってくれよな~~たのむよ~~

拡張ユークリッド互除法について

f:id:daruma3940:20160520223745p:plain
ユークリッド互除法について書くのじぇ
ユークリッドの互除法について書くといったけどユークリッドの互除法についてはさらっと流して
拡張ユークリッド互除法について書きたいのじぇ
蟻本の拡張ユークリッドについての記述がよくわからなかったので自分でわかるようにした備忘録なのじぇ
f:id:daruma3940:20160521003616p:plain
まりちゃがこういうこといいだすのは大体よくわかってないことについてなので鵜呑みにしないほうがいいわよ
f:id:daruma3940:20160520223745p:plain

まずはユークリッド互除法について
自然数a,bに対して最大公約数を求める方法をユークリッド互除法というのじぇ
自然数a,bの最大公約数と毎回書くのは面倒なのでgcd(a,b)と書くことにするのじぇ。
ここでaをbで割った商をp余りをqとすると
a=b*p+qと書けることは明らかなのじぇ
b,qの最大公約数をgcd(b,q)とすると
b=x*gcd(b,q),q=y*gcd(b,q)のようにかけるので
a=gcd(b,q)*(x*p+y)と変形できるののじぇ。
gcd(b,q)は明らかにbを割り切り、上の式からaも割り切るのでgcd(b,q)はgcd(a,b)も割り切ることができなければならないのじぇ

ここで逆に
q=a-b*qとすると先ほどと同様にgcd(a,b)はgcd(b,q)を割り切れなければならないのじぇ
つまり
gcd(a,b)=X*gcd(b,q)
gcd(b,q)=Y*gcd(a,b)
となり、最大公約数は少なくとも1(>0)より大きいはずなのでX=Y=1となり
gcd(a,b)=gcd(b,q)となるのじぇ

aをbで割った余りをqとしたのでC言語的にq=a%bと書くとするのじぇ.
上の処理を何回も繰り返すと第二引数はどんどん減っていき最終的にゼロとなるのじぇ
gcd(a,b)=gcd(b,a%b)=gcd(a%b,b%(a%b))=.......=gcd(x,0)
0とxの最小公倍数はxなのでこれでgcd(a,b)は求まるのじぇこれがユークリッドの互除法なのじぇ

コードにするとこんな感じ(from 蟻本)
aをbにbをa%bにどんどん変えていく感じ

int gcd(int a,int b){
    if(b==0){return a;}
    return gcd(b,a%b);
}

f:id:daruma3940:20160520223745p:plain
これを拡張した拡張ユークリッドの互除法を利用することで
ax+by=gcd(a,b)となる整数解x,yをを求めることができるのじぇ

ちなみに右辺がなんでgcd(a,b)じゃないといけないのかについてはここを見てほしいのじぇ
mathtrain.jp

f:id:daruma3940:20160520223745p:plain
ユークリッド互除法のaをbに、bをa%bにどんどん変えていってどんどん係数を簡単にしていってgcd(a,b)を求めるのに似た方法を用いて
ax+by=gcd(a,b)もどんどん簡単にしていくのじぇ
f:id:daruma3940:20160520223745p:plain
最大公約数を求める方法と似てるなーーと思いながら読み進めてほしいのじぇ
f:id:daruma3940:20160520223745p:plain
とりあえず一回簡単にしてみるのじぇ
bx'+(a%b)y'=gcd(a,b)

aがbになってbがa%bになってxがx'にyがy'になっているのじぇ
gcd(b,a%b)=gcd(a,b)なので右辺は変わらんのじぇ

さっきのと似てるなーーーー


そしてこれをどんどん繰り返していくのじぇ
そうするとyの係数bはどんどん減少していって最終的にゼロになるのじぇ

ax'+0y'=gcd(a,b)

gcd(a,b)を求めるときbがゼロになったらgcd(a,b)=aだったのを思い出してほしいのじぇ
ということはx'の値は1でないといけないのじぇ

a*1+0*y=gcd(a,b)

0がかかっているのでy'の値はなんでもいいのじぇ。
ここではy'=0ということにしておくのじぇ
f:id:daruma3940:20160520223745p:plain
さてここまでやってきたことはax+by=gcd(a,b)をどんどん簡単にしていくことだったのじぇ

(a,b,x,y)

(b,a%b,x',y')
.
.
.

(a,0,1,0)
f:id:daruma3940:20160520223745p:plain
ax+by=gcd(a,b)となる整数解x,yを求めるには今度はこの逆をたどって最初の式に戻らなければならないのじぇ

(a,b,x,y)
↓  
(b,a%b,x',y')
.
.
.
↓  
(a,0,1,0)

f:id:daruma3940:20160520223745p:plain
その方法をここで説明するのじぇ

ax+by=gcd(a,b)

bx'+(a%b)y'=gcd(a,b)
にした場合のもとに戻し方を考えるのじぇ。

a%bはプログラミングの整数の割り算の扱い的に
a%b=a-(a/b)*b
とかけるので
bx'+(a%b)y'=gcd(a,b)は
bx'+(a-(a/b)*b)y'=gcd(a,b)となって
aとbについてまとめると
ay'+b(x'-(a/b)*y')=gcd(a,b)
となるのじぇ。
これを
ax+by=gcd(a,b)
と比較すると
x=y'
y=x'-(a/b)*y'

となるのでこのようにしてx,yの値を順にもとに戻していくことで
ax+by=gcd(a,b)となるx,yを求めることができるのじぇ
f:id:daruma3940:20160520223745p:plain
コードにするとこんな感じ(from 蟻本)

int extgcd(int a,int b,int& x,int& y){
    int d=a;//b==0のときgcd=a
    if(b!=0){
        //ここでgcdが求まる(b=0になる)まで潜る
        //さかのぼりながらx yを更新してゆく
        d=extgcd(b,a%b,y,x);

        //y=x'-(a/b)y'
        //xは1つ前のy yは一つ前のx
        y-=(a/b)*x;
    }
    else{
        //gcdが求まるとき(b==0のとき) x=1,y=0にする
        x=1;y=0;
    }
    return d;//gcdを返す
}


f:id:daruma3940:20160520223745p:plain

yの値をなんで0にしたのかについてはこれはほんとに何でもよくて他の値にすれば
他のax+by=gcd(a,b)を満たすx,yが求まるだけなのじぇ。
そもそもax+by=gcd(a,b)の解は一意ではないので...

f:id:daruma3940:20160520223745p:plain

ここからは完全に余談なのじぇ

さっき紹介したサイトを見ると
ax+by=gcd(a,b)の解き方は

この方程式の解を一つ求めてそれを(x0,y0)として

a(x1)+b(x2)=0
(x1:=x-x0 , y1:=y-y0 とした)

の形にすると
これは線形方程式なのでこの一般解はすぐに求まる。
そこからx=x1+x0,y=y1+y2の関係を用いて
x,yの一般式を求めることができると書いてあるのじぇ。
ここで(x0,y0)はさっき述べた拡張ユークリッド互除法で求めることができるのでこれで一般解を求めることができるのじぇ


f:id:daruma3940:20160520223745p:plain
これは非同次微分方程式を解くときと同じ考え方なのじぇ。

y''+ay'+y=f(x)
のを満たす解を一つ見つけてy_1とすると
(y-y_1)''+a(y-y_1)'+(y-y_1)=0 という形にできて
これは線形微分方程式なので簡単にとけて、その解をy_2とすると
y=y1+y2の関係から非同次微分方程式の一般解も求まるのじぇ

f:id:daruma3940:20160520223745p:plain
なるほど~~~~って感じなのじぇ

SDT5個人的注目ソフト

  • きふわらべ 

今年はどのような伝説を残してくれるのか期待wお話をするのが楽しみ。

  • messiah

個人的にすごく応援してます。

  • CGP

WCSC27でSquirrelにかなりヘイトを溜めているはずなのでガチで殺りに来てそう。怖い。

  • 海底

私がかなり苦労したbitboardを軽く実装するつよい人学習部も軽く実装しててガクブル。競プロもつよい。2駒勢。怖い。

  • shogi686

コンピュータ将棋将棋についていろいろなことを教えてくださった師匠。競プロもめっちゃつよい。2駒勢。怖い。お話をするのが楽しみ。

Unityのml-agents使ってみた

github.com
f:id:daruma3940:20160520223745p:plain

Unityのml-agents使ってみたので忘れる前に記録しておこうじぇ??


.......


......


...
f:id:daruma3940:20160520223745p:plain

もう忘れたwww





とりあえず記憶を呼び起こしながらでも今書いておかないとあとからまた思い出すのに手間取ってしまうのでがんばれ自分!!!

とりあえず上のリンクにあるwikiを見ながらやっていけば基本いける

Unityで作った環境(実行ファイルは ml-agents\pythonに置くのが安全(exeファイルだけじゃなくてDataフォルダも一緒に置かないと動かないので注意な)
env_name
には
env_name = "TestBall2"
みたいな感じで実行ファイルから.exeを除いた名前を置いておく。ちょっとわかりにかった(ソースコード読んで理解した。)
この指定の仕方が嫌ならソースコードを改変すればよろし

学習するときはBrain(Script)のType of BrainをExternalにしておく。
Type of brainのtypeにはPlayer, Hulistic, External, Internalの4タイプがある。
PlayerはPlayer Actionを決めることで人間が動かすモード
Hulisticは手で書いたコードを動かすモード
Externalは学習するときのモード
Internalは学習させたモデルを使って動かすときのモード
InternalモードはTensorFlow# AssetをUnity Projectに入れて PlayerSettingのScripting Define SymbolにENABLE_TENSORFLOWを追加して Scripting Runtime VersionをExperimental(.Net4.6)にしないといけないはず


自分の作った環境で
学習させたりするのは大変そうに見えるかもしれないけどそんなに難しくないはず
AcademyとAgentとDecisionを作成した環境に合わせて弄ってやるだけでいい。
ML-AgentをgithubからダウンロードしてきたところについてるTemplateフォルダに雛形があるのでそれを弄れば基本出来るはず


CollectState()関数でStateにおける特徴をリストにして返す{{Brain(Script) InspectorのStateSizeはここの特徴の数を入力する}}
AgentStep(float[] act),AcademyStep()はFixedUpdate()で毎回呼ばれる処理を記述する
Agent.done=trueで今回の試行が終わったフラグを立てる
AgentReset(),AcademyReset() はdone=trueになった後次の試行へ移るときに呼び出される
Agent.rewardは報酬
Decisionについてはよくわかってない

学習時にはml-agents\python内のppo.pyかPPO.ipynbを使えばいい
exeファイルのレスポンスがタイムアウトすることがよくあった。私の場合jupyter notebookを再起動したらタイムアウトしなくなったのでよくわからない。


wikiにはSplashImage(ゲームを開始したときに出てくるロゴ)の表示はOFFにしておけと書いているが、nonProUnityではオフにできないしONにしてても起動に時間がかかるようになるだけでちゃんと動くから大丈夫。

AndoridでInternalmodeを動かすことに成功したのでこれを利用したAndroidアプリも作れる

f:id:daruma3940:20160520223745p:plain
まあwikiに書いてないことで特記しておくことはこれぐらいかだじぇ...????

行列 ベクトル 微分

f:id:daruma3940:20160520223745p:plain
これであってるはずなのだけれど間違ってたら教えてほしいのじぇ。




 A = \left(
    \begin{array}{ccc}
      a_{11} & a_{21}  \\
      a_{12} & a_{22} 
    \end{array}
  \right)
\\
 \vec{x }= \left(
    \begin{array}{c}
      x_1 \\
      x_2 
    \end{array}
  \right)
\\とすると\\
\vec{y}=A \vec{x}=\left(
    \begin{array}{c}
      a_{11} x_1+a_{21} x_2 \\
      a_{12} x_1+a_{22} x_2 
    \end{array}
  \right)
≡ \left(
    \begin{array}{c}
      y_1 \\
      y_2 
    \end{array}
  \right)
\\となる\\
\\ここで\\
\frac{ \partial \vec{y} }{ \partial \vec{x} } = 
\left(
\begin{array}{ccccc}
\frac{ \partial {y}_{1} }{ \partial {x}_{1} } & \cdots & 
\frac{ \partial {y}_{j} }{ \partial {x}_{1} } & \cdots & 
\frac{ \partial {y}_{m} }{ \partial {x}_{1} }\\
\vdots & \ddots & & & \vdots \\
\frac{ \partial {y}_{1} }{ \partial {x}_{i} } & & 
\frac{ \partial {y}_{j} }{ \partial {x}_{i} } & & 
\frac{ \partial {y}_{m} }{ \partial {x}_{i} } \\
\vdots & & & \ddots & \vdots \\
\frac{ \partial {y}_{1} }{ \partial {x}_{n} } & \cdots & 
\frac{ \partial {y}_{j} }{ \partial {x}_{n} } & \cdots & 
\frac{ \partial {y}_{m} }{ \partial {x}_{n} }\end{array}
\right)
\\というベクトルに対するベクトル微分の定義
http://www.r.dl.itc.u-tokyo.ac.jp/~nakagawa/SML1/math1.pdf

qiita.com

から

\frac{ \partial \vec{y} }{ \partial \vec{x} } 
= 
\left(
    \begin{array}{ccc}
      \frac{ \partial {y}_{1} }{ \partial {x}_{1} } & \frac{ \partial {y}_{2} }{ \partial {x}_{1} }  \\
      \frac{ \partial {y}_{1} }{ \partial {x}_{2} } & \frac{ \partial {y}_{2} }{ \partial {x}_{2} } 
    \end{array}
  \right)
= 
\left(
    \begin{array}{ccc}
        a_{11} & a_{12}  \\
      a_{21} & a_{22} 
    \end{array}
  \right)
=
 {}^t A
\\
\\
\\





\\また\\

\vec{a}_{1}=
\left(
    \begin{array}{ccc}
      a_{11} \\ a_{12}
    \end{array}
  \right)
,
\vec{a}_{2}=
\left(
    \begin{array}{ccc}
      a_{21} \\ a_{22}
    \end{array}
  \right)
\\とすると\\
A=\left(
    \begin{array}{ccc}
      \vec{a}_{1} & \vec{a}_{2}
    \end{array}
  \right)
\\であるので\\


\frac{ \partial f( \vec{x}) }{ \partial \vec{x} } 
= 
\left(
    \begin{array}{ccc}
      \frac{ \partial  f(\vec{x}) }{ \partial {x}_{1} } \\ \frac{ \partial f(\vec{x}) }{ \partial {x}_{2} }
    \end{array}
  \right)

\\より\\


\frac{ \partial \vec{y} }{ \partial A } 
= 
\left(
    \begin{array}{ccc}
      \frac{ \partial \vec{y} }{ \partial \vec{a}_{1} } & \frac{ \partial \vec{y} }{ \partial \vec{a}_{2} }
    \end{array}
  \right)
= 
\left(
    \begin{array}{ccc}
        \left(
    \begin{array}{ccc}
      \frac{ \partial {y}_{1} }{ \partial {a}_{11} } & \frac{ \partial {y}_{2} }{ \partial {a}_{11} } \\
       \frac{ \partial {y}_{1} }{ \partial {a}_{12} } & \frac{ \partial {y}_{2} }{ \partial {a}_{12} } 
    \end{array}
  \right)
&&
\left(
    \begin{array}{ccc}
      \frac{ \partial {y}_{1} }{ \partial {a}_{21} } & \frac{ \partial {y}_{2} }{ \partial {a}_{21} } \\
       \frac{ \partial {y}_{1} }{ \partial {a}_{22} } & \frac{ \partial {y}_{2} }{ \partial {a}_{22} } 
    \end{array}
  \right)
     
    \end{array}
  \right)
\\
= 
\left(
    \begin{array}{ccc}
        \left(

    \begin{array}{ccc}
      x_1& 0 \\
      0 & x_1
    \end{array}
  \right)

&
\left(
    \begin{array}{ccc}
     x_2& 0 \\
      0 & x_2
    \end{array}
  \right)
     
    \end{array}
  \right)
\\

= 
\left(
    \begin{array}{ccc}
      x_1 I & x_2 I
    \end{array}
  \right)
= 
\left(
    \begin{array}{ccc}
      x_1 & x_2 
    \end{array}
  \right)
I
= 
\left(
    \begin{array}{ccc}
      x_1 & x_2 
    \end{array}
  \right)
=

 {}^t \vec{x}

電王トーナメント申し込み始まったみたいなのじぇ

f:id:daruma3940:20160520223745p:plain
電王トーナメントの申し込みが始まったみたいなのじぇ

denou.jp


f:id:daruma3940:20160521003616p:plain
他人事みたいに言うわね....
f:id:daruma3940:20160520223745p:plain
出ようかなぁどうしようかなぁ....

なんかほかにもいろいろやりたいことあるし,研究もあるし、以前みたいに熱中して開発する気は起らないし、
でも研究ややりたいことがいっぱいあるなんて言っても、どうせごろごろして無駄な時間を過ごしそうだし。出なかったら出なかったで後悔もしそうだし。
でも参加申し込みしたら大会のことで頭のメモリが常に何割か占拠されるし、
何日もウィンウィンうなるPCの前で寝るのもつらいし、
学習にi7 4770Kを占拠されて 開発はクッソ遅いi5でしないといけないのもつらいし....

もし出るとしてもまだWCSC27後からRootStrapを成功させたことによるR上昇100しかないし
おそらくトップ勢はR200ぐらい上げてくる気がするので今のままだと相対的に後ろに下がっていることになるし...

まあうちもR200ぐらい上げようと思えばKPPを取り入れてRootStrapで学習させればそれぐらいは行く気がする(簡単だとは言っていない)んだけれど
どこまで巨人の肩に乗り続けるんだよって感じもするし....
一応次の大会でやりたいことはあるけれど、モチベーションを保てるかどうか....
f:id:daruma3940:20160520223745p:plain
とりあえず申し込み期間はまだあるのでもう少しなやもうじぇ??

カオスについて2

daruma3940.hatenablog.com
f:id:daruma3940:20160520223745p:plain
院試も終わったし前回の続きをかくのじぇ
前回2次元調和振動子は摂動を加えられなければトーラスになり、
摂動が弱い場合でも元のトーラスは変形して新しいトーラスができ、摂動が強くなるとトーラスが壊れてカオスが発生するらしいことについて書いたのじぇ。

f:id:daruma3940:20160520223745p:plain
ここでポアンカレマップについて説明するのじぇ。
ポアンカレマップとはトーラスの断面を切り取ったものなのじぇ。
f:id:daruma3940:20170824112349p:plain

2次元調和振動子の場合は変数はx,y,px,pyの4つであるけれど、
ハミルトニアンが一定という条件
{ \displaystyle
H_0 = \sum_{i=1,2} (p_i^2+\omega_i^2q_i^2)/2 
}
から変数(例えばpy)を一つ消せ、
y=0の点となった時に点をプロットするということをすれば
変数はx,pxだけになるので2次元平面に点をプロットできるようになるというわけのじぇ。
↓こんな感じ
f:id:daruma3940:20170824114046p:plain
f:id:daruma3940:20160520223745p:plain
ここで非加積分系について考えようじぇ。
{ \displaystyle
H_0 = (p_x^2+p_y^2+x^2+y^2)/2 -kxy
}

f:id:daruma3940:20170824142059p:plain
周期軌道から少しずらしたものがポアンカレマップのどこに来るのかを考えると

{ \displaystyle
\delta x(T)=M \delta x(0)
}

となるのじぇ。
このMは安定性行列と呼ばれ、保存系ではdet(M)=1となるのらしいのじぇ(よくわかってない)
これはリウビルの定理と呼ばれるものなのじぇ(リウビルの定理って正準変換のところで出てきたけど...)
{ \displaystyle
 M = \left(
    \begin{array}{ccc}
      a & b  \\
      c & d 
    \end{array}
  \right)
}
とすると
{ \displaystyle
M \left(
    \begin{array}{c}
      1 \\
      0 
    \end{array}
  \right)
=
 \left(
    \begin{array}{c}
      a \\
      c 
    \end{array}
  \right)
}
{ \displaystyle
M \left(
    \begin{array}{c}
      0 \\
      1 
    \end{array}
  \right)
=
 \left(
    \begin{array}{c}
      b \\
      d 
    \end{array}
  \right)
}
となり、
f:id:daruma3940:20170824143433p:plain
            ↓
f:id:daruma3940:20170824143442p:plain

面積S=ab-bc=det(M)=1
になるのでdet(M)は面積が何倍されたかを表すのじぇ。
ここで行列Mの固有値について考えると
{ \displaystyle
det(M)=\lambda_1 * \lambda_2=1               

\\
\lambda_1 = 1/\lambda_2                
}

ここでもし固有値が実数で重解ではない場合、どちらか片方は固有値が1より大きく片方が固有値が1より小さくなるのじぇ
ということは片方は写像を繰り返すたびに指数関数的な拡大をし、
片方は写像を繰り返すたびに指数関数的に小さくなっていくのじぇ。
f:id:daruma3940:20170824163326p:plain
f:id:daruma3940:20170824163339p:plain
安定方向の直線上にない、つまり安定成分と不安定成分を持っているものは不安定性の指数関数的増加に引っ張られて不安定方向に流されていくのじぇ
f:id:daruma3940:20170824164353p:plain
もし固有値複素数の場合、
λ1はλ2の複素共役かつλ1*λ2=|λ1|^2=1なので
{ \displaystyle
\lambda_1 =exp(i \theta) \\
\lambda_2 =exp(-i \theta)
}
になりこのような行列Mは
{ \displaystyle
 M = \left(
    \begin{array}{ccc}
      cosθ & sinθ  \\
      -sinθ & cosθ 
    \end{array}
  \right)
}
であり、これは回転行列なのじぇ。
これは指数関数的に拡大しないため安定なのじぇ。
f:id:daruma3940:20170824191827p:plain
もし重解を持った場合、
λ1=λ2=1であり、このような固有値を持つ行列は
{ \displaystyle
 M = \left(
    \begin{array}{ccc}
      1 & a  \\
      0 & 1 
    \end{array}
  \right)
}
であり
{ \displaystyle
M \left(
    \begin{array}{c}
      1 \\
      0 
    \end{array}
  \right)
=
 \left(
    \begin{array}{c}
      1 \\
      0 
    \end{array}
  \right)
}
{ \displaystyle
M \left(
    \begin{array}{c}
      0 \\
      1 
    \end{array}
  \right)
=
 \left(
    \begin{array}{c}
      a \\
      1 
    \end{array}
  \right)
}

となるのじぇ。
図に書くとこんな感じ。
f:id:daruma3940:20170824191749p:plain

これも指数関数的に大きくなっていくわけではないので安定(中立安定)なのじぇ(この図あんまり理解できていない)
f:id:daruma3940:20170824191805p:plain


もし上のハミルトニアンでk=0のときこれは加積分系であり、周期軌道から少し初期値をずらしても別の周期軌道に乗るだけなのでこれは指数関数的に大きくなっていかず、中立安定的なのじぇ
しかしkが大きくなっていくとトーラスが壊れ、もしそれが複素数のほうへ壊れたのなら、安定固定点が発生し、実数のほうへ壊れたのなら不安定固定点が発生するのじぇ

このようすをkicked rotatormodelで描いてみるとこんな感じになったのじぇ
kicked rotatormodelとは
周期的に撃力を受ける回転子のことで....って書くのはめんどくさいので
http://a-phys.eng.osaka-cu.ac.jp/suri-g/phys8.html
このページを見てくれだじぇ
f:id:daruma3940:20170824160407p:plain

これがK(キックの強さ)=0の時
キックのじゃまを受けないのでひたすらθは大きくなっていくのじぇ
f:id:daruma3940:20170824160307p:plain
こんな感じ まさしく中立安定
K=0.01のとき
f:id:daruma3940:20170824160435p:plain
K=1.2の時
f:id:daruma3940:20170824160514p:plain
島みたいなのが現れているのじぇ。これが複素数の方向に壊れたやつなのじぇ。
島と島の間には実数の方向に壊れたやつが存在しているけどカオス的になってしまっているので構造はみてとれないのじぇ
ちなみにこの島を拡大すると島の中にも島が存在するのじぇ。これはいわゆるフラクタル的構造なのじぇ
f:id:daruma3940:20170824160921p:plain
さらにキックを強くするとこんな感じ
f:id:daruma3940:20170824160713p:plain


これらの図は範囲が違ったりするので注意なのじぇ
Kが大きくなると完全にカオス的になるのはすべて実数方向にずれるから??
安定固定点と不安定固定点が交互に出てくるのはなぜ???
さっきの
{ \displaystyle
H_0 = (p_x^2+p_y^2+x^2+y^2)/2 -kxy
}

についてもやってみたけどそんなにきれいなのは見られなかったのじぇ
でもまあ一応張っておくか
f:id:daruma3940:20170824183931p:plain
f:id:daruma3940:20170824172501p:plain
f:id:daruma3940:20170824195423p:plain
f:id:daruma3940:20170824223802p:plain
f:id:daruma3940:20170824223818p:plain
f:id:daruma3940:20170825113033p:plain
f:id:daruma3940:20170824223831p:plain
f:id:daruma3940:20170826195420p:plain
f:id:daruma3940:20170826184006p:plain

なんかうまくいかないなーー
なんでだろうだじぇ?


う~~んう~~ん......
f:id:daruma3940:20160520223745p:plain
分かったのじぇ
ハミルトニアン
{ \displaystyle
H_0 = (p_x^2+p_y^2+x^2+y^2)/2 -kxy
}
で表されるので運動方程式
{ \displaystyle
d x / dt=p_x / m \\
d p_x / dt =-x +ky \\
d y / dt  =p_y /m \\
d p_y / dt = -y +kx
}

これは線形なので非線形性がないためカオスは現れないのじぇ。
f:id:daruma3940:20160520223745p:plain
そしてk=0の時トーラスっぽいものが現れていないのは
これは単純な調和振動子なので一度ポアンカレマップに乗った後運動をして、もう一度ポアンカレマップに乗ったときに完全に同じ点に乗ってしまう、
つまりトーラス上の一部しか運動をしないためトーラスが現れていないのじぇ

それに変数変換をすれば、2つの2次元調和振動子に分離もできるのじぇ
{ \displaystyle

f(x)=(x^2+y^2)/2-kxy \\


\\
{\bf x}=\left(
    \begin{array}{c}
      x \\
      y 
    \end{array}
  \right)
\\
\\
 A = \left(
    \begin{array}{ccc}
      1/2 & -k/2  \\
      -k/2 & 1/2 
    \end{array}
  \right)
\\

とすると
\\
固有ベクトルから
\\
P= 1/\sqrt{2} \left(
    \begin{array}{ccc}
      -1 & 1  \\
      1 & 1 
    \end{array}
  \right)
\\
\\
f(x)= ^t{\bf xAx} \\
\\
^t{\bf xAx}={\bf (^txP) P^{-1}AP (P^{-1}x)\\=^t(^tPx)P^{-1}AP (P^{-1}x)\\=^t(^tPx)P^{-1}AP (^tPx)\\=^tx^{'}A^{'}x^{'}  } \\

\\
\\
P^{-1}=^tP= 1/\sqrt{2} \left(
    \begin{array}{ccc}
      -1 & 1  \\
      1 & 1 
    \end{array}
  \right)
\\
\\
A^{'}=P^{-1}AP=\left(
    \begin{array}{ccc}
      (1+k)/2 & 0  \\
      0 & (1-k)/2 
    \end{array}
  \right)
\\
\\

{\bf x^{'} = ^tPx}=1/\sqrt{2}
\left(
    \begin{array}{c}
      -x+y \\
      x+y 
    \end{array}
  \right)
\\

運動量は\\
P_{x'}=1/ \sqrt{2}(-P_x+P_y)\\
P_{y'}=1/ \sqrt{2}(P_x+P_y)\\

}

A'は体格化されているのでx',y'のように変数変換をすれば2つの調和振動子に分けることができるのじぇ。
つまりこれはkの値がなんであれトーラスでカオス性は持たないのじぇ

でもkが0.6の時にはトーラスが出てきてないよね???これはなんで???
そしてトーラスの中心が2つに分かれているのはなんでなんだじぇ???
~~~追記~~~
f:id:daruma3940:20160520223745p:plain
k=0.6のときトーラスが現れていないのは新しい座標では振動数がそれぞれsqrt(1+k), sqrt(1-k) になっているのでk=0.6のとき周期の比が2になるので
片方が1周する間にもう片方がちょうど2周するので同じ点に帰ってきてしまい、トーラスが見えにくくなっているのじぇ

トーラスの中心は片方の振幅がゼロになる箇所に対応するのでどの方向の振幅がゼロになるかによって2種類あるので中心は2つあるのじぇ