周波数特性　続き　大熊

 投稿日：2008年11月27日(木)15時55分50秒
  山中　様　　大熊です。
確認のつもりが間違って投稿を押してしまいました。



（３）突然ですが、同じ行数と列数のMATには、
       固有値なるものが、あるそうですが、
　　　 この命令語は単独で在るのでしょうか。
　　　または、此れを解いて表示するプログラム
　　　があれば御教えください。



　敬具 		 

Re: 周波数特性　続き　大熊

 投稿者：山中和義  投稿日：2008年11月28日(金)20時05分5秒
  > No.124[元記事へ]

大熊　正さんへのお返事です。

> （３）突然ですが、同じ行数と列数のMATには、固有値なるものが、あるそうですが、
>　　　　略
> 　　　または、此れを解いて表示するプログラムがあれば御教えください。

!直接法による行列の固有値を求める

OPTION ARITHMETIC COMPLEX

LET i=SQR(-1) !虚数単位

LET cEps=1e-8 !誤差　※単精度


LET N=3 !N次正方行列

FUNCTION tr(A(,)) !行列Aのトレース
   LET t=0
   FOR j=1 TO N
      LET t=t+A(j,j)
   NEXT j
   LET tr=t
END FUNCTION


DIM c(N) !多項式 X^N+c(1)*X^(N-1)+c(2)*X^(N-2)+ … +c(N-1)*X+c(N) の係数
SUB DKA_00(A(),Xr()) !DKA法（Durand Kerner Aberth）
   LET r=1 !初期値を仮定する
   FOR j=2 TO N
      LET rn=ABS(A(j))^(1/j)
      if r<rn then LET r=rn
   NEXT j
   FOR j=1 TO N !半径rの円に等間隔に配置する
      LET Xr(j)=-A(1)/N+r*EXP( 2*PI*i/N *(j-3/4) ) !アーバスの初期値
   NEXT j

   FOR m=0 TO 100 !反復　※調整要
      LET mfx=0
      LET maj=0
      FOR j=1 TO N
         LET Xk=1
         LET fx=1
         FOR w=1 TO N
            LET fx=fx*Xr(j)+A(w)
            IF w<>j THEN LET Xk=Xk*(Xr(j)-Xr(w))
         NEXT w
         LET Xr(j)=Xr(j)-fx/Xk
         IF mfx<ABS(fx) THEN LET mfx=ABS(fx)
         IF maj<ABS(fx/Xk) THEN LET maj=ABS(fx/Xk)
      NEXT j
      IF mfx<cEps AND maj<cEps THEN EXIT FOR !収束したら
   NEXT m
END SUB
!-------------------- ここまでがサブルーチン


DIM A(N,N) !行列A
!DATA 1,0,0 !λ=1（３重根）
!DATA 0,1,1
!DATA 0,0,1

!DATA 0,1,1 !λ=2,-1（重根）
!DATA 1,0,1
!DATA 1,1,0

!DATA 3,0,0 !λ=3,±i
!DATA 0,2,-5
!DATA 0,1,-2

DATA 2,1,-1 !λ=3,2,1
DATA 0,3,0
DATA 0,2,1

MAT READ A

MAT PRINT A;


!n次正方行列Aの固有多項式 det(tE-A)=t^n+c1*t^(n-1)+ … + cn を求める。
DIM X(N,N),cE(N,N)
MAT X=IDN !frame法
FOR k=1 TO N
   MAT X=A*X
   LET c(k)=-tr(X)/k
   MAT cE=(c(k))*IDN
   MAT X=X+cE
NEXT k
MAT PRINT c;

!ニュートン法などで解く。解が固有値になる。
DIM lmd(N)
CALL DKA_00(c,lmd)

FOR k=1 TO N
   PRINT "固有値=";lmd(k)
NEXT k



!※N個求まった場合の検算
LET s=1
FOR k=1 TO N
   LET s=s*lmd(k)
NEXT k
PRINT s, DET(A) !固有値の積=行列の行列式　Πλi=|A|

LET s=0
FOR k=1 TO N
   LET s=s+lmd(k)
NEXT k
PRINT s, tr(A) !固有値の和=行列のトレース　Σλi=trA


END
 

戻る