[基]連立一次方程式（第六章）

連立一次方程式を行列で解くプログラム


線形代数でやったのをプログラムにしようか！っていうだけ。
だから解説はいらないと思うので関数だけさらしておきます。
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