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

/****************** 図7.3の数値解を求めるプログラム *********************

	・∂u/∂t = -u(∂u/∂x)+μ(∂^2u/∂x^2)
		上記バーガース方程式を，アダムス・バッシュフォース法と
			クランク・ニコルソン法を用いて解きます．

	・係数行列は3重対角行列なので，メモリを有効に使うため，
		A[M-1][M-1] → A[M-1][3+1]
			として値のあるところのみを配列に代入します．

	・初期条件を変えて
		グラフがゆがむような場合は，
			Nを大きくするなどしてみて下さい．

	・結果は，標準出力のみでファイルへの出力はありません．
		グラフにする場合は，ファイルにリダイレクトするか，
			ファイル出力のプログラムを書き加えて下さい．

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

  #define eps pow(10,-15)
  #define Pi acos(-1.)
  #define Mu 0.01

/******************************/
/********** 初期条件 **********/
/******************************/
  /******  t=[t0, tN]  ******/
  #define t0 0
  #define tN 2
  #define N  200			//区間[t0, tN]の分割数
  #define k  (tN-t0)/(double)N		//時間ステップ
/******************************/
  /******  x=[x0, xM]  ******/
  #define x0 0
  #define xM 1
  #define M  128			//区間[x0, xM]の分割数
  #define h  (xM-x0)/(double)M		//差分間隔
/******************************/
/******************************/

/******************************/
/********** 境界条件 **********/
/******************************/
  #define u0_t 0			//u(x0, t)
  #define uM_t 0			//u(xM, t)
/******************************/
/******************************/

double u_0(double x){			//u(x, 0)
	double y;
	y = sin(2*Pi*x);
	return y;
}

/***************************************************/
/************ 係数行列を初期化する関数 *************/
/***************************************************/

void initialize_(double A[M][3+1], double p, double q){
	int j;

//	A[1][1] = 0;
	A[1][2] = q;
	A[1][3] = p;
	A[M-1][1] = p;
	A[M-1][2] = q;
//	A[M-1][3] = 0;

	for(j=2; j<M-1; j++){
		A[j][1] = p;
		A[j][2] = q;
		A[j][3] = p;
	}
}
/***************************************************/


/****************************************************************/
/******************* 三重対角行列を解く関数 *********************/
/****************************************************************/
/**************       m : xの分割区間数        ******************/
/**************       A : 係数行列             ******************/
/**************       u : 次ステップの解       ******************/
/**************       b : 右辺                 ******************/
/****************************************************************/

int TRI(int m, double A[M][3+1], double u[M+1], double b[M]){

	int j;
	double temp;

	/***** 前進消去により，上三角行列に直す *****/
	for(j=2; j<m; j++){
		temp = A[j][1]/A[j-1][2];
		A[j][2] -= A[j-1][3]*temp;
		b[j] -= b[j-1]*temp;
	}

	/***** 後退代入により，次ステップのu(x, t)を求める *****/
	/***** 配列のA[j][2]にj行目のピボットが入っている  *****/
	if(fabs(A[m-1][2]) < eps){
		printf("divided by zero!!\n");
		return 2;
	}
	u[m-1] = b[m-1]/A[m-1][2];		//n+1:次ステップ

	for(j=m-2; 0<j; j--){
		if(fabs(A[j][2]) < eps){
			printf("divided by zero!!\n");
			return 3;
		}
		u[j] = (b[j]-A[j][3]*u[j+1])/A[j][2];
	}

	return 1;		//正常に終了したときのみ，1(真)を返す
}
/****************************************************************/
/****************************************************************/
/****************************************************************/

int main(void){
	int j, n, i;		
	double A[M][3+1], x[M+1], u_o[M+1], u_n[M+1], b[M];
	double p, q;

	p = -Mu*k/(2*h*h);
	q = 1+Mu*k/(h*h);

	/***** 境界条件を代入 *****/
	u_o[0] = u0_t;
	u_o[M] = uM_t;
	/**************************/

	/******* 初期条件を代入 *******/
	x[0] = x0;
	x[M] = xM;
	for(j=1; j<M; j++){
		x[j] = x0+(double)j*h;
		u_o[j] = u_0(x[j]);
	}
	/******************************/

	for(j=0; j<M+1; j++)
		printf("%5.3lf, %12.7lf\n", x[j], u_o[j]);

	for(j=1; j<M; j++){
		b[j] = u_o[j]-k*(u_o[j]*(u_o[j+1]-u_o[j-1]))/(2*h)					-p*(u_o[j+1]-2*u_o[j]+u_o[j-1]);
					//式(7.29)を中心差分に直した式

	}

	initialize_(A, p, q);		//係数行列の初期化
	if( TRI(M, 0, A, u_n, b) != 1 )	return 1;

//	for(j=0; j<M+1; j++)
//		printf("%5.3lf, %12.7lf\n", x[j], u_n[j]);

	for(n=1; n<N; n++){			//n: 現在の時間ステップ
		for(j=1; j<M; j++){
			b[j] = u_n[j]										-k*(3*u_n[j]*(u_n[j+1]-u_n[j-1])						-u_o[j]*(u_o[j+1]-u_o[j-1]))/(4*h)						-p*(u_n[j+1]-2*u_n[j]+u_n[j-1]);
						//式(7.29)
		}

		for(j=0; j<M+1; j++)
			u_o[j] = u_n[j];

		initialize_(A, p, q);		//係数行列の初期化
		if( TRI(M, n, A, u_n, b) != 1 )	return n+1;

		if( fmod(n+1,20)==0 )
			for(j=0; j<M+1; j++)
				printf("%5.3lf, %12.7lf\n", x[j], u_n[j]);
	}

	return 0;
}