update:对一处表述错误进行了修正,求过qwq
常规的思路:高斯消元
高斯消元
这是一种用来求解线性方程组的算法,在方程数较多的时候尤其省时。
对于一个 (n) 元线性方程组:
[egin{cases}
a_{1_1}x_1+a_{1_2}x_2+ cdots +a_{1_n}x_n=b_1\
a_{2_1}x_1+a_{2_2}x_2+ cdots +a_{2_n}x_n=b_2\
vdots\
a_{n_1}x_1+a_{n_2}x_2+ cdots +a_{n_n}x_n=b_n\
end{cases}
]
要解它,我们的第一反应显然是消元。
这里给出对于方程组中方程变换最常用的三个法则:
-
两个方程位置互换,解不变。
-
方程两边同时乘一个非零实数 (k) ,解不变
-
方程乘一非零实数 (k) ,再加上另一个方程,解不变
对应到系数矩阵中去,就是矩阵的初等行变换。
我们还知道,对于任意矩阵都能通过初等行变换转化为一个上三角矩阵
所以这个线性方程组的系数矩阵:
[egin{bmatrix}
a_{1_1}&a_{1_2}&dots&a_{1_n}\
a_{2_1}&a_{2_2}&dots&a_{2_n}\
vdots&vdots&&vdots&\
a_{n_1}&a_{n_2}&dots&a_{n_n}\
end{bmatrix}
]
可以被转化成上三角形式:
[egin{bmatrix}
a_{1_1}&a_{1_2}&dots&a_{1_n}\
&a_{2_2}&dots&a_{2_n}\
&&&vdots\
& & &a_{n_n}\
end{bmatrix}
]
对应到方程就是:
[egin{cases}
a_{1_1}x_1+a_{1_2}x_2+ cdots +a_{1_n}x_n=b_1\
a_{2_2}x_2+ cdots +a_{2_n}x_n=b_2\
vdots\
a_{n_n}x_n=b_n\
end{cases}
]
显然最后方程能被直接解出来,解出的根又可以向后代入上一个方程里面。
具体步骤:
枚举每一列c:(相当于把 (x_c) 作为当前处理的列主元)
- 找到绝对值最大的一行,将这一行与顶行交换位置,这行不再移动。
这一步是为了将绝对值最大的 (a_{r,c}) 做除数,减少计算时的误差
- 等式两边同除,将该行第一个数变为 (1)。
系数化一,方便后面求解
- 通过等式之间加减,将当前列下面的所有数消为(0)。
要构造上三角矩阵,这是必要
- 重复至最后一列,然后回代求解。
如果求解过程中发现有方程出现 ( ext{非零常数}=0) 的情况,则方程组无解。
出现了(0=0)的情况,则有无数多组解。
code:
#include <bits/stdc++.h>
using namespace std;
const int N=110;
const double eps=1e-6;//double精度
int n;
double a[N][N];
int gauss()
{
int c,r;//c表示在枚举哪一列,r表示在枚举哪一行
for(c=0,r=0;c<n;c++)
{
int t=r;
for(int i=r;i<n;i++)
{
if(fabs(a[i][c])>fabs(a[t][c]))//定位到当前列主元的最大系数所在行
t=i;
}
if(fabs(a[t][c])<eps) continue;
for(int i=c;i<=n;i++) swap(a[t][i],a[r][i]);//交换这两行
for(int i=n;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;j>=c;j--)
a[i][j]-=a[r][j]*a[i][c];//相加减,消去列下面的系数
}
}
r++;
}
for(int i=n-1;i>=0;i--)//向后回代求解
{
for(int j=i+1;j<n;j++)
{
a[i][n]-=a[i][j]*a[j][n];
}
}
if(r<n)
{
for(int i=r;i<n;i++)
{
if(fabs(a[i][n])>eps)
return 2;//无解
}
return 1;//有无穷多组解
}
return 0;//有唯一解
}
int main()
{
cin>>n;
for(int i=0;i<n;i++)
{
for(int j=0;j<n+1;j++)
{
cin>>a[i][j];
}
}
int t=gauss();
if(t==0)//有唯一解
{
for(int i=0;i<n;i++) printf("x%d=%.2lf
",i+1,a[i][n]);
}
else if(t==1)puts("0");
else puts("-1");
return 0;
}