というわけで、自分の大学のサーバーにファイル転送したりするのが最近はセキュリティの関係で面倒になってきているので、とりあえずこちらに続きを作っておきます。
誰が読んでも分かるページを目指しているわけではありませんので、専門家以外の方には分かりにくくて申し訳ありません。一応、
http://www.kusm.kyoto-u.ac.jp/~sugawa/bers.html
からの続きではあるのですが、そちらを読んでもらってもたぶん何のことか分からないと思います。一応
http://www.kusm.kyoto-u.ac.jp/complex/Bers/
に関連するページがあって、定義を含めて少し説明してありますが、英語で書かれています。では、以下進展に従って、日記風に書いていく予定です。（進展がなければ続きません(^_^;)）
少し前から、関連する論文の草稿を書いていて、色々と必要な計算チェックなどをしてはいたのですが、そういうことが一段落したのでようやく本格的に再開する気になったわけです。（その予備的な段階の考察はあえてこのページでは紹介しません。Heun's
differential equationに置き換えた場合に、計算の精度が上がるかどうかなど、調べていたのですが結局あまり変わらないことが分かったので、そちらの方を熱心に考えることは最近は休止しています。）
この日と、その前の日辺りからやっていたのは、1次元タイヒミュラー空間のBers埋め込みにおけるbending locusをきれいに描くということでした。ただ、そのlocusに対応する単純閉曲線を表すwordが複雑になるにつれて、精度が悪くなっていくのですがある段階からプログラムが止まってしまう（というか、止まらなくなって終了しない）という現象が起こるようになり、これはもう少しプログラムを改善する必要があるのではないか？と思い始めたわけです。
そのために、複素平面から抜く(0,1に続く)3つめの点を今までは -1にしていたのを、1/2に変更してみました。こうすると、色々なことが結構きれいに見えてくるようなので、こちらの方がいいかな？と最近は思い始めています。（山下さんも、最近の絵ではそちらを選ばれているようです。） 次に、プログラムをもう少し見直してみて、例えば誤差と扱っている数字のスケールが違った場合の処理や、点のプロットの間隔をこれまで等差数列的に取っていたのを、等比数列的に変えたりなどしてみました。
以下がモノドロミーを計算させるMathematica 3.0でのプログラムです。ただし、1/2の場合に特化させたものです。適当にいじれば、もちろん一般の場合にも計算することが可能だと思います。（ただし、アクセサリパラメータの計算はこの中には含まれていません。）
phi[x_,t_,c_,lam_]:=1/x^2/(x-1)^2/2+1/(x-lam)^2/2+(t+c)/x/(x-1)/(x-lam)
psi[x_,t_]:=phi[x,t,0,1/2]
hor1[a_,b_,t_,h_]:={z[b], y[b]}/.NDSolve[{2z'[x]+psi[x+h I,t]y[x]==0,y'[x]==z[x],y[a]==0,
z[a]==1},{y,z},{x,b}]
hor0[a_,b_,t_,h_]:={z[b], y[b]}/.NDSolve[{2z'[x]+psi[x+h I,t]y[x]==0,y'[x]==z[x],y[a]==1,
z[a]==0},{y,z},{x,b}]
hm[a_,b_,t_,h_]:={hor1[a,b,t,h][[1]],hor0[a,b,t,h][[1]]}
ver1[a_,b_,t_,h_]:={-z[b]I, y[b]}/.NDSolve[{-2z'[x]+psi[x I+h,t]y[x]==0,y'[x]==z[x],y[a]==0,z[a]==I},{y,z},{x,b}]
ver0[a_,b_,t_,h_]:={-z[b] I,y[b]}/.NDSolve[{-2z'[x]+psi[x I+h,t]y[x]==0,y'[x]==z[x],y[a]==1,z[a]==0},{y,z},{x,b}]
vm[a_,b_,t_,h_]:={ver1[a,b,t,h][[1]],ver0[a,b,t,h][[1]]}
mma={{-1-I,0},{-2,1-I}}*I/Sqrt[2] //N
mmb={{1-I,0},{-2,-1-I}}*I/Sqrt[2] //N
a[t_]:=hm[1/2,-1/2,t,1/2].vm[1/2,-1/2,t,-1/2].hm[-1/2,1/2,t,-1/2].mma
b[t_]:=hm[1/2,3/2,t,1/2].vm[1/2,-1/2,t,3/2].hm[3/2,1/2,t,-1/2].mmb
trace[m_]:=Sum[m[[i,i]],{i,Length[m]}]
ta[t_]:=trace[a[t]];tb[t_]:=trace[b[t]]; tab[t_]:=trace[a[t].b[t]]
t1:=x*y-z; t01:=x*t1-y;t11:=y*t1-x;
t001:=x*t01-t1; t011:=t1*t01-x; t101:=t1*t11-y; t111:=y*t11-t1;
t0001:=x*t001-t001;t0011:=t001*t01-x; t0101:=t01*t011-t1; t0111:=t1*t011-t01;
t1001:=t1*t101-t11;t1011:=t11*t101-t1; t1101:=t11*t111-y;t1111:=y*t111-t11
tf[x_,y_,z_]:=x
tr[t_]:=tf[ta[t],tb[t],tab[t]];
npart=50;eta=0.00001;delta=0.0001;delta2=0.000001;
locus[0]=0;coord[0]={0,0};conj[0]={0,0};
a0=Re[tr[0]];sold=0;aold=0;
If[Abs[Im[tr[0]]/a0]>delta,
Print["imaginary part of Trace(0) is too large:" tr[0]];Abort[]];
If[Abs[a0]<2,Print["|Trace(0)|<2; Error"];Abort[]];
inc=Exp[-Log[Abs[a0]/2]/npart];
snew=delta*a0*(1-inc)/(a0-tr[delta]);anew=tr[snew];
Do[
a0=a0*inc;
For[j=1,Abs[(anew-a0)/a0]>eta,j++,
s0=If[Abs[snew-sold]>delta,
snew-delta(anew-a0)/(tr[snew+delta]-anew),
snew-(snew-sold)(anew-a0)/(anew-aold)];
sold=snew;snew=s0; aold=anew; anew=tr[snew];Print["*"];
If [Abs[snew-sold]<delta2,
Print["Warning: error term
is still large but we continue the process"];Break[]
]
If [j>10,Print["Newton's method diverges!"];Abort[]
]
];
locus[i]=snew;sold=snew;aold=anew;value[i]=anew;
coord[i]={Re[locus[i]],Im[locus[i]]};
conj[i]={-Re[locus[i]],Im[locus[i]]};
snew=2locus[i]-locus[i-1];
anew=tr[snew];
Print[StringForm["i=``, a0=``, ideal value=``, error
term=``",i,aold, a0,
Abs[(aold-a0)/a0]]
],
{i,npart}];
points=Table[coord[i-1],{i,npart+1}];
conjpoints=Table[conj[i-1],{i,npart+1}];
minuspoints=Table[-coord[i-1],{i,npart+1}];
minusconjpoints=Table[-conj[i-1],{i,npart+1}];
values=Table[{Abs[locus[i-1]],Abs[value[i-1]]},{i,npart+1}];
logvalues=Table[{Abs[locus[i-1]],Log[Abs[value[i-1]]]},{i,npart+1}];
segments=Line[points];
semiallsegments={Line[points],Line[minuspoints]};
allsegments={Line[points],Line[conjpoints],Line[minuspoints],
Line[minusconjpoints]};
tempinclist=Table[Abs[locus[i-1]],{i,npart+1}]
あえて詳しい説明はしませんが、上のプログラムだとlocusを50個のセグメントに区切って計算させています。計算法はNewton法に基づいていますが、微分係数を理論的に求めるのは困難なので、そこを差分に置き換えて計算しています。
後半は対応するword（上の場合だと単に x=trace(A) )の hyperbolic locusを描くプログラムです。（最後の方のは色々なグラフィックスを出力するためのもので、多分に試験的なものを含んでいます。）一般のwordの場合は、例えば
t101とかだったら、
Expand[t101]
として、その結果をコピーして tf[x,y,z]の定義式に貼り付ければいいわけです。この101とかいうのは、実際には0から1までの2進有理数の2進展開を表すもので、0.101の意味です。こういうのも再帰的に定義できる多項式なので、本当はプログラムを組めば簡単に生成できるのでしょうが、Mathematicaの場合は引数に函数を入れることが出来ないようなので、そういうプログラムを組むのは当面は断念しています。（どなたか良いアイデアがあればお教え下さい。）
まあ、プログラム自体は比較的改善されたと思いますので、あとはひたすら計算させてデータを集めるだけですね。そのうち良い絵が出来上がったらここにも載せたいと思います。
次にさせていたのは、さらに具体的に正則二次微分 1/z(z-1)(z-1/2) の双曲ノルムを求めることで、これについては次の論文(Israel Math. J.から出版予定だそうです）
Estimates of the hyperbolic metric of the punctured plane and applications, A. Yu. Solynin and M. Vuorinen
に出ている3点穴あき球面の双曲密度を数値計算する方法を用いました。Mathematicaでその双曲密度を函数として実現するプログラムを書けば次のようになります。
G[z_]:=-1+2(z+Sqrt[z^2-1])^2 //N
Density[z_]:=(
If[Re[z]<0,w=-z,w=z];
delta=0.0001;new=-Log[2];
mnew=-Log[Abs[z]Log[Abs[z]]];old=0;mold=0;w0=w;
While [Abs[new+mnew-old-mold]>delta,
w=G[w];If[Re[w]<0,w=-Conjugate[w]];old=new;
new=new+Log[2Abs[G[w0]+1]/Sqrt[Abs[w0^2-1]]];
mold=mnew;mnew=-Log[Abs[w]Log[Abs[w]]] ;w0=w
];
Return[Re[new+mnew]])
ここで、Density[z]というのが函数になります。これは複素平面から1,-1を抜いた領域の曲率-4の双曲密度を与えます。（実は、なんでもないようですが、このプログラムのバグ取りに多大な時間を要しました。単に私が少し数式のいじり方を間違っていたために起こったミスでしたが、どこで間違えているのか分からないので、その分析に時間がかかったんですね。）
ただ、このプログラムの中の変数 newとmnewがともに非常に大きくなってしまう数で、和を取るとそれが収束するという性質を持つものなので、実は誤差評価に用いているdeltaを小さくして精度を上げようとすればするほど、実は打ち切り誤差の影響が出てきて、実際には精度が悪くなるというdilemmaをかかえています。その辺のバランスが非常に難しいですね。ともあれ、これを使って函数
r(z)^{-2}/|z(z-1)(z-1/2)|
の上限を計算するのですが（ここでr(z)は複素平面から 0,1,1/2を除いた領域の双曲密度）、これは函数としては次のようなものになります。
もう少しtechnicalなことを述べておくと、どうして3点穴あき球面の結果が4点穴あき球面の結果に使えるかですが、実は函数
f(z)=8z(z-1)+1 を考えるとこれは、我々の4点穴あき球面から上の3点穴あき球面への不分岐2重被覆になっていることが分かります。従って、上のプログラムで計算される計量を函数
f で引き戻せば4点穴あき球面の計量 r(z)が計算できるというわけです。
ここで見えている3つの穴は 0,1/2,1に対応しています。最大値を取りそうなのは従って、見えている2つのピークということになりますが、どうもこれは見たところ z=(1+i)/2, (1-i)/2のようです。実際、その部分を拡大してみると次のようになります。
これは (1+i)/2 を中心に幅 1/1000の正方形でスケールを変えて図示したものです。もし、これが正しいとすると、このノルムの値はHartogsによる古典的な結果を使えば正確に
|| 1/z(z-1)(z-1/2) ||=Γ(1/4)^8/16π^4 = 19.157071....
と表示できます。（ここで、Γは通常のガンマ函数） 理論的にこれが示されば面白いと思うのですが、まだそこまでは考えていません。いずれにせよ、このデータとあとは一番大きなinner cusp（wordで言えば、1, i に対応する基本群の元が表す曲線のhyperbolic locusの端点）の位置を計れば、タイヒミュラー空間の内半径が計算できます。手許にあるデータによれば、これは大きさが0.112889くらいの値なので、square torusに対応するタイヒミュラー空間のBers埋め込みの内半径は
0.112889×19.1570717=2.16263
という感じになります。Nehariの結果などから、2と6の間であることは分かっているのですが、こういう具体的な値が分かったのは今回が初めてではないかと思います。（結構小さい数字ですが、ある意味でlogに近いということなのかもしれません。） もちろん、同様に外半径を計算することもできるはずですが、外半径はこのように具体的にcuspの位置などで計ることが出来そうにないので（たとえば、http://www.kusm.kyoto-u.ac.jp/complex/Bers/ にあるBers sliceの絵をご覧になればおわかりでしょう）、なかなか難しい問題かもしれません。この辺は山下さんに数値的に計算結果を出して頂ければと考えています。
とりあえず、今回はこの辺にしておきましょう。
ということで、昨日までの進展についてもう少し報告しておきたいと思います。昨日は、結局なんだかんだ言いながらも、Farey三角形の頂点に対応する単純閉曲線のトレースを与える函数として、A,
B, ABのトレースだけで表す多項式を再帰的に生成するプログラムを作成しました。Farey三角形は、ある意味でtreeの構造を持っており、分岐は2本ずつですから、普通に考えて2進法による記述をするのが便利です。詳しい規則は書きませんが、それを第N世代まで対応するトレース函数を計算させるのが次のプログラムです。コンピュータの内部演算ではもちろん本質的に二進法が用いられているはずですが、Mathematicaのcommandでは、例えば2進法の足し算すらなく、いちいちコンピュータの内部演算に相当することを明示的にプログラミングする必要があるなど、非常に馬鹿馬鹿しい思いをしながら作っておりました。やはり、C++とかでプログラムを組んだ方が圧倒的に楽そうな気はするのですが、いまさらC++を勉強するのもなんだし、というところです。（とか言いながら、C++の日本語の本を2冊ばかり持ってきてはいるのですが。まあ、最後の手段として残しておきましょう。）
plus1[m_]:=(
local=Rest[m];If[local=={},Return[{}]];
While[local[[1]]==1,local=Rest[local];If[local=={},Return[{}]]];
Return[ReplacePart[local,1,1]]
)
FD[m_]:=(
If[m=={},Return[Infinity],Return[FromDigits[m,2]]]
)
mean[m_,n_]:=(
If[FD[m]>FD[n],local1=m;local2=n,local1=n;local2=m];
If[local1=={},Return[ReplacePart[local2,0,1]]];
If[FD[Rest[local1]]==FD[local2],Return[plus1[local1]],Return[Rest[local1]]]
)
partner[m_]:=mean[Rest[m],plus1[m]]
Clear[t];lev=5;
For[t[Infinity]=z;t[0]=x;t[1]=x z-y;i=2,i<2^lev+1,i++,
m=IntegerDigits[i,2];
t[i]=t[FD[Rest[m]]]*t[FD[plus1[m]]]-t[FD[partner[m]]];
]
このプログラムでは深さのパラメータがlev=5に設定されています。これを動かせばもっと多くの多項式を生成することも可能です。なお、変数に対するメモリーを節約するため(?)多項式 t[i]のパラメータ i は整数(10進法)で表しています。例えば第5世代までなら、2^5=32までの数字ですみます。なお、実際にこれをプログラムの中で使うには、t[i]は単に代数的な意味しか持っていませんから、これをevaluateする必要があります。例えば、トレース函数を tf[x,y,z]で表すとすれば、
tf[x_,y_,z_]:=Evaluate[t[i]]
として実行させればいいわけです。
例えば、lev=5で計算させてみた結果が次のページです。（多項式が2^5=32個生成されています。ちなみに、A,
B, ABのトレースをそれぞれ x, y, zと置いたときの多項式です。）
これを前に紹介したbending locusを計算させるプログラムに組み込んで、自動的にある深さまでのbending
locusを全て描かせるプログラムを作り、計算させてみました。次の絵は深さ6まで計算させた結果です。（2^6=64個の多項式のreal
locusを計算させ対称性でばらまいたので、延べ256本のbending locusが見えているはずです。）
一目見ておわかりのように、深さを一定にするとlocusは均一には分布しないようなのですね。これは、対応するwordの複雑さにも関係しているかもしれません。深さでは計れない何か別の量が支配しているようにも思えます。
★注意：できれば画面を最大にしてご覧下さい。
なお、上下・左右に対称なのはまあ納得が出来るでしょうが、例えば90度回転についてはあまり対称には見えない、ひょっとして絵が歪んでいるのでは？とお思いの方もいらっしゃるかもしれません。実際に、多少絵は歪んでいるとは思いますが、実は本質的に90度回転については対称ではありません。例えば、実軸に対応するlocusの長さと虚軸に対応するそれとでは、実際に違っています。ですから、単射性内半径を測るにはこの実軸に対応する方を採用する必要があります。（虚軸の方は少し長い。）これは、対応する曲線を考えてもらえれば自然に納得が出来る話で、実軸に対応する二つのlocusはそれぞれ
A, Bに対応しているのに対して、虚軸に対応する二つのlocusはそれぞれ AB, AB-1に対応しています。
次回は、この軸に近づいていく曲線に対応するlocusに絞って計算をしてみたいと思います。
では、予告通り軸に近づいていくようなlocusを描いてみましょう。正確には次のような曲線族を考えます。例えば、正の実軸に対応する曲線はAに対応します。これに近づいていくbending
locusは、もう一つの生成元BによるDehn twistによって与えることができます。wordで書けば、ABnに対応します。（n=1,-1の場合がちょうど虚軸になっていて、nを増やすと
どんどん実軸に近づいていくというわけです。次の絵はこの nの絶対値が20までに対応する曲線のbending
locusおよび、それらを対称性でばらまいたlocusを描いてあります。wordで書けば、AnBに対応します。（分かりやすいように、A,Bに対応するlocusも入れてあります。）
面白いことに、境界点自体は結構近づいていくようですが、中身がなかなか詰まりません。同じことを、AのABによるDehn twistの列 A(AB)nについても計算させてみると次のような感じになります。
これらだけでは分かりにくいので、これら2つの画像を前回与えた深さ6までのlocusの絵に重ね合わせてみます。少しlocusの密度に当然ながら偏りができますが、徐々にタイヒミュラー空間の姿が浮かび上がってきました。
続いて、アクセサリパラメータについて計算をしたので、その結果を報告しておきます。とりあえず、アクセサリパラメータとは何かを復習しておきます。（文脈によって定義はまちまちなので、きちんと定義を与えておいた方がいいでしょう。） aを0,1とは異なる複素数とし、Yを複素平面から3点 0,1,aを抜いて得られる領域とします。（無限遠点もリーマン球面から抜けていると思えば、これは4点穴あき球面とも思えます。） この領域の単位円板からの正則な普遍被覆写像を p:D→Y とします。そこで、この写像の局所的な逆写像を考えてそのSchwarz微分を取ります。すると、Schwarz微分の基本的な性質から簡単に分かるように、このSchwarz微分 {p-1,z}はこの逆写像の分枝や被覆写像の取り方にはよらずに定まり、Y上の正則函数になることが分かります。しかもYの境界における挙動も局所的には 1/2z2 を最低次数とする有理型になることが分かりますので、これらのことからこの函数は
の形に表せることが分かります。ここで c は定数です。この定数は a のみによって定まる定数ですので、c(a)と書くことができます。この定数をアクセサリパラメータと呼ぶわけです。このパラメータは一般には計算するのが非常に難しいとされており、少し文献を探してみましたが計算法について書かれているものは見あたりませんでした。ただ、0,1,∞の置換を引き起こすメビウス変換の作用に関してはcの挙動が完全に分かるので、そのことからいくつか特別なaに対する値は知られています。少しそれらの式を紹介しておきましょう。
たとえば、最初の式で a=1/2を代入すると、c(1/2)=0が得られます。また、aが実数なら c(a)も実数、aの実部が1/2ならば、c(a)は純虚数になることなども分かります。なお、Hempel, Kraらの結果によりc(a)はaについて実解析的ではあるものの、複素解析的でないことが知られています。複素数値函数のグラフを考えるのは我々凡人にとっては非常に難しいことなので、ここでは上に述べた非常に特殊な場合、aが実数であるときと、aの実部が1/2になる場合に限ってグラフを書いてみましょう。その場合は本質的に実函数と思えるわけです。
まず易しい方からということで、aが実の場合に考えます。上の変換公式を使えば実質的に0<a≦1/2の場合にのみ考えればよいので、その場合に限定して計算してみました。なお、この場合というのは対応する1点穴あきトーラスを表現する格子が長方形になることに対応しています。我々が考えているのは次の微分方程式のモノドロミーです。
aを固定すればある特定のwordに対応するモノドロミーのトレースはパラメータ
t についての整函数になります。t=c(a)の場合にちょうどモノドロミーが面を一意化するFuchs群となるわけです。従ってその場合は全てのトレース函数の値は実となります。逆にこの性質が（少なくとも局所的には）t=c(a)の位置を特徴づけることになります。実際には全て必要というわけではなく、常に使っているように例えば一般のトレース函数はA,
B, ABのトレース函数の整係数多項式で表せますから、要はこれらのトレースが実である点を計算すればc(a)が計算できるはずです。
さて、今の場合はA,Bに対応するトレースは対称性から、tが実なら実となりますので、要はABに対応するトレースが実となるような実数tを探せばよいということになります。少なくとも
a=1/2の時はその値は0であることが分かっているわけですから、そこから出発して徐々にずらして計算を続けてみました。計算法はいつものようにニュートン法をmodifyしたセカント法によりました。
基本的にはあまりにも単調すぎてあまり面白くもありませんが、よくみると a=0の近くでは少しグラフが尖ってきているのが分かります。実はこの点ではこの函数は極限値 1/2を取ることが分かってはいますが、そこでは微分可能ではなく、傾きは無限大になっていきます。実際、Hempelは次のような評価式を証明しています。（特にこの式から、c(a)が正則でないことも分かります。）
実際、この評価式が正しいことを見てみましょう。そのために、意地悪く （1/2-c(a))×2(log a)2 /π2のグラフを書いてみました。これがa=0の近くで値1に近づけばいいわけです。ただ、a=0の近傍ではグラフが込み入っていて変化が見づらいので、さらに意地悪くaの代わりに log aをx-軸として取ってあります。（いわゆる片対数グラフです。）
ということで、これでは全然1にまで到達していませんね。もう少し頑張れば1に届く勢いがあるようには見えますが、実は普通に計算するのでは精度的にこれが限界で、a=3.2×10-13を越えたところで計算が行き詰まったのでそこまでしかデータはありません。そこでもまだc(a)の値は0.495781くらいで、何も知らなければこれが 0.5に収束すると判断は出来ないでしょう。改めてHempelの理論計算の的確さを実感します。
続いてRe a=1/2におけるアクセサリパラメータの計算結果です。この部分というのは、対応する1点穴あきトーラスの格子が菱形になることに対応しています。従ってこの時はABのトレースがt-平面の虚軸上で実数値を取り、対称性からはさらにAとBのトレースの値が等しいことも分かります。そこで虚軸に沿って Aのトレース函数が実数値を取る点を探せばよいということになります。下に、y=i c(1/2+xi) のグラフを与えました。（ここで i は虚数単位）
なんだか、直線みたいで全然面白くないですね。もう少し面白いグラフが出てくれば良かった(?)のですが。なんとなくしまりませんが、今回はこの辺にしておきましょう。