連立一次方程式を行列で解くプログラム
線形代数でやったのをプログラムにしようか!っていうだけ。
だから解説はいらないと思うので関数だけさらしておきます。
jgp(a):行列aを解く本関数
maxrow(a,k):k列目の係数の絶対値が最大となる行をk行目以降から探す補助関数
abs(x):xの絶対値を出す補助関数
swap(a,i,j):配列aのi項目とj項目を入れ替える補助関数
ただ、誤差が生じる※のを防ぐために、k行目使ってk列を消去する前に、|ai,k|が最も大きい行i=pを選び、k行目とp行目を入れ替える必要があるのに注意。
※ak,kが0だったら計算できない。また0に近いと、計算途中に非常に大きな値が大量に発生し、情報落ち誤差が生まれてしまう。
def gjp(a)
row = a.length()
col = a[0].length()
for k in 0..(col-2)
max = maxrow(a,k)
swap(a,k,max)
akk = a[k][k]
for i in 0..(col-1)
a[k][i] = a[k][i]*1.0 / akk
end
for i in 0..(row-1)
if i != k
aik = a[i][k]
for j in k..(col-1)
a[i][j] = a[i][j] - aik*a[k][j]
end
end
end
end
a
end
def maxrow(a,k)
s = k
t = abs(a[k][k])
if k != a.length()-1
for y in (k+1)..(a.length()-1)
if t < abs(a[y][k])
s = y
t = abs(a[y][k])
end
end
end
s
end
def abs(x)
if x<0
-x
else
x
end
end
def swap(a,i,j)
t = a[i]
a[i] = a[j]
a[j] =t
a
end
0 件のコメント:
コメントを投稿