一.线性方程组的概念.
线性方程组,即方程中每一项次数都为一次的方程组.
回忆初中学的二元一次方程组,就是一种特殊的线性方程组.而线性方程组的解法与二元一次方程组的解法类似,通常只有代入法和加减消元法.
但是计算机的求解过程讲究程序化,所以我们要有一种特定的求解方式,于是高斯消元应运而生.
二.高斯消元.
高斯消元的过程就是在模拟加减消元得到一个未知量的解后,通过把这个解回代得到方程组的解.
具体来说,若我们现在有一个由n个未知量和n个方程组成的方程组,表现形式如下:
那么高斯消元会将这个方程组对应成一个
的矩阵,其中前
列是方程中每一项的系数,第
列是方程中的常数项,具体为:
得到了这么一个矩阵后,高斯消元就要把这个矩阵消成一个上三角矩阵,即对角线以下的所有值均为0的矩阵.
具体来说,我们模拟加减消元,第i步把第i行以下的所有矩阵的第i列的元素删完(消成0),然后我们就会得到一个上三角矩阵.实现时,通常会先把第i行的第i项变为1,然后用第行去消除它下面的每一行的第i项的系数.
之后我们考虑把考虑把从第n行开始到第2行为止的所有值回代到上一行中得到一个未知量的解即可.
分析一下时间复杂度,发现求上三角矩阵时间复杂度为 ,回代的过程是 的,所以总时间复杂度为 .
这里有一个避免精度误差的小技巧,就是求上三角矩阵时,每次消掉第i项系数绝对值最大的行,感性理解一下这可以尽量避免精度误差.
代码如下:
void gauss(matrix a){
int r,n=a.n;
double d;
for (int i=1;i<=n;++i){
r=i;
for (int j=i+1;j<=n;++j)
if (fabs(a.a[r][i])<fabs(a.a[j][i])) r=j;
if (i^r)
for (int j=1;j<=n+1;++j)
swap(a.a[i][j],a.a[r][j]);
d=a.a[i][i];
for (int j=i;j<=n+1;++j)
a.a[i][j]/=d;
for (int j=i+1;j<=n;++j){
d=a.a[j][i];
for (int k=i;k<=n+1;++k)
a.a[j][k]-=a.a[i][k]*d;
}
}
ans[n]=a.a[n][n+1];
for (int i=n-1;i>=1;--i){
ans[i]=a.a[i][n+1];
for (int j=i+1;j<=n;++j)
ans[i]-=a.a[i][j]*ans[j];
}
}
三.线性方程组解的情况判定.
我们发现这个线性方程组的解不一定有唯一解,还有无解与无穷解的情况,为了判定解的情况,我们可以讨论一下.
首先在高斯消元过程中,若出现了方程
这类情况(
为常数),方程无组解.然后,在高斯消元结束后,若出现了系数均为0,常数不为0的行,说明方程组无解;若出现了系数不全为0的行,则说明方程有无穷解;进一步的,若系数不为0的行共有
行,则最少需要确定
个未知量的值,才能确定所有解.
四.例题与代码.
题目:luogu3389.
代码如下:
#include<bits/stdc++.h>
using namespace std;
#define Abigail inline void
typedef long long LL;
const double eps=0.0000001;
const int N=100;
struct matrix{
double a[N+9][N+9];
int n;
}a;
double ans[N+9];
int gauss(matrix a){ //唯一解返回0,其它情况返回1,答案存在ans数组中
int r,n=a.n;
double d;
for (int i=1;i<=n;++i){
r=i;
for (int j=i+1;j<=n;++j)
if (fabs(a.a[r][i])<fabs(a.a[j][i])) r=j;
if (fabs(a.a[r][i])<eps) return 1;
if (i^r)
for (int j=1;j<=n+1;++j)
swap(a.a[i][j],a.a[r][j]);
d=a.a[i][i];
for (int j=i;j<=n+1;++j)
a.a[i][j]/=d;
for (int j=i+1;j<=n;++j){
d=a.a[j][i];
for (int k=i;k<=n+1;++k)
a.a[j][k]-=a.a[i][k]*d;
}
}
ans[n]=a.a[n][n+1];
for (int i=n-1;i>=1;--i){
ans[i]=a.a[i][n+1];
for (int j=i+1;j<=n;++j)
ans[i]-=a.a[i][j]*ans[j];
}
return 0;
}
Abigail into(){
scanf("%d",&a.n);
for (int i=1;i<=a.n;++i)
for (int j=1;j<=a.n+1;++j)
scanf("%lf",&a.a[i][j]);
}
Abigail outo(){
if (gauss(a)) puts("No Solution");
else
for (int i=1;i<=a.n;++i)
printf("%.2lf\n",ans[i]);
}
int main(){
into();
outo();
return 0;
}
题目:BZOJ1013.
题目大意:给定在n维空间中n+1个点的坐标,求一个球的球心位置满足这个球经过这n+1个点.
很显然这n+1个点到球心的距离r相等,也就是说设球心为
,每一个点
都满足:
发现这个式子中的未知数只有a数组和r,由于题目中给出了n+1个坐标,考虑把相邻两个式子相减去掉
,式子变为:
发现这样子就只剩下一个包含n个方程和n个未知量的方程组了,高斯消元一下即可.
代码如下:
#include<bits/stdc++.h>
using namespace std;
#define Abigail inline void
typedef long long LL;
const int N=10;
struct matrix{
double a[N+9][N+9];
int n;
}a;
double ans[N+9],e[N+9][N+9];
int n;
void gauss(matrix a){
int r,n=a.n;
double d;
for (int i=1;i<=n;++i){
r=i;
for (int j=i+1;j<=n;++j)
if (fabs(a.a[j][i])>fabs(a.a[r][i])) r=j;
if (r^i)
for (int j=1;j<=n+1;++j)
swap(a.a[i][j],a.a[r][j]);
d=a.a[i][i];
for (int j=i;j<=n+1;++j)
a.a[i][j]/=d;
for (int j=i+1;j<=n;++j){
d=a.a[j][i];
for (int k=i;k<=n+1;++k)
a.a[j][k]-=a.a[i][k]*d;
}
}
ans[n]=a.a[n][n+1];
for (int i=n-1;i>=1;--i){
ans[i]=a.a[i][n+1];
for (int j=i+1;j<=n;++j)
ans[i]-=a.a[i][j]*ans[j];
}
}
Abigail into(){
scanf("%d",&n);
for (int i=1;i<=n+1;++i)
for (int j=1;j<=n;++j)
scanf("%lf",&e[i][j]);
}
Abigail work(){
a.n=n;
for (int i=1;i<=n;++i)
for (int j=1;j<=n;++j)
a.a[i][j]=2.0*(e[i+1][j]-e[i][j]),a.a[i][n+1]+=e[i+1][j]*e[i+1][j]-e[i][j]*e[i][j];
gauss(a);
}
Abigail outo(){
for (int i=1;i<n;++i)
printf("%.3lf ",ans[i]);
printf("%.3lf",ans[n]);
}
int main(){
into();
work();
outo();
return 0;
}