2014-05-16

連立一次方程式

計算統計学の講義の中では順を追って作って行きますが、完成したものをすぐ使う場合のために載せておきます。
解が一組だけ存在する方程式用です。
最小二乗法の最適解を求める場合は未知数の個数と式の個数が一致します。
また最適値が一組でない場合は無数にあり、最適値がないということはありません。
そこでこのプログラムでは、一意に定まらない未知数は0とおくことで一組だけ求めるようにしています。
実際には行列式が0に近い時に一般逆行列を使ったほうが良いですが、そのような処理は入れていません。

Sub 連立一次方程式(係数(), 定数(), 解(), r)
Rem (c) 2014 Kaoru Fueda
Dim 係数行列() As Double, 最大値() As Double
Dim i As Integer, j As Integer, k As Integer, 最大行 As Integer, 最小行 As Integer
Dim temp As Double, 最大比 As Double, 最小値 As Double
ReDim 係数行列(r, r + 1), 最大値(r)
For i = 1 To r
  For j = 1 To r
    係数行列(i, j) = 係数(i, j)
  Next j
  係数行列(i, r + 1) = 定数(i)
Next i
Rem iはピボット
For i = 1 To r
Rem ピボット選択
  If i < r Then
    最大比 = 0
    最大行 = i
    For j = i To r
      最大値(j) = 0
      For k = i To r + 1
        If Abs(係数行列(j, k)) > 最大値(j) Then
          最大値(j) = Abs(係数行列(j, k))
        End If
      Next k
      If 最大値(j) > 0 Then
        If 最大比 < Abs(係数行列(j, i) / 最大値(j)) Then
          最大比 = Abs(係数行列(j, i) / 最大値(j))
          最大行 = j
        End If
      End If
    Next j
    If 最大行 > i Then
      For k = 1 To r + 1
        temp = 係数行列(i, k)
        係数行列(i, k) = 係数行列(最大行, k)
        係数行列(最大行, k) = temp
      Next k
    End If
Rem ピボットが0の場合は最も0=0に近い行を第i行に移動して、後に無視する
    If 最大比 = 0 Then
      最小値 = 最大値(i)
      最小行 = i
      For j = i + 1 To r
        If 最小値 > 最大値(j) Then
          最小値 = 最大値(j)
          最小行 = j
        End If
      Next j
      If 最小行 > i Then
        For k = 1 To r + 1
          temp = 係数行列(i, k)
          係数行列(i, k) = 係数行列(最小行, k)
          係数行列(最小行, k) = temp
        Next k
      End If
    End If
  End If
Rem ピボットを1に
  If 係数行列(i, i) <> 0 Then
    For k = i + 1 To r + 1
      係数行列(i, k) = 係数行列(i, k) / 係数行列(i, i)
    Next k
  Else
    MsgBox ("解けないので解の一つを0とします。")
Rem 本来は解ベクトルの長さが最小とすべき。
    For k = i + 1 To r + 1
      係数行列(i, k) = 0
    Next k
  End If
  係数行列(i, i) = 1
Rem ピボット以外を0に
  For j = 1 To r
    If j <> i Then
      For k = i + 1 To r + 1
        係数行列(j, k) = 係数行列(j, k) - 係数行列(i, k) * 係数行列(j, i)
      Next k
      係数行列(j, i) = 0
    End If
  Next j
Next i
For i = 1 To r
  解(i) = 係数行列(i, r + 1)
Next i
End Sub

0 件のコメント: