Gauss列主元消去法（封装函数）

àì夳堔傛蜴生んèń 2022-12-28 00:47 125阅读 0赞

//高斯列主元消去法（封装函数一）
    #include<bits/stdc++.h>
    #include "windows.h"
    using namespace std;
    int SolverEqGuass(double **A,double *b,int n,double eps,double *x)
    { 
    	int m,i,j,k;
    	double **a =new double *[n];	
    	double *btemp= new double[n];
    	double d,temp;
    	for(i=0;i<n;i++)
    	{ 
    		a[i]=new double[n+1];	//系数矩阵+常数项矩阵
    	}
    	m=n;
    	//赋值
    	for(i=0;i<n;i++)
    	{ 
    		for(j=0;j<n;j++)
    		{ 
    			a[i][j]=A[i][j];
    		}
    		a[i][n]=b[i];
    	} 
    	if(m<=0||n<=0)
    	{ 
    		cout<<"Input Mistake"<<endl;
    	}
    	for(k=0;k<n-1;k++)	//找列主元最大值
    	{ 
    		double max=0;
    		int row=k,num=0;
    		for(i=k;i<n;i++)
    		{ 
    			if(fabs(a[i][k])>max)
    			{ 
    				max=fabs(a[i][k]);
    				row=i;	//记录换行行标
    			}
    		}
    		if(a[row][i]==0)
    		{ 
    			cout<<"无法计算"<<endl; 
    		}
    		if(k!=row)//换行
    		{ 
    		 	for(i=0;i<m+1;i++)//该行全部元素交换
    		 	{ 
    		 		temp=a[k][i];
    		 		a[k][i]=a[row][i];
    		 		a[row][i]=temp;
    			 }
    		}
    		for(i=k+1;i<m;i++)
    		{ 
    			d=a[i][k]/a[k][k];	//定义解法
    			for(j=0;j<n+1;j++)
    			{ 
    				a[i][j]=a[i][j]-d*a[k][j];
    			}
    		}
    	} 
    	for(i=0;i<n;i++)
    	{ 
    		btemp[i]=0.0;
    	}
    	for(i=m-1;i>=0;i--)	//求解向量
    	{ 
    		d=0;
    		for(k=0;k<n;k++)
    		{ 
    			d=d+btemp[k]*a[i][k];
    		}
    		btemp[i]=(a[i][n]-d)/a[i][i]; 
    	} 
    	for(i=0;i<m;i++)
    	{ 
    		x[i]=btemp[i];
    	}
    	for(i=0;i<n;i++)
    	{ 
    		delete[]a[i];
    	}
    	delete[]a;
    	delete[]btemp;
    }
    int main()
    { 
    	LARGE_INTEGER lFrequency, lEndCount, lBeginCount;//LARGE_INTEGER为数据类型
    	QueryPerformanceFrequency(&lFrequency);//获取CPU的时钟频率
    	QueryPerformanceCounter(&lBeginCount);//获取CPU计数器数字，放在计时开头
    	double CompuTime;
    	double **A,*b,*x;
    	int kk=4;	
    	A=new double *[kk];
    	b=new double[kk];
    	x=new double[kk];
    	for(int i=0;i<kk;i++)
    	{ 
    		A[i]=new double[kk];
    	}
    	int Loop_num=1e5;
    	int n=4;	//根据题目选择赋值或输入
    	A[0][0]=1;A[0][1]=-1;A[0][2]=2;A[0][3]=-1;
    	A[1][0]=2;A[1][1]=-2;A[1][2]=3;A[1][3]=-3;
    	A[2][0]=1;A[2][1]=1;A[2][2]=1;A[2][3]=0;
    	A[3][0]=1;A[3][1]=-1;A[3][2]=4;A[3][3]=3;
    	b[0]=-8;b[1]=-20;b[2]=-2;b[3]=4;
    	SolverEqGuass(A,b,n,1.0e-10,x);
    	for(int i=0;i<n;i++)
    	{ 
    		cout<<"x"<<i+1<<"=";
    		printf("%.6f\n",x[i]);	//输出六位小数
    	}
    	QueryPerformanceCounter(&lEndCount);//获得CPU计数器数字，放在计时结尾
    	CompuTime = (double)(lEndCount.QuadPart - lBeginCount.QuadPart) / (double)lFrequency.QuadPart;//得到运行时间，单位微秒
    	cout<<"Gauss法计算时间为"<<CompuTime; 
    }

//封装函数二
    void SolverEqGauss(double** A, double* b, int n, double* x)
    { 
        double aef = 0, bt = 0, gm = 0, delta = 0, omiga = 0;
        int lamda = 0;
        bool bo = false;
        for (int I = 1;I < n;I++)
        { 
            delta = abs(A[I][I]);
            lamda = I;
            for (int k = I+1;k < n + 1;k++)
            { 
                if (abs(A[k][I]) > delta)
                { 
                    delta = abs(A[k][I]);
                    lamda = k;
                }
            }
            if (delta < pow(0.1, 8))
                bo = true;
            if (bo)
            { 
                cout << "此方程无解";
                break;
            }
            else
            { 
                for (int l = I;l < n + 1;l++)
                { 
                    omiga = A[I][l];
                    A[I][l] = A[lamda][l];
                    A[lamda][l] = omiga;
                }
                omiga = b[I];
                b[I] = b[lamda];
                b[lamda] = omiga;
            }
            for (int i = I+1;i < n + 1;i++)
            { 
                aef = A[i][I] / A[I][I];
    
                for (int j = I;j < n + 1;j++)
                { 
                    A[i][j] = A[i][j] - A[I][j]*aef;                
                }
                b[i] = b[i] - b[I]*aef;
            }
        }
        
        for (int u = n;u > 0;u--)
        { 
            for (int v = n;v > u ;v--)
            { 
                b[u] -= A[u][v] * x[v];
            }
            x[u] = b[u] / A[u][u];
        }
    }

//高斯列主元消去法（封装函数三）
    void SolverEqGauss(double **A,double *b,int n, double eps,double *x)
    { 
    	int i,j,k;
    	double **a = new double *[n];
    	double d,temp;
    	for(i=0;i<n;i++)
    	{ 
    		a[i]=new double[n+1];
    	}
    	//赋值
    	for(i=0;i<n;i++)
    	{ 
    		for(j=0;j<n;j++)
    		{ 
    			a[i][j]=A[i][j];
    		}
    		a[i][n]=b[i];
    	}
    	if(n<=0)
    	{ 
    		cout<<"输入格式有误！\n";
    	}
    	for(k=0;k<n-1;k++)
    	{ 
    		double max=a[k][k];
    		int row=0,num=0;
    		for(i=k+1;i<n;i++)
    		{ 
    			if(fabs(a[i][k]>max)
    			{ 
    				max=fabs(a[i][k]);
    				row=i;
    			}
    		}
    		if (fabs(a[row][k]) < eps) 
    		{  
    			cout << "主元素绝对值太小，无法计算" << endl; 
    			return; 
    		}
    		if (k != row) //换行 
    		{  
    			for (i = 0; i < n + 1; i++) 
    			{  
    				temp = a[k][i]; 
    				a[k][i] = a[row][i]; 
    				a[row][i] = temp; 
    			} 
    		}
    		for (i = k + 1; i < n; i++) //消元 
    		{  
    			d = a[i][k] / a[k][k];
    			for (j = 0; j < n + 1; j++) 
    			{  
    				a[i][j] = a[i][j] - d * a[k][j]; 
    			} 
    		} 
    	}
    	if (fabs(a[n - 1][n - 1]) < eps) 
    	{  
    		cout << "a[n - 1][n - 1]绝对值太小，无法计算" << endl; 
    		return; 
    	}
    	for (i = n - 1; i >= 0; i--) //回代 
    	{  
    		d = 0; 
    		for (k = i; k < n; k++)) 
    		{  
    			d = d + x[k] * a[i][k]; 
    		}
    		x[i] = (a[i][n] - d) / a[i][i]; 
    	}
    	for ( i = 0; i < n; i++) 
    	{  
    		delete[]a[i]; 
    	}
    	delete[]a; 
    }