なんとなく手探りで作成。
function fn_gaussianElimination(dm)
# ax+by+cz = d の連立方程式を解く 項数はいくつでも可能
# ただし0チェックや同じ値などによるチェックをしていないので解ナシの場合いろいろなエラーがでる可能性あり
# (a1)x+(b1)y+(c1)z=d1・・・のような連立1次方程式を
# dm=[a1 b1 c1 d1 ; a2 b2 c2 d2 ; a3 b3 c3 d3]のように係数として入力する
#戻り値 [x y z]の値
dmSize_tup = size(dm)
#次元サイズチェック
if dmSize_tup[1] != dmSize_tup[2]-1
throw(DomainError(dm, "サイズが正しくありません"))
return
end
#---前進消去--- エラーチェックしていないがよいのか??
for rowParent_num in 1:dmSize_tup[1]-1
# rowParent_num ベースとなる行(分母となる方)番号
dmChild =[]
for rowChild_num in rowParent_num+1:dmSize_tup[1]
# rowChild_num 分子となる 行番号
dm_x = dm[rowChild_num,:]-dm[rowParent_num,:] .* dm[rowChild_num,rowParent_num]/dm[rowParent_num,rowParent_num]
if size(dmChild)[1] == 0
dmChild = dm_x'
else
dmChild = vcat(dmChild,dm_x') #行列と行ベクトルの結合
end
end
#dmを更新する
dm = vcat(dm[1:rowParent_num,:],dmChild)
end
# この時点で消去後の行列がdmに格納
#---後退代入---
answerValue_ary = zeros(1,dmSize_tup[2]-1) #ここに変数の正解値を入れる
for checkRow_num in dmSize_tup[1]:-1:1
checkColumn_num = checkRow_num #現在解いている変数の位置
totalVariablesValues_num = 0
for checkColumnAnswer_num in size(answerValue_ary)[2]:-1:checkColumn_num+1
# checkColumnAnswer_num すでに解き終わった変数の値の位置
# 変数合計値 += 係数 * 変数値
totalVariablesValues_num += dm[checkRow_num,checkColumnAnswer_num]*answerValue_ary[checkColumnAnswer_num]
end
#今解いている変数の値 = (数値項+変数合計値)÷変数の係数
variableValue_num = (dm[checkRow_num,dmSize_tup[2]]-totalVariablesValues_num)/dm[checkRow_num,checkColumn_num]
answerValue_ary[checkColumn_num] = variableValue_num
end
return answerValue_ary
end
確認
dm1 = [2 1 1 15; 4 6 3 41; 8 8 9 83]
println(fn_gaussianElimination(dm1)) # [5.0 2.0 3.0]
dm2=[2 3 7; 4 -5 3]
println(fn_gaussianElimination(dm2)) #[2.0 1.0]
dm3=[2 3 1 7; 4 -5 -1 3; 6 -7 1 5]
println(fn_gaussianElimination(dm3)) #[2.0 1.0]
確認 解無しなど
dm4=[2 1 1 15; 4 6 3 41; 2 9 3 37] #解なしの場合
println(fn_gaussianElimination(dm4)) #[NaN NaN NaN]
dm5=[2 3 7;4 6 14] #解が無数にある場合
println(fn_gaussianElimination(dm5)) #[NaN NaN]
dm6=[1 3 -2 2;2 5 1 4;3 8 -1 0] #こちらも解無し
println(fn_gaussianElimination(dm6)) #[NaN -Inf -Inf]
これだけで済むならばコンパクトかもですが。。
方程式と正解は以下を参考にしました。