计算方法实验(五):高斯列主元消去法

Gauss列主元消去法数学原理

高斯(Gauss)列主元消去法:对给定的 n n n阶线性方程组 A x = b Ax = b Ax=b,首先进行列主元消元过程,然后进行回代过程,最后得到解或确定该线性方程组是奇异的。

如果系数矩阵的元素按绝对值在数量级方面相差很大,那么,在进行列主元消元过程前,先把系数矩阵的元素进行行平衡:系数矩阵的每行元素和相应的右端向量元素同除以该行元素绝对值最大的元素。这就是所谓的平衡技术。然后再进行列主元消元过程。

如果真正进行运算去确定相对主元,则称为显式相对Gauss列主元消去法;如果不进行运算,也能确定相对主元,则称为隐式相对Gauss列主元消去法。

显式相对Gauss列主元消去法:对给定的 n n n阶线性方程组 A x = b Ax = b Ax=b,首先进行列主元消元过程,在消元过程中利用显式平衡技术,然后进行回代过程,最后得到解或确定该线性方程组是奇异的。

隐式相对Gauss列主元消去法:对给定的 n n n阶线性方程组 A x = b Ax = b Ax=b,首先进行列主元消元过程,在消元过程中利用隐式平衡技术,然后进行回代过程,最后得到解或确定该线性方程组是奇异的。

Gauss列主元消去法

1对 k = 1 , 2 , ⋯   , n − 1 k = 1,2,\cdots,n - 1 k=1,2,,n1,做1.1—1.3,消元过程

1.1 寻找最小的正整数 p p p k ≤ p ≤ n k \leq p \leq n kpn ∣ a pk ∣ = max ⁡ k ≤ j ≤ n ∣ a jk ∣ \left| a_{\text{pk}} \right| = \max_{k \leq j \leq n}\left| a_{\text{jk}} \right| apk=maxkjnajk。如果 a pk = 0 a_{\text{pk}} = 0 apk=0,输出奇异标志,停机;

1.2 如果 p ≠ k p \neq k p=k,那么交换 p , k p,k p,k两行;

1.3 对 i = k + 1 , ⋯   , n i = k + 1,\cdots,n i=k+1,,n , 记 ,记 m ik = a ik / a kk m_{\text{ik}} = a_{\text{ik}}/a_{\text{kk}} mik=aik/akk$,计算

{   a ij = a ij − a kj m ik   i = k + 1 , ⋯   , n   j = k + 1 , ⋯   , n   b i = b i − b k m ik   i = k + 1 , ⋯   , n   \left\{ \begin{matrix} \ a_{\text{ij}} = a_{\text{ij}} - a_{\text{kj}}m_{\text{ik}} \\ \ i = k + 1,\cdots,n \\ \ j = k + 1,\cdots,n \\ \ b_{i} = b_{i} - b_{k}m_{\text{ik}} \\ \ i = k + 1,\cdots,n \\ \end{matrix} \right.\  aij=aijakjmik i=k+1,,n j=k+1,,n bi=bibkmik i=k+1,,n 

  1. 如果 a nn = 0 a_{\text{nn}} = 0 ann=0输出奇异标志,停机;

  2. x n = b n / a nn x_{n} = b_{n}/a_{\text{nn}} xn=bn/ann,回代过程

  3. k = n − 1 , ⋯   , 2 , 1 k = n - 1,\cdots,2,1 k=n1,,2,1,置 x k = ( b k − ∑ j = k + 1 n a kj x j ) / a kk x_{k} = (b_{k} - \sum_{j = k + 1}^{n}{a_{\text{kj}}x_{j}})/a_{\text{kk}} xk=(bkj=k+1nakjxj)/akk

程序流程

代码

#include <cmath>
#include <cstdio>
#include <iostream>
using namespace std;
#define N 10

int n;
double a[N][N], b[N], x[N];

int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n; j++) scanf("%lf", &a[i][j]);
    for (int i = 1; i <= n; i++) scanf("%lf", &b[i]);
    for (int k = 1; k < n; k++) {
        int p = k;
        double maxabs = fabs(a[k][k]);
        for (int j = k + 1; j <= n; j++)
            if (fabs(a[j][k]) - maxabs > 0) {
                p = j;
                maxabs = fabs(a[j][k]);
            }
        if (a[p][k] == 0) {
            printf("Singular");
            return 0;
        }
        if (p != k) {
            double tmp;
            for (int j = 1; j <= n; j++) {
                tmp = a[p][j];
                a[p][j] = a[k][j];
                a[k][j] = tmp;
            }
            tmp = b[p];
            b[p] = b[k];
            b[k] = tmp;
        }
        for (int i = k + 1; i <= n; i++) {
            double m_ik = a[i][k] / a[k][k];
            for (int j = k + 1; j <= n; j++) a[i][j] -= a[k][j] * m_ik;
            b[i] -= b[k] * m_ik;
        }
    }
    if (a[n][n] == 0) {
        printf("Singular");
        return 0;
    }
    x[n] = b[n] / a[n][n];
    for (int k = n - 1; k >= 1; k--) {
        double sigma = 0.0;
        for (int j = k + 1; j <= n; j++) sigma += a[k][j] * x[j];
        x[k] = (b[k] - sigma) / a[k][k];
    }
    for (int i = 1; i <= n; i++) printf("%lf\t", x[i]);
    return 0;
}
#include #include #define N 100 #define epsilon 1e-6 float a[N][N+1]; void menu( ) { printf("\t\t%c%c%c^_^Gauss列主元消去法求解线性方程组^_^%c%c%c\n\n",1,1,1,1,1,1); printf("强烈建议您先阅读以下几点后在运行:\n"); printf("1.这是用Gauus列主元消去法求解线性方程组的应用程序\n"); printf(" (Gauus全主元消去法类似可做,读者有兴趣的话可自行而做)\n"); printf("2.请您先了解Gauus列主元消去法的主要思想\n"); } void main( ) { int i,j,k,n; float t,s=0; char choice; menu( ); loop: printf("\n请输入系数方阵的阶数:"); scanf("%d",&n); while(n>0) { printf("\n"); printf("请输入增广阵矩:\n"); for(i=0;i<n;i++) for(j=0;j<n+1;j++) scanf("%f",&a[i][j]);/*存储增广阵矩*/ for(k=0;k<n-1;k++) { for(i=k+1;i fabs(a[k][k]) ) for(j=k;j<n+1;j++) { t=a[k][j]; a[k][j]=a[i][j]; a[i][j]=t; } if( fabs(a[k][k]) < epsilon)/*最大数小于很小数时退出*/ { printf("\n错误,Gauss列主元消去法无法忍受,在%d步退出!\n",k+1); printf("还要再计算其他的么(Y/N)?"); scanf("%c",&choice); if(choice=='Y' || choice=='y')/*判断用户输入*/ goto loop; else return; } for(i=k+1;i<n;i++) { a[i][k]=a[i][k] / a[k][k]; for(j=k+1;j=0;k--) { s=0; for(j=k+1;j<n;j++) s+=a[k][j]*a[j][n]; a[k][n]=( a[k][n]-s ) / a[k][k]; } printf("\n*****运行结果*****\n"); for(i=0;i<n;i++) printf(" x[%d]=%.4f\n",i+1,a[i][n]); printf(" 谢谢使用!\n"); printf("还要再计算其他的么(Y/N)?"); getchar( ); scanf("%c",&choice); if(choice=='Y' || choice=='y')/*判断用户输入*/ goto loop; else return; } }
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值