小弟写了一个解线性方程组的雅克比迭代算法，但是发现有时候会一直迭代，不收敛

[richardzrc]

附上代码

#include <iostream>
using namespace std;
bool Converge(int A[4][4],double x[4],int b[4],int n);
void Show(int M[4][4],int b[4],double x[4],int n);
void Jacob(int A[4][4], double x[4],int b[4],int n)
{
	int i,j,iterative_count=0;
	double *newx=new double[n+1];
	double sigma;
	double delta;
	//while(!Converge(A,x,b,n))
	while(1)
	{
		for(i=1;i<=n;i++)
		{
			sigma=0;
			for(j=1;j<=n;j++)
			{
				if(i!=j)
				{
					sigma+=A[i][j]*x[j];
				}
			}
			newx[i]=double(double(b[i])-sigma)/A[i][i];
		}
		//max|x[i]-newx[i]|最大值作为误差分析
		delta=abs(x[1]-newx[1]);
		for(i=2;i<=n;i++)
		{
			if(delta<abs(x[i]-newx[i]))
				delta=abs(x[i]-newx[i]);
		}
		if(delta<1e-7)
			break;
		for(i=1;i<=n;i++)
			x[i]=newx[i];
		iterative_count++;
		if(iterative_count==6||iterative_count==100)
		{
			cout<<endl<<iterative_count<<endl;
			Show(A,b,x,n);
		}
	}
	delete []newx;
}
bool Converge(int A[4][4],double x[4],int b[4],int n)
{
	int i,j;
	double delta=0,tmp;
	for(i=1;i<=n;i++)
	{
		tmp=0;
		for(j=1;j<=n;j++)
		{
			tmp+=double(A[i][j]*x[j]);
		}
		//delta+=(tmp-b[i])*(tmp-b[i]);
		delta+=abs(tmp-b[i]);
	}
	//delta=sqrt(delta);
	if(delta<1e-4)
		return true;
	else
		return false;
}
void initial(int M[4][4],double x[4],int b[4],int n)
{
	int i,j;
	for(i=1;i<=n;i++)
	{
		x[i]=0;
		/*b[i]=rand()%10;
		for(j=1;j<=n;j++)
			M[i][j]=rand()%10;*/
	}
	

}
void Show(int M[][4],int b[4],double x[4],int n)
{
	int i,j;
	cout<<"M:"<<endl;
	for(i=1;i<=n;i++)
	{
		for(j=1;j<=n;j++)
			cout<<M[i][j]<<" ";
		cout<<endl;
	}
	cout<<"b:"<<endl;
	for(i=1;i<=n;i++)
		cout<<b[i]<<endl;
	cout<<"x:"<<endl;
	for(i=1;i<=n;i++)
		cout<<x[i]<<endl;
}
int main()
{
	const int n=3;
	int M[n+1][n+1]={{0,0,0,0},{0,16,4,78},{0,4,5,-4},{0,8,-4,22},/*{0,3,-7,2},{0,9,12,4,7,8}*/},b[n+1]={0,-4,3,10};
	double x[n+1];
	int i;
	/*for(i=0;i<=n;i++)
		M[i]=new int[n+1];*/
    //initila M and b,x;
	initial(M,x,b,n);
	Show(M,b,x,n);
	Jacob(M,x,b,n);
	cout<<"x:"<<endl;
	for(i=1;i<=n;i++)
		cout<<x[i]<<endl;
	//delete []x;
	system("pause");

}

另外有一个问题，如果方程组的系数都是随机生成的，那么是不是一定能算出一个结果呢？这是我一直在考虑的问题，方程组无解，也能计算出一个近似解的吧，但是为啥我的程序就是算不出呢？

[kerbcurb]

迭代法需要一定的条件，比如AX = B;A = D + C,D是对角上矩阵；X1 = B/D - (C/D )X，要求：||C/D||<1

[richardzrc]

[Quote=引用 1 楼 的回复:]

迭代法需要一定的条件，比如AX = B;A = D + C,D是对角上矩阵；X1 = B/D - (C/D )X，要求：||C/D||<1
[/Quote]

但是好像计算出近似解的吧，这个精确的数学方法计算不一样吧，程序如果满足设定的条件 应该就可以收敛的吧。

[horris]

请参考一下矩阵应用解线性方程组方面的数学资料，不是什么样的系数都能迭代收敛的，随机的系数肯定会有发散的情况出现，而且还会出现病态的系数矩阵，就是说系数发生很小的变化，引起解的很大的变化。

[kerbcurb]

看来你对“近似解”的理解不明确。
如果一个线性方程组有解，不管用什么算法，计算机求出来都是近似的！所谓“近似解”是指这个方程组有解，求出来的结果与真实解很接近或者说有一定的误差，但是误差很小。

而“无解”是指这个方程组没有解！
如果不满足那个条件，数值求解的结果发散。

[horris]

[Quote=引用 4 楼 的回复:]
看来你对“近似解”的理解不明确。
如果一个线性方程组有解，不管用什么算法，计算机求出来都是近似的！所谓“近似解”是指这个方程组有解，求出来的结果与真实解很接近或者说有一定的误差，但是误差很小。

而“无解”是指这个方程组没有解！
如果不满足那个条件，数值求解的结果发散。
[/Quote]
部分不同意楼上的观点，几乎所有的迭代方法，都会有发散的情况，主要取决于方程组的系数（矩阵），和方程组有没有解好象也有关系，但大多数情况下，不收敛主要是由于方程组的系数不适宜迭代,还有某些系数会造成病态方程组，这时可以考虑对方程组做一些变换。

[cnmhx]

同学，象这样的问题，使用已有的程序最好。

[horris]

查了下数学书（Jacobi迭代法解线性方程组，在随便哪本矩阵应用或解线性方程组的书中都会讲到），总结如下：
首先要把线性方程组表示成矩阵形式，设其系数矩阵为A，右端项为b。
1。方程组是否能用Jacobi解，前提是A非奇异。线性代数的东西忘的差不多了，只记得A可逆则A是非奇异的。A非奇异，则线性方程组有唯一解，这时才可以用Jacobi迭代。
2。Jacobi迭代的矩阵表示，是把A分解为一个对角阵D、上三角阵U、下三角阵L，A=D-L-U，则M=D^-1(L+U)称为Jacobi迭代矩阵（D^-1表示D的逆）。Jacobi迭代是否收敛，取决于迭代矩阵M。M只依赖于系数矩阵A。
1）迭代矩阵M的所有特征值的模都小于1，是Jacobi迭代收敛的充要条件。
2）M的1范数、2范数、无穷范数，这三者中的任何一个，如果小于1，则Jacobi迭代收敛。范数越小，则收敛越快。（矩阵的范数是什么东东，请看矩阵方面的资料，矩阵的范数相当于实数的绝对值）
3）若系数矩阵A按行（列）严格对角占优，且非奇异，则Jacobi迭代收敛。（如果A的对角线元素的模大于同行（列）中其他元素的模的和，则称A按行（列）严格对角占优）。

由上可见，线性方程组是否可用Jacobi迭代、是否收敛，只取决于系数矩阵A，右端项都没有影响。所以楼主的系数随机产生，肯定不行，可能会出现：1）无解或有无穷组解，不适用Jacobi迭代；2）Jacobi迭代不收敛。
那么，在迭代前，应检查系数矩阵A，不适用的话要终止计算并提示用户。
这个检查涉及到：
1）判断系数矩阵A是否可逆，这不难（其实我也忘了怎么判断，但应该有成熟算法，记得对称阵必定可逆）。
2）若系数矩阵A按行（列）严格对角占优，则Jacobi迭代收敛。判断对角占优很直观，编程很简单。若A严格对角占优，则不必进行下面的判断了。
3）将系数矩阵A分解为对角阵D、上三角阵U、下三角阵L，构造Jacobi迭代矩阵M=D^-1(L+U)，这很简单。
4）计算M的任意一个算子范数（1范数、2范数、无穷范数），只要三个范数有一个小于1，则Jacobi迭代收敛。不必进行下面的判断了。范数听起来很深奥，其实计算很简单，我记得1、2范数好象是矩阵A的行或列的元素的模的最大值，无穷范数好象是所有元素的模的最大值。
5）计算M的所有特征值，如果都小于1，则Jacobi迭代收敛。矩阵特征值的数值计算，有算法，可能这个稍复杂点。

就介绍这些，其中的概念请看线性代数的书或google，算法请google。

[horris]

其实在工程上，方程组的系数矩阵严格对角占优是很常见的，即使不满足，也可以将方程组变换一下，比如变换一下未知变量的顺序。你在随机产生参数时，加上对角占优的限制也很简单。
保证系数矩阵非奇异，即可逆，稍有复杂。但是工程实践上，系数矩阵经常是对称的。比如在建立数学模型时，或测量数据时，刻意做成对称的系数矩阵就可以了。你在随机产生参数时，保证系数对称也很简单。
总之，要写一个实用的迭代算法，可能应用场合就保证了你的迭代是适用的。但是如果你要写一个教学或实验算法，就必须验证用户给的系数。

[kerbcurb]

horris:你不同意我哪一部分？1楼和4楼是我的回答。
我看你的回复里面有很多“几乎”，“一般”，“好像”，“大多数情况下”等修饰语。这样很不严谨，容易误导别人，切记！

[FancyMouse]

你输入的矩阵不满足收敛条件。当然不收敛。只有收敛了才能求近似解。

[richardzrc]

谢谢各位大侠，小弟找到有解的数据了！

[kerbcurb]

horris：
我觉得你比较自以为是，你从头到尾看看楼主说的，我说的和你说的再回想一下自己说的什么。不过到此为止，我不再与你说了！

[horris]

[Quote=引用 9 楼 的回复:]
horris:你不同意我哪一部分？1楼和4楼是我的回答。
我看你的回复里面有很多“几乎”，“一般”，“好像”，“大多数情况下”等修饰语。这样很不严谨，容易误导别人，切记！
[/Quote]
不要看到别人反对你的意见就急着吐槽。
看了你的楼层，你没有把问题说清楚，并非无解时迭代不收敛，系数矩阵不满足一定条件时，即使有解，也必定不收敛。楼主的问题是为什么不收敛。
我用了那些模糊的词，是因为对数学细节记不清楚了，也懒得查书去，但关键地方自认为讲的很严谨。“一般”、“大多数”这种词，对于讲解数学问题总是得用到，因为你一条结论，总是有例外的地方，提问者如果想了解细节，给他一个指路牌，让他自己查详细资料去，我不是老师，没有时间没有义务把每个细节都讲到。

[horris]

比如楼主做的算法如果是要在工程上实用的，那么大多数这种工程实践的系数矩阵都是对角占优的甚至是对称的，不然数学上研究对角占优和对称阵干什么？这种表达怎么能算是模糊呢？

[horris]

[Quote=引用 4 楼 的回复:]
看来你对“近似解”的理解不明确。
如果一个线性方程组有解，不管用什么算法，计算机求出来都是近似的！所谓“近似解”是指这个方程组有解，求出来的结果与真实解很接近或者说有一定的误差，但是误差很小。

而“无解”是指这个方程组没有解！
如果不满足那个条件，数值求解的结果发散。
[/Quote]
最后两句话我可能误解了，我理解成你要说“无解即迭代发散”。怪不得你看起来火气这么大。这也不能全怪我，你文字上没表达清楚，应该把无解时根本不能用迭代法讲出来。不管有意无意，这真有不严谨和误导人之嫌。

[Quote=引用 1 楼 的回复:]
迭代法需要一定的条件，比如AX = B;A = D + C,D是对角上矩阵；X1 = B/D - (C/D )X，要求：||C/D||<1
[/Quote]
这里，你给的公式是通用的迭代格式，而不是具体的Jacob迭代。

也就是说，从技术层面上，你说的没错，我说的也没错。
说完技术，再讲为人：
你说我不严谨，误导人，这就强词夺理了。我给的比你的具体得多，就差我直接写出代码了。我在算法关键点上没有模棱两可。“几乎”，“一般”，“好像”，“大多数情况下”等修饰语，即使在数学和工程实践中也不能避免使用这些词，而且我基本上用在了和编程有关的地方。我没有义务去讲清楚怎样写具体代码。楼主完全可以凭我给的信息，去网上搜索的关键词。授人以渔而不是授人以鱼。而你这样一遇到点不同意见就蹦高的人，才是自以为是，自找别扭。把你评价我的词放到你身上想想吧。

[horris]

忘了说一个重要条件：Jacobi迭代要求线性方程组有唯一解，所以系数矩阵要非奇异，而矩阵可逆是非奇异的充分条件。