#include<stdio.h>
#include<math.h>

/*************** 表5.3の数値解を求めるプログラム ****************

	・式(5.39)を用いて解を求めます．

	・ただし，i=0,1 とします．
		  l=0,1

*****************************************************************/

#define eps pow(10,-8)		//収束判定条件
#define w 1.2			//緩和係数
#define kmax 1000		//収束判定条件との比較回数の上限

int main(void){
	int i, l, k;
	double L2n, sum1, sum2;
	double a[2][2] = {{5,4},{2,3}};
	double b[2] = {13,8};
	double x1[2] = {0,0}, x2[2];

	for(k=0; k<kmax+1; k++){
		for(i=0; i<2; i++){
			sum1 = 0.;
			sum2 = 0.;
			for(l=0; l<i; l++)
				sum1 += -a[i][l]*x2[l];
			for(l=i+1; l<2; l++)
				sum2 += -a[i][l]*x1[l];
			x2[i] = (1-w)*x1[i]+w*(sum1+sum2+b[i])/a[i][i];
		}

		for(i=0; i<2; i++){
			printf(" k=%2d, x[%2d]=%13.6e\n"							, k+1, i+1, x2[i]);
		}

		L2n = 0.;
		for(i=0; i<2; i++)
			L2n += pow((x2[i]-x1[i]),2);

		if(fabs(L2n) < eps)		return 0;

		for(i=0; i<2; i++)
			x1[i] = x2[i];
	}
	return 1;
}