忍者ブログ
20080511~ 13と7と11の倍数の論理積は13と7と11の積の倍数である。 和ァ・・・
[4361] [4360] [4359] [4358] [4357] [4356] [4355] [4354] [4353] [4352] [4351]
昨日のブログに書いた、9次の特殊ユニタリSU(9)の実装を、scilabで行ってみます。
scilabは、信号処理に特化したオープンソースのプログラミング言語なので
複素行列の扱いが得意です。


まず80個の自由度を持つので、0から1までの間の一様(?)乱数を80個生成しましょう。
1行80列のベクトルsとして生成します。

s=rand(1,80)

次に、このベクトルのノルム(長さ)を1にします。
s=s/norm(s)

できあがったら一応、長さが1になっていることを確認しましょう。

norm(s)


9つの対角成分の入力をします。
x1~x9の9つです。

x1=s(3)+s(8)/sqrt(3)+s(15)/sqrt(6)+s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x2=-s(3)+s(8)/sqrt(3)+s(15)/sqrt(6)+s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x3=-2*s(8)/sqrt(3)+s(15)/sqrt(6)+s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x4=-3*s(15)/sqrt(6)+s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x5=-4*s(24)/sqrt(10)+s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x6=-5*s(35)/sqrt(15)+s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x7=-6*s(48)/sqrt(21)+s(63)/sqrt(28)+s(80)/sqrt(36)

x8=-7*s(63)/sqrt(28)+s(80)/sqrt(36)

x9=-8*s(80)/sqrt(36)

ここで、計算ミスがないように、検算をしてみましょう。
x1からs9までの和はゼロになるはずです。



次に、非対角成分のうちの上三角部分を定義しましょう。

a12=s(1)+%i*s(2)

a13=s(4)+%i*s(5)
a23=s(6)+%i*s(7)

a14=s(9)+%i*s(10)
a24=s(11)+%i*s(12)
a34=s(13)+%i*s(14)

a15=s(16)+%i*s(17)
a25=s(18)+%i*s(19)
a35=s(20)+%i*s(21)
a45=s(22)+%i*s(23)

a16=s(25)+%i*s(26)
a26=s(27)+%i*s(28)
a36=s(29)+%i*s(30)
a46=s(31)+%i*s(32)
a56=s(33)+%i*s(34)

a17=s(36)+%i*s(37)
a27=s(38)+%i*s(39)
a37=s(40)+%i*s(41)
a47=s(42)+%i*s(43)
a57=s(44)+%i*s(45)
a67=s(46)+%i*s(47)

a18=s(49)+%i*s(50)
a28=s(51)+%i*s(52)
a38=s(53)+%i*s(54)
a48=s(55)+%i*s(56)
a58=s(57)+%i*s(58)
a68=s(59)+%i*s(60)
a78=s(61)+%i*s(62)

a19=s(64)+%i*s(65)
a29=s(66)+%i*s(67)
a39=s(68)+%i*s(69)
a49=s(70)+%i*s(71)
a59=s(72)+%i*s(73)
a69=s(74)+%i*s(75)
a79=s(76)+%i*s(77)
a89=s(78)+%i*s(79)

要素を書きだしたら、これを行列に入れます。

いっぺんにやると凡ミスをしやすいので、1行ずつ入れることにします。
B1=[0,a12,a13,a14,a15,a16,a17,a18,a19]
B2=[0,0,a23,a24,a25,a26,a27,a28,a29]
B3=[0,0,0,a34,a35,a36,a37,a38,a39]
B4=[0,0,0,0,a45,a46,a47,a48,a49]
B5=[0,0,0,0,0,a56,a57,a58,a59]
B6=[0,0,0,0,0,0,a67,a68,a69]
B7=[0,0,0,0,0,0,0,a78,a79]
B8=[0,0,0,0,0,0,0,0,a89]
B9=[0,0,0,0,0,0,0,0,0]

A=[B1;B2;B3;B4;B5;B6;B7;B8;B9]

こうして、上三角部分は埋まりました。

エルミート行列なので、下三角は上三角の複素共役です。

A=A+A'

ここまでで、非対角成分すべてが埋まりました。

念のため、ここでエルミートであることを確認してみましょう。
A-A'=0

対角要素を入れます。
A(1,1)=x1
A(2,2)=x2
A(3,3)=x3
A(4,4)=x4
A(5,5)=x5
A(6,6)=x6
A(7,7)=x7
A(8,8)=x8
A(9,9)=x9


これでエルミート行列Aが完成したので
今一度A-A'=0が成立するのを確かめておきましょう。


ここで、いきなりexpm(%i*A)とやってもいいのですが、せっかくなので寄り道して

Aの固有値方程式と固有値を見てみましょう。

変数xを定義します。
x=poly(0,"x")

次に、固有値方程式Bを見てみましょう。

B=det(A-x*eye(9,9))
detは行列式、eye(n,m)はn行m列の単位行列を意味します。

 0.0000013 + 0.0000969x + 0.0011569x^2 - 0.0025147x^3 - 0.0637862x^4 - 0.1498150x^5 + 0.2886920x^6 + x^7 - x^9

しっかりと、7次の符号が9次の符号の逆相になっていて、8次の係数がゼロになっていますね。
係数がちゃんと実数だけでできていることも確認しておきましょう。

この9次方程式の解はすべて実数です。エルミート行列であることからそういえます。

Aの固有値がこの9次方程式の解なのですが

spec(A)を計算してみるとわかります。
  - 0.6196769  
  - 0.4362829  
  - 0.2728734  
  - 0.2080910  
  - 0.0761649  
  - 0.0165342  
    0.1421505  
    0.4476562  
    1.0398167 

エルミート行列をただの数に例えると実数
歪エルミート行列は純虚数に対応するため
expm(%i*A)は、複素平面上の単位円に相当します。

これがユニタリ=ユニットな行列と呼ばれる所以です。

C=expm(%i*A)
を計算してみて(行列指数関数などの場合は、関数名の後ろに行列(matrix)を意味するmを付ける必要があります)

spec(C)を計算してみますと
    0.5063784 + 0.8623114i  
    0.9014641 + 0.4328538i  
    0.8140662 - 0.5807721i  
    0.9899136 + 0.1416723i  
    0.9063287 - 0.4225735i  
    0.9630005 - 0.2694997i  
    0.9784271 - 0.2065925i  
    0.9998633 - 0.0165335i  
    0.9971009 - 0.0760913i  
このようになり、固有値がすべて複素平面における単位円周上にあることが確認できます。


また、このユニタリ行列は「特殊」ユニタリ行列といって、
abs(det(C))とやるまでもなく
det(C)そのものが1になるという性質があります。

つまり個の場合は9つすべての固有値の積が、ちょうど実正数の1になるということです。



なお、今回乱数にすべて正の数を用いたのは、
複素共役を見やすくするためでした。
Aの上三角の虚部がすべて正の数で、下三角の虚部がすべて負の数になります。

拍手[0回]

PR

コメント


コメントフォーム
お名前
タイトル
文字色
メールアドレス
URL
コメント
パスワード
  Vodafone絵文字 i-mode絵文字 Ezweb絵文字


忍者ブログ [PR]
カレンダー
10 2018/11 12
S M T W T F S
2 3
7 8 9
12 14 16 17
18 19 20 21 22 23 24
25 26 27 28 29 30
ブログランキング
ブログランキング参戦中
にほんブログ村 アニメブログ 深夜アニメへ
にほんブログ村 漫画ブログ SF・ファンタジー漫画へ
にほんブログ村 科学ブログ 自然科学へ
よかったらポチッとお願いします^^
最新CM
[02/28 GIC結晶ファン]
[12/01 量子きのこ]
[11/06 NONAME]
[10/12 量子きのこ]
[10/11 にしもり]
プロフィール
HN:
量子きのこ
年齢:
37
性別:
男性
誕生日:
1981/04/04
職業:
WinDOS.N臣T
趣味:
妄想・計算・測定・アニメ
自己紹介:
日記タイトルの頭についてるアルファベットは日記の番号です
26進数を右から読みます
例:H→7番目、XP→15(P)×26+23(X)=413番目。
A=0とする仕様につき一番右の桁はAにできませんのでご了承くださいズコー
バーコード
ブログ内検索
アクセス解析