zoukankan      html  css  js  c++  java
  • [笔记] 浅谈高斯消元

    高斯消元。。。当初以为自己学会了。。。后来。。。


    话说这个东西好像最早出现于《九章算术》诶(古代人就是强)

    废话不说,进入正题。。。

    前置知识

    高斯消元法是解线性方程组的方法之一 

    首先,线性方程组了解一下:

      可认为线性方程组就是一次方程组。如图:

    如果存在常数c1,c2,c3,...,cn代替x1,x2,x3,...,xn,使上图的每个方程成立,则称(c1,c2,c3,...,cn)为方程组的一个解. 解方程就是求其全部解

    显然对线性方程组进行以下三种变换,所得方程组解集不变(这三种变换被称为同解变换)

      (1) 对换两个方程(换法变换)(其实就是交换两个方程的位置,上下交换) ;

      (2) 用非零数 c 乘以某一个方程(倍法变换);

      (3) 把某一个方程的 k 倍加到另一个方程上去( 消法变换) .

    的确显然吧...

    再来了解一下矩阵(=_=):

    由 s*n 个数 aij ( i = 1, 2, … , s; j = 1, 2, … , n) 排成的矩形表(可理解为长是n,宽是s的二维数组。。),称为 s 行 n 列矩阵 .aij称为矩阵的( i , j) 元 .通常用大写字母或 ( aij )sn 表示矩阵 .如果 s = n(长等于宽)称A是 n 阶方阵或 n 阶矩阵 .

    全零矩阵记做Osn,或O;

    n阶单位矩阵记做Inn,或I(单位矩阵是个方阵,从左上角到右下角的对角线(称为主对角线)上的元素均为1。除此以外全都为0。)

    对于矩阵(aij)sn、(bij)ns,若aij=bji(1<=i<=s,1<=j<=n)则称B为A的转置矩阵,记做B=AT

    另:系数矩阵和增广矩阵

    把线性方程组的系数aij列入矩阵,称为系数矩阵。

    如图:

    而在系数矩阵右侧加入一列常数项,称其为增广矩阵。

    如图:

    (好吧增广矩阵就是类似把未知数扣掉了再扔到矩阵里。。)

    对矩阵的行(列)做 线性方程组 的三种同解变换(上面的),称为矩阵的初等变换。

    高斯消元法

    对方程组的增广矩阵做初等变换,将系数部分化为对角线,上三角形,或阶梯形

    思路:找到一个不为零的系数,设其所在行为pos,列为i,然后拿这个数去消掉其它行,但同列的数;但是消元的时候要有一定顺序

    步骤:

    先把一个个方程扔到矩阵中,构建增广矩阵

    然后从第i行(i从1到n,按顺序)开始,在i+1行到n中的第i列(上面的方程已经符合上三角性质了)选取一个不为零的系数,记其所在行为pos,

    交换第pos行与第i行,然后pos行变成了第i行(注:后面的第i行为交换完后的,即交换前的pos行)

    至于为什么要找最大的,这样有利于提高精度。。(学长说的。。)

    于是拿第i行的第i项去消i+1行到n行的第i项

    这样可以消成上三角矩阵;

    然后从最后一行开始,向上消去不需要的系数(最后一行只有一个系数,可以直接计算答案)

    代码(Luogu 模板):

    #include<cstdio>
    #include<iostream>
    #define R register int
    using namespace std;
    const double eps=1E-11;
    inline int g() {
        R ret=0,fix=1; register char ch; while(!isdigit(ch=getchar())) fix=ch=='-'?-1:fix;
        do ret=ret*10+(ch^48); while(isdigit(ch=getchar())); return ret*fix;
    }
    int n;
    double a[110][110];
    inline bool ck0(double x) {return x<eps&&x>-eps;}
    inline void Gauss() {
        for(R i=1,pos;i<=n;++i) {
            pos=i; for(R j=i+1;j<=n;++j) if(ck0(a[pos][i])||a[pos][i]<a[j][i]) pos=j;
            if(pos!=i) for(R j=1;j<=n+1;++j) swap(a[i][j],a[pos][j]);
            if(ck0(a[i][i])) {printf("No Solution
    "); return ;}
            for(R j=i+1;j<=n;++j) {
                register double tmp=a[j][i]/a[i][i];
                for(R k=1;k<=n+1;++k) a[j][k]-=tmp*a[i][k];
            }
        } for(R i=n;i>=1;--i) {
            for(R j=n;j>i;--j) a[i][n+1]-=a[i][j]*a[j][n+1];//把这行系数代表的答案消掉
            a[i][n+1]/=a[i][i];
        } for(R i=1;i<=n;++i) printf("%.2lf
    ",a[i][n+1]);
    }
    signed main() {
        n=g(); for(R i=1;i<=n;++i) for(R j=1;j<=n+1;++j) a[i][j]=g();
        Gauss();
    }

     

    (例题)

    Luogu P3389 【模板】高斯消元法

    Luogu P4035 [JSOI2008]球形空间产生器

    POJ1830   

    etc. 


    2019.05.10

     
     
  • 相关阅读:
    django集成django-xadmin
    Django设置 DEBUG=False后静态文件无法加载解决
    Django ORM必会的查询方法
    django admin-过滤器
    Django settings.py 中设置访问 MySQL 数据库【一种是直接在 settings.py 文件中直接写数据库信息,另一种是读文件获取数据库信息】
    django-admin之ModelAdmin最全解释
    SPL(Standard PHP Library 标准PHP类库)
    rsync 数据同步
    PHP 安装memcache.so 和memcached.so
    linux 安装memcached
  • 原文地址:https://www.cnblogs.com/Jackpei/p/10845807.html
Copyright © 2011-2022 走看看