ヘキサボナッチ数列のn番目のnを自然数から正負の実数まで拡張すると(前日譚)

数列自体は自然数から複素数まで一気に広がる。


その様子を、scilabとExcelを使って見てみましょう。

まずscilabに、ヘキサボナッチ数列を象徴する行列あるいはその特性方程式を多項式として入力します。

＝＝＝＝＝＝＝＝
方法その１～行列の固有値を求める方法～

scinotesに以下のように入力して実行します。

A=[ones(1,m);eye(m-1,m-1),zeros(m-1,1)]
x=poly(0,'x');
z=det(x*eye(A)-A)
w=roots(z)


＝＝＝＝＝＝＝＝
方法その２～多項式を直接作って根を求める方法～

scinotesに以下のように入力して実行します。


x=poly(0,'x');
z=x^6-x^5-x^4-x^3-x^2-x-1
w=roots(z)

＝＝＝＝＝＝＝＝

 
どちらの方法でもだいたいこんな感じになるので、赤矢印で示したwという6行2列の行列に格納された左と右が、それぞれ固有値の実部と虚部です。
 
全部選択してExcelにコピペしましょう。


これを、m=1から10まで、m=6以外にも適用してみましょう。
モノボナッチからデカボナッチまでの固有値を、異なるシートの同じセルに貼り付けます。

配列の大きさが異なりますが、先頭をE4セルに揃えておきます。



ご覧の状態で、10個のシートの中の10というシートを選択し、シフトを押しながら1というシートを選択すると、モノボナッチからデカボナッチまでの10個のシート全部が選択されるので
この状態で作業すると、選択したすべてのシートに対して同じ操作が行えます。

たとえば、D4セルにa1と入力し、セルの右下のポッチをダブルクリックすると、「オートフィル」(オートフィルタではない)が自動認識し、a10まで連れてってくれます。


ただし、10のシートでオートフィルをしたので、固有値が1個しかなくても、a10まで連れていかれてしまいます。


それから、G4からG13まで、complex(E列,F列)を行って、複素数をG列に入れる作業もオートフィルで行ってしまいましょう。

なぜ、scilabで一旦実部と虚部にわけたのかと言われますと
そのままコピペすると、scilabの予約語である虚数単位%iが、この仕様のまま吐き出されてしまうのを防ぐためでした。

G列をコピーし、値だけ貼り付け、E,F列を列ごと抹消します。

今度は、D列のa1と書かれている部分から、F列の10番目の固有値までをコピーし、
「行と列を入れ替えて」、F2を先頭にして貼り付けます。



何がしたいのかというと、このwikiの分母を作りたいのです。

F4セルにIMSUB($E4,F$3)
と入力すると、以下のような表が出来上がります。


この表を、横に掛け算して、先ほどのwikiの分母を算出したいのですが、ゼロが邪魔です。
なので、同じnのan同士が引き算する場合は、1を返すようにif文を構築します。

IF(F$3<>$E4,IMSUB($E4,F$3),1)
この式で大丈夫です。
Excelで複素数を扱う場合、文字型として格納されてしまうので、算術演算は本来できないのですが
数値型としてではなく文字型として「等しい」か「等しくないか」については判別可能なのです。

もちろん、ifの条件の中身をF$3-$E4<>0にしてはいけません。引き算という算術演算が入ってしまうからです。


10個全部掛け算してみましょう。
P4=improduct(F4:O4)
を、P4からP13まで繰り返します。

これで各分母は計算できました。今度は分子を計算しましょう。

Q1にn、R1に1などと入れておいて
Q4=impower(E4,$R$1)としてから、Q4からQ13まで繰り返すと、分子の出来上がりです。

それから、分子を分母で割ります。

R4=imdiv(Q4,P4)として、R4からR13まで繰り返します。
 
さらに、R3=imsum(R4:R13)とすると、nが0から8まで、実部も虚部も、ほぼゼロになって、9番目の数列がいきなり約1+i0になったりしないでしょうか。

n=9から20くらいまで入力しても、nが整数であれば、数列の値は実整数に近似した値が出るはずです。

デカボナッチ数列の漸化式と一般式を載せておきます。


さて、ほかのシートはどうなったでしょうか。
たとえば、「9」のシート、ノナボナッチ数列のところは
固有値が9つしかないのに枠を10個設けているので、おかしなことになっています。

9から1までのシートを選択して、一旦罫線を外しましょう。

 
13行目は残ってても割りとどうでもいいのですが、悪さをするのはO列です。
ここをゼロにしてしまうとゼロで割り算することになってしまうので、O列のO4からO13までにはすべて1を入れることにします。


同様に、シート1から8までのすべてに関してはN4からN13まで
シート1から7までのすべてに関してはM4からM13まで
シート1から6まではL4からL13までといった風に、1を代入することで、問題は解決します。

(追記：n≦0でQ13とかのQ列の「見かけゼロ」が悪さすることが判明したため、不必要な式impowerを取り除いて「0」と定義してください)


ここまでのExcelファイルをここに置いておきます。

にほんブログ村

ボナッチ表題変形アニメ

ボナッチ系(なんとかボナッチ)数列の一般式の公式を導くあの手順な
 
これ↓

をこれ↓

みたいに具体化できて、なんか集大成っぽい変形アニメができそうな気がするんだけどな
なんつうか規模がな・・・一人で大丈夫か？おっぱいない。って規模な気もするし
新しいパワポがな、久々か初めての運用だから自信がなくてな
めっちゃごちゃごちゃしそうじゃん、アニメを作るさまが。

テトラナッチ数列で類推する、ボナッチ系数列の一般式の一般化～その２～

前回の続きです。

テトラボナッチ数列の一般式を出すために計算していて、|P|までは計算が終わりました。

|P|Q11を計算してみましょう。


2列目から3列目を引いて2列目に代入しますと

1列目から3列目を引いて1列目に代入しますと

行列式を2次に縮めます。

1列目と2列目をそれぞれ、r2-r4とr3-r4で約分すると

1列目から2列目を引いて1列目に代入


|P|で割ると




同様に、|P|Q21も計算しますと


ここまでは同じ方法なので、よかったらやってみてください。
それから、また同じように|P|で割ります。









あとは、これのQ11～Q41に代入すればいいだけなので




こういう風になるわけです。
ちゃんとwikiとも一致していますね。つづく

テトラナッチ数列で類推する、ボナッチ系数列の一般式の一般化

こう書かれるわけですが(あ、rは固有値です)


Pというのが以下のようにあらわされ

Pの逆行列をQとすると


このように計算されますが
 
n番目のFしか求める必要がないのと、初期値がF3以外全部ゼロなことから

必要な計算は青の網掛け以外のところだけとなり

ここまで省略することができます。

また、Q11～Q41は

このうちの

こんな風になります。

まず|P|を計算してみましょう。
 
2列目から1列目を引いたものを、2列目に代入しますと

さらに、3列目から4列目を引いたものを3列目に代入しますと
 
4列目から1列目を引いたものを4列目に代入しますと

ここで、4行目のほとんどがゼロになったのでくくりだして3次の行列に縮めたいのですが
1行1列から数えて偶数番目に1があるので、縮める際にマイナスを付け忘れないでください。

1列目、2列目、3列目がそれぞれ、r2-r1、r3-r4、r4-r1で約分できるので

2列目から3列目を引いて2列目に代入しますと

1列目から3列目を引いて1列目に代入しますと

こうなるので3行目ががら空きになりました。
そこでまた行列式を2次に縮めます。
今回の1は1行1列から数えて奇数番目にあるので符号はプラスです。


ここでも1列目と2列目で、それぞれr2-r4とr3-r1で約分できるので、外に出してしまいますと

こうなるので、1列目から2列目を引いて1列目に代入しますと

あとはサラスの方法を使っても使わなくても迷わないと思いますので

いちおう、添え字の小さいほうから大きいほうを引く約束にして、符号を調整しておきましょう

|P|はこのようになりました。


とりあえず今回はこの辺で。つづく。

トリボナッチ数列の一般式の、分子が実数であることは確かめる

続きです



分母が実数であることはわかりました。

じゃあ分子も実数でないと、この数列は整数はおろか、実数にはなりませんよね。


そういうことで、分子も実数であることを確かめます。


分子はこちらになります。

第一項はBが実数なので、その2乗ももちろん実数。

2項目と3項目を変形してみましょう。


はい、これで、2項目と3項目の和も実数であることがわかりましたね。はい。

しかし、分母の整数倍かどうかを確かめるには程遠い様子・・・

任意のnで確かめるのはちょっと荷が重そうです。

nが3か4ぐらいでしたら、まだいけそうなんですけどねー

二項定理使って直交座標系でいくべきか、ド・モアブル使って極座標系でいくべきか
おそらく前者・・・あわよくば数学的帰納法みたいな感じで任意のnに対応できればうれしい


無理数の無理性から、約分するほかに無理数を有理数にする道はないはずですからねえ

トリボナッチ数列の一般式の分母を整理する

昨日の続きです。


トリボナッチ数列の一般式はこれでした。


この分母を整理したいと思います。

まずはBの2乗を求めましょう。

だったので、その2乗は

こうなります。

続いて、|b|^2を求めましょう。bは

このような複素数だったので、その絶対値の2乗は

こうなります。

最後に
-2BRe(b)を求めましょう。Bとbは先ほど出したのを使って

こうなります。

ではこの3つを合わせて分母にしますと

このようになることがわかりました。無理数ですが、実数ですね。
ということは、分子も実数で、nがいくつであっても、この分母の整数倍であることが予想されます。

トリボナッチ数列の一般式を行列で求める

 
トリボナッチ数列はこのようになっていて、
フィボナッチ数列では「前の2つの合計」だったのが、「前の3つの合計」に変わっています。

これを行列で表すとこのような漸化式から一般式に変わります。



ここで、行列のべき乗が現れます。これを効率的に求める方法として、対角化がありますが、
対角化の準備段階として、固有値と固有ベクトルを求める必要があります。

ここに、以下に示すAという行列があって、Aのべき乗を求めたいとしましょう。


Aの固有値が3つバラバラで、対角成分だけが固有値のみでできているJという行列があり


行列JとAとの間に

AP=PJ

という変換法則があったとしたら 

Aのn乗は次のように簡単に求められます。


このPは、固有ベクトルである縦ベクトルを、横に並べたものとなります。


まず、固有値を求めてみましょう。

固有値はこのように求めることができます。

この3次方程式の解は、カルダノの方法で


このように求まります。
ここで、λ1は実数、λ2とλ3は複素数で、λ3はλ2と複素共役であることがわかるので
B=λ1、b=λ2とおくことにします。

また、固有ベクトルは、この行列の場合
固有値がBのとき、bのとき、b*(複素共役)のときでそれぞれ

とわかっているので、AとJの間の変換のための行列Pは

であることがわかりますし、Pの逆行列は
行列式|P|と余因子adjPを用いて

こう書けるので

行列式は

こうなり

余因子は

こうなります。


一般式は

この行列の、一番下の列だけあればよいですし、縦ベクトル(1,0,0)の最初の1以外はゼロなので

青く塗った部分の計算はいりません。

結局、一般式Tnは
 
このようになります。

複素数でなおかつ、無理数なのに、
不思議なことに、トリボナッチ数列は自然数として算出されるのです。

nが自然数ならトリボナッチ数列も自然数で合ってるのですが、nは実は負の整数まで拡張でき
トリボナッチ数列自体も、負の整数まで拡張されます。


もうちょっと綺麗にできそうな気がしますね。


テトラナッチ数列を行列の固有値から求めることもできなくはないんですが
4次方程式の解である固有値そのものが、解析的に化け物じみた狂気をはらんでいるため
テトラナッチは諦めましたｗｗｗｗ
狂気の様を見たいなら、ウルフラムαに「行列の固有値として」代入してみると見れます。
なぜかx^4-x^3-x^2-x-1=0という高次方程式にすると近似値だけが返ってくるんですよね。
3次方程式x^3-x^2-x-1=0も同様です。

僕がやりたいのは、解析的に整数に落ち着くまでの姿を見たいのです；ω；

トリボナッチ数列～複素無理数が正の整数に化ける不思議～

トリボナッチ数列というのは以下のように

n番目の数とn+1番目の数を足してn+2番目の数列の値にするというもので
式にするとこのような漸化式となるのですが



これを行列で表現すると、このようになります。


この漸化式を、一般式に直すのは、行列では簡単です。

F1=1,F0=0とすると

こうすればいいだけです。

しかし、この式の中の行列のべき乗がややこしくなってくるため、
固有値・固有ベクトル・対角化の技術が必要になってくるのです。

この行列はフィボナッチ数列とは違って、対称ではないため、固有値は複素数になりえます。
ただし、非負なので、ペロンフロベニウスの定理には従います。


このような3次方程式を解くことになります。


フリーのプログラミング言語、scilabで解いてみましょう。

scilab上で、
A=[1,1,1;1,0,0;0,1,0]
p=spec(A)

と書くと、解いてくれます。

ただし、scilabはあくまでフリーの言語なので
Tを対角化された行列
Dを対角化するための行列として
[D,T]=bdiag(A)

というコマンドを指定しても、非対称(非エルミート)行列では望み通りの対角化をしてくれないかもしれません
(ただ、T*D*inv(T)を計算させたら矛盾なくAになったので、道理は通っているみたいです)


ここでは、spec(A)で得られた固有値を用いることにします。

p=spec(A)
というコマンドで、
λ1=1.8392868、λ2,3=− 0.4196434 ± 0.6062907i

という固有値を得ました。
p=
このように、列ベクトルの形で入ってます。

マセマティカによると、この3つの固有値で対角化された行列
J=
を算出する、対角化のための行列は
P=
だそうなので、

scilabに

P=[p(1)^2,p(2)^2,p(3)^2;p(1),p(2),p(3);1,1,1]
Q=inv(P)
J=clean(Q*A*P)

というコマンドを打ちこむと、対角化された行列Jを算出できます。
クリーンってコマンドで、誤差範囲の純虚数成分をカットしてくれます。

そして、行列Aのべき乗は

P*J^n*Qで算出可能なので

forループを使って、以下のように算出してみました。
Cという列ベクトルに、逐次トリボナッチ数列を吐き出してます。
forでnを増やすたびに、これまでのCの値が逐一出てしまうので、コマンドの後ろにセミコロンを入力して、計算結果の表示を省略してます。
最後の行でCとだけ書いているのは、「列ベクトルCの値を表示せよ」という意味です。

roundは、実数の範囲で出てきてしまった誤差を取り除くためのものです。
床関数floorじゃダメです。0.99999とかいうのが0になってしまいます。
型変換int32もダメです。Cという入れ物が一度整数型になってしまうと、動作テストのために「やっぱりint32なしで計算させてみたい」ってなっても「型の不一致」というエラーが出ます。融通が利くように、roundにしたほうがいいようです。

for n=1:31
B=round(clean(P*J^n*Q));
C(n)=B(1,1);
end
C

算出してみると

1.
2.
4.
7.
13.
24.
44.
81.
149.
274.
504.
927.
1705.
3136.
5768.
10609.
19513.
35890.
66012.
121415.
223317.
410744.
755476.
1389537.
2555757.
4700770.
8646064.
15902591.
29249425.
53798080.
98950096.


0,1に続く3番目からの数列ですが、ちゃんと出ているのがわかるかと思います。
1e+10とかって指数表示になるような大きな数になる手前までにしておきました。

少々無理強いしましたが、しっかりと複素無理数が打ち消されて、整数に化けています。


そう考えると、デカボナッチ数列を実装する際には、最初9つの初期値が必要で
もはや数列の体をなしていないともいえそうです。
実際、フィボナッチの初期値を変えたものとして、リュカ数が有名？ですよね


じゃあ501次方程式なんか500個の初期値を持つわけだから、いったい数列として何の意味を成すんだと思っちゃいます。

が、これは、数列と行列と高次方程式がつながっているからこそ
意味があるのではないかと思いました。


scilabは最初見たときはコマンドプロンプトのように見えました。
しかし、scipad(今のscinotes)があることを知って、ちゃんと閉じた格好して部品として保存できることがわかりました。
また、scinotesではないscilab本体にfor文を入れても、ちゃんとendと打つまで演算を待ってくれます。

ただ、新たな疑問がわいてきました。
こいつは実行ファイルを作ることができるのかと。

たとえばC言語だったら実行ファイルを作ることができて
コード次第ではいろんな拡張子のファイルを作ってそこに有用なデータを入れることができます。
bmp仕様のbmpファイルを作る、実行ファイルも作れることでしょう。

しかし、このscilabはそんなことをするように設計されてはいないのかもしれません

なんとかボナッチ行列の固有値をより早く、より正確に

 
 


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

(mを変えるだけでいくらでも：だいたい500くらいまでだったら有限時間内に終わるんじゃないすかね)

scilabの力を借りると、こんなにもあっさりと解けてしまうんですねえ
フリー言語なのにすごいなぁ
解析的な代数方程式をなかば仮想的にでも定義できるし
それを複素数の範囲で解くこともできるなんて。


まあ、超越方程式の複素根とかはさすがに無茶ブリですよね＾＾；

こはるん「低みを目指して！」

トリボナッチ以上のボナッチ系数列は非対称なんですが

2次行列以下のフィボナッチなどは実対称行列(エルミート)なので、固有値は実数になります。


特に、

の漸化式に従うモノボナッチ数列の場合は


このようになるため、固有値を求めるための特性方程式は
λ-1=0となり、λは正の実数λ=1ひとつだけになります。

固有ベクトルは
λ=1を代入した

 
1行の固有ベクトルvが1列だけあればいいので、規格化してv=1となります。

この逆行列は と、ユニタリ性(直交行列)もあるので
の対角化は


となって、モノボッチ数列の一般式は常に1だということがわかります。

もちろん、たった1つの初期値F0が1の場合の話です。
初期値が任意の値F0だった場合はすべてFn=F0となります。

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

グラム・シュミットの正規直交化・・・だと！？

昨日の、

ケットベクトルを規格化してユニタリにしようと思ったんです。実数行列のユニタリってそんな数なくね？とか思いながら

 
1列目から3列目までをそれぞれ√((a-1)^2+x0^2)、√((b-1)^2+y0^2)、√((c-1)^2+z0^2)で割ればいいのかなって安直に考えてたんです。

そしたら、なんか変なんですよ。逆行列が見るからに転置行列（エルミート共役）にならないんです。

ん？この式どっかで見たな～って思って
よく見たらこのPって、この未規格化ブラ・ケットもまたアフィン変換そのものじゃないですか。

しかもPって

この積だから
それぞれの逆行列は、逆変換(拡大→縮小、平行移動→逆方向に平行移動)と考えればいいから、Pの逆行列は積の順序を入れ替えて

こうなって当たり前で、じゃあ転置なんかしたらめちゃくちゃじゃないですか。



何か僕の認識同士が明らかに矛盾してる・・・
まあ、逆行列は昨日のと矛盾しないので、これ自体は合ってるはず。


そう思って、規格化の方法のほうを疑ってぐぐってみると
ただ単にノルムを1にするだけじゃだめみたいなことが判明しましてね
「グラム・シュミットの正規直交化」って知識が必要みたいなんですわ。
wiki見るとこれがさっぱりなんですわｗｗｗｗ
おそらく、ゆっくり時間をかけて読めば理解はできるんだろうけど
なんか「俺の知ってる規格化とチガウ！」感がすごくて、まだ読む気になれんのです；ω；

3次元アフィン変換行列のべき乗

ついつい慣例に従って、アルゴリズムの定まっているほうに堕落してしまいます。

ほんとならpixivにもアップしたほうが、というかpixivのほうを重点的にアップしたほうが宣伝効果がずっと高いとは思うのですが

ブログにアップするのとちょっとアルゴリズムが違うんですよね。


ブログだったらテキストの合間に図を挟むことが可能なのに対し
pixivのほうは、上にキャプションがあって、メインは画像という形を取るので
テキストを画像に埋め込む形も取らなければならず

ブログをはいアップしました、といってもそこからそのまんま流用してpixivに張り付けるという単純な作業ではないわけですよ。


なんかね、ブログのほうでさえ時々アップするアルゴリズムを忘れそうになるんです。

「文字を作りながらその都度計算して、その様子をキャプチャーしてアップするだけの簡単なお仕事」ってアルゴリズムに落として考えないと、最近はろくにブログも更新できなくって・・・


精神的余裕がないと、そのアルゴリズムの存在をすっかり忘れてどんよりしてしまうので困るんですよね。


pixivのアップ方法のアルゴリズムも、さっさと確立させて慣れていかないといけませんね。
といってももうpixivもだいぶ老年期に入ったようなコンテンツのようなので
何か別のSNSを探したほうがもしかしたら有意義なのかもしれませんが
僕はいつもいつも後手に回ってしまうのです。



＝＝＝＝＝＝＝＝＝
さて、今回はただの自己満足以外の何物でもないんですが

3次元のアフィン変換(でいいんでしたっけ？)行列のべき乗についての話です。


ジョルダン標準形以前に、まずこれを整えておかないとと思いまして。



こんな変換行列Aがあるとします。
x軸方向にa、y軸方向にb、z軸方向にcだけ拡大しながら、
x軸方向にx0、y軸方向にy0、z軸方向にz0だけ平行移動させる変換行列です。
この行列Aのべき乗を考えてみましょう。


平行移動を含んでいるので、3次元の変換といっても、ダミー次元のせいで4次の行列になっています。


べき乗を計算したいので、この4次の行列の固有値λをまず求めてみましょう。
(a-λ)(b-λ)(c-λ)(1-λ)=0なので、λ=a,b,c,1となりますね。
4次の行列式なのでサラスの方法が使えませんが、すでに行ごとに掃き出し終えている形なのですんなり計算できるのがわかるかと思います。




では、λにa,b,c,1を順に代入していった際の固有ベクトルを求めてみましょう。


まずはλ=a
この連立方程式(ただし永年方程式)を解けばいいので

v2=v3=0とわかります。
また、x0v1=(a-1)v4であることもわかるので
求めたい固有ベクトルvは以下のようになります。
 


次にλ=bについて解きますと
 
λ=aのときと同様に
v1=v3=0と
yv2=(b-1)v4が条件なので、固有ベクトルは
 


λ=cでは

v1=v2=0
zv3=(c-1)v4



λ=1では

v1=v2=v3=0



これらの固有ベクトル(縦ベクトル)を、λがa,b,c,1の順に横に並べると
対角化のための行列の片割れPが以下のように求められます。
 
ちなみに、固有ベクトルを4つとも規格化できるときはPはユニタリ行列になるはずです。
量子力学でいうところのケットベクトルに相当しますね。


この「対角化のための行列P」のもう片方は、Pがユニタリならエルミート共役(転置して複素共役)をとればすぐに求まるのですが

規格化してないときは、Pの逆行列を求めることになります。

その際に、Pの行列式も求めることになりますが、これは以下のようになります。



Pの逆行列は、次のように余因子展開をしてから転置して求めます。

 


ここでは、invP×A×P=Jという対角化された行列になる確認は割愛して

変換行列Aのべき乗を直接求めてみます。

A^n=P×J^n×invPとなるので、以下のようになります。



実は、この計算にはa-1とb-1とc-1という部分が分母に含まれています。
ですので、a=1またはb=1またはc=1で計算不能になってしまいます。
この「または」は排他的XORではない論理和ORなので、a=b=1やc=a=1、a=b=c=1のときも
もちろん計算不能になります。


しかし、この式をよく見ると
たとえば をa-1で割る計算だったら以下のように多項式の割り算を行って


と、aの0乗からaの2乗までの級数として表現しなおすことが可能で
ｎを一般化しても


とすることが可能です。

特に、これを使うときはa=1の場合に限られるため

とすることが可能で

の例外処理は、最悪、a=b=c=1だった場合でも

このように処理することが可能です。
もちろんa,b,cどれか1つか2つが1でないときの例外処理も大丈夫です。

まあ、等倍の平行移動というイメージを忘れていなければ、当然といえば当然の帰結です。

また、平行移動がなくそれぞれの軸方向にスカラー倍だけしたいときの計算も問題ありません。


解析的に計算すれば、これらの例外はこのままでも簡単に取り除くことができます。


しかし、もし行列のまま計算結果を得たい場合は少々異なってきます。

たとえばa=bだった場合やb=c=1だった場合などは、固有値がカブるので、
それなりの措置を取らなくてはいけません。

ランク落ちがあるかもしれないので、必ずしもジョルダン標準形になるとは限らないようなのですが。

ただ、解析的に計算した場合の結果のみを用いる際は
このように 横ベクトルの後ろに変換行列を掛け算する方式でも
このように 縦ベクトルの前に変換行列を掛け算する方式でも、議論は同じということになります。(平行移動の項の位置に注意)


このような、数値計算ではない解析計算でありがちな困ることは、汎用性がなく、「計算が閉じている」ことですね。


以下に、アフィン変換の様子の例をgifであげておきます。

 
x軸倍率→y軸倍率→z軸倍率→x軸平行移動→y軸平行移動→z軸平行移動


x軸a倍しながら平行移動→y軸b倍しながら平行移動→z軸c倍しながら平行移動


3軸同時拡大→3軸同時平行移動→3軸同時拡大しながら平行移動



あ、そうだ。忘れてましたがExcelソースファイルを置いておきますね。
ダウンロード

応援したいと思ってくれたら押してね＾＾

にほんブログ村

アフィン変換で理解し始める「ジョルダン標準形」例題

一緒に考えてみようや＼PON／

まずは、2次元(平面)同士の図形の変換「拡大・縮小」に、「平行移動」を合わせるためのダミー次元を加えた、3行3列の行列で考えてみることにしよう。



Aというアフィン変換が、以下のようにｘ、ｙ軸共通の縮尺aと、ｘ、ｙ軸それぞれの平行移動の値ｘ、ｙで定義されているとすると、



固有値λはλ=aが2つと、λ=1が1つとなるので、

対角化された行列すなわちジョルダン細胞は以下のようになる。

僕はまだ不慣れなので、nとmに0が入るか1が入るかは試行錯誤で判断している。ほかの数が入ることはない。

対角化のための行列(ブラ・ケット)Pを以下のように定義すると

AP=PJが恒等的に成り立つためには
n=m=0で、
d=g=0である必要があり(a≠1である限り)

bx+ey=(a-1)h
cx+fy=(a-1)i

を満たす必要がある。

簡単のため、h=i=e=1、f=2とすると(f=1にしなかったのは、1列目と2列目を線形従属にしないため)

bx=cx=(a-1-y)
であることが要求されるため、
Pは
 
または、全体にxをかけて
 
あるいは、簡単のためj=2を代入して(j=1にしなかったのは、一番下の行を線形従属にしないため)
 
このように定義して、対角化を行ってもいい。



なお、nかmどちらかもしくは両方が1だった場合、自己矛盾が起こることを確かめてみよう。


ここで、Pの行列式と逆行列を求めてみることにすると

行列式
 
逆行列


対角化を確かめてみると



対角化の逆算で、アフィン変換Aをn回連続で行った結果を求めると


このようになることがわかる。
a、ｘ、ｙに適当な数値か文字を入れて、n=2やn=-1とかで試してみるといい。


しかしながら、この対角化は、a=1の結果については教えてくれない。


と思うかもしれないが、実は
このような級数関係があって
こうなっちゃったりするので
a=1でも問題なくなってしまって
こんなんだったりする。ただの等倍平行移動の繰り返しだ。
 

ともあれ、結果は出せたが、途中経過でバグるので、そのための対策として
a=1の場合のジョルダン標準形を考える、という言い訳を作る。

固有値と固有ベクトルで暗号、復号ごっこ

復号が容易であることを想定した暗号wwwwww遊びかwwwwww



ごっこで終わらない可能性が素粒子レベルで存在する？