数値解析(2)演習

Gauss-Jordan法

以下のプログラム(gj.rb)は係数行列aを与えて連立一次方程式を解くプログラムである．

#
# Linear Equation (without pivoting)
#
def gj(x)
  row=x.size 
  col=row+1
  a=Array.new(x.size){|i| x[i].dup}
  for k in 0..(row-1)
    # normalize a[k][i] s.t. a[k][k]=1
    akk = a[k][k]
    for i in 0..(col-1)
      a[k][i]=a[k][i]/akk
    end
    # adjust i-th row s.t. a[i][k]=0,  skipping k-th row
    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
  Array.new(row){|i| a[i][row]}
end

このファイルを gj.rb というファイルに保存して，

load "./gj.rb"
a = [[1.0,1.0,-1.0,2.0],
     [3.0,5.0,-7.0,0.0],
     [2.0,-3.0,1.0,5.0]]
gj(a)

のようにして使える．

演習問題

N変数の連立一次方程式をGauss-Jordan法で解く(係数行列が N*N)プログラムの (時間)計算量のオーダはどうなるのか，O(N)?, O(N log N)?, O(N**2)?, O(N**2 log N)?, O(N**3)?．オーダを答えると共に，どうしてそうなるのか説明をしなさい．

以下のプログラム(gj_pivot.rb )はpivotingを行い誤差が大きくならないように工夫をしたプログラムである(配布した紙のスライドでは一部の式が間違っているので注意が必要)．

#
# Linear Equation (without pivoting)
#
def gj_pivot(x)
  row=x.size 
  col=row+1
  pivot=Array.new(row){|i| i}
  a=Array.new(row){|i| x[i].dup}
  for k in 0..(row-1)
    # find absolute-maximal coefficient 
    # among a[k+1][k]..a[row-1][k]
    max = k
    for i in k+1..row-1
      if a[pivot[i]][k].abs > a[pivot[max]][k].abs
        max = i
      end
    end
    # swap pivot[k] and pivot[max]
    tmp=pivot[max]
    pivot[max]=pivot[k]
    pivot[k]=tmp
    # normalize a[k][i] s.t. a[k][k]=1
    akk = a[pivot[k]][k]
    for i in 0..(col-1)
      a[pivot[k]][i]=a[pivot[k]][i]/akk
    end
    # adjust i-th row s.t. a[i][k]=0,  skipping k-th row
    for i in 0..(row-1)
      if i != k
        aik = a[pivot[i]][k]
        for j in k..(col-1)
          a[pivot[i]][j] = a[pivot[i]][j] - aik * a[pivot[k]][j]
        end
      end
    end
  end
  b=Array.new(row)
  for i in 0..(row-1)
#   b[pivot[i]]=a[i][row] から修正
    b[i]=a[pivot[i]][row]
  end
  b
end

このファイルを gj_pivot.rb というファイルに保存して，

load "./gj_pivot.rb"
a =  [[1.0,  -50.0, -3.0,  -90.0],
     [-85.0,  2.0,  -25.0, -6.0],
     [79.0,  5.0,  30.0, -1.0]]
gj_pivot(a)
gj(a)

のようにして結果を比較してみなさい． gj_pivot.rbのプログラムに誤りがあったので12/3に修正しました．

演習問題

Gauss-Jordan法でPivotingを行わないと誤差が大きくなるような入力で， [[1.0, -50.0, -3.0, -90.0], [-85.0, 2.0, -25.0, -6.0], [79.0, 5.0, 30.0, -1.0]] 以外の例を作成し，両方の実行結果を比較しなさい．

演習問題(発展)

直接法で(反復によらずに)連立方程式を解く場合は，Gauss-Jordan法ではなく LU分解を使う方が効率的である．LU分解について調べてみなさい．また，自分でプログラムを組まなくてもLAPACK, BLAS等のライブラリを使うのが一般的である．線形代数を扱うためにRubyでどのようなライブラリが利用できるかを調べてみなさい．

モンテカルロ積分

以下のプログラムはモンテカルロ法で円周率を求めるプログラム(pi.rb)である．

include Math
def pi1(n)
  m=0
  for i in 1..n
    x=rand
    y=rand
    if(x*x+y*y)<1.0
      m+=1
    end
  end
  4.0*m/n
end
def pi2(n)
  r=0.0
  for i in 1..n
    x=rand
    r+=Math.sqrt(1-x*x)
  end
  4.0*r/n
end

注: Math.sqrt(x) はxの平方根(square root)を計算する関数 x**0.5と同じ
実行例は以下のようになる．

>> pi1(100)
pi1(100)
=> 2.92
>> pi2(100)
pi2(100)
=> 3.20524531832742

演習問題

プログラム pi1とpi2がそれぞれどのようにして円周率を求めているかを説明し， nを変えた時の誤差に関して観察しなさい．