注意c语言中,abs()为求整型的绝对值的函数。对于float型需用fabs(),代码被这个不起眼的错误坑了好久,最后得出结论不熟悉的系统函数一定要慎用,只有在确切知道该函数的顶以后才取使用它,否则自己写函数都比使用不明函数来的好。
c语言版
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <ctime>
const int maxn=1001,maxl=100;
double x[maxn],y[maxn],b[maxn],A[maxn][maxn];
int flag;
int gdsl(int n){
int k=0;
for(int i=0;i<n;i++){
x[i]=0;
}
while(1){
for(int i=0;i<n;i++)
y[i]=x[i];
for(int i=0;i<n;i++){
double z=b[i];
for(int j=0;j<n;j++){
if(j!=i)z-=A[i][j]*x[j];
}
if(abs(A[i][i])<1e-10||k==maxl){
flag=0;return 0;
}
z/=A[i][i];x[i]=z;
printf("%lf ",x[i]);
}
puts("");
for(int i=0;i<n;i++)
printf("%lf ",y[i]);
puts("@");
double norm=0;
for(int i=0;i<n;i++)
if(fabs(y[i]-x[i])>norm)norm=fabs(y[i]-x[i]);
printf("%lf#\n",norm);
if(norm<1e-5)break;
++k;
}
return k;
}
int main(){
int n;
scanf("%d",&n);
for(int i=0;i<n;i++){
for(int j=0;j<n;j++)
scanf("%lf",&A[i][j]);
}
for(int i=0;i<n;i++)
scanf("%lf",&b[i]);
int tmp=gdsl(n);
printf("%d\n",tmp);
for(int i=0;i<n;i++)
printf("%lf\n",x[i]);
return 0;
}
matlab版
function[x,k,flag,nor]=gdsl(A,b)
% A:方程组的系数矩阵
% b:方程组的右边项
% delata:精度1e-5
% maxl:最大迭代次数
n=length(A);k=0;
x=zeros(n,1);y=zeros(n,1);flag='OK!';
delta=1e-5;maxl=100;
while 1
y=x;
for i=1:n
z=b(i);
for j=1:n
if j~=i
z=z-A(i,j)*x(j);
end
end
if abs(A(i,i))<1e-10|k==maxl
flag='Fail!';return;
end
z=z/A(i,i);x(i)=z;
end
if norm(y-x,inf)<delta
break;
end
k=k+1;
end
直接调用gdsl函数即可。