ONLY DO WHAT ONLY YOU CAN DO

こけたら立ちなはれ 立ったら歩きなはれ

Go 言語で固有値を求める (LR分解)

f:id:fornext1119:20150711180736p:plain
対称行列 f:id:fornext1119:20150805175030p:plain を 左三角行列 f:id:fornext1119:20150806120047p:plain と 右三角行列 f:id:fornext1119:20150806115724p:plain
積に分解し,
f:id:fornext1119:20150806120116p:plain
逆順に掛けたものを新たな f:id:fornext1119:20150805175030p:plain として
f:id:fornext1119:20150806120132p:plain
処理を反復すると, 元の行列の固有値が対角成分に並んだ右三角行列に収束する.

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
参考文献