Go 言語で固有値を求める (LR分解)
対称行列 を 左三角行列 と 右三角行列 の
積に分解し,
逆順に掛けたものを新たな として
処理を反復すると, 元の行列の固有値が対角成分に並んだ右三角行列に収束する.
package main import "fmt" import "math" const N = 4 // LR分解で固有値を求める func main() { var a [N][N]float64 = [N][N]float64{{5.0, 4.0, 1.0, 1.0}, {4.0, 5.0, 1.0, 1.0}, {1.0, 1.0, 4.0, 2.0}, {1.0, 1.0, 2.0, 4.0}} var l [N][N]float64 = [N][N]float64{{0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}} var u [N][N]float64 = [N][N]float64{{0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}, {0.0 ,0.0 ,0.0 ,0.0}} for k := 1; k < 200; k++ { // LU分解 decomp(a, &l, &u) // 行列の積 multiply(u, l, &a) // 対角要素を表示 fmt.Printf("%3d\t", k) disp_eigenvalue(a) // 収束判定 var e float64 = 0.0 for i := 1; i < N; i++ { for j := 0; j < i; j++ { e += math.Fabs(a[i][j]) } } if (e < 0.00000000001) { break } } fmt.Println("\neigenvector") disp_eigenvalue(a) } // LU分解 func decomp(a[N][N]float64, l *[N][N]float64, u *[N][N]float64) { for i := 0; i < N; i++ { for j := 0; j < N; j++ { l[i][j] = 0.0 u[i][j] = 0.0 } } l[0][0] = 1.0 for j := 0; j < N; j++ { u[0][j] = a[0][j] } for i := 1; i < N; i++ { u[i][0] = 0.0 l[0][i] = 0.0 l[i][0] = a[i][0] / u[0][0] } for i := 1; i < N; i++ { l[i][i] = 1.0 var t float64 = a[i][i] for k := 0; k <= i; k++ { t -= l[i][k] * u[k][i] } u[i][i] = t for j := i + 1; j < N; j++ { u[j][i] = 0.0 l[i][j] = 0.0 t = a[j][i] for k := 0; k <= i; k++ { t -= l[j][k] * u[k][i] } l[j][i] = t / u[i][i] t = a[i][j] for k := 0; k <= i; k++ { t -= l[i][k] * u[k][j] } u[i][j] = t } } } // 行列の積 func multiply(a[N][N]float64, b[N][N]float64, c *[N][N]float64) { for i := 0; i < N; i++ { for j := 0; j < N; j++ { var s float64 = 0.0 for k := 0; k < N; k++ { s += a[i][k] * b[k][j] } c[i][j] = s } } } // 対角要素を表示 func disp_eigenvalue(a[N][N]float64) { for i := 0; i < N; i++ { fmt.Printf("%14.10f\t", a[i][i]) } fmt.Println("") }
z:\go>8g GO1103.go z:\go>8l -o GO1103.exe GO1103.8 z:\go>GO1103 1 8.6000000000 1.8444444444 4.6143790850 2.9411764706 2 9.6046511628 1.1012311902 4.8997514499 2.3943661972 3 9.8377723971 1.0107124514 4.9934604403 2.1580547112 4 9.9219788334 1.0010980897 5.0142271522 2.0626959248 5 9.9611291643 1.0001111458 5.0138771204 2.0248825695 6 9.9805335651 1.0000111820 5.0095550207 2.0099002322 7 9.9902522897 1.0000011216 5.0057991328 2.0039474559 8 9.9951218389 1.0000001123 5.0033019159 2.0015761328 9 9.9975597740 1.0000000112 5.0018103821 2.0006298327 10 9.9987795937 1.0000000011 5.0009686042 2.0002518010 (中略) 140 10.0000000000 5.0000000000 1.9999999973 1.0000000027 141 10.0000000000 5.0000000000 1.9999999986 1.0000000014 142 10.0000000000 5.0000000000 1.9999999993 1.0000000007 143 10.0000000000 5.0000000000 1.9999999997 1.0000000003 144 10.0000000000 5.0000000000 1.9999999998 1.0000000002 145 10.0000000000 5.0000000000 1.9999999999 1.0000000001 146 10.0000000000 5.0000000000 2.0000000000 1.0000000000 147 10.0000000000 5.0000000000 2.0000000000 1.0000000000 148 10.0000000000 5.0000000000 2.0000000000 1.0000000000 149 10.0000000000 5.0000000000 2.0000000000 1.0000000000 eigenvector 10.0000000000 5.0000000000 2.0000000000 1.0000000000