忍者ブログ
20080511~ 13と7と11の倍数の論理積は13と7と11の積の倍数である。 和ァ・・・
[1] [2] [3]
これの続きです。


 


この余因子展開の具体例は、q31とq41だけ計算しておきますので
残りは読者さんへのプレゼントにしますね。




  


 





しいて似てるとすれば、テトラナッチ数列の一般式の、βとγの項の分母がこれに相当しますね



約分のイメージはこんな感じです

6角形から5角形の辺(と対角線)を抜き取った、補グラフ?が分母をとったら、ちょうど手の指のような感じになった。そんなイメージです

拍手[0回]

PR


前回の続きです。


かばん「ヘキサボナッチ数列を具体的に言ってみよう!」

サーバル「・・・うーん」

かばん「・・・あ、それじゃあサーバルちゃん、フィボナッチ数列の具体例から始めてみようか」

サーバル「3番目が1つ目と2つ目の合計だから・・・1つ目と2つ目っていくつだっけ?」

かばん「0と1だよ」

サーバル「なんで?」

かばん「なんていったらいいかなあ。1つ目が0、2つ目が1っていうのは「初期値」っていうんだけどこの組み合わせ以外だと、3番目以降が全然違ってきちゃうでしょ。こういう数列は別の名前があって、「リュカ数」って呼ばれてるんだ」

サーバル「そうなんだ。」

かばん「「まずはなるべくシンプルに」っていう理念が数学の根っこにあるんだよ」

サーバル「じゃあ、トリボナッチだと初期値は0,0,1で、テトラナッチだと0,0,0,1ってこと?」

かばん「その通り!やったねサーバルちゃん!」

サーバル「やったぁ!」

かばん「じゃあヘキサボナッチの場合はどうなると思う?」

サーバル「ヘキサが6だから、6つ合わせた合計を7つ目にするから、0,0,0,0,0,1が初期値で、その6つの合計1が7つ目?」

かばん「すごいねサーバルちゃん!」

サーバル「えっへん!サーバルの技だよ!」

かばん(そうなのかなあ^^;)「じゃあ、100個目はわかる?」

サーバル「え・えー!?」

かばん「結構大変だよね。前回「漸化式」って言ったのはこういうことで、100個目の芋を掘りだすには、それまでの99個の芋を掘り当ててないとたどり着けないっていう弱点が、漸化式にはあるんだ。これを、いきなり「100個目はいくつですって言えるようにしたのが「一般式」って呼ばれるんだよ」

サーバル「一般式、すごいねー!」

かばん「ところで、「行列」って知ってる?」

サーバル「知ってる!PPPのチケットをもらうために並ぶフレンズのことでしょ?」

かばん「ま、まあ、それも行列、というか、行だけの行列だね^^;数学には別の意味の「行列」があって、1行や1列とは限らないんだ。たとえば、PPPを見てるフレンズさんたちが、縦4列、横10列に規則正しく並んでるとしても、それも行列なんだよ」



サーバル「なんだかラジオ体操みたいだね」

かばん「フレンズさんたちもラジオ体操って、するの?」

サーバル「夜起きてすぐのラジオ体操は気持ちいよね!ウェルカムtoようこそジャパリパーク~!今日もどったんばったんおおさわぎ!♪元気元気!

かばん(やっぱりその歌なんだ^^)

サーバル「ラジオ体操とヘキサボナッチ数列って関係あるの?そういえば数列と行列って似てるかも!」

かばん「ヘキサボナッチ数列を行列表現すると、まるで漸化式が一般式で表されてるかのような錯覚を受けるんだ。今地面に数式を書くから、見ててね」

 


サーバル「なにこれー!?」

かばん「これが、「行列」っていう数学だよ。数が行と列になって並んでるでしょ。元々は連立方程式を解くために考案されたもの、みたい。僕の宿主の記憶ではそういうことになってる。」

サーバル「れんりつほーてーしき?」

かばん「この行列方程式は、3本の連立方程式からできているんだ。」

x=1x+2y+3z
y=5x+4y+6z
z=9x+7y+8z
 
サーバル「これを行列で表すと、こうなるの?」





かばん「そういうこと」

サーバル「でも、なんか変なの~」

かばん「やっぱりサーバルちゃんもなんか変だと思うでしょ^^実は、行列っていうのは行列Aと行列Bの掛け算を入れ替えると、同じ行列cになるとは限らないんだ」

サーバル「そっか!私そこが引っかかってたんだ!積の交換法則っていうんだよね?」

かばん「サーバルちゃんよく覚えてたね!よしよし、いい子いい子」

サーバル「かばんちゃんの毛づくろいは気持ちいいなぁ。ぐるぐるぐるぐる・・・もっと撫でて~」

かばん「え、そこも撫でていいの?」

サーバル「かばんちゃんにだったらどこを撫でられても安心だよ」

かばん「えー?もう、サーバルちゃんったら~〃∇〃


ご起立ください
(中略)
座れ


かばん「はっ!それはさておき」

サーバル「積の交換法則が乱れちゃうなんて、そんなの数として大丈夫なの?」

かばん「やっぱりそこ、気になるよね。実は大丈夫なんだ。数の3大法則、覚えてる?」

サーバル「交換法則、分配法則、結合法則、だったっけ?」

かばん「そうそう。8元数っていう数学では、この3つのうちまともに機能してるのは1つだけらしいしね」

サーバル「あ、そういえばベクトル!あれも交換法則がいまいちだったよね」

かばん「ベクトル積のことだね。大きさは同じだったけど、向きが逆だったね。そうそう。数っていうのはね、必ずしもこの3大法則を守らなくてもいいんだよ」




かばん「話をヘキサボナッチに戻すと、この行列は結局、こういう風に表すこともできるよね」

サーバル「うんうん。」

かばん「その先はどうなると思う?」

サーバル「うーん・・・あ!そうか!漸化式か!芋づる式にずーっと連なって・・・!!初期値までさかのぼれるね!」




かばん「そう。n番目の数列を、行列の積で、しかも初期値までさかのぼって表現することができたね。」

サーバル「この、「同じ行列をn回掛け算する」っていう表現はなんとかならないのかなあ?」



かばん「サーバルちゃん、いいところに気が付いたね!掛け算が足し算の重ね合わせだったように、掛け算を重ね合わせた何かがあってもおかしくないって思うのは自然なことで、これを「べき乗」って呼んで、たとえばaって行列をn回掛け算するっていうのは、aのn乗って言って、a^nって表現するんだよ。nをaの右肩に乗せて表現する方法もあるんだ。」

サーバル「ほんとだぁ、漸化式だと思ってた式が、行列のべき乗のおかげで、いつの間にか一般式みたいになってる。すごーい!」

かばん「これはダジャレって言ってね」

サーバル「そうなの?」

パシャ

かばん「いい顔いただき!」

サーバル「え!?冗談!?ひどいよかばんちゃーん。ボスも私の変な顔の写真消して?」

かばん「きゃー消さないでくださーいv

ボス「(やれやれ)サーバルはお客様じゃないので要求には答えられないよ





かばん「行列のべき乗は、さっき言った通り、積の交換法則が通用するとは限らないことから、複雑になることが多いんだ。そこで、「対角化」という方法を使うんだけどね

A、P、Jっていう行列があったとして

A=P×A×inv(P)

っていう変換をすることで、Aのべき乗を計算しやすくするんだ。」

サーバル「inv(P)ってなに?」

かばん「Pの逆行列って意味だよ。P×inv(P)もinv(P)×Pも、単位行列Iになる、そんな行列のことを言うんだ。」

サーバル「単位行列?」

かばん「数字で言うところの「1」みたいなものだね。1にaを掛け算すると、aそのものになるでしょ?そういうのを「単位なんとか」って呼ぶんだ。行列の場合は



これが単位行列」

サーバル「こんなのが単位行列なの?なんか変なのー」

かばん「連立方程式に戻して考えてみよう。



この式は、
a=1×a+0×b+0×c+0×x+0×y+0×z
b=0×a+1×b+0×c+0×x+0×y+0×z
c=0×a+0×b+1×c+0×x+0×y+0×z
x=0×a+0×b+0×c+1×x+0×y+0×z
y=0×a+0×b+0×c+0×x+1×y+0×z
z=0×a+0×b+0×c+0×x+0×y+1×z

こういうことになる。」

サーバル「当たり前の式だね!」

かばん「つまりA=IA」

サーバル「あ、ほんとだ。Aが行列じゃなくて数字だと考えると、確かに当たり前だし、単位行列Iにaをかけると、Aそのものになってるね。」

かばん「同様に、Pに右側からinv(P)をかけても、Pの左からinv(P)をかけても、単位行列Iになるものを、逆行列って呼ぶんだよ。」

サーバル「そっか。逆数みたいなものだね!」

かばん「そゆこと!」

サーバル「でも逆行列の中の人はどうやって計算して出せばいいの?」

かばん「inv(P)=Q/|P|って方法が知られてるよ」

サーバル「|P|ってなに?」

かばん「|P|はdet(P)とも書くんだけど、Pの行列式。行列式は計算していくと、ただの数になるから、割り算ができるんだよ。あとで詳しく教えるね。」

サーバル「じゃあQは?」

かばん「QってのはPの余因子だね。たとえばPの中身がこうなってたら

 



こんな風に表現できるよ。」

サーバル「じゃあここのtは?」

かばん「転置行列、つまり、行と列を入れ替える意味の記号だね。」

サーバル「余因子を求めるときにも、行列式を使うの?」

かばん「うん。元の行列Pより1行1列だけ小さくした行列の行列式を計算するんだよ」

サーバル「大変なんだねー」

かばん「そうだね^^;ちょっと大変だね。でもこのひと手間のおかげで、行列のべき乗がすっごくおいしく料理されるんだ。さっきのA=P×J×inv(P)って式に、Jってのが出てきたでしょ?」

サーバル「うんうん。」

かばん「Jっていうのは対角化された行列のことなんだ。



たとえばこんな感じ。」

サーバル「ちょっとだけ単位行列に似てるね!」

かばん「当然です。対角化された行列は、普通の行列が単位行列化したものと言われているのですよ

サーバル「マーゲイ!?

かばん「ちょっとやってみたかったんだ^^;まあ今のは冗談として」

サーバル「え?」

かばん「対角化された行列っていうのは単位行列みたいに、すごくべき乗がしやすい構造になっているんだ。」

サーバル「うみゃみゃみゃみゃみゃー!

行列のべき乗(右クリックで等身大)

ほんとだ!簡単になってる!すごーい!」


かばん「こんなことのために、JとPの中身をわざわざ考えるんだよ。これが理性開放の技の1つ。「固有値(J)」と「固有ベクトル(P)」だね。これを使って、Aのべき乗を考えてみよう」

サーバル「A^n=PJinv(P)×PJinv(P)×・・・×PJinv(P)×PJinv(P)。あ!隣り合うPとinv(P)が単位行列になって消えるよ!」

かばん「よく気が付いたね!結局、A^n=P×J^n×inv(P)ってなって、Jのべき乗は楽にできるから、行列Aのn乗A^nも楽に計算できることになるんだ。」

拍手[0回]





 

じかんねえええええ!




追記

項は全部で15本。奇数だから、15個の引き算を全部入れ替えると、全体の符号は反転する


こういう、解けることと解き方をあらかじめ知っておいてから
完全に予定調和なパズルゲームに挑む姿勢が僕は好きだ。
人はそれを「ゲーム性がない」という・・・なんでや。こんなに楽しいのに

パズルゲーム大会で何度もリハーサルしたら楽しいだろうが・・・!
マリオだってシーケンサー痛快でたのしかろう!?(よくしらんけども)

拍手[0回]


行列のべき乗 固有ベクトル 固有値 対角化 複素行列 非負行列 ペロン・フロベニウスの定理 グラフ理論 有向 重みなし 隣接行列 それが大事 ノルム 規格化 エルミート共役 ユニタリ行列 複素共役 転置


これの


この行列の、固有値と固有ベクトルを計算してみましょう。

scilabでも求まるんですが、勝手に固有値の順番を 指定されてしまうので
手計算でやってみました。


そんな計算を海外旅行のホテルや機内でも淡々と行う、相変わらずのきのこですこんにちは


det(A-λE)=0 (Eは単位行列)
から、固有値λを求めてみます。


このように、
1行目-(2行目×λのべき乗)
の繰り返しで、どんどん行列の次数が小さくなってくれますので
最終的には

λ^6=1

という式になって

 λ=λ1^n
という、1の原始複素6乗根のべき乗で、すべての固有値を表すことができます。


では、固有ベクトルはどうなるかといいますと

nを0から5までいっぺんに計算させますと、以下のような関係を導出することができます。


また、6乗したら1に戻る性質を利用して、見やすくした上
すべての固有値の絶対値が1であることを利用して、√6で割り算することで、規格化も可能です。


ということは、対角化のためのブラ・ケットの片方はこのようになります。

規格化しているのでこのPはユニタリ行列です。
expを何度も描くのが煩わしいので、肩の部分だけ取って、便宜的に下部分のように表現することにします。


Pがブラだとすると、ケットはブラのエルミート共役
つまり
複素共役を取って転置したものなので、以下のようになります。

6乗して元に戻るので、冗長な書き方かもしれませんが、
規則性を表すためにあえてこのように書いています。

複素共役を取ったものが対称行列なので、その転置は同じものとなり
エルミート共役は以下で表されます。
 



ケットがこのようにシンプルに表せたので、ついでにブラのほうもシンプルに表してみましょう。


このようになりました。
まあつまり、expの指数部分においての複素共役は、マイナスを取るだけでいいのです。



まとめると、対角化のためのブラとケットと、対角化された行列は、このように表されます。

具体例は以下のようになりますが


実は、1の複素6乗根はいずれも、1の複素原始3乗根
ω={-1+j√(3)}/2
とそのべき乗、あるいは複素共役か符号を反転したものだけで構成できることがわかります。

また、対角化された行列のべき乗は、以下のように簡潔に表すことができます。


これらを使って、元の行列Aのべき乗を表してみましょう。

A^L=P×(λ^L)×P†

なので

その過程をExcelにぶち込んで表現してみたのがこちらです。

原寸大はこっち
ソースファイルはこっちからー!

黄色部分(左上)が対角化のための行列ブラP、青部分(右上)が逆行列=エルミート共役のケットP†、
緑部分(左下)が対角化された行列λのべき乗
それらの複素行列を掛け算したのが赤(右下)の部分で、一番左には計算して得られた行列Aの実整数を出力させてます。

行列の中身が複素数になるため、Excelの行列ライブラリは使えません。
そこで、複合参照と複素数ライブラリを用いて計算してます。


空白セルを指定して、そこでdel押しまくるなり再計算しまくるなりすると、動くと思います。
少なくともPCでは動きますが、タブレットではどうなるかわかりませふんせふん

拍手[0回]

大変遅くなってしまい申し訳ありません

2017/02/17の日記のDL用Excelファイルを改造します。


シート1から10までを選択した状態にして、セルも全選択し、

枠線を外します。

それから、右側の表も消します。


準備が整ったところで、5行目のa2の行を、
A5からO13まで選択し、Q4にコピペではなくカットアンドペーストします。
コピペでは絶対参照にしてない場合に、コピー元がずれるからです。


同様に、Q5からAE12までをAG4にカッペーします。


あとは、aなんとかとついているシリーズの2段目を、一連の表の最後尾まで選択し、
表の右側の、さらに1列空けてカッペーします。

これをa10の行まで繰り返します。

最後尾がFC列になってれば作業はOKだと思います。

それから、A列を選択して、列の挿入を行い、P1セルの42とか書いてる数字を
A4セルにカッペーします。

O1のnと書いてあるラベルも、A3セルにカッペーしておきましょう。

今度は、imdivのあるP列を切り取って
FF列目指して、切り取ったセルの挿入を行いましょう。
これと同様のことを、a2ならAE列をFF列に列カッペーといったように、a10まで繰り返します。

これで、
EV列からFE列まで、imdivの式が揃いました。


ここで、

EV3の数式に注目してください。
コロンの指示で区別している複数のセルは、コピペーではもちろん、カッペーでもうまく移動されません。
これまでの一連の作業は、グラフにしやすいようにと見やすさのほかに、imsumのコロン参照をうまくいかせるためでもあったのです。


縦参照合計ではなく、横参照合計にしてやると、ちゃんと数列の値が出ることがわかります。

×4244727807.02364-2.37723717217999E-07i
○4244727806.99991-2.37723717222155E-07i

それでは、EV3セルの値をFH4セルにカッペーし
FG4セルにA4セルの値を参照して代入しましょう。


これで1つの番号が、1つの行にまとめられました。

あとは、A列を好きなようにいじるだけです。

とりあえずA4セルを0にして、
A5をA4+1にしてみましょう。

それから、FHセルの字の色を黒に戻し太字でもなくして、
imreal関数で数列の実部だけ取りだして、文字型から数値型に変えましょう。


A1にdn、A2に1を代入して

A5=A4+$A$1

としておきます。
そして、結構行コピペしますと

あ、忘れ物がありました。
O、AD、AS、BH、BW、CL、DA、DP、EE、ET列でのnの参照が絶対参照でした。
列だけ絶対参照の複合参照に直します。

O列にある$A$4を、$A5に変更し、これを残りの9列にコピペーします。

それから、D4~FE4の行を下までコピペーします。

そうした上で、デカボナッチ数列を見てみると、以下のようになります。


・・・あーなんか煩雑になってきましたね。
とりあえず出来上がったExcelファイルを置いておきます。
モノボナッチからデカボナッチまでの。


とにかく、n番目1つを1行にまとめたかったんです。
それで、nの刻み幅を整数から実数にしてみたり、nの初期値をゼロ未満にしてみたりすると
たとえばヘキサボナッチだったらこんな風に、

4番目まで整数だけ見るとゼロのまま静まり返ってますが、nを整数から実数に拡張してみると、数列は実は複素数の範囲で結構ダイナミックにうごめいているし
nが負の場合にも拡張できて、負の実数のnだとこんな風になってるよーってのを
ここ何日かずっと、言いたかったんだけど余裕がなかったんです!今もない!


なんつーか、真空の量子ゆらぎやビッグバウンスを彷彿とさせるかなーなんつって!

拍手[0回]

昨晩、寝る前になんとなく放送大学をつけたら「入門線形代数」の最終回をやっていて

ちょうど対角化で締めるところだったらしいです。

2番目の例題で、具体的な数値はいまいち忘れてしまったのですが

実対称行列で、

A=[a,zeros(1,2);zeros(2,1),B]
かつ
B=[b,c;c,b]のような構造をしていて

固有値が1重根だったのは覚えてました。

たとえば

A=[-1,0,0;0,2,3;0,3,2]

それで、ジョルダン標準形にはならなかったんです。

たぶんこれが「ランク落ち」とかいうやつで

1つの固有値から2つの線形独立な固有ベクトルが作れることがわかりました。

その根拠として「次元定理」が用いられていたようで

僕の一番知りたかった情報はここにあったようです。

rankA+dim(kerA)=dim(ImA)+dim(kerA)=n


rank:階数、ランク
dim:次元
ker:カーネル、核
Im:像

だそうです。


scilabにとりあえずrankを計算する関数だけは見つけたので遊んでみたところ
少しだけわかってきました。

あと、ウルフラムαでも少しだけ、3次元アフィン変換について遊んでみてました。

今はだるいので画像を作る気力がないのですが

[a,0,0,0;0,b,0,0;0,0,c,0;x,y,z,1]
かその転置によるアフィン変換の


xとyとzを有限に保ちながらaかbかcを1にすると(排他的ではない)、2次のジョルダン標準形になるみたいで、3次以上のジョルダン細胞にはならないみたいです


また、aを1にしながらxを0にするみたいにすると、ジョルダン標準形からただの対角化に戻るみたいです。

おそらく、これはディラック行列の中にパウリ行列を入れ子にしたのと同様の理屈が通るのではないかと類推しました。


[a,0,0,0;0,b,0,0;0,0,1,0;x,y,0,1]というのは、2次行列が2次行列になって入っているマトリョーシカにできて
[A,O;X,E]
とできるのではないか。
A=[a,0;0,b]、O:2×2ゼロ行列、E:2×2単位行列、X=[0,0;x,y]

また、


[a,0,0,0;0,1,0,0;0,0,1,0;x,0,0,1]

というのも、
[A,O;X,E]
ただし、今度は
A=a(スカラー)、O:1行3列のゼロ行列(横ベクトル)、E:3次の単位行列、X=[0;0;x](縦ベクトル)

と考えられるのではないかと類推してみたのです。

拍手[0回]

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


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

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

========
方法その1~行列の固有値を求める方法~

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)


========
方法その2~多項式を直接作って根を求める方法~

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ファイルをここに置いておきます。
ブログランキング・にほんブログ村へ
にほんブログ村

拍手[0回]

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

をこれ↓

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

拍手[0回]

前回の続きです。

テトラボナッチ数列の一般式を出すために計算していて、|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とも一致していますね。つづく

拍手[0回]



こう書かれるわけですが(あ、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|はこのようになりました。


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

拍手[0回]

続きです



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

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


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


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

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

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


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

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

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

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

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


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

拍手[0回]

昨日の続きです。


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


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

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

だったので、その2乗は

こうなります。

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

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

こうなります。

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

こうなります。

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

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

拍手[0回]

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

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

拍手[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はそんなことをするように設計されてはいないのかもしれません

拍手[0回]

 
 


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


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

拍手[0回]



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