3元連立方程式をガウスの消去法で解いてみた

タイトルまんまです（ガウスの消去法についてはwiki参照のこと）。逆行列を求めるのに似たような方法を使うのですが、微妙に忘れかけてたので復習を兼ねて解いてみた。本当は精度の問題とかもあって、もっとちゃんとしたアルゴリズムがあるらしいが、そんなのは知らん。あまりにひどいハードコーディングなので、明日はN元まで拡張するつもり。

view sourceprint?

01	#include <iostream>

02	#include <iomanip>

03

04	void PrintArray(const double *pArray, int col, int row)

05

{

06	for(int i=0; i<col; ++i)

07

{

08	for(int j=0; j<row; ++j)

09

{

10	std::cout << std::setw(10) << pArray[i*row+j] << " ";

11

}

12	std::cout << std::endl;

13

}

14	std::cout << "------------------------------" << std::endl;

15

16

}

17

18

// 入力用データ

19	// x,y,z = 2, 3, 4

20	// 2x + 3y + 2z = 21 (1)

21	// 3x + 4y + 3z = 30 (2)

22	// 5x + 2y + 4z = 32 (3)

23	double mat[4*4] =

24

{

25	2, 3, 2, 21,

26	3, 4, 3, 30,

27	5, 2, 4, 32,

28	1, 1, 1, 1,

29

};

30

31	int main(int argc, char **argv)

32

{

33

//計算用配列

34	double tmp[4*4] = {1};

35

36	//まずは(1)式のxの係数を1にする

37	for(int i=0; i<4; ++i)

38

{

39	tmp[i] = mat[i] / mat[0];

40

}

41	PrintArray(tmp, 4, 4);

42

43	//(2),(3)式のxの係数を１に合わせる

44	for(int i=0; i<4; ++i)

45

{

46	tmp[1*4+i] = mat[1*4+i] / mat[1*4];

47	tmp[2*4+i] = mat[2*4+i] / mat[2*4];

48

}

49	PrintArray(tmp, 4, 4);

50

51	//(1)式から(2),(3)式を引いたものをtmp[2,3][0-3]に格納

52	for(int i=0; i<4; ++i)

53

{

54	tmp[1*4+i] = tmp[0*4+i] - tmp[1*4+i];

55	tmp[2*4+i] = tmp[0*4+i] - tmp[2*4+i];

56

}

57	PrintArray(tmp, 4, 4);

58

59	//次に(2)式のyの係数と(3)式のyの係数を合わせる

60	double yk = tmp[1*4+1] / tmp[2*4+1];

61	for(int i=1; i<4; ++i)

62

{

63	tmp[2*4+i] = tmp[2*4+i] * yk;

64

}

65	PrintArray(tmp, 4, 4);

66

67	//(2)式から(3)式を引く

68	for(int i=1; i<4; ++i)

69

{

70	tmp[2*4+i] = tmp[1*4+i] - tmp[2*4+i];

71

}

72	PrintArray(tmp, 4, 4);

73

74

//z成分を求める

75	double z = tmp[3*4+2] = tmp[2*4+3] / tmp[2*4+2];

76	PrintArray(tmp, 4, 4);

77

78

//y成分を求める

79	double y = tmp[3*4+1] = (tmp[1*4+3] - (z*tmp[1*4+2])) / tmp[1*4+1];

80	PrintArray(tmp, 4, 4);

81

82

//xを求める

83	double x = tmp[3*4+0] = (tmp[0*4+3] - (z*tmp[0*4+2] + y*tmp[0*4+1])) / tmp[0*4+0];

84	PrintArray(tmp, 4, 4);

85

86	std::cout <<"(x,y,z) = (" << x << ", " << y << ", " << z << ")" << std::endl;

87

88

return 0;

89

}

どうでもいいがなんで2次元配列使わなかったんだろう俺…。あと所々でPrintArrayを呼び出してるのはデバッグのためです。あしからず。