Gauess列主元消去法C++写法

mac2022-06-30  18

高斯列主元消去法C++代码

一、思想原理:

算法:

第1步  输入系数矩阵A,右端项b,置k:=1;

第2步  对k=1,…,n-1进行如下操作:  

(1)选列主元,确定第一行第一个元素,保证它在第一列中绝对值是最大的,即为主元素;

(2)若第二行第二例的元素比第三行第二例的绝对值小,则交换的这两行;以此类推每一行的,直至保证每一行的对应位置的元素皆为绝对值最大的元素即可;

(3)消元,将矩阵运算化为上三角形即可

第3步  回代

将线性方程组做成增广矩阵,对增广矩阵进行行变换。

在第列中,第行及以下的元素选取绝对值最大的元素,将该元素最大的行与第行交换,然后采用高斯消元法将新得到的消去第行以下的元素,一次进行直到。从而得到上三角矩阵。再对得到的上三角矩阵进行回代操作,即可以得到方程组的解。

二、C++程序代码

#include<iostream> #include<cmath> using namespace std; #define N 10                //定义矩阵的维数,可自定义 static double A[N][N];       //系数矩阵 static double b[N];             //右端项 static double Aug[N][N + 1]; //增广矩阵 //高斯列主元消去法解线性方程组 void gaussin_LZY(double A[N][N],double b[N],double x[N],int n) {     int i,i1, j, k;     double Aug[N][N + 1], maxele, Temp, M;     for (i = 0; i < n; i++) {            //构造增广矩阵Aug         for (j = 0; j < n; j++)             Aug[i][j] = A[i][j];             Aug[i][n] = b[i];     }     for (i1 = 0; i1 < n - 1; i1++) {  //列主元         maxele = fabs(Aug[i1][i1]);  //取主元,fabs是取绝对值函数         k = i1;         for (i = i1; i < n; i++)             if (maxele < fabs(Aug[i1][i]))//比较i1列                 k = i1;         //换行         for (j = 0; j < n + 1; j++)         {             Temp = Aug[i1][j];             Aug[i1][j] = Aug[k][j];             Aug[k][j] = Temp;         }         //消元,以第i1行为工具行处理下一行的元素         for (k = i1 + 1; k < n; k++)         {             M = -Aug[k][i1] / Aug[i1][i1];             for (j = i1; j < n + 1; j++)                 Aug[k][j] = Aug[k][j] + M * Aug[i1][j];         }     }     //回代求解     double s;     x[n - 1] = Aug[n - 1][n] / Aug[n - 1][n - 1];     for (i = n - 2; i >= 0; i--)     {         s = 0;         for (j = i + 1; j < n; j++)             s = s + Aug[i][j] * x[j];         x[i] = (Aug[i][n] - s) / Aug[i][i];     }

三、实例验证结果

//Aug=[1 1 1 1 10;-1 2 -3 1 -2;3 -3 6 -2 7;-4 5 2 -3 0] //x[i]={1 2 3 4}

 

四、MATLAB程序验证

function x=magauss2(A,b,flag) if nargin<3,flag=0;end n=length(b); for k=1:(n-1) %选主元 [ap,p]=max(abs(A(k:n,k))); p=p+k-1; if p>k t=A(k,:);A(k,:)=A(p,:);A(p,:)=t; t=b(k);b(k)=b(p);b(p)=t; end m=A(k+1:n,k)/A(k,k); %消元 A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n); b(k+1:n)=b(k+1:n)-m*b(k); A(k+1:n,k)=zeros(n-k,1); if flag~=0,Ab=[A,b],end end x=zeros(n,1); %回代 x(n)=b(n)/A(n,n); for k=n-1:-1:1 x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k); end

验证结果:结果一致,代码正确

最新回复(0)