数値解析課題-ヤコビ法，ガウス・ザイテル法

3次の行列しか対応せず
解が収束しない場合の対応は考えていない
一応動くけど細かい動作は未検証
プログラム的にはLU分解よりなんぼか簡単

ちょっと修正．

// 数値解析-ヤコビ・ガウスザイテル法.cpp 
// 解が収束しなかった場合の処理は特に考えていない

#include "stdafx.h"
#include "iostream"
#include "cstdio"
#include "cmath"

#define N 3
#define EPSI 0.001	//ε
#define LMT 100000	//最大反復回数

bool converge_judge(void);	//収束(事前)判定
int yakobi(void);
int gaussdel(void);

double mat[N][N] = { { 9.0, 3.0, 5.0 },  //簡単のためグローバル変数
						{ 2.0, 9.0, 6.0 },
						{ 4.0, 3.0, 9.0 }
	};
double b[N] = { 30.0, 38.0, 37.0 };
double x[N] = { 0.0, 0.0, 0.0 };			//解を格納
double x2[N];								//解を一時的に格納(ヤコビ法で使用)

int main(int argc, _TCHAR* argv[])
{
	int c;

	if (converge_judge() == true){
		printf("収束判定:真\n");
	}

	//c = gaussdel();		//ガウスザイテル法
	c = yakobi();			//ヤコビ法

	printf("x={ %1.3lf, %1.3lf, %1.3lf }\n",x[0],x[1],x[2]);	//解の表示
	printf("反復回数：%d\n",c);
	return 0;
}

bool converge_judge(void){	//収束判定(解が収束するかどうか) 
	for(int i = 0;i < N ;i++){
		for(int j = 0;j < N;j++){
			if(i !=j){
				if(mat[i][i] < mat[i][j]){
					return false;
				}
			}
		}
	}
	return true;
}

int yakobi(void){		//ヤコビ法
	double s = 0.0, emax = 0.0;
	for(int c = 0;c < LMT;c++){
		emax = 0.0;
		for(int i = 0;i < N;i++){
			s = b[i];
			for(int j = 0;j < N;j++){
				s -= mat[i][j] * x[j];
			}
			s = s / mat[i][i];
			x2[i] = x[i] + s;
                        if(x[i] != 0){
			        if(emax < fabs(1 - (x2[i] / x[i]) )){	//相対誤差の最大を求める
				        emax = fabs(1 - (x2[i] / x[i]) );
			        }
                        }
                        else{                   //ループの１回目で終了しないように
                         emax = 1.0;
                        }
		}

		for(int k = 0;k < N;k++){
			x[k] = x2[k];	
		}
		if(emax <= EPSI){
			return c;
		}
		
	}
	return 0;		//解が収束しなかった
}

int gaussdel(void){		//ガウス・ザイテル法
	double s = 0.0, e = 0.0, emax = 0.0;
	for(int c = 0;c < LMT;c++){
		emax = 0.0;
		for(int i = 0;i < N;i++){
			s = b[i];
			for(int j = 0;j < N;j++){
				s -= mat[i][j] * x[j];
			}
			s = s / mat[i][i];
			e = x[i] + s;				//上書きする前に相対誤差を求める
                        if(x[i] != 0){
			        if(emax < fabs(1 - (e / x[i]) )){
				        emax = fabs(1 - (e / x[i]) );
			        }
                        }
                        else{                   //ループの１回目で終了しないように
                         emax = 1.0;
                        }
                        x[i] = e;
		}

		if(emax <= EPSI){
			return c;
		}
		
	}
	return 0;		//解が収束しなかった