数学とか語学とか楽しいよね

ドイツ語、ロシア語、オランダ語、英語、フランス語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。

【数値計算】Gauss-Seidel法(ガウス・ザイデル法)のC++コード

現在「【数値計算】Gauss-Seidel法(ガウス・ザイデル法)の説明 実践編」を準備中ですが、それに先立ってGauss-Seidel法のC++コードを公開します。私はコーディングの専門家ではないのであまり綺麗ではありませんが悪しからず。皆様の勉強に役立てて頂ければ幸いです。解いている問題は

 \left( \begin{array}{ccc} 3 & 2 & 1 \\ 1 & 4 & 1 \\ 2 & 2 & 5 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} 10 \\ 12 \\ 21 \end{array} \right)

で、答えは
 \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \end{array} \right) = \left( \begin{array}{c} 1 \\ 2 \\ 3 \end{array} \right)

となります。実際に動かして遊んでみて下さい。

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <vector>

using namespace std;

//inner product//
inline double ip(vector<double> &a, vector<double> &b, int &size)
{
	double value=0.;
	int i;
	for(i=0;i<size;i++) value+=a[i]*b[i];
	return value;
}

//Gauss-Seidel
inline void GS(vector<vector<double> > &A,vector<double> &x,vector<double> &b, int &size, double &eps)
{
	int i,k;
	double in,sum,max;
	
	for(k=0;k<1000000;k++)
	{
		max=0.0;
		for(i=0;i<size;i++)
		{
			in=x[i];
			x[i]=(b[i]-(ip(A[i],x,size)-A[i][i]*x[i]))/A[i][i];
			if(max<abs(in-x[i])) max=abs(in-x[i]);
		}
		if(eps>max) break;
	}
}

int main()
{
	int size=3;
	int i;
	double eps=pow(2.0,-50);
	double in,max;
	vector<double> b(size,0.);
	vector<double> x(size,0.);
	vector<vector<double> > A(size,(vector<double>(size,0.)));
	
	A[0][0]=3;A[0][1]=2;A[0][2]=1;
	A[1][0]=1;A[1][1]=4;A[1][2]=1;
	A[2][0]=2;A[2][1]=2;A[2][2]=5;
	b[0]=10;b[1]=12;b[2]=21;
	
	GS(A,x,b,size,eps);

	for(i=0;i<size;i++) cout<<x[i]<<endl;

        return 0;
}