c++でガウスの消去法(部分ピボッティング付き)

programming c/c++

前のエントリ - c++でガウスの消去法 - で
http://d.hatena.ne.jp/yatt/20080214/1202987642

ピボッティングは0、あるいは0に非常に近い値で除算しないため、また、精度を高めるためにやると追記しました。ただ、その内容に関しては自分で実際に確かめてはいませんでした。

そこで理解のために、部分ピボッティングするガウスの消去法のプログラムを作って見ました。前に作った、ピボッティングを行わないプログラムと比較してみようという試みです。

// gauss_piv.cpp
#include <iostream>
#include <algorithm> // swap
using namespace std;

#define N (3)
#define FMT "%5.7f "
float A[N][N], x[N], B[N]; // Ax = Bの関係


////
void initialize()
{
  for (int i = 0; i < N; i++){
    for (int j = 0; j < N; j++)
      cin >> A[i][j];
    cin >> B[i];
  }
}

void print_A()
{
  for (int i = 0; i < N; i++){
    for (int j = 0; j < N; j++)
      printf(FMT, A[i][j]);
    printf("\n");
  }
}

void print_B()
{
  for (int i = 0; i < N; i++)
    printf(FMT, B[i]);
  printf("\n");
}
void print_x()
{
  for (int i = 0; i < N; i++)
    printf(FMT, x[i]);
  printf("\n");
}
////


void partial_pivoting(int row)
// 部分ピボッティング
{
  int   max_line = row; // 絶対値が最大の要素を持つ行
  float max_elem = fabsf(A[row][row]);

  // 絶対値が最大の要素を持つ行を探す
  for (int i = row+1; i < N; i++)
    if (max_elem < fabsf(A[i][row])){
      max_elem = fabsf(A[i][row]);
      max_line = i;
    }
  // row行とmax_line行を交換
  for (int i = 0; i < N; i++)
    swap(A[row][i], A[max_line][i]);
  // Bも交換
  swap(B[row], B[max_line]);
  ///printf("partial_piv row=%d\n",row);  print_A();
}

void forward_elimination()
{
  // 前進消去
  for (int i = 0; i < N-1; i++){
    partial_pivoting(i);
    for (int j = i+1; j < N; j++){
      float ratio = A[j][i] / A[i][i];
      for (int k = i+1; k < N; k++) // 消去後の値を綺麗に出力させるならk=0で初期化
	A[j][k] -= ratio * A[i][k];
      B[j] -= ratio * B[i];
    }
  }
}

void backward_substitution()
{
  // 後退代入
  x[N-1] = B[N-1] / A[N-1][N-1];
  for (int i = N-2; i >= 0; i--){
    for (int j = i+1; j < N; j++)
      B[i] -= A[i][j] *x[j];
    x[i] = B[i] / A[i][i];
  }
}

void solve()
{
  forward_elimination();
  backward_substitution();
}


int main()
{
  initialize();

  solve();
  // print x
  puts("x:");
  print_x();
  return 0;
}


まずは消去中に対角要素が0になる場合。前のエントリで書いたプログラムは、解を小数点第7位まで表示させるように変更しました。

問：
 A          x     = B
 ( 1  1  1)(x1) = (6)
 ( 2  2 -1)(x2) = (3)
 (-1  3  1)(x3) = (8)
解：
x1 = 1
x2 = 2
x3 = 3

まずは、前回のピボッティングなしのプログラムに入力してみました。

$ ./a.out
input n: 3
1 1 1 6
2 2 -1 3
-1 3 1 8
x:
  nan   nan   nan

解が求まっていません。次に部分ピボッティングするプログラムで実行。

$ ./a.out
1 1 1 6
2 2 -1 3
-1 3 1 8
x:
1.0000000 2.0000000 3.0000000

部分ピボッティングすることで対角要素が最大の行と交換され、解が求まっています。


今度は、精度が劣化する問題の例。

問：
 A              x    = B
(0.0003  1 1)(x1) = 5
(     1 -1  2)(x2) = 6
(    -1  3  1)(x3) = 8
解：
x1 = 10000/9979 = 1.0021044192804891...
x2 = 19970/9979 = 2.0012025253031367...
x3 = 29922/9979 = 2.9984968433710795...

前回のプログラムで実行

$ ./a.out
input n: 3
0.0003 1 1 5
2 -1 2 6
-1 3 1 8
x:
1.0013580 2.0004263 2.9992733

小数点以下第2位まで解に一致。次に今回作ったプログラムで実行。

$ ./a.out
0.0003 1 1 5
2 -1 2 6
-1 3 1 8
x:
1.0021031 2.0012019 2.9984977

解の精度が3桁改善されました。


ピボッティングに効果があることを確認できました。部分ピボッティングよりも精度を高められる完全ピボッティングというのもあるようです。ググってそれっぽいページも発見。
http://www.eng.auburn.edu/department/me/courses/nmadsen/me231/gauss.html