scilabには複素行列などの構造的な変数が最初から入れられるため

for文を使う用事にあまり出くわさない気がする。


でもこれならどうだ

この解をピンポイントで求めて、nを上げてってn次方程式の解を次々表示するなんてのはどうだ！？


いまんとこscilabのグラフ機能はいまいちわかってないから、Excelにデータを移そうと思ってんだけど
その際にfor文を使う用事ができそうな気がする！



たとえば、そうだな。疑似言語で書くと

for(n++){

a1=ones(1,n)
a2=eye(n-1,n-1)
a3=zeros(n,1)
A=[a1;a2,a3]
x=poly(0,'x')
q=det(x*eye(A)-A)
p=roots(q)
re(n,:)=real(p)
im(n,:)=imag(p)
}


みたいなことをやりたい

reとimはベクトルじゃなく行列になるけど、領域の半分ぐらいが空になる感じ。

もう複素二分法のためのスツルム列とか考える意味なくなっちゃうよね・・・そこなんだよな、プログラミングのちょっと虚しいところ

せっかくこの図のソースになるExcelファイルを改良しようと思ってたのに、意欲が削がれる削がれる

scilabで計算するフィボナッチ系501次方程式

a1=ones(1,500)
a2=eye(500,500)←2重ループなので多少時間かかります
a3=zeros(500,1)
A=[a1,1;a2,a3]←2重
x=poly(0,'x')
q=det(x*eye(A)-A)←ここで一番時間食いそうです
p=roots(q)
re=real(p)
im=imag(p)

これを慣れたExcelに張り付けて、
re
im
complex(re,im)
imargument←これでソートして、re対imでグラフ化


さすがにちょっと時間かかりましたね。パソコンが

ベンチマークとしていかがでしょうｗｗｗｗｗｗ


もうね、解が、ほぼほぼ線です。吹き出しみたくなってます

1001次は終わる気しなかったので止めました。


実数だけども、非対称行列だとspecとbdiagの結果が食い違いますねえ。
specのほうがもっともらしく見えるけど、
もしかしてどっちも合ってるなんてことあるかな？？
今度、実際に対角化して確かめてみましょうか。

二分法の複素版

こないだ、高次方程式の複素解を求めてる最中にふと思ったんだけど
ニュートン法や二分法にも複素版がほしいよなぁって。


どうも、重解とかガウス平面をてさぐる範囲とかの問題で、一般的な方法論はないみたいなんだけど
スツルム列とかいうのを使って、暫定的に使える方法論とかはあるみたい。

スツルム・・・どっかで聞いたよな。スツルム・リウビル問題とか。なんだっけあれ



やっぱ効率よく精度ほしいもんな。まあ手動でできなくもないけどね


10次方程式
x^10-x^9-x^8-x^7-x^6-x^5-x^4-x^3-x^2-x-1=0

の解を求めた際に、全部足しても0.5にしかならないんだもんな。ホントは1になるはずなのに。
10個の解の積のほうがマイナス1に近かった気がする

そりゃあメッシュの刻み幅0.05固定だったら限度あるって。

あのExcelファイル、要改良だよなぁ色々と

備忘録：ペンローズ図について

これはシュバルツシルトブラックホールにおけるペンローズダイヤグラムらしいです。
ほかのブラックホールではそれぞれ少しずつ図が違うみたいです。

右側の青いひし形の中身を我々の宇宙全体とすると、左側の緑色のひし形はもう1つの宇宙全体です。

白い領域はホワイトホール、白から青や緑へは一方通行で、白に戻ることはできないという意味で
矢印つきの線で描かれています。


同様に、黒い領域はブラックホールの事象の地平面の内部です。

ギザギザはそれぞれ、ホワイトホールと ブラックホールの特異点らしいです。

青の領域の中心、(t,r)=(0,0)が原点(T,R)=(0,0)です。

t,rとT,Rの変換則はt±r=tan(T±R)

(t±r)/2=tan(T±R)という流派もあるみたいですが、式がより綺麗になるほうを選びました。
2で割らないとと、端っこが±πではなく±π/2となります。


t=r=0だったらT=R=0なのはいうまでもないでしょう。っていうかさっき言いましたね。
イマ・ココです。

無限の未来t=∞、ココr=0だったら
t+r=tan(T+R)
から
t-r=tan(T-R)
を引いた
2r=0=tan(T+R)-tan(T-R)なので
tan(T+R)=∞
tan(T-R)=∞

ということは

T+R=T-R=π/2つまりT=π/2、R=0となります。

これが無限の過去t=-∞、ココr=0だったら同様に、T=-π/2、R=0となります。



イマt=0、無限右r=∞だったら
t+r=tan(T+R)
と
t-r=tan(T-R)
を足した
2t=0=tan(T+R)+tan(T-R)なので
tan(T+R)=∞
tan(T-R)=-∞

T+R=π/2
T-R=-π/2

ということで、T=0、R=π/2

同様に、イマt=0、無限左r=-∞だったらT=0、R=-π/2となります。



上下左右の点が求まったので、四辺の様子も見てみましょう。


右上にある右肩下がりの線の式は
T=-R+π/2なので、移行して

T+R=π/2となります。

T+R=atan(t+r)なので

atan(t+r)=π/2ということは

t+r=∞ということになります。



右下にある左肩上がりの線の式は
T=R-π/2なので、移行して

T-R=atan(t-r)=-π/2なので
t-r=-∞です。


左下にある右肩上がりの線は
T=-R-π/2なので、移行して
T+R=atan(t+r)=-π/2なので
t+r=-∞


左上にある右肩上がりの線は
T=R+π/2なので、移行して
T-R=atan(t-r)=π/2なので
t-r=∞


となります。


ブラックホール、ホワイトホール内部と特異点に関しては、もう少し経ったら計算してみたいです。



とりあえず、これで借りてた本を返却する準備ができました。

カー解、ライスナーノルドシュトルム解、カーニューマン解についても解いてみたいですねえ。
もう少し将来かなあ、(一般)相対論自体あまりよくわかってないからなあ

備忘録：ペンローズダイヤグラムにびびるな

時間をt、距離3つのうち1つをrとし、3次元空間を極座標で表すとすると、緯度とか経度をpやsで表すとする。

ペンローズ図(ペンローズダイヤグラム)では特に、tとrに関して、無限を有限に押し込める。

atan(t+r)=(T+R)/2
atan(t-r)=(T-R)/2

と定義して、Tを縦軸、Rを横軸にプロットすることで、時間・空間ともにゼロを含んでマイナス無限大からプラス無限大までを
-πから+πに押し込める。
正方形を斜め45度に傾けたひし形の中に、すべての観測可能な時間と空間がすっぽりと収まる。
この正方形(ひし形)の面積は(2π)^2/2=2*π^2かな。だから正方形でいうところの一辺の長さはπ√(2)ってことになる。

目盛変換はただこれだけだ。
初等関数で、ゼロを含む無限の変数値を有限の関数値にするのはアークタンジェントかハイパボリックタンジェントくらいで、ハイパボリックのほうはクルスカルスゼッケルのほうにすでに使われている。


プロならともかく、素人にはこれで十分。
思いついてしまえば、先人から学ぶ分にはなんてことはない。
びびることはなにもないから、みんなもネットにももっと軽々と書いてほしい。

全然見当たらなくて困ってるんだ。


実はこのペンローズ図、観測不可能な時空領域まで描ける。
ブラックホールやホワイトホール、アインシュタイン・ローゼンの橋の向こう側にあるとされるもう1つの宇宙まで描画可能だったりする。



気になるのは、このペンローズダイヤグラムは平面限定なのかということだ。
ミンコフスキー時空を、次元を落として3次元「時空」としているのは理解できる。

しかし、ペンローズ図の場合はどうなのだろう？
極座標で表すのは必須なのだろうか？緯度経度のsやpではなく、直交座標系のyやzとすることは可能なのか？
ブラックホールの場合はそうだろうが、もしこれがブラックホールではなく宇宙の「観測可能なほうの」果て(やホワイトホールに相当する？ビッグバン)だったら？

ペンローズダイヤグラムを(4次元は無理でも)3次元に起こすことは物理的に意味があることなのだろうか？？それとも意味をなさないのだろうか？？

気がかりなのは、「向こう側の宇宙」が隣り合って1つしか描かれないことだ。いや2つかもしれない
なぜ1つか2つなのか。

だからといって、ペンローズ図を3次元に起こしたところで、1つや2つがいくつかに増えるだけで、無数になるわけではない。(数珠つなぎにという意味では2次元のままでも無数なんだが)

向こう側の宇宙はどういう基準で選ばれてこっちの宇宙とER橋でつながるのか。


そもそもブラックホールごとにつながる宇宙は異なるのか


ああそういえば、BHに生えている「毛」次第で、ER橋の向こう側との連絡状況が異なるんだっけ
3本のアホ毛、回転・質量・電荷(もしくは色荷も追加されるかもしれない)


シュバルツシルトの場合と
ライスナーノルドシュトルムの場合と
カーの場合と
カーニューマンの場合

それぞれでペンローズ図は、やっぱ違ってくるのかな？そうなのか？？




やっぱりこの辺の分野での日本人は貧弱だな
日本語のペンローズ図と、英語のペンローズ図ではこれだけ違う。
グルーオンのときもそうだった。
日本語のグルーオン　英語のグルーオン。

ノナッチ行列とか

といきなり方程式の形で出されると、この高次方程式の解の法則性はよく見えないように思えますが、
この方程式が

こういう規則性のうえになり立った行列の固有値を求める特性方程式であるとわかると
得られる知識は一気に広がります。


まず、この行列のトレースは1なので、高次方程式の解の和は1になることがわかりますし

中身が全部正の実数であることから、ペロンフロベニウスの定理が活きて、ガウス平面における解の分布が右寄り横長であることもわかります。

det(デターミナントつまり行列式)も±1ということから、解を全部掛け算したものも±1であることがわかります。
行列式は、行列の大きさが偶数次だとマイナス1、奇数次だとプラス1だということもわかります。
行列式は、たとえば5次行列のペンタボナッチ数列なら、このような再帰的アルゴリズムで求めることができます。



また、ペロンフロベニウスの定理から、絶対値の一番大きな固有値は正の実数であることがわかるため、


この式のλに1や2を入れると、
解の1つが1以上2未満であることがわかるかと思います。

λ=1を代入した時

λ=2を代入した時

λが1と2のときで符号が変わっていますね。

このことから、固有値λをてさぐる範囲は、最大でも実部および虚部が-2から+2の間までであることがわかります。


それでは実際にこの多項式の零点を探してみましょう。

固有値の和も積も1あるいは±1ということは
おおよそ、ガウス平面で上下左右対称気味の配置になりそうな気がしますね。



実際、以上のように、2付近の正の実数の固有値が1つあり
1付近にはなく
あとは大方、複素平面上で半径1の円より少し小さい半径のところを等しく分割したような感じの配置になります。
半径1より少し小さいのは、2に漸近していく実数の固有値との兼ね合いだと思います。

多項式の次数が大きくなるにつれて、真円に近付いていくでしょうし、配置もほぼ円周上びっしりに漸近していくはずです。
が、ちょっと直感からズレるかもしれませんが、どれだけ次数が高くなっても、この固有値の和は1なんです。

どういうことかというと、正の実数1のところには解がなく、その代わりに正の実数2のところにあるからなのです。
単位円の円周すべての点を足し合わせたら打ち消しあってゼロになるはずです。
そこから、正の実数1を取り除くと、マイナス1になります。
どれだけ円周上に固有値を敷き詰めた円でも、正の実数1を取り除くとポッカリと空いてすべての和が-1になるのです。
そこに、正の実数1の代わりに正の実数2を付け足すと、-1+2で1になるのです。これがトレース＝1の由来になります。


DL用Excelファイル を載せておきます。


モノボナッチ　フィボナッチ　トリボナッチ　テトラボナッチ　ペンタボナッチ　ヘキサボナッチ　ヘプタボナッチ　オクタボナッチ　ノナボナッチ　デカボナッチ

モノボッチ　ジボナッチ　テトラナッチ　ノナッチ　ノナボッチ　ノナッチ





にほんブログ村

スーパー中之条ロゼ「ザマスくん、この髪型は校則違反ザマス。」

「横から生えてこないんです！」

「あなたもか」

「私だ」

「私だったのか。まったく気がつかなかった。」

ぼくのscilab冒険記

読み進めていって、ようやくscilabがプログラミング言語に見えてきた。

scipad(今名前違う？)とかいうメモ帳みたいなところに、命令文をまとめて書いて保存することができるらしい。まるでコマンドプロンプトのバッチ処理みたいに。


ようやくコードが閉じてくれたって感じがした。
そう、「開いたコード」に違和感があったんだ。


とはいえ、今日はもう、疲れた。
なんもする気出ない
参考の本は読んだけど、PCでの作業はほとんどなにもしてない。まじ疲れた

17000歩歩いたのに対して、昼飯が100円マックシリーズだからなあ
チキンクリスプ、コンビニのおにぎり、ファミチキみたいの、コンビニの100円サラダ、アップルパイ、これだけだからね

そういえば寝不足でもあったんだった。

(1対)^2応

あーよかった。scilabちゃんはちゃんと複素行列の逆行列は普通に出してくれた。
ようやく、Excelより手間が省ける実感がわいてきたというか

縦ベクトルを横に並べる入れ子もうまくいった。

気を付けなきゃいけないのは、行列AのA'ってやると自動的にエルミート共役になってしまうことだね。ただの転置にしたければA.'とするらしい。(ダッシュとの間にピリオドをいれた)


で、なんとかナッチ数列の行列モデルは
有向グラフのごく狭い部分集合なんだな。
いや、こじつけ以上の意味があるかどうかはわかんないけど。

非負バイナリ行列というか非負ブーリアン行列というか
有向グラフの他のバージョンには、それぞれなにかしらの数列が1対1対応してるんだろうか


あ、でもあれだな、等比数列だったらもはやバイナリとかブーリアンじゃなくなるな。
等差数列だったらなんか非線形っぽいっつーかアフィン変換みたいなことになりそう

きのこは新しい「物理的な」高次方程式のネタを得た！

トリボナッチ以上のテトラナッチとかの数列の一般式を求めるための行列って
非負だけど対称じゃないんですよね。


そうすると、エルミート行列の性質「固有値がすべて実数」が当てはまらなくて
でもペロンフロベニウスの定理は生きるんで、大きなLボナッチ数列の概算はほぼ実数の指数関数に近似するんですけど

途中計算に複素の無理数を含みながら、数列としての結果は整数っていう奇妙な状態になりましてね
フィボナッチの無理数に飽き足らず複素数まで拡張するんですよｗｗｗ


それで、ちょっと計算させてみようと思ってscilabを立ち上げたら
さすが無償！どうもエルミート行列以外の行列に対する固有値の求め方が、数値解析でもscilabちゃんはできないみたいでして＾＾；
結局ウルフラムアルファとExcelで対処しましたｗｗｗｗｗあほかｗｗｗｗｗ
せめてガウス平面を固有値てさぐったりしてなんとかならないんですかねえ・・・まあてさぐる範囲があらかじめわかんなかったら即手詰まりな感じはありますけどね・・・？



今回はほぼ確認のためだったので、サラスの公式が使える3次の行列にとどめまして
対角化のための行列Pの逆行列も、有効数字3桁でウルフラムアルファに求めてもらいました＞＜
だって複素数なんだもんｗｗｗｗｗ

あ、でもそうですね、Excelをscilabちゃんに取りこむ方法が見つかれば
逆行列ぐらいは出してくれますよね！そのくらいできてもらわないと困る！だって信号処理のために複素行列扱うプログラミング言語ですもんね！


ウルフラムアルファの弱点は、無償バージョンではデータがコピペで移動できない点です。まーた狡猾なお手前で。＾＾


あれ？でも結局サラス使わなかったじゃん。今さっきウルフラムアルファに頼ったって言ったよね・・・？
じゃあテトラナッチでもいけたやん。手間めっちゃかかるけど。


あ、それでな
ユニタリ化できたら楽なのになーって思ってサラス使ったんでした。
でもユニタリ化失敗ｗｗｗｗｗウルフラム先生がユニタリ化して出してくれなかったら不可能なんだよたぶん俺たちの中ではなｗｗｗｗｗｗ
よく見たら固有ベクトルがなんか、線形従属っぽいんですよ。だめだこりゃー



あーアホみたいなことしたｗたのしかったー







でもこの場合、固有ベクトルには規則性があって
たとえば3次方程式だったらa1～a3までの複素根があるとすると
n番目の固有ベクトルが、以下みたいにかけるってわかったのはすごいね。

して、それを即座に見抜くウルフラムさんも化けもんですわ。

ネタ切れから逃亡！Lボナッチ数列の一般式の一般化

うおおおおお！！！ネタが切れているときに友達からネタが降ってきた！親方ぁ！

一般化したLボナッチ数列の一般式を行列で求める際には
このような行列の固有値と固有ベクトルを解けばいいらしい！

特性方程式はこうなるっぽい！

ちなみにEは単位行列、nは整数、Lは自然数

3Dアフィン変換とジョルダン標準形の例題

どうもうまく理解できないので、とりあえず結果から書こうと思った。

ウルフラムアルファにぶっこんでみたところ


このどちらの形でa=b、a=c、b=c、a=b=cのいずれでもジョルダン標準形にはならなかった。

a=1かb=1かc=1のいずれかのときにジョルダン標準形になるようで

ジョルダン細胞の大きさは最大で2次のようだ。

つまりa=b=1とかa=b=c=1とやっても、細胞は3次や4次にならない。

あまりいい例題ではないのかもしれない。じゃあどうしろと・・・orz



ところで、scilabはジョルダン標準形に対応してないらしい・・・。見るからにおかしな値が出る。
matlabは対応してるらしい。さすが有償・・・(まああくまで別の言語らしいが)

scilabはグラムシュミットの正規直交化も特にやらないみたいだ。そういうのが必要になってくると、固有ベクトルのユニタリ化を諦めるっぽい。なかなか堅実な判断じゃないか。

scilabで特殊ユニタリのお勉強

 
8つのゲルマン行列σ1～σ8と任意の8つの実数s1～s8(ただし8つの2乗和は1)について

8つの生成子の線形結合A=s1σ1+s2σ2+s3σ3+s4σ4+s5σ5+s6σ6+s7σ7+s8σ8
exp(iA)が3次の特殊ユニタリSU(3)。


生成子の固有ベクトルを出す過程で、ただのユニタリU(3)も現れる。
生成子Aの固有値と固有ベクトルは行列指数関数で特殊ユニタリを解析的に求める際に必要不可欠。


4次でやると何か途中で間違えるので、3次で我慢した。


固有値をゆらして多項式をなぞったり、その多項式自体のパラメータをゆらしたりなんてことはまだscilabではできていませんTдT

っていうか歳明けてscilabの本借りなおしたばっかりなんです＞＜


「変数」ってやつの意味とか、グラフへのプロットとか、ちょっと何言ってるかよくわかんないですorz


これまでの記述をきれいなまま保存したり、データを出したり入れたりしたいです
っていうかまだforが出てこなくて、伏線として出てるだけなんですが
プログラム言語として大丈夫なんでしょうか・・・
たぶん本を読み進めると安心するんでしょうが、なんか挫折気味です

時々気分転換に、アフィン変換のジョルダン標準形とかやってみましょうか

直継「なあシロエ、エルダーテイルの世界に俺みたいな」

「やつしかいなかったら、どうする？」

シロエ「経済の概念も生命の概念もない、やることがおパンツの研究くらいしかない、ということ？」

直継「そう。何が残ると思う？」

シロエ「ノーイベント、グッドライフあたりが残るんじゃないかな」

直継「ゆゆしき問題だな」

シロエ「そうだろうか？たとえば3人の会話の途中で誰か一人が大火傷したとする。これは罪に問えると思う？」

＼PON／

直継「念話がどこまで公開されているかにもよるんじゃないかな・・・」

シロエ「じゃあ、ここに僕という直継、そっちに直継という直継がいて、すべての直継の区別がつかなかったら？」

直継「そりゃお前とアカツキの子供がマスコットキャラみたいになった契約詐欺のやつじゃん」

シロエ「じゃあ温度下げてみる？手始めに摂氏マイナス269度くらいでどう？」

直継「やめやめやめ！俺らはどうがんばってもボーズアインシュタイン凝縮なんて起きないから。」

シロエ「それがありえるかも？僕たちはヘリウムアリにはなれるかもしれないよ？」

直継「うん・・・変なとこだけゲームだからな、この世界は。」



直継「ところで、アカツキはアサシンなんだろうか？」

シロエ「必殺技で、あれは殺してるのか殺してないのか、結局よくわからないからね」

直継「どうしてこうなったんだ？(ﾃﾞﾗｯﾍﾞｽﾃﾞﾗｯﾍﾞｽ)」

シロエ「さあ？過去の人類による数知れないクレームのせいじゃないかな」




直継「ところで、図書館の司書ロボットは、下半身を概念にするべきだと思うか？」

シロエ「下半身どころか全身が形而上のものでいいような気がするんだけど」

直継「そうしたら俺は誰を愛すればいいんだ？ラブ・アカツキん？それともラブ・ライブラリ？」

シロエ「これもうわかんないよね＾＾；ん？ちょっと待って？これって三角関数？僕も直継もアカツキのことが好きってことになるような」

直継「自覚あったの！？」

シロエ「いや、ないけど。設定資料集読んでたらなんかそういうことになってた。」

直継「手をつなげば腋から孕むとかいう設定を流布したらいいんじゃないか」

シロエ「ああ～先週のモーガンフリーマンね」

べ、別にペンローズダイヤグラムを2次元に押し込めたいわけじゃないんだからね

ああそういえば、3Dでもいいんだ。4Dにさえならなければ描写可能なんだ。

3Dだったらどうなるんだろう？

ライトコーンは3Dに描けたんだから、それに対応するペンローズダイヤグラムも見てみたい