と行列を定義してその固有値を計算させたとしても，実行結果は
実行結果の例 [ 1 1/2 1/3 1/4 1/5]
[1/2 1/3 1/4 1/5 1/6]
[1/3 1/4 1/5 1/6 1/7]
[1/4 1/5 1/6 1/7 1/8]
[1/5 1/6 1/7 1/8 1/9]
Traceback (click to the left of this block for traceback) ...
ArithmeticError: could not determine eigenvalues exactly using symbolic matrices; try using a different type of matrix via self.change_ring(), if possible
となって，固有値を具体的に決定できなかったことが分かります*17。この行列の固有値は５次方程式であ り，一般解を具体的に表すことはできません。しかし，数値的には固有値を求めることはできます*18。
そこで，例えば上で作った2×2行列の係数体をCDFにして計算します：
1 mat = m a t r i x ( CDF , mat ) # b a s e _ r i n gをC D Fに 変 更 2 p r i n t mat . e i g e n v a l u e s ()
実行結果の例 [-0.372281323269, 5.37228132327]
固有値の数値のリストが表示されました。はじめから数値結果だけ得られればよい場合にはSRを使わずに RDFやCDFを使いましょう。そのほうがSRで固有値を計算するよりもずっと高速に計算できます。
固有ベクトルは行列に対するメソッドeigenvectors_right()で求めます。
1 m a t 2 = m a t r i x ( SR , 2 ,2 ,[1 ,2 ,3 ,4]) # b a s e _ r i n gはS Rに 設 定 2 e i g v = m a t 2 . e i g e n v e c t o r s _ r i g h t () # 固 有 ベ ク ト ル の 組 を 定 義 3 p r i n t ei g v
実行結果の例 [(-1/2*sqrt(33) + 5/2, [(1, -1/4*sqrt(33) + 3/4)], 1), (1/2*sqrt(33) +
5/2, [(1, 1/4*sqrt(33) + 3/4)], 1)]
これは，次を意味しています：
[ (第1固有値,[対応する固有ベクトル],多重度), (第2固有値,[対応する固有ベクトル],多重度) ] なので，少し面倒ですが，固有ベクトルを取り出すには，上のプログラムに続けて次のようにします。
1 p r i n t ei g v [ 0 ] [ 1 ] [ 0 ] # 一 つ 目 の 固 有 ベ ク ト ル 2 p r i n t ei g v [ 1 ] [ 1 ] [ 0 ] # 二 つ 目 の 固 有 ベ ク ト ル
実行結果の例 (1, -1/4*sqrt(33) + 3/4)
(1, 1/4*sqrt(33) + 3/4)
49.7 対角化とその応用
49.7.1 行列の対角化
行列の対角化は，行列の冪や指数関数を計算するために必要な手順の一つです。行列Aを対角化するとは，
ある正則行列Pを見つけて
P−1AP = diag[λ1,· · ·, λn] : 対角行列
を計算することです。通常，対角化行列Pを作るためには，Aの固有ベクトルを計算すればよいのですが，
Sageでは，eigenmatrix_right()で対角化を簡単に行うことができます。
*175×5行列であるところを4×4行列にすると厳密な固有値が得られます。
*18もちろん結果は誤差を含む近似値になりますが。
行列 (1 4
2 3 )
を対角化してみましょう。
1 A = m a t r i x ( SR , [[1 ,4] ,[2 ,3]] ) # 行 列 の 定 義
2 A . e i g e n m a t r i x _ r i g h t () # 固 有 値 と 対 角 化 行 列 を 求 め る
実行結果の例 (
[ 5 0] [ 1 1]
[ 0 -1], [ 1 -1/2]
)
これはAの固有値が5,−1で
P−1AP = (5 0
0 −1 )
, P =
(1 1 1 −1/2
)
となることを意味しています。実際 (1 1
1 −1/2 )−1(
1 4 2 3
) (1 1 1 −1/2
)
=
(1/3 2/3 2/3 −2/3
) (1 4 2 3
) (1 1 1 −1/2
)
= (5 0
0 −1 )
となり，ちゃんと対角化されているのが分かります。
49.7.2 行列の指数関数と線形微分方程式
行列Aに対して，その指数関数exp(tA)を exp(tA) =
∑∞ n=0
tn
n!An, t∈C (29)
によって定義します。Sageで行列の指数関数を計算することができますが，以下ではこれを応用して微分方 程式を解きます。
最も簡単な例に対して，具体的な計算を行ってみましょう。ばねにつながれた質量mからなる1次元の力 学系を考えます。ばねのつり合いの位置を原点としましょう。このとき運動方程式から
md2x(t)
dt2 =−kx(t)
が成り立ちます。ここでkはバネ定数です。粒子は時刻0に位置X にいて，速度V であるとします。こ こでの目標は時刻t での位置x(t)を求めることです。この微分方程式は定数変化法などにより解くこと ができますが，いまは行列を用いて解いてみます。まず方程式を線形化します。√
k/m = wと置きます。
v(t) =dx(t)/dtとすると上の微分方程式は d
dt (x(t)
v(t) )
=
( v(t)
−w2x(t) )
=
( 0 1
−w2 0 ) (x(t)
v(t) )
と同値です。つまり
A=
( 0 1
−w2 0 )
と置けば，
d dt
(x(t) v(t) )
=A (x(t)
v(t) )
です，これはベクトル(x(t), v(t))tを一回微分するごとに行列Aが前に出てくることを意味しているので，
(x(0), v(0)) = (X, V)を思い出すと
(x(t) v(t) )
= exp(tA) (X
V )
49.7 対角化とその応用 49 行列計算
となります。微分方程式を解く問題は行列の指数関数exp(tA)を求める問題に置き換わりました。Sageで行 列T の指数関数exp(T)を求める命令がT.exp()です。具体的にSageに行列を計算させて微分方程式を解 いてみましょう：
1 var (’ t w ’) # t wを 変 数 と す る
2 A = m a t r i x ([[0 ,1] ,[ - w ^2 ,0]]) # 行 列 Aを 定 義 す る
3 B = t * A # Bをt Aと す る
4 e x p o = B . exp () # e x p oをe x p ( B )と す る
これで，B= exp(tA)が求まりました。つぎに
1 var (" X V ") # X , Vを 変 数 と す る
2 y = v e c t o r ([ X , V ]) # y を ベ ク ト ル( X , V )と す る
3 sol = e x p o * y # sol をexp ( tA ) を( X , V )に 作 用 さ せ た も の と す る
とします。最後得られた量sol[0]は(x(t), v(t))なのでその最初の成分（第０成分）が時刻tの位置を与え ます。つまり微分方程式の解x(t)はsol[0]です。得られる関数は実なので実部を取って（これは何の変化 ももたらしません）full_simplifyを用いて整理します。上のプログラムに続けて
1 p r i n t r e a l _ p a r t ( sol [ 0 ] ) . f u l l _ s i m p l i f y () と入力すると最終的な答え
1 ( X * w * cos ( t * w ) + V * sin ( t * w ))/ w
が得られます。これが微分方程式の解x(t)です。通常の数学の記法で書けば，微分方程式の解は x(t) =Xcos(wt) + V
W sin(wt), x(0) =X, x′(0) =V
となるという事を意味しています。
49.7.3 行列の対角化とその応用：フィボナッチ数列の一般項
フィボナッチ数列a1= 1, a2= 1,
an+2=an+1+an, n= 1,2,3,· · · (30) を考えます。まず，上の関係式は
(an+2 an+1 )
= (1 1
1 0
) (an+1 an
)
となる事に注意します。右辺の行列をAとします。すると (an+1
an
)
=A ( an
an−1
)
=A2 (an−1
an−2
)
=· · ·=An−1 (a2
a1
)
=An−1 (1
1 )
(31) なのでAnが具体的に計算できれば，数列の一般項anが求まる事になります。Anを求めるために，Aを対 角化してみましょう。
1 A = m a t r i x ( SR , [ [1 ,1] ,[1 ,0] ]) 2 sol = A . e i g e n m a t r i x _ r i g h t () 3 P = sol [1]
4 D = P ^( -1)* A * P
5 p r i n t ’ P = ’, P , ’ ’ 6 p r i n t ’ D = P ^( -1) AP = ’ 7 p r i n t D . e x p a n d ()
実行結果の例 P =
[ 1 1]
[-1/2*sqrt(5) - 1/2 1/2*sqrt(5) - 1/2]
D = P^(-1)AP =
[-1/2*sqrt(5) + 1/2 0]
[ 0 1/2*sqrt(5) + 1/2]
ここでDを単純化するためにexpand()を行いました。この事からAを対角化する行列は P =
( 1 1
−(√
5 + 1)/2 (√
5−1)/2 )
であり，P−1AP = diag[(−√
5 + 1)/2,(√
5 + 1)/2]となることが分かります。この式の両辺をn乗すると
(
(−1+2√5)n 0 0 (√5+12 )n
)
=P−1AnP となります。したがってAnは
An=P (
(−1+2√5)n 0 0 (√5+12 )n
)
P−1 (32)
と計算する事ができます。これをつかってAn−1を計算させてみましょう：
1 var (’ n ’)
2 Dn = m a t r i x ([ [ D [ 0 ] [ 0 ] ^ ( n -1) , 0] , [0 , D [ 1 ] [ 1 ] ^ ( n - 1 ) ] ] ) # Dのn -1乗 を 作 る
3 An = P * Dn * P ^( -1) # こ れ がAのn -1乗
4 p r i n t e x p a n d ( An ) # A nを 少 し 単 純 化 し て プ リ ン ト
実行結果の例 [ 1/5*sqrt(5)*(1/2*sqrt(5) + 1/2)^n/(sqrt(5) + 1) +
1/5*sqrt(5)*(-1/2*sqrt(5) + 1/2)^n/(sqrt(5) - 1) + (1/2*sqrt(5) + 1/2)^n/(sqrt(5) + 1) - (-1/2*sqrt(5) + 1/2)^n/(sqrt(5) - 1) 2/5*sqrt(5)*(1/2*sqrt(5) + 1/2)^n/(sqrt(5) + 1) +
2/5*sqrt(5)*(-1/2*sqrt(5) + 1/2)^n/(sqrt(5) - 1)]
[
2/5*sqrt(5)*(1/2*sqrt(5) + 1/2)^n/(sqrt(5) + 1) + 2/5*sqrt(5)*(-1/2*sqrt(5) + 1/2)^n/(sqrt(5) - 1) 1/5*sqrt(5)*(1/2*sqrt(5) + 1/2)^n/(sqrt(5) + 1)
-1/5*sqrt(5)*(-1/2*sqrt(5) + 1/2)^n/(sqrt(5) - 1) + (1/2*sqrt(5) + 1/2)^n/(sqrt(5) + 1) - (-1/2*sqrt(5) + 1/2)^n/(sqrt(5) - 1)]
Dn−1, An−1をDn, Anと定義しました。複雑ですが正しく計算できています。上の行列にベクトル
(1 1 )
(33) を掛けて得られるベクトルの第2成分が一般項anとなります。上のプログラムに続けて次を実行します。
1 a s s u m e ( n ,’ i n t e g e r ’) # nを 整 数 と 仮 定 す る
2 v = v e c t o r ([1 ,1]) # ベ ク ト ル( 1 , 1 )をvと 定 義
3 aa = An * v # A ^ nをvに 掛 け る と
4 aa [ 1 ] . s i m p l i f y _ f u l l () # そ の 第2成 分 が 一 般 項 と な る
実行結果の例 -1/5*((sqrt(5) - 1)^n*(-1)^n - (sqrt(5) + 1)^n)*sqrt(5)/2^n
最後の結果が，フィボナッチ数列の一般項anとなります。これを整理すれば，
an =− (
512n(−1)n(√
5−1)n
−(√
5 + 5)n)
5−12n−12
2n = 1
√5
{(1 +√ 5 2
)n
−(1−√ 5 2
)n} となります。
実際，最後の出力の結果をn=1,2,...に対して表示させると