Thomas 求解三对角矩阵（追赶法）

最新推荐文章于 2022-10-12 14:27:36 发布

转载最新推荐文章于 2022-10-12 14:27:36 发布 · 1.4k 阅读

4 ·

CC 4.0 BY-SA版权

原文链接：https://my.oschina.net/u/929434/blog/98603

文章标签：

#python #matlab #ux

本文介绍了一个使用LU分解求解线性方程组的C++函数实现。该函数采用追赶法进行求解，并详细展示了从矩阵分解到最终求解的过程。通过实例解释了如何将原始矩阵转换为上三角和下三角矩阵，再逐步求得未知数。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

2019独角兽企业重金招聘Python工程师标准>>>

此处只给定一个函数，传入的数组计算是从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;
}

转载于:https://my.oschina.net/u/929434/blog/98603