碼上快樂

相關內容    簡體    繁體

高斯消元學習筆記及算法實現與運用

本文轉載自   查看原文   2021-11-11 15:28   182 

高斯消元學習筆記及算法實現與運用

目錄

0.前言

上大學了，好久沒寫博客了，剛好近期學習了高斯消元，打算在這篇文章中將知識點進行整理， 同時對以前的OI知識進行一個復習。

1.高斯消元

階梯形線性方程組

\[\left\{ \begin{aligned} a_ {1i_{1}}x_{i_{1}}+……+a_ {1i_{2}}x_{i_{2}}+……+a_ {1r_{1}}x_{1_{r}}+……a_ {1n}x_{n}=b_1\\ a_ {1i_{2}}x_{i_{2}}+……+a_ {2r_{1}}x_{2_{r}}+……a_ {2n}x_{n}=b_2 \\ …… \\ a_ {ri_{r}}x_{i_{1}}+……+ a_ {rn}x_{n}=b_r\\……\\0=b_{r+1}\\0=0\\……\\0=0\end{aligned} \right. \]

其中1=i1<i2<……<ir\(\leq\)n, \(a_{ki_k}\neq0,k=1,2,……,r\),則稱該線性方程組是階梯形線性方程組。

Gauss消元法的基本思想可以描述為在不改變線性方程組解的前提下，將一般的線性方程組化為”階梯形“線性方程組求解。其過程可以分為兩步

  • 消元過程
  • 即將原線性方程組等價地化為接替形線性方程組；
  • 回代過程
  • 即從階梯形線性方程組求得原線性方程組的解

線性方程組的解有3種情況

  • 有唯一解
  • 無解
  • 有無窮解

線性方程組的初等變換（同解變換）

  1. 將線性方程組中某兩個方程互換
  2. 將線性方程組中某個方程兩邊同時乘一個非零常數
  3. 將線性方程組中某個方程的k倍與另一個方程相加

兩個定理

  • 齊次線性方程組中，若方程個數s小於未知量個數n，則其必有非零解。
  • n元齊次線性方程組，只有零解的充分必要條件為線性方程組系數行列式不為零。

這里對第二個定理做證明

證明： 若線性方程組的系數行列式不為零，則由Cramer法則可知，線性方程組只有零解。

若線性方程組只有零解，則由Gauss消元法可知，消元過程結束時的階梯形線性方程組為

image

注意到階梯形線性方程組是由線性方程組經過一系列初等變換所得，並且對方程組做初等變換時，發生變化的只是相應的方程中未知量前面的系數所以行列式D1可以視為由行列式D經過一系列運算化為上三角形行列式而得，因此兩者只差一個非零因子。由於D1\(\neq0\)，故D\(\neq0\)，即線性方程組的系數行列式不為0.

階梯形矩陣

若一個矩陣的零行（如果有的話）均在非零行的下方，並且每個非零行左起第一個非零元素所在的下標隨着行標的增大而嚴格增大，則稱次矩陣為階梯形矩陣。如下：

\[\begin{bmatrix} 2&1&-1&8\\ 0&-1&2&11\\ 0&0&2&-1\\0&0&0&3\\0&0&0&0 \end{bmatrix}\]

這里給出一種更通俗的解釋

\[\left\{ \begin{aligned} 非零行首個非零元素下方都是0\\ 零行在最下方 \end{aligned} \right. \]

2.算法實現

算法分析

如果給定一個形如以下式子得多元方程式

\[\left\{ \begin{aligned} 2x+y-z & = & 8 \\ -3x-y+2z& = &11 \\ -2x+y+2z & = & -1 \end{aligned} \right. \]

我們首先提出各項的系數（因為我們知道，高斯消元其實只跟系數有關）

我們可以寫成一下矩陣形式

其中左邊是各項的系數，最右邊一列是等式右邊的常數列

\[\begin{bmatrix} 2&1&-1&8\\ -3&-1&2&11\\ -2&1&2&-1\\ \end{bmatrix}\]

接下來我們進行加減消元

首先我們要用第i個方程來消去第k個方程的第k列，那么第k行所有元素A [ k ] [ j ]都應減去A [ i ] [ j ]的A [ k ] [ i ]/A [ i ] [ i ]倍。（我們約定，A [ i ] [ j ]是第i行第j列的系數）

上面的矩陣經過加減消元得

\[\begin{bmatrix} -3&-1&2&-11\\ 0&1/3&1/3&2/3\\ 0&5/3&2/3&13/3\\ \end{bmatrix}\]

接下來就是不斷重復這個過程，再代入。

根據我們的算法，可以直接從最后一行推出Xn的值，然后倒數第二行可以推出Xn-1的值，以此類推。最后我們可以求出所有X的唯一解。

各部分代碼詳解

1.經過r行和第i行交換和加減消元

for(int i = 0; i < n; i ++) {
        r = i;
        for(int j = i + 1; j < n; j ++)
            if(fabs(A[j][i]) > fabs(A[r][i])) r = j; //fabs()表示計算浮點數的絕對值
        if(r != i) for(int j = 0; j <= n; j ++) std :: swap(A[r][j], A[i][j]);
       				/*↑↑      經過r行和第i行交換      ↑↑*/
                    
        for(int k = i + 1; k < n; k ++) {
			double f = A[k][i] / A[i][i];
			for(int j = i; j <= n; j ++) A[k][j] -= f * A[i][j];
		}            
        				/*↑↑     加減消元（低精度）    ↑↑*/
        
        for(int j = n; j >= i; j --) {
            for(int k = i + 1; k < n; k ++)
                A[k][j] -= A[k][i] / A[i][i] * A[i][j];
        }
        				/*↑↑     加減消元（保證精度）    ↑↑*/
    }

2.回代過程

	for(int i = n - 1; i >= 0; i --) {
		for(int j = i + 1; j < n; j ++)
			A[i][n] -= A[j][n] * A[i][j];
		A[i][n] /= A[i][i];
	}

總代碼

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define ll long long
using namespace std;

int n,m;
double a[110][110];


int main()
{
	scanf("%d",&n);
	m=n;
	for(int i=1;i<=m;i++)
	{
		for(int j=1;j<=n+1;j++)
		{
			scanf("%lf",&a[i][j]);
		}
	}
	int flag=0;
	for(int i=1;i<=n;i++)
	{
		int temp=i;
		while(a[temp][i]==0&&temp<=n)
		{
			temp++; 
		}
		if(temp==n+1)
		{
			flag=1;
			continue;
		}
		for(int j=1;j<=n+1;j++)
		{
			swap(a[i][j],a[temp][j]);
		}
		double k=a[i][i];
		for(int j=1;j<=n+1;j++)
		{
			a[i][j]/=k;
		}
		for(int j=1;j<=m;j++)
		{
			if(i!=j)
			{
				double k=a[j][i];
				for(int l=1;l<=n+1;l++)
				{
					a[j][l]-=(k*a[i][l]);
				}
			}
		}
	}
	if(flag)
	{
		printf("No Solution\n");
		return 0;
	}
	for(int j=1;j<=m;j++)
	{
		printf("%.2lf\n",a[j][n+1]/a[j][j]);
	}
	return 0;
}

3.高斯消元的實例運用

球形空間產生器（BZOJ1013）

有一個球形空間產生器能夠在 n 維空間中產生一個堅硬的球體。現在，你被困在了這個 n 維球體中，你只知道球面上 n+1 個點的坐標，你需要以最快的速度確定這個 n 維球體的球心坐標，以便於摧毀這個球形空間產生器。數據保證有解，答案精確到小數點后三位。

一個球體上的所有點到球心的距離相等，因此只需求出一個點（x1，x2，……，xn），使得：

\[\sum_{j=0}^n(a_i,_j-x_j)^2=C \]

其中C為常數，i屬於[1,n+1],球面上第i個點的坐標為（ai，1，ai，2，……，ai，n）。該方程組由n+1個n元二次方程構成，不是線性方程組。但是我們可以通過相鄰兩個方程做差，把它變成n個n元一次方程，同時消去常數C：

\[\sum_{j=1}^n(a_i,_j^2-a_i+_1,_j^2-2x_j(a_i,_j-a_i+_1,_j))=0 (i=1,2,…，n) \]

把變量放左邊，把常數放右邊：

\[\sum_{j=1}^n(a_i,_j^2-a_i+_1,_j^2)=\sum_{j=1}^n2x_j(a_i,_j-a_i+_1,_j) (i=1,2,…，n) \]

這就是一個線性方程組了。題目保證方程組有唯一解，我們直接對下面的增廣矩陣進行高斯消元，變為簡化階梯形矩陣，即可得到每個xj的值。

\[\begin{bmatrix} 2（a_1,_1-a_2,_1）&2（a_1,_2-a_2,_2）&……&2（a_1,_n-a_2,_n）&\sum_{j=1}^n(a_1,_j^2-a_2,_j^2)\\ 2（a_2,_1-a_3,_1）&2（a_2,_2-a_3,_2）&……&2（a_1,_n-a_2,_n）&\sum_{j=1}^n(a_1,_j^2-a_2,_j^2)\\ ……&……&……&……&……\\2（a_n,_1-a_n+_1,_1）&2（a_n,_2-a_n+_1,_2）&……&2（a_n,_n-a_n+_1,_n）&\sum_{j=1}^n(a_n,_j^2-a_n+_1,_j^2)\\\end{bmatrix}\] 

#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define ll long long
using namespace std;

int n;
double a[20][20],b[20],c[20][20];



int main()
{
	scanf("%d",&n);
	for(int i=1;i<=n+1;i++)
	{
		for(int j=1;j<=n;j++)
		{
			scanf("%lf",&a[i][j]);
		}
	 } 
	 for(int i=1;i<=n;i++)
	 {
	 	for(int j=1;j<=n;j++)
	 	{
	 		c[i][j]=2*(a[i][j]-a[i+1][j]);
	 		b[i]+=a[i][j]*a[i][j]-a[i+1][j]*a[i+1][j];
		 }
	 }
	 for(int i=1;i<=n;i++)
	 {
	 	for(int j=i;j<=n;j++)
	 	{
	 		if(fabs(c[i][j])>1e-8)
	 		{
	 			for(int k=1;k<=n;k++)
	 			{
	 				swap(c[i][k],c[j][k]);
				 }
				swap(b[i],b[j]);
			 }
		 }
		 for(int j=1;j<=n;j++)
		 {
		 	if(i==j) continue;
		 	double t=c[j][i]/c[i][i];
		 	for(int k=i;k<=n;k++)
		 	{
		 		c[j][k]-=t*c[i][k];
			}
			b[j]-=b[i]*t;
		 }
	 }
	 for(int i=1;i<=n;i++)
	 {
	 	printf("%.3lf ",b[i]/c[i][i]);
	 }
	return 0;
}

Barracuda（luoguP5027）

題目描述

大三角形給小正方形講起自己的過去：過去的它是一個挖寶工，后來被黑暗之主污染才會落到此番境地。

它也希望小正方形去戰勝黑暗之主，不過限於黑暗之主的眼線密布，因此必須給小正方形設置障礙才能騙過那些“眼線”。

他給小正方形的問題是：它有 n 個小三角形，每個小三角形有一定的質量，它對這些三角形進行了 n+1 次稱量，然而由於托盤天平(?)的問題，有一次稱量的結果是有誤的。

現在，大三角形想要知道最重的小三角形的 編號。

一組輸入是合法的，當且僅當輸入滿足以下條件：

不存在一組 i,j，使得當我們假定第 i 條稱量數據有誤時能求出一種合法方案且我們假定第 j 條稱量數據有誤時也能求出一種合法方案。

合法方案定義如下：

1、最重的三角形只有一個。

2、不存在重量不確定的三角形。

3、所有三角形的重量均為正整數。

題目還是很好理解的，每次給定幾個三角形和重量，就可以建立n+1個方程，轉化為高斯消元模板

想一想高斯消元的模板長什么樣子，一個n∗(n+1)的矩陣，但是按照上面的思路來建立的話，會搞出來一個(n+1)∗(n+1)的矩陣，那么就會遇到一些坑

思路就是假設每一行都是錯誤的情況

  1. 並不是每一個編號的系數為1，最開始我這里想了很久，可以抓多次同個三角形
  2. 如何假設是x*行出錯，對於高斯消元時的處理，並不能直接continue當前行，還是應該存儲一個模板矩陣

然后就繼續說如何判斷題目中的條件

  1. 最重的三角形只有一個，我們只需要在回代過后，對max x的個數記錄一下即可
  2. 不存在重量不確定的三角形，說明該方程在轉化為“上三角形”之后，出現了0=0這類恆等式，表示無數組解，當然無解也是這樣
  3. 三角形重量是正整數，我們用double來存儲答案，如果int強制轉換后的值比double類型小，說明不是正整數
  4. 對於不是合法方案，其實就是相當於當i錯誤時，我們能求出來一組解,當jj錯誤時，我們也能求出來一組解，說明illegal
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
using namespace std;
const double eps=1e-6;
double a[150][150];
double b[150][150];
int n,m,x;
int maxx,pre;
int gauss() 
{
	int c,r;
	for(c=1,r=1; c<=n; c++) 
    {
		int t=r;
		for(int i=r; i<=n; i++) 
        {
			if(fabs(a[t][c])<fabs(a[i][c])) t=i;
		}
		if(fabs(a[t][c])<eps) return 0;
		for(int i=c; i<=n+1; i++) swap(a[r][i],a[t][i]);
		for(int i=n+1; i>=c; i--) a[r][i]/=a[r][c];
		for(int i=r+1; i<=n; i++) 
        {
			if(fabs(a[i][c])>eps) 
            {
				for(int j=n+1; j>=c; j--) a[i][j]-=a[r][j]*a[i][c];
			}
		}
		r++;
	} //高斯消元模板 
	for(int i=n; i>=1; i--) 
    {
		for(int j=i+1; j<=n; j++) 
        {
			a[i][n+1]-=a[j][n+1]*a[i][j];
		}
		if(a[i][n+1]<eps||((int)a[i][n+1])<a[i][n+1]) return 0; //無解或是重量為小數的情況 
	} //回代 
	int sum;
	maxx=-1;
	for(int i=1; i<=n; i++) 
	{
		if((int)a[i][n+1]>maxx) maxx=a[i][n+1],sum=i;
	} //記錄最值 
	int num=0;
	for(int i=1; i<=n; i++)
	{
		if((int)a[i][n+1]==maxx) num++; //最值數量 
	}
	if(num>1) return 0;
	return sum;
}
int tot;
int main() 
{
	scanf("%d",&n);
	for(int i=1; i<=n+1; i++) 
	{
		scanf("%d",&m);
		for(int j=1; j<=m; j++) 
		{
			scanf("%d",&x);
			b[i][x]++;
		}
		scanf("%d",&x);
		b[i][n+1]=x*1.0;
	} //存儲初始的(n+1)*(n+1)的矩陣 
	for(int i=1; i<=n+1; i++) 
	{
		int cnt=0;
		for(int h=1; h<=n+1; h++) 
		{
			if(h==i) continue;
			cnt++;
			for(int k=1; k<=n+1; k++) 
			{
				a[cnt][k]=b[h][k];
			}
		} //存儲一個n*(n+1)的矩陣 
		int p=gauss();
		if(p==0) continue; //無解或者無數組解 
		else tot++,pre=p; //tot表示有幾組解，pre是答案 
	}
	if(tot>1||tot==0) puts("illegal"); //全部無解或者出現多組解 
	else cout<<pre;  
	return 0;
}

4.如何理解Gauss消元法解線性方程組的正確性

線性方程組Ax=b求解的Gauss消元法就是將方程組的增廣矩陣 \(\overline{A}\)=(A,b)作初等行變換化為階梯形\(\overline{C}\)=(C,b)，然后再求解方程Cx=d，從而得到原方程組解的方法。在這個過程中我們既沒驗證方程組Cx=d的解就是Ax=b的解，又沒證明求解過程中方程組的解會不會丟失。換句話說，為什么說方程組Cx=d和Ax=b一定同解？這個問題用方程組的向量形式來描述就會看得清楚。

首先將方程組Ax=b寫成向量形式，記A=（\(\alpha\)1,\(\alpha\)2,……,\(\alpha\)n）,其中\(\alpha\)i=（a1i,a2i,……,ami）T (i=1,2,……，n),則方程組可寫成

image

容易看出，對\(\overline{A}\)=(A,b)作初等變換等價於對方程組中的方程做相應的運算，這些運算不會改變上式中各列向量的線性關系。也就是各列向量的系數x1，x2,……,xn與初等變換無關。另一方面，初等變換是可逆的，從而方程組Ax=b與Cx=b在初等行變換下可以互化，所以它們是同解的。

5.后記

這篇文章只是本人對於高斯消元的一些小小的思考總結，很多知識我也只是”知其然而不知其所以然“。我也再次意識到，在前人的偉大發現和思考前，我是多么渺小。如果本文有什么錯誤請正在閱讀的您指出並告知，便於本人修改，也希望您能願意和我交流學習相關知識。同時萬分感激您能讀到這里。

PS.

十分感謝我的同學，徐同學，願意將她的筆記給我作為參考，她的幫助是我能完成此文的重要因素。


免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。



猜您在找 高斯消元算法 高斯消元 & 線性基【學習筆記】 「筆記」高斯消元 順序高斯消元法（Python實現） 淺談高斯消元的實現和簡單應用 高斯消元法 高斯消元法 高斯消元 高斯消元 高斯消元總結
 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM