解方程有高斯消元和高斯-约当消元两种了,前者是先化为上三角形,然后化为最简形,后者是直接化为最简形的。不过前者效率稍稍高一点,然后不知道为什么我就写了个前者啊……
方程要自己实现构造好放在a数组里,a数组最边上存增广矩阵,var存变量数(矩阵的列数),equ存方程数(矩阵的行数)
自己试着码了一个模板,然后放在这里保存下,慢慢更改
然后一些说明:
- 这个模板应该是自己码的考虑范围比较全面的一个了,但是,这个不是效率最高的
- 如果我们明确知道要求的未知数是哪一个,那么我们在最后把行阶梯型化成最简形的过程中,一解出那个未知数就可以结束函数了(毕竟整个化简过程中交换的是行,不是列,未知数没有受到影响)
- 这个函数可以对无解和无穷组解做出反应(无解返回-1,无穷组解返回0),但是这只是个模板而已,具体问题肯定需要不同的处理才行
- 重复一下,具体问题需要不同的处理,每题说不定都有些可以化简的部分,要好好思考
#include<cstdio> #include<cstring> #include<cmath> #include<algorithm> using namespace std; const int MAXN=100+5; const double EPS=1e-7; double a[MAXN][MAXN]; //数组的有效范围是[0,equ),[0,val]因为想了想还是把x放在里面比较方便 int Gauss(int equ,int val) { int maxrow,r,c,x_num; //每一次对一个未知数进行消元 for(r=c=0;r<equ&&c<val;++c) { maxrow=r; for(int i=r+1;i<equ;++i) if(fabs(a[i][c])>fabs(a[maxrow][c])) maxrow=i; //当前未知数每一个未处理的行都为0,产生了自由未知数,换下一个未知数 if(fabs(a[maxrow][c])<EPS) continue; if(maxrow!=r) for(int i=c;i<=val;++i) swap(a[r][i],a[maxrow][i]); //目标还算明确,化成行阶梯形(上三角形),所以这里我们不对自己进行化简 for(int i=r+1;i<equ;++i) if(fabs(a[i][c])>EPS) for(int j=val;j>=c;--j) A[i][j]-=A[i][c]/A[r][c]*A[r][j]; //一行消完再消下一行,前面的continue就是为了躲避这句话 ++r; } //存在(0 0 0 0 0 0 d[j][c]) (d[j][c] != 0) 的情况,方程无解 for(int i=r;i<equ;++i) if(fabs(a[i][val])>EPS) return -1; //这部分是自己加的,要是线性无关的方程数没有未知数那么多,貌似就是无穷组解了啊,先写在这里 if(r<val) return 0; //接下来枚举每一列(每一个未知数),c>=1就够了,我们没必要处理第一列,因为当其它列处理完的时候,第一列就不需要处理了 for(c=val-1;c>=1;--c) { for(r=c-1;r>=0;--r) if(fabs(a[r][c])>EPS) { a[r][val]-=a[r][c]/a[c][c]*a[c][equ]; // a[r][c]=0; } a[c][val]/=a[c][c]; // a[c][c]=1; } return 1; } int main() { int n,m; while(~scanf("%d%d",&n,&m)) { for(int i=0;i<n;++i) for(int j=0;j<=m;++j) scanf("%lf",&a[i][j]); Gauss(n,m); for(int i=0;i<n;++i) { for(int j=0;j<=m;++j) printf("%f ",a[i][j]); putchar(10); } } }
main中的作用是读入一个方程,然后调用函数求解,大家自己看看也能看懂的吧。
然后然后就是要有题目练练手了。
第一道:HDU 3359 Kind of a Blur,构建方程的时候要处理好,其实就是裸的高斯消元,代码就不放了,前面的模板可以ac的。
第二道:HDU 3364 Lanterns,判断方程自由未知数的个数,若方程的自由未知数个数为k,答案就是2^k,如果无解,答案就是0,模板后面求解的部分完全不需要了,只要化成上三角形即可。