NB. 7.固有値、固有ベクトル計算のPower法とDeflation法−Jプログラム NB. QR method =============================== qr =: 3 : 0 A1 =. (128!:0) y. ((>@{:) (+/ . *) (>@{.)) A1 ) diag =: i.@# {"_1 ] eigen_value =: 3 : 0 40 eigen_value y. : diag qr^:x. y. ) NB. Roger Hui / M. Shimura 2010/8/3 rheval =: +/ . * & >/@ |. &(128!:0) rh =: [: (<0 1)&|: rheval ^:_ rhn =: [: (<0 1)&|: rheval ^:(50) NB. 50 repeat safely ! by TN NB. Eigen Calc. / Power Iteration and Deflation ============ NB. 戸田、小野「入門数値計算」p.170-2 DT0 =: 2 2$1 8 2 1 NB. p.170 NB. Usage: NB. DT0 pow^:(4) 1 1;1 => x; r; el DT1 =: 3 3$6 1 1 1 6 1 3 3 6 wr =: 1!:2&2 amax =: >./ @: | ip =: +/ . * cleanz =: * | >: 1e_4"_ pow0 =: 3 : 0 : A =. x. x0 =. > {. y. el =. amax x0 y0 =. x0 % el x1 =. A ip y0 r =. (y0 ip x1) % (y0 ip y0) wr 'x: ', ": x0 wr 'l: ', ": el wr 'y: ', ": y0 wr 'x: ', ": x1 wr 'r: ', ": r wr '-------------------' x1; r; el ) pow1 =: 3 : 0 : A =. x. x0 =. > {. y. el =. amax x0 y0 =. x0 % el x1 =. A ip y0 r =. (y0 ip x1) % (y0 ip y0) xx =. x1 % ({. x1) xx; r; el ) pow =: 3 : 0 A =. y. 9!:11 (5) NB. floating 5 digits cleanz L:0 A pow1^:(40) (({. $A)#1);1 ) max =: ({.@\:) { ] defpow =: 3 : 0 0 defpow y. : j =. x. wr 'j: ', ": j A =. y. N =. {. $ A wr 'A' wr A 'v1 lam1 one' =: pow A wr 'v1' wr V1 =: (N, 1)$v1 wr 'lambda1 = ', ": lam1 wr 'X' lam1 * (0{v1) AT =: (N, 1) $ {. A wr X =. AT % (lam1 * (0{v1)) wr 'lam1*v1*Xt' wr C =. lam1 * ((N, 1)$ , v1) ip ((1, N)$ , X) wr 'B' wr cleanz B =. A - C wr 'B1' wr BB =. }. }."(1) B wr pow BB 'uu2 lam2 one' =. pow BB wr 'u2x' u2x =. uu2 wr u2x =. cleanz u2x % (max u2x) wr 'lambda2 = ', ": lam2 wr '---' wr u2 =: (N, 1)$0, u2x wr 'v2 ----' wr 'v2a' wr v2a =. (lam2 - lam1) * u2 wr 'v2b' wr ((1, N)$ , X) wr '*** ' wr lam1 * ((1, N)$ , X) wr ((N, 1)$,u2) wr v2b1 =. lam1 * ((1, N)$ , X) ip ((N, 1)$,u2) wr '---' wr v2b =. ({. , v2b1) * (N, 1)$v1 wr 'v2(result:)' wr v2 =: cleanz v2a + v2b '** end **' return. ) dpow =: 3 : 0 0 dpow y. : j =. x. A =. y. N =. {. $ A 'v1 lam1 one' =: pow A V1 =: (N, 1)$v1 wr 'lambda1 = ', ": lam1 wr v1;lam1 lam1 * (0{v1) AT =: (N, 1) $ {. A X =. AT % (lam1 * (0{v1)) C =. lam1 * ((N, 1)$ , v1) ip ((1, N)$ , X) B =. A - C BB =. }. }."(1) B 'uu2 lam2 one' =. pow BB u2x =. uu2 u2x =. cleanz u2x % (max u2x) wr 'lambda2 = ', ": lam2 u2 =: (N, 1)$0, u2x v2a =. (lam2 - lam1) * u2 v2b1 =. lam1 * ((1, N)$ , X) ip ((N, 1)$,u2) v2b =. ({. , v2b1) * (N, 1)$v1 v2 =: , cleanz v2a + v2b v2;lam2 )