Thomas 求解三對角矩陣 （追趕法）

此處只給定一個函數，傳入的數組計算是從1開始，而不是從0開始

這個版本比較複雜，可是便於理解求解過程

int TM(double a[], double b[], double c[], double r[],  double x[],int n)
{
	//原矩陣形式 要求 對角佔優 a(1)=0 c(n)=0
	//b1 c1
	//a2 b2 c2
	//   a3 b3 c3
	//    ···
	//         an bn
	//分解後的矩陣 LU分解 上三角爲單位陣 
	//alpha                           1 beta
	//gamma alpha                         1  beta
	//      gamma alpha                        1  beta
	//            `````` ```                        ``  ``
	//                   gamma alpha                     1
	// Ax=b  LUx=d  Ux=y  Ly=d
	//僅對A分解
	double *gamma, *alpha, *beta, *y;
	gamma=new double[n+1];
	alpha=new double[n+1];
	beta=new double[n+1];
	y=new double[n+1];
	gamma[1]=0;//gamma沒有1,這裏設爲0
	alpha[1]=b[1];
	beta[1]=c[1]/alpha[1];
	for (int i=2;i<=n-1;i++){
		gamma[i]=a[i];
		alpha[i]=b[i]-beta[i-1]*gamma[i];
		beta[i]=c[i]/alpha[i];//beta只到n-1
	}
	gamma[n]=a[n];
	alpha[n]=b[n]-beta[n-1]*gamma[n];
	//由Ly=d求y 追
	y[1]=r[1]/alpha[1];
	for (int i=2;i<=n;i++){
		y[i]=(r[i]-gamma[i]*y[i-1])/alpha[i];
	}
	//由Ux=y求x 趕
	x[n]=y[n];
	for (int i=n-1;i>=1;i--){
		x[i]=y[i]-beta[i]*x[i+1];
	}
	delete[] gamma;
	delete[] alpha;
	delete[] beta;
	delete[] y;
	return 0;
}