Gauss-Seidel迭代法求方程组的解

今天药忘吃喽~ 2022-07-28 04:06 87阅读 0赞

注意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函数即可。