to site top page

Ruby Programing : Inverse Matrix[加筆修正:2008/11/6]

注)Rubyにはちゃんと逆行列はき出す関数あるので,間違っても以下のcodeを実装なんかしないように.

Rubyで逆行列演算のプログラム.ソースコードは下の方に.

原理は単純に,元行列の各成分の余因子を行列式で割り,その転置行列を求める,というだけのもの. 念頭に“次元の制約をつけない,構造は出来るだけ単純に”というのがあったので, 計算速度はあまり速くないと思われる.もう少し綺麗なコードがかけないか,後でまた色々弄る予定.

で,どのくらいの速度で計算できるのか

day1 = Time.now
[main block]
day2 = Time.now
print day2 - day1

以上のコードをプログラムの頭と末尾に追加,Check部分をコメントアウトして計算時間を計ってみた.以下それぞれの次数の計算速度.

DimTime(s)DimTime(s)DimTime(s)
CPU Spec:Intel Celeron(Northwood-128K) 2.00GHz
40.01660.375826.453
50.04772.9699256.672

で,これをプロット.

plot graph

ごくごく単純に,アバウトに推定すれば10次で40分位かかるので,10次以上は計算しなかった. 果たして手計算で40分で10次の逆行列計算できるのか.誰かやってみて結果教えて\(>∀<)/

以下,逆行列演算プログラム.

#inverse.rb(Shift_JIS:左クリックよりも右クリック→保存の方が良いかも)

def det(mt)
    row = mt.size - 1
    col = mt[0].size - 1
    if(row == col)
        val = 0
        case row
            when 0
                val = mt[0][0]
            when 1
                val = mt[0][0] * mt[1][1] - mt[0][1] * mt[1][0]
            else
                sub_mtx = []
                sub_mtx_col = []
                sft = 0
                for l in 0..row
                    for m in 0..(row - 1)
                        if l == m
                            sft = 1
                        end
                        for n in 0..(col - 1)
                            sub_mtx_col << mt[m + sft][n]
                        end
                        sub_mtx << sub_mtx_col
                        sub_mtx_col = []
                    end
                    sft = 0
                    val += (((-1)**(l + row)) * mt[l][row] * det(sub_mtx))
                    sub_mtx = []
                end
        end
    else
        val = 0
    end
    return(val)
end


def cof(mtb,m,n)
    sub_mtx = []
    sub_mtx_col = []
    sub_r = mtb.size - 2
    sub_c = mtb[0].size - 2
    sft_r = 0
    sft_c = 0

    for x in 0..sub_r
        if x == m
            sft_r = 1
        end
        for y in 0..sub_c
            if y == n
                sft_c = 1
            end
            sub_mtx_col << mtb[x + sft_r][y + sft_c]
        end
        sft_c = 0
        sub_mtx << sub_mtx_col
        sub_mtx_col = []
    end
    return(((-1)**(m + n + 2)) * det(sub_mtx))
end


def transpose(mtt)
    trs_mtx = []
    trs_mtx_col = []
    for a in 0..(mtt.size - 1)
        for b in 0..(mtt[0].size - 1)
            trs_mtx_col << mtt[b][a]
        end
        trs_mtx << trs_mtx_col
        trs_mtx_col = []
    end
    return(trs_mtx)
end

#※行列サンプル(Encoding : Shift_JIS)
mtx = <適当な行列>

inv_mtx = []
inv_mtx_buf = []
inv_mtx_col = []

inv_r = mtx.size - 1
inv_c = mtx[0].size - 1

for m in 0..inv_r
    for n in 0..inv_c
        inv_mtx_col << cof(mtx,m,n)
    end
    inv_mtx_buf << inv_mtx_col
    inv_mtx_col = []
end

inv_mtx = transpose(inv_mtx_buf)

######################## display part ##########################################

print("matrix:\\n")
for m in 0..(mtx.size - 1)
    for n in 0..(mtx[0].size - 1)
        print(mtx[m][n].to_s(10).rjust(8))
    end
    print "\\n"
end

print("\\ndet(matrix) : ")
print det(mtx)

print("\\n\\ninverse matrix(Elements isn't divided by determinant.) : \\n")
for m in 0..(inv_mtx.size - 1)
    for n in 0..(inv_mtx[0].size - 1)
        print(inv_mtx[m][n].to_s(10).rjust(8))
    end
    print "\\n"
end

# check
print("\\nCheck\\n")
for m in 0..(mtx.size - 1)
    for n in 0..(inv_mtx[0].size - 1)
        calc_buf = 0
        for k in 0..(mtx[0].size - 1)
            calc_buf += mtx[m][k] * inv_mtx[k][n]
        end
        print (calc_buf / det(mtx)).to_s(10).rjust(8)
    end
    print "\\n"
end

can't load my result

前後の記事

最近の記事(5件分)

する事

そのうち記事にするかもリスト

欲しい本