4次の特性方程式の実装

前回の続きなんですが、実際にSU(4)の生成子の中の3つの実数と6対の共役複素数から、固有値を求めるための4次の特性方程式を実装し、生成子のエルミート性から4次方程式の解が4つとも実数であることを確かめたのが以下の図です。

対角成分の実数をx,y,zとし、それ以外の共役複素数をa,b,c,d,f,gとして
a,b,c,d,f,の実部をa1,b1,c1,d1,f1,g1、虚部をa2,b2,c2,d2,f2,g2とそれぞれ定義し、
15個のパラメータを-0.5～+0.5の範囲で一様乱数として、左から2列目に出現させてみました。

左から3列目は左から2列目の2乗で、一番下に15個の2乗和のルートを算出しています。
15個のランダムの値それぞれをこの2乗和のルートで割ったのを、規格化後として左から4列目に配置しました。特性方程式にはこの4列目のデータを用います。

念のため5列目には4列目の2乗を算出して一番下で和のルートを取り、1になることを確かめています。

そうすると特性多項式のグラフは右のようになり、
x軸と交わる点つまり4次方程式の解は必ず4つあることがわかるかと思います
4次方程式の解はフェラーリの方法で算出しました。

SU(4)の特性方程式の実数解の解析

前回の続きですが、SU(4)に使ったa1,a2,b1,b2,c1,c2,d1,d2,f1,f2,g1,g2,x,y,zにrand()-1/2つまり-1/2～1/2までの一様乱数を入れて様子を見てみました。








4次方程式におけるAは純虚数、Bは実数、B+AとB-Aは複素共役、その3乗根同士も複素共役の関係にいつもなるようで、足すと常に正の実数になるようでした。
A>B≧0ってことですね

つまりそれに2/3を足したCも常に正の実数で、Cのルートも同様でした。
C≧2/3
√C≧√(2/3)

一度複素数になってから実数の範囲に留まるのは興味深いですね。

SU(4)は変数が多すぎて解析できる気がしないのですが、傾向を見るとほとんど1に近い確率でこうなるようです

この傾向まで出せると、ここから先の解析は比較的楽かもしれません

Excelでマイクラ風フルル

珈琲ととんがりこーんの容積

昨日かおとといの朝、こんな入れ物のカフェオレを飲んでいたら回転体の体積のことが気になり始めて
久しぶりに回転体の体積の公式を帰宅後にググっていました。



xを変数としたx=aからx=bまでの関数f(x)をx軸周りに回転させた回転体の体積Vは

と出ましてね


あれ？こんな感じだったっけ？って思ったんですがまあ道具が見つかったので使うまでです。
(二重積分で習った気もするんですけどね)


x=0でf(0)=y1、x=Lでf(L)=y2となる直線をy(x)とすると
y=(y2-y1)/*x/L+y1
となるので、このyをVの式に代入すれば体積は求まります。


なので、まずはyの2乗を計算しますと


なので、



このようになります。


台形のときのような、上底とか下底の代わりになる上底面積や下底面積などの概念はありませんね？
これで本当に合っているのでしょうか？




＝＝＝＝＝＝＝＝＝＝＝
ということで、別のアプローチを拾ってきました。
大きな円錐の体積から、小さな円錐の体積を引き算するという方法です。


大きな円錐の高さはA2で、底面積は半径y2の円の面積。
小さな円錐の高さはA1で、底面積は半径y1の円の面積です。


とんがりこーんの先っぽを少しかじって捨てるようなイメージです。
A1とA2はいくつでしょうか？


A1は、y(x)=0となるxの絶対値です。
A2は、A2-A1=Lを満たします。



を満たすxなので

よって、A1=|x|なので

A2は



とわかります。

大きな円錐の体積はV2=πy2^2A2/3
小さな円錐の体積はV1=πy1^2A1/3なので、求めたい体積Vはその差である

です。


2つのアプローチから算出された体積Vの式を見比べてみてください。
全然違うように見えるかもしれませんが、これは同じ式です。

因数分解を考えてみましょう。

このような公式がありましたよね？同じ式なのです。


また、y1をそのままにしてy2=0にしたり、y2をそのままにしてy1=0にしてみると、
ちゃんとただの円錐の体積の式になります。

y1=y2のときはどうなるでしょうか？
これが何を意味するのかというと、円筒の体積です。
片方は円筒の式になりますね？

もう片方は0/0になってしまいますね？
ここでロピタルの定理の出番です。分母と分子をy1で微分し、y1→y2の極限をとってみますと
ちゃんと円筒の体積になります。

分母と分子をy2でも微分してy2→y1の極限でも円筒の体積になることを確かめてみましょう

構造定数、お前だったのか

先月の通院くらいのころにはっと思いついたやつ

パウリ行列同士の積あるいは(反)交換関係がクォータニオンの基底同士の積みたいになるんだったら
ゲルマン行列でできねえ道理はねえじゃねえか

っつって
e1=-iσ1
e2=-iσ2
e3=-iσ3
e4=-iσ4
e5=-iσ5
e6=-iσ6
e7=-iσ7
e8=-iσ8

[en,em]
とか
{en,em}

とかって綺麗な形に…ならねええええええ！

あのクソ忌々しい構造定数にたどり着いたぞ！！！！！ちくしょおおおおおお！

綺麗なのは実数ってことだけじゃねえか！！！

でもこれはこれで添え字3つだし、レヴィチビタ記号のときはグリッド3つだけだったから
グリッド8つの3次元配列とかやってみてもいいかもな…

ユークリッドの互除法の需要と供給が揃った

６５６３というナンバーを見かけた。

64の2乗から1を引くわけだけど

65には13、63には7が素因数として含まれているので

65*63と1001=11*13*7の最大公約数は91になるはずだ。

ところで、63*65はいくつか

64は2の3乗の2乗なので2の6乗
64の2乗は2の12乗

2の10乗が1024なので、それを4倍して4096

4096から1引いて、4095

なるほどちゃんと5の倍数になっている。

また、4+5=9なので、9の倍数でもある。

5を2倍した10を409から引くと、399

9を2倍した18を39から引くと21でこれは7の倍数なので、4095も7の倍数だし

5を4倍した20を409に足したら429で

9を4倍した36を42に足すと78、

8を4倍した32を7に足すと39になってこれは13の倍数なので、4095も13の倍数

よしOKだ。


4095を1001で割ると、4あまり91

倍数判定式の作り方

ある数10m+nがあります。m,nは整数です。
これが任意の整数Aの倍数かどうかを判定する式を作りたいとき、どうするか。



たとえば7の倍数判定だったらA=7で

-m+2nが7の倍数だったら10m+nも7の倍数


13の倍数判定だったらA=13で
m+4nが13の倍数だったら10m+nも13の倍数なんですが、
この4とか-2とか、Aを一般化したいわけです。



10m+n=A*z1
B2*m+B1*n=A*z2


としてみましょう。
(z1,z2、B1,B2も整数です)


1つ目の式を変形して、
n=A*z1-10m
このnを2つ目の式に代入して

B2*m+B1*(A*z1-10m)=A*z2
展開して
B2*m+A*B1*z1-10m*B1=A*z2

mを含む項と含まない項に分け
B2*m-10m*B1=A*z2-A*B1*z1

(B2-10B1)*m=A*(z2-B1*z1)


こうした上で、A=7だったら
B2=-1、B1=2とすると、左辺のカッコの中が-1-10*2=-21になって、両辺が7で約分できるので、判定式として成立します。

同様にA=13だったら
B2=1、B1=4とすると、左辺のカッコの中が1-10*4=-39になるため、両辺が13で約分できて、これも判定式として成立します。


つまり、以下の形に変形することにより、任意のAに対して

(B2-10B1)*m=A*(z2-B1*z1)

B2とB1を定めやすくなる、というわけです。

ためしにA=4だと
B1=1、B2=6とか
B2=2、B1=-1とかで成立します。

32が4の倍数だと判別するために、後者を採用すると
下1桁2から、それ以外3に2をかけて6にして引き算
-4になってこれは4の倍数なので、元の数32も4の倍数と、ちゃんとできています。



A=3だったら

B1=1、B2=1とすると
たとえば36だったら下1桁6とそれ以外3を足した9が3の倍数だったら、元の数36も3の倍数


A=11だったら
B2=-1、B1=1
121を相手にしてみると
下1桁1からそれ以外12を引いて11なので、元の121も11の倍数


とできます。



任意の数の倍数判定に使うことができます。
が、速度は遅いし、あまり役に立たないです。
2のべき乗と3のべき乗と5のべき乗と11のべき乗などを区別しません。

ですので、2や5のべき乗の倍数判定(十進数)では必ず割り切れるのでほとんど意味がありませんし

3や9の倍数判定、「各桁を足す」に対して、完全下位互換となりますし(おそろしく遅い)
11の倍数判定「偶数桁の合計と奇数桁の合計を引く」に対しても完全下位互換になります。



たとえば1053が3の倍数だということを確かめたい場合、
1+5+3=9なのですぐにわかりますが、今回の方法では

105と3を足した108が3の倍数か知るために8に10を足した18が3の倍数か知るために
1+8=9と、数珠つなぎのような地獄を味わうことになり、大変判定が遅くなります。



なお、ある数Aの倍数判定のための式は1つではなく、無数に作りだすことができます



100m3+10m2+m1とかにしたらどうなるんでしょうね？

今日のできもの

カメラ「おい磁場おまえ、俺らどっから出入りするよｗｗｗｗ」
無限長ソレノイドコイルにて。

先生「ここに無限長のソレノイドコイルがあります単位長さあたりのインダクタンスを求めなさい」
学生「？？？無限！？単位長さあたり！？」


あくまでこの描写は近似ですが、無限に続く螺旋を表現してみました。
一定間隔で巻かれているコイルの中を、一定の速度でカメラが移動するとき
遠近法として無矛盾になっているかを確認したくて、螺旋なら数式1つで表現できるから手っ取り早いなって思ってこれにしました。たぶん大丈夫そうですね

(たぶんあるある)整数論舐めプ

昨日までの、35を13でわったあまりの判別方法

証明しようにも証明するべき式がうまく確立できなくてずっとモヤモヤしてた。

今朝さっきトイレに入って、ああそうか両辺を4倍すれば割り算入らないからいけるやん！

ってことにようやく気づいた。


m,nを整数として、
10m+n
って数があったとすると、これを13でわったあまりd1は
mod(10m+n,13)=d1

mod(4n+m,13)=4d2

のd2と等しいので、d1=d2

これを証明すればいい。
たぶんこれでかつる。


modが邪魔なので
勝手な整数L1、L2を用意して

10m+n=13*L1+d1
4n+m=13*L2+4*d2
のd1=d2である(L1=L2とは限らない)

という問題に変換できそう。


2本目の式から
m=13L2+4d2-4n
を1本目の式に代入して
10(13L2+4d2-4n)+n=13L1+d1

130L2+40d2-40n+n=13L1+d1
130L2-13L1+40d2-d1-39n=0

13(10L2-L1-3n)=d1-40d2

もしd1=d2=だったら、
13(10L2-L1-3n)=-13*3d
両辺ともに13の倍数になる


というのはやり方が強引だろうか？？


もしこの証明ができたら、13や7以外の、4とか5とか17とか11とか3とか
必要だったりまったく不要な下位互換だったりの任意の数で割ったあまりの判別方法も
作りだせるんじゃないかって期待してたんだけど
うーん…これだけだとどうなのかなあ？いまいちピンとこない



今日は起きたばっかりでストレスもほとんどなかったから
整数を扱っても頭の端っこがキューっとならなかった。
どうも大学を出てから整数を扱おうとすると時々不調になる。
元々整数論みたいの好きだったと思うんだけどね、
意味もなく3や9で割ったあまりとか計算しまくってたし
でもなんかこう、あらかじめ得意だった分野をあとから学業で学ぶときの僕って
舐めプする癖があるみたいで、そういうのに限ってダメになる傾向があるような気がする

整数論とか確率とか、あるいは英語全般とかね


整数論は特に、「ああ、デジタルね～」って感じで舐めプする人が多いんじゃないかって
数学ガール読んだときに思った

昨日の別アプローチ(35を13で割ったあまりの判別)

昨日は、「下1桁を4倍してそれ以外と足して13で割る」っていうのは
本来は「下1桁以外を4分の1して下1桁と足して13で割る」から派生したオマケなんだよ

っていうのを強調したくて、つい忘れていたんだけど、別アプローチがあるんだ。


というのは、「下1桁を4倍してそれ以外と足す」までは同じなんだ。

35を13で割ったあまりの例で見てみよう。

下1桁5を4倍して20
20と3を足して23

ここから事情が変わる。
ここで23を4で割ってさらに、それを13で割ったあまりを考えるんだけど

23は4の倍数じゃない。

そこで、今度は23に13の整数倍を足し引きして、4の倍数に持ってくる。
ちょうど、23+13=36が4の倍数になっている。
36÷4は9だ。

だから、元の数である35を13で割ったあまりも9
昨日の話とつじつまが合っている。


計算尺的な楽しみはこのアプローチでも相変わらず残ってくれる
4等分するときに物差しの折り紙をしたいのだ。
都合のいいことに、折り紙をした際の誤差は、整数を扱っているからという理由で
キズモノのCDをコピーするみたいにノイズを消し去ることができる。

たとえば35を13で割ったあまりを判別する計算尺

下1桁5を除いた3を4分の1したいが、3はそもそも偶数ですらない。

なので、3に13の整数倍を足したり引いたりして、4の倍数になるまで続ける。
3+13=16
4の倍数になった。

16を4等分して4

これに下1桁5を足して9

つまり、35を13で割ったあまりは9

35÷13=2あまり9
合ってる。


本来はこうなんだ。
でも、「あまり」じゃなくて「割り切れるかどうかの判定」だけ知りたい場合は
下1桁以外を4等分するより、下1桁を4倍するほうがはるかにたやすいので、みんな（？）こっちを使うんだ

今北産業ローレンツアトラクタ

chouchoが舞えば気象兵器が儲かる
ローレンツ方程式
ほれ3行

'(ダッシュ・プライム)は微分を意味するので、この3行の方程式は連立微分方程式です。時間tで微分します。
なお、非線形ですので行列は帰って、どうぞ。

x,y,zが求められるべき関数で、x(t)、y(t)、z(t)らしいです。

p,r,bはパラメータです。
p=10、r=28、b=8/3っていうのが割りと注目されるらしいです。

今回数値計算でシミュレーションしましたが、初期値はx=10、y=12、z=15を参考にしました。

一度オイラー法での差分方程式によるシミュレーションを試しましたがすぐに変な値に収束してしまうので、ルンゲ・クッタ法でやり直したらExcelでもうまくいきました。陰的とかブッチャー配列とか知りません。古典的でやりました。



ルンゲクッタ法は以上のような手続きを取るそうです。
積分における短冊・台形近似をシンプソンの公式にしたやつの微分方程式版と言われて言葉ではなんとなく把握したけども、いまいち理解していません。

ですので、応用が利かず、たとえばy'=-xz+rx-yの場合任意の関数fはf=-xz+rx-yなのでf(x,z,y,t)なんでしょうが今はx,zのことは考えなくていいとして
yにxとかzとか掛け算されてたらどうしよう！？ｇｋｂｒとか思ってたんですが3本の式とも杞憂に終わってしまって助かりました。
たとえば
y'=xyとかだったらどうしよう？って話でした！


(PC版ではクリックすると原寸大が見れます、スマホ？ピンチアウト(スワイプみたいなあれ)でなんとかなるんじゃないすか(適当))
 

xの構成方法と、xに関するk1～k4の計算方法を図に記述しました。
y,zも同様です

この状態で、h=0.01、n=2000くらいまで計算させました。

それからですね、

こんな感じの3次行列を用意します。
「x軸回転」のすぐ下にdegreeの角度θ、そのさらに下にradの角度θを格納し、行列の意味は

こんな感じです。これをy軸回転、z軸回転で同様に繰り返し

このような感じに並べます。
一番下の3次行列はx,y,z、3軸回転行列の積で

このような意味合いです。
=MMULT(MMULT(x軸:x軸,y軸:y軸),z軸:z軸)
という演算をやっています。配列全体に染み渡らせたいときは、演算させたい範囲を囲ってｃｔｒｌ＋shift＋エンターです。
なお、sinの符号は、片方がもう片方の逆符号であればどっちがプラスでもいいです。どうせ逆回転になるだけなんで。


そうして出来上がった3D回転行列を、
元のx,y,zに作用させましょう。
=MMULT(x:z,$合成行列)

そうして出来上がった新しいx,y,z座標のxとyだけをグラフに反映させると

このようになります。

x→y→z軸の順番で、180度ずつ回転させてますので、全部の回転が済めば元の回転角に戻ります


この図が動く仕組みを説明しましょう

まず
now()-today()
を入力します。
関数now()は時刻込み今のの日付のシリアル値を出力します。
関数today()は、時刻抜きの今日の日付のシリアル値を出力します。

シリアル値というのは、1900年1月1日の0時0分0秒を「0」とし、24時間を「1」と定義した表現方法なので、
たとえば同日の朝6時というのは0.25という数値になります。

６：００と打って、書式をシリアル値から数値に戻すと、0.25になっているはずです。
また、1900/2/1と「イコールなしで」打てば、1月1日から数えて32日目なので、書式を数値に戻すと32になるはずです。
これらを踏まえて、日付に半角スペースを入れてから時刻を入力すると、たとえば
1900/2/1 6:00:00
でしたら、32.25という数値になります。


now()はたとえば今が2019/6/4の10:57だったら
2019/6/4 10:57のシリアル値を算出するので

43620.45758

ぐらいの数値になります。

実はこの数値、一番小さい桁を見ると、ミリ秒単位ぐらいで動いています。
F4を連打するとうごめいているのがわかるかと思います。

これを動画の媒介変数に用いることができます。

たとえばnow()の値を700万倍に拡大してグラフに反映させることができれば、目に見えて動きます。
が、整数部が5桁もあるので、700万倍したらオーバーフローしてしまいます。

これを解決する手っ取り早い方法として、(色んな方法がありますが)「関数today()で関数now()を引き算する」というのがあげられると思います。

必ず1未満の数値になりますし、引き算という演算を行っているため、Excelはシリアル値ではなく数値と勝手に解釈して出力してくれます。


今回は、180度の回転を3回繰り返したいので、180×3=540°の周期が必要です。
ですので、now()-today()を700万倍した数値を540で割り算した剰余(あまり)が役に立ってくれます。
剰余の関数はモジュロ演算なので、modという関数です。
A=mod(700万*(now()-today()),540)
という演算を行います。


x,y,z軸の回転角θ1、θ2、θ3のdegreeにそれぞれ
θ1=IF(A<=180,A,180)
θ2=IF(A<=180,0,IF(A<=180*2,A-180,180))
θ3=IF(A<=180*2,0,A-180*2)

と入力すると、上のように動くグラフが作れます




バタフライ効果という名前の由来そのものが
・「蝶が舞えば気象兵器が儲かる」
と
・ローレンツアトラクタのグラフが蝶のように見える
の2種類あるというのがなんとも、
2つの過去が1つの未来に収束したみたいで興味深いです

ｼｭﾀｯL・ψ・卍=3ｹﾞｰ

GWのハイライトは、複素解析の同じ例題が2通りできたら楽しかったです

この複素積分が楽しかった。

の線Cに沿って


を計算するんですけどね、まずはコーシーの積分公式を使った方法から。

 が の範囲内でごっそり正則なので、これをf(z)と置きます。

コーシーの積分公式では以下のように

積分が微分で表されてしまうので、今回の例題の場合はn=0、α=-iとすると、簡単に

解を導出できてしまいます。




＝＝＝＝＝＝＝＝
それでは、もう一つの方法、コーシーの積分定理と部分分数分解を用いた方法もやってみることにしましょう。

2つではなく3つの積になっているのですが、まずは(z+3)と(z-3)を部分分数分解して、(z+i)に分配して掛け算されるとみなします。


この式を恒等的に満たすAとBを見つけます。

A(z+3)+B(z-3)=z
(A+B)z=z
なので、A+B=1
もう1つの式が
3(B-A)=0
なのでA=B
これを1つ目の式に代入すると、A=B=1/2が出てきます。


となったので、分配して、再度部分分数分解します。



つまり元の積分は以下のようになりますが


線Cで囲まれた領域の中にz=±3はないので、1/(z±3)の積分は正則によりゼロになります。

よって元の式は以下のようになり、

結局、さきほどの積分公式(積分定理ではなく)を使った計算結果と見事一致します。



コンピュータの技術の積層のボトムアップ感と違って
複素解析のまほう！はどんどん略されていくトップダウン感がありますね

複素解析つれづれ

今日は仕事が疲れたので、ファミマで複素解析の勉強しませんでした。

得た知識が曖昧なので、今日は割りと定性的なネタで日記をつぶします。

昨日の寝る前に、対数で何が得られるのか？対数という演算で何をやっているのか？

ということを考えていたら、

指数(とかオイラーの公式とか)：極座標→直交座標変換

対数：直交座標→極座標変換

という意味合いを見つけました。(ちょう今更なので、デジャブってるかもしれませｎ)

そして、対数では、極座標に変換したデータを、実部に絶対値、虚部に偏角として出力し

足し算をしている！

ここちょっと驚きました。絶対値と偏角は、足し算できる間柄だったんですよ！！

じゃあこれを物理の現象に翻訳したらどうなる！？

僕は電気工学出身なので、一番身近なものだと、インピーダンスとかアンプのゲイン・位相とかそっちになります。量子力学はいまいちよくわかりません。

対数を取ると、電気工学的な意味合いでは、ゲインと位相を足し算していることになります。
まあ、無次元量同士なら足し算しても問題ないですし、実際無次元量として足し算しています。

ゲインとかもデシベルにしてから位相と足し算しますしね。
デシベルにする時点で、基本となる増幅率で割り算してから対数を取っているので、もちろん無次元量になってます。

興味深いのは、ゲインも位相も、周波数特性としてボード線図にすると
ゲインは両対数、位相は片対数を取って、必然的に対数(無次元)にしている点です。

おそらく、縦軸にゲイン・横軸に位相を取ったような作図があるはず…
ニコルス線図かナイキスト線図か、それともスミスチャートか。なんだったかな





ところで、複素積分するとよく出てくる2πiとかいうやつ。
なんで積分した結果に出てくる語尾みたいのが、純虚数なんすかね？？

僕の中のイメージではまだ、複素平面をぐるりんと回して右ネジの法則～みたいなふわっとしたイメージしかないので
もしそれなら複素平面の紙面に垂直な成分～？って稚拙なイメージしか出てきません。

しかし、あくまで元が2つの複素数を扱っているわけだから、複素平面の外側を考えることに意味がないのもなんとなくわかって、モヤっとしますね。


この辺、ベクトル解析とかだったらどうなるんだろう？
アンペールの法則とかに相当するよね？たぶん。あるいはストークスの定理かなあ？
本にはグリーンの定理の複素版って書いてた気がするけど、グリーンの定理がいまいちわからない


線積分や面積分についてもまだまだ知らないことばかりだし、回転楕円体の積分(二重積分？)とかとの関連もまださっぱり。


で、おとといの日記だったか、1/zの積分は有限なのに、1/z^2とかのそれ以外のやつはゼロになるってやつが妙に気になりだしてて

じゃあ、周回積分で有限になるのって、逆1乗則と密接な関係でもあるの！？とか考え始めたりして。

1/zを真四角経路で積分したったｗｗｗｗ

複素解析


この複素積分をｗｗｗｗｗ


この経路でｗｗｗｗ

経路CをC1,C2,C3,C4に分けます。

C1：z=1+it
C2：z=-t+i
C3：z=-1-it
C4：z=t-i

媒介変数tは実数で、範囲は4つとも-1≦t≦1です。


2項目、3項目、4項目の分母分子にそれぞれ、-i、-1、iを掛け算すると、
4項とも同じ数式になるので、1項目の4倍になります。



定積分を行うと、このような対数になります。
しかし中身が複素数に拡張された対数なので、気を付けて計算しなければなりません。


指数やオイラーの公式に立ち戻って考えますと、以上のようになります。
複素数の場合、わりと常に多価関数を意識する必要がありますが、今回は特に意識しなくて大丈夫でした。

与式はA-Bを計算したいので、最終的に√2の対数の項が打ち消しあって、
iπ/4が2倍され、さらに4倍されるので、



積分の答えは2πiと、昨日の定理と同じ値になることが確認できました。