C言語で連立一次方程式を解く

C言語で掃き出し法を用いて連立一次方程式を解きます。

掃き出し法はガウスの消去法とも呼ばれるもので、線形代数学を学んだことがある人はご存知かと思いますが、簡単に説明すると、連立一次方程式について係数拡大行列を作り、対角成分のみを1とし、他の成分を0とすると対応する変数について解が求まるというものですね。

要するに、

これを行列で表して、

左辺の行列と右辺のベクトルを合わせて拡大係数行列をつくると、

これを変形していきます。今回は行列のrankは3である(x, y, zはそれぞれ独立した解をもつ)場合にのみ限った連立方程式を解いていくと前提にします。すると、基本変形を駆使して最終的には

の形になり、x = l, y = m, z = nが求まるという仕組みです。これをC言語で表現してみましょう。

#define N 4

int i, j, k, n;

/*
  例として、
   x +  y +  z +  w = 10
  2x +  y + 2z +  w = 14
   x + 2y + 3z - 4w = -2
   x -  y -  z +  w =  0
  を解く。
*/

// 拡大係数行列 M
double M[N][N + 1] =
{
	{1,	 1,  1,  1, 10},
	{2,	 1,  2,  1, 14},
	{1,	 2,  3, -4, -2},
	{1,	-1, -1,  1,  0}
};

double pivot, mul;

// 対角成分が1で正規化された階段行列を作る(前進消去)
for (i = 0; i < N; ++i)
{
	// 対角成分の選択、この値で行成分を正規化
	pivot = M[i][i];
	for (j = 0; j < N + 1; ++j)
	{
		M[i][j] = (1 / pivot) * M[i][j];
	}

	// 階段行列を作る為に、現在の行より下の行について
	// i列目の成分が0になるような基本変形をする
	for (k = i + 1; k < N; ++k)
	{
		mul = M[k][i];
		for (n = i; n < N + 1; ++n)
		{
			M[k][n] = M[k][n] -  mul * M[i][n];
		}
	}
}

// 下から上に向かって変数に代入して、独立した解の形にする(後進代入)
// このとき一番下の行はすでに独立した解を得ている
for (i = N - 1; i > 0; --i)
{
	for (k = i - 1; k >= 0; --k)
	{
		mul = M[k][i];
		for (n = i; n < N + 1; ++n)
		{
			M[k][n] = M[k][n] - mul * M[i][n];
		}
	}
}

/* 実行結果:
   M =
   [1 0 0 0 1]
   [0 1 0 0 2]
   [0 0 1 0 3]
   [0 0 0 1 4]
   これより、x = 1, y = 2, z = 3, w = 4
*/

 

 

“C言語で連立一次方程式を解く” への1件のフィードバック

  1. for文の再設定形式のところは「++k」などにしておくと配列の最初の成分は扱われないのでは?

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です