読者です 読者をやめる 読者になる 読者になる

ibaibabaibaiのサイエンスブログ

サイエンス中心の予定ですが,何を書くかわかりません.統計とかの話はこっちに書くつもり. https://sites.google.com/site/iwanamidatascience/memberspages/ibayukito  ツイッターは@ibaibabaibai

平方根の日付変更線

複素数をなめたらあかんで

複素関数論のテキストを持って歩いていたら,情報系の研究者に「複素数って何に役立つんですか」みたいなことを聞かれた.心外である.2重に間違っているではないか. 

まず,複素数は役に立つ.たとえば,量子力学以降の自然記述では必須だし,これがないと電気関係の技術者もたいへん不便なことになる.次に,何に役立つかとばかり聞くのは立派な数学に対して失礼である.音楽や美術についてはそんなことはいわないのに.

複素数平方根

複素数の範囲で考えると,実数だけ見ていても思いつかないようなことがいろいろ出てくる.

たとえば,平方根の符号である.

2の平方根{ \sqrt{2} }{ -\sqrt{2} }があるのは誰でも知っているが,これらは数直線上でゼロの両側の隔離された場所にあり,決して混じったり入れ替わったりはしないように思える.しかし,複素数の範囲で考えると話が違ってくる.複素数の実部を横軸,虚部を縦軸にとった平面を考えて,2からはじめて,原点0のまわりの半径2の円の上を動かしてみよう(下の図の左).すると,ちょうど一周したときに,平方根は半周しかしないので,{ \sqrt{2} }{ -\sqrt{2} }になってしまう(下の図の右).2周すると元に戻る.

f:id:ibaibabaibai_h:20151118111147p:plain

これは「複素数の掛け算では偏角は和に絶対値は積になる」ということの結果である.{a}平方根{z}は方程式{ z^2=a}の2つの解のひとつとして定義される.いま,{z}偏角{\theta}とすると,その2乗である{a}偏角{ \theta+\theta=2\theta}だから,{a}を連続的に動かすと,{z}偏角{a}のそれの半分の角度だけ変わることになる.{a}が半回転で{z}が90度回り,{a}が1周すると{z}が180度回転して,ちょうど { \sqrt{2} }{ -\sqrt{2} }と入れ替わる.

f:id:ibaibabaibai_h:20151118001906p:plain

平方根の日付変更線

このことから,ある{a}についての{z}をひとつに決めるためには,{a}の値を連続的に変えたときに,どこかで{z}の値が不連続に飛ぶ必要があることがわかる

すなわち「日付変更線」に相当するものが必要になる.これは約束事だから,原理的にはどこに入れてもよい.

たとえば,データサイエンティスト御用達のR言語複素数型では

>(-2+0.01i)^0.5
[1] 0.003536+1.414218i
> (-2-0.01i)^0.5
[1] 0.003536-1.414218i

となって,実部が負のときに,虚部の符号が変わるところで値が飛ぶ.R言語では下の図の右側で緑色の線を入れたところ(実軸の負の部分)に日付変更線があるのだ.

f:id:ibaibabaibai_h:20151118111445p:plain

いくつかの円の上の点の平方根がどうなるかをみると,次のようになる.

f:id:ibaibabaibai_h:20151117173147p:plain

ほかの言語では確かめていないが,どれも同じなのだろうか,

もっと複雑な場合は

平方根の場合は簡単に様子がわかったが,もっと複雑な対象について「入力を徐々に動かすと結果がどう動くか」みたいなことを考えるにはどうすればよいか.たとえば「5次方程式のある係数を複素数の平面で動かすとき,方程式の解(これも一般には複素数)はどう動くか」とか.

数学の本をみると,たとえば,べき級数展開(テーラー展開)を利用して,そういう議論がされている.複素数の平面の上のある点の「近く」で収束するべき級数展開を使って「その点での値」から「ちょっと離れた点での値」を計算する.次に,新しい点のまわりで展開し・・みたいなことを繰り返す.いわゆる解析接続である.こうして少しずつずらしていって,ぐるっと回って戻ってくると,関数の値はもとに戻るときもあるし,戻らないときもある.

ここでは,ちょっと違うやり方として,方程式を数値的に解くいちばん簡単な方法であるニュートン法を使う方法を考えてみよう.すなわち

ひとつ前の{a}の値に対して求めた{z}を初期値としたニュートン法によって,新しい{a}の値についての式{ z^2=a}を解き,新しい{z}の値を求める

とする.各ステップで初期値を介して前の結果が伝わることになる.ベキ級数展開の場合と同じように,これを繰り返してつないでゆく.普通の本で定義されている解析接続との関係をきちんと示すのは簡単でないかもしれないが,実際に試すのはすぐである.

どうせなら,マウスでクリックして動かしてみたい.というので,マウスでクリックした座標で決まる{a}についての{ z^2=a}の解を,その前のときの解を初期値にしてニュートン法で計算するプログラムをR言語で書いて,カチカチやって動かしてみた.

黒い点で描いたのがカチカチやって入力したほう,2色の点で描いたのが対応する「ニュートン法による解析接続もどき」の結果である(平方根の初期値が{ +1 }{ -1 }の場合についての結果を同時並行的に計算している).図は「1の平方根」の場合の結果である(上に合わせるには2の平方根にすればよかったorz.まあ同じだよね).ちゃんと半分の角度だけ回っている.

f:id:ibaibabaibai_h:20151117190313p:plain:w350f:id:ibaibabaibai_h:20151117190322p:plain:w350

ただ,これは結果の図だけを見てもあんまり面白くない.私はプログラマーとしては超最低レベルなので,普段はソースコードは決して人に見せないのだが,恥を忍んで掲載するので,Rを入れている人はカットペーストして,実際にカチカチやって試してみてほしい.ついでにもっとマシなプログラムができたら教えてほしい.

なに,カチカチやっても別に面白くないって? うーん,それでは,次回は続編として,5次方程式の解をカチカチやって動かしてみよう.

おまけ

{ z^2=a}のかわりに{ z^4-1=a}だったらどういう風に日付変更線を入れるか」「変更線の位置によらない性質を直観的に理解するにはどうしたらよいか」という話から,リーマン面が生まれて,三角関数(実数の周期をひとつ持つ)の一般化としての楕円関数(複素数の周期をふたつ持つ)の理論と結びつき,それは20世紀の数学の出発点になった.しかし,その話はきっといろんな人が書いているから,とりあえず,ここではやらない.

ソースコードと使い方

「1の平方根」の場合についてのR言語のコードで,平方根の初期値が{ +1 }{ -1 }の場合について,並列に{ z^2=a}の解を計算するようになっている(もっと汎用性があるように書けよ>自分).

実行するとふたつのウィンドウが重なって出る.マウスで上のを移動させるともうひとつ下から出てくる(技術力がないので重ならないようにできなかったorz).最初に上にあるのが,
マウスでクリックして{a}を入れるウィンドウで,下にあるのが結果{z}が出るほう.どちらも縦軸が虚軸,横軸が実軸.{a=1}(2次元座標でいえば(1,0))のあたりをまずクリックして,クリックする点を徐々にずらしていくと,本文中の最後の図のように,1回転→半回転が実演できる.

ニュートン法が100ステップでも収束しないとき(許容誤差は1.0e-04とかに決め打ち)は,そう印刷してあきらめてその値を出すようにしている(そこで終了はしない).

終了は実部か虚部の絶対値が2より大きい点をクリックすると終わるはず.終わらないうちにウィンドウを消すとRごと落ちてしまうが,落ちると終わるわな(何言ってんだか・・)

f=function(z,a)
{
w=z^2-a
return(w)
}

df=function(z,a)
{
w=2*z
return(w)
}

dev.new()
a=0+0i
plot(a,xlim=c(-2,2),ylim=c(-2,2),col=1,pch=20,cex=2)

dev.new()
th=c(0,pi)
z=cos(th)+sin(th)*(0+1i)
plot(z,xlim=c(-2,2),ylim=c(-2,2),col=1,cex=3,lwd=3)

zz=z
for(j in 1:1000)
{
  dev.set(dev.prev())
  aa=locator(1)
  points(aa,xlim=c(-2,2),ylim=c(-2,2),pch=20,col=1)
  if(abs(aa$x>2)|abs(aa$y)>2) break
  a=aa$x+aa$y*(0+1i)
  icount=0
  repeat
  {
    zzold=zz
    zz=zz-f(zz,a)/df(zz,a)
    if(sum(abs(zz-zzold))<1.0e-04){print(paste(i," converge")); break}
    if(icount>100) {print("not converge"); break}
    icount=icount+1 
  }
  dev.set(dev.prev())
  points(zz,xlim=c(-2,2),ylim=c(-2,2),pch=20,col=c(2:6))
}