$Gauss$消元

$Gauss$消元

  今天金牌爷来问我一个高消的题目,我才想起来忘了学高消...

  高斯消元用于解线性方程组,也就是形如:

  $left{egin{matrix}a_{11}x_1+a_{12}x_2+...+a_{1n}x_n=b_1\ a_{21}x_1+a_{22}x_2+...+a_{2n}x_n=b_2\a_{31}x_1+a_{32}x_2+...+a_{3n}x_n=b_3end{matrix} ight.$

  好像也可以写成这样:

  $AX=B$

  其实就是小学学的加减消元...

  举个栗子:

  $left{egin{matrix}3x_1+2x_2=5\ 2x_1+3x_2=10\end{matrix} ight.$

  首先从第一列开始,找到第一项系数的绝对值最大的一行放到第一行,把这个系数除去(系数化为$1$):

  $left{egin{matrix}x_1+frac{2}{3}x_2=frac{5}{3}\ 2x_1+3x_2=10\end{matrix} ight.$

  发现第二行减去两个第一行就可以消掉第一个未知数,那么就减掉两个好咯:

  $left{egin{matrix}x_1+frac{2}{3}x_2=frac{5}{3}\ quad space space frac{5}{3}x_2=frac{20}{3}\end{matrix} ight.$

  再消掉第二行的系数:

  $left{egin{matrix}x_1+frac{2}{3}x_2=frac{5}{3}\ quad space space x_2=4\end{matrix} ight.$

  现在就解出了第二个未知数,再从底往上代回去:

  $left{egin{matrix}x_1+frac{2}{3} imes 4=frac{5}{3}\ quad space space x_2=4\end{matrix} ight.$

  依次解出所有的未知数即可:

  $left{egin{matrix}x_1=-1\ x_2=4\end{matrix} ight.$

  为什么要将系数绝对值最大的一项作为主元进行消元呢?因为实际做题中用的不是分数而是浮点数,有误差的问题,如果用于消元的主元太接近零,那么下一行就需要减掉非常多倍的上一行,导致精度大量损失.

  怎样判断无解:如果消元结束后发现有一行的系数全为$0$,但是此行的$b$不为$0$,那么未知数取任何值都不满足要求;

  怎样判断无穷组解:如果消元结束后有一行系数全为$0$,$b$也是$0$,那么可以随意取值.

  也就是平时数学中做题可能会碰到的“两个方程左边的本质是相同的,但是答案却不相同,就无解,如果答案相同,说明有$n$个未知数,限制条件却少于$n$个,此时有多组解”。

  注意一点:如果有的方程无解,有的有多组解,总方程组还是无解,所以先判断无解再判断多解.

  

  线性方程组:https://www.luogu.org/problemnew/show/P2455

  题意概述:完全的模板题.

  
 1 # include <cstdio>
 2 # include <iostream>
 3 # include <queue>
 4 # include <cstring>
 5 # include <cmath>
 6 # include <string>
 7 # define R register int
 8 # define ll long long
 9 
10 using namespace std;
11 
12 const double eps=0.0000001;
13 const int maxn=101;
14 int n,cnt,None,Endless;
15 bool fre[maxn];
16 double A[maxn][maxn],x,ans[maxn];
17 
18 void Gauss()
19 {
20     for (R i=1;i<=n;++i)
21     {
22         int id=i;
23         for (R j=i+1;j<=n;++j)
24             if(fabs(A[id][i])<fabs(A[j][i])) id=j;
25         for (R j=1;j<=n+1;++j)
26             swap(A[i][j],A[id][j]);
27         if(fabs(A[i][i])<eps) continue;
28         x=A[i][i];
29         for (R j=1;j<=n+1;++j) A[i][j]/=x;
30         for (R j=1;j<=n;++j)
31         {
32             if(i==j) continue;
33             x=A[j][i];
34             for (R k=1;k<=n+1;++k)
35                 A[j][k]-=x*A[i][k];
36         }
37     }
38 }
39 
40 int main()
41 {
42     scanf("%d",&n);
43     for (R i=1;i<=n;++i)
44         for (R j=1;j<=n+1;++j)
45             scanf("%lf",&A[i][j]);
46     Gauss();
47     for (R i=1;i<=n;++i)
48     {
49         R j=1;
50         while (fabs(A[i][j])<eps&&j<=n+1) j++;
51         if(j>n+1) Endless=1;
52         else if(j==n+1) None=1; 
53     }
54     if(None) { printf("-1"); return 0; }
55     if(Endless) { printf("0"); return 0; }
56     for (R i=n;i>=1;--i)
57     {
58         ans[i]=A[i][n+1];
59         for (R j=i-1;j>=1;--j)
60         {
61             A[j][n+1]-=ans[i]*A[j][i];
62             A[j][i]=0;
63         }
64     }
65     for (R i=1;i<=n;++i)
66         printf("x%d=%.2lf
",i,ans[i]);
67     return 0;
68 }
线性方程组

   

  游走:https://www.lydsy.com/JudgeOnline/problem.php?id=3143

  题意概述:给定一张$n$个点$m$条边的无向简单连通图,从一号点出发,每次从这个点发出的所有边中随机选择一条走过去,到达$n$之后就不再走了,要求给每一条边赋一个独特的,$[1,m]$的权值,使得整条路径上期望的权值总和最小.

  为什么要做这道题呢?$shzr:$我学了高斯消元;$asuldb$:高斯消元有什么用啊,又不能做题;$shzr$:点开"线性方程组";$asuldb$:那有什么用啊,你除了会做模板题还是什么都不会啊,你会用高斯消元做期望吗?$shzr$:...

  也许这题应该放到期望的标签下? 首先根据期望的线性可加性,权值和的期望等于权值的期望和,所以可以对于每一条边分别计算贡献.每一条边的贡献就是经过这一条边的期望次数乘上它的权值;这时贪心策略就很明显了,首先求出每条边的期望经过次数,按照这个次数进行排序,出现次数多的优先赋值小权值即可.那么怎样计算一条边的贡献呢?一条边有两个端点,经过它必然是从某一个端点走过来的,假设我们已经计算出了每个点的期望经过次数记作$E_i$,每个点的度数记为$d_i$,那么对于一条端点为$x,y$的边,它的期望经过次数就是$frac{E_x}{d_x}+frac{E_y}{d_y}$.

  如何计算每个点的期望经过次数?枚举每一个与它有边相连的点,它的经过次数就是$E_u=sum_{v->u}frac{E_v}{d_v}$,到这里问题就完美解决了。

  然而你可能突然发现这道题是无向图,事情并没有那么简单·_· 每个点之间的相互关系不仅无法进行拓扑排序而且错综复杂,问题进行到这里我们好像陷入了知识盲区。但是仔细整理思路可以发现虽然每个点互相关联,但是却正好组成了一个$N$元一次方程组,$N$个方程,于是可以使用高斯消元.因为走到第$n$个点就不能再走了,这个点不能计算贡献,所以可以直接在矩阵中删去这个点.因为第一个点是开始点,所以除了从别的点到这里的贡献之外还另外有一个$1$,不要忘了.初始化矩阵时要注意:第$x$行的方程是用于计算第$x$个未知数的,第$y$列的系数是对于第$y$个未知数给出去的贡献的,不要写反了.

1 for (R i=1;i<=m;++i)
2     {
3         if(x[i]==n||y[i]==n) continue;
4         a[ x[i] ][ y[i] ]-=1.0/d[ y[i] ]; //注意这里
5         a[ y[i] ][ x[i] ]-=1.0/d[ x[i] ]; //
6     }
  
 1 # include <cstdio>
 2 # include <iostream>
 3 # include <algorithm>
 4 # define R register int
 5 
 6 using namespace std;
 7 
 8 const int maxn=502;
 9 int u,v,n,m,h,firs[maxn],d[maxn];
10 double a[maxn][maxn],ans[maxn];
11 int x[maxn*maxn],y[maxn*maxn];
12 double co[maxn*maxn];
13 
14 double ab (double a) { if(a<0) return -a; return a; }
15 
16 int main()
17 {
18     scanf("%d%d",&n,&m);
19     for (R i=1;i<=m;++i)
20     {
21         scanf("%d%d",&u,&v);
22         x[++h]=u;
23         y[h]=v;
24         d[u]++,d[v]++;
25     }
26     for (R i=1;i<=m;++i)
27     {
28         if(x[i]==n||y[i]==n) continue;
29         a[ x[i] ][ y[i] ]-=1.0/d[ y[i] ];
30         a[ y[i] ][ x[i] ]-=1.0/d[ x[i] ];
31     }
32     for (R i=1;i<n;++i)
33         a[i][i]=1.0;
34     a[1][n]=1.0;
35     for (R i=1;i<n;++i)
36     {
37         int maxx=i;
38         for (R j=i+1;j<n;++j)
39             if(ab(a[maxx][i])<ab(a[j][i])) maxx=j;
40         swap(a[i],a[maxx]);
41         double x;
42         for (R j=i+1;j<n;++j)
43         {
44             x=a[j][i]/a[i][i];
45             for (R k=i;k<=n;++k)
46                 a[j][k]-=a[i][k]*x;
47         }
48     }
49     for (R i=n-1;i>=1;--i)
50     {
51         for (R j=i+1;j<n;++j) a[i][n]-=a[i][j]*ans[j];
52         ans[i]=a[i][n]/a[i][i];
53     }
54     for (R i=1;i<=m;++i)
55     {
56         co[i]+=ans[ x[i] ]/d[ x[i] ];
57         co[i]+=ans[ y[i] ]/d[ y[i] ];
58     }
59     double fans=0;
60     sort(co+1,co+1+m);
61     for (R i=1;i<=m;++i)
62         fans+=co[i]*(m-i+1);
63     printf("%.3lf",fans);
64     return 0;
65 }
游走

---shzr

原文地址:https://www.cnblogs.com/shzr/p/9845799.html