SECONDさんへのお返事です。
> どなたか、固有ベクトルを高速に算出する方法を、ご指導ください。
> CADなどは、速いですが、どんなアルゴリズムが、使われているのでしょうか。
・対称行列
ヤコビ法
アルゴリズムの本に掲載されている。手元にコードなし。
・最大の固有値・固有ベクトルの算出
べき乗法(パワー法)
アルゴリズムの本に掲載されている。
数値計算の専門書には、複素数への拡張がされている。
!べき乗法による行列の固有値と固有ベクトルを求める
!※固有値が0、重複する場合は適用できない。
!※実数の固有値のみ。虚数を含む解は得られない。
!Ax=λIx、λ:固有値、x:固有ベクトル
LET N=3 !N次正方行列
DATA 2,1,-1 !λ=3,2,1
DATA 0,3,0
DATA 0,2,1
!DATA 0,1,1 !λ=2,1,0
!DATA -4,4,2
!DATA 4,-3,-1
!DATA 1,0,0 !λ=1(3重根)
!DATA 0,1,1
!DATA 0,0,1
DIM A(N,N) !行列A
MAT READ A
MAT PRINT A;
LET cEps=1e-6 !誤差 ※調整要、単精度
DIM u(N) !固有ベクトル
DIM AA(N,N) !作業用
MAT AA=A
FOR s=1 TO N !s番目
CALL EigenPower(N,AA, lambda,u)
PRINT "固有値=";lambda
PRINT "固有ベクトル"
MAT PRINT u;
FOR i=1 TO N !残差行列を求めて、次へ
FOR j=1 TO N
LET AA(i,j)=AA(i,j)-lambda*u(i)*u(j)
NEXT j
NEXT i
NEXT s
DEF norm(v())=SQR(DOT(v,v)) !ノルム
SUB EigenPower(N,A(,), lambda,u()) !固有値(絶対値最大)、固有ベクトルを求める
DIM u0(100),u2(100) !※最大100次
MAT u=CON !初期値 ※ノルムが1
MAT u=(1/norm(u))*u
LET cMax=100
FOR i=1 TO cMax !最大回数まで繰り返す
MAT u0=u !直前のu
MAT u=A*u0
WHEN EXCEPTION IN
MAT u=(1/norm(u))*u !正規化する
USE
PRINT "0ベクトルになりました。"
STOP
END WHEN
MAT u2=u-u0 !収束したか確認する
IF norm(u2)<cEps THEN EXIT FOR
MAT u2=u+u0
IF norm(u2)<cEps THEN EXIT FOR
NEXT i
IF i>cMax THEN
PRINT "収束しません。"
STOP
END IF
MAT u2=A*u
LET lambda=DOT(u2,u)/DOT(u,u) !固有値
END SUB
END