【问题标题】:Calculate matrix determinant with partial pivoting Gauss in C用C中的部分旋转高斯计算矩阵行列式
【发布时间】:2011-10-19 16:29:16
【问题描述】:

我正在尝试在 C 中创建一个简单的控制台应用程序,它将使用高斯部分旋转消除方法计算矩阵的行列式。我遇到的两个问题是: - 有人告诉我,有些矩阵不适用于这种方法(从数学上讲),在阅读了谷歌上的文章后,我找不到那个特殊情况 - 经过大量测试后,我发现我的程序不适用于某些矩阵,经过 2 天的“浪费”时间编辑和撤消后,我找不到问题。

我们非常欢迎任何类型的改进。我只是从 C 开始。

#include<stdio.h>
#include<cstdlib>
#include<math.h>
#include<conio.h>
#include<windows.h>

// calculate biggest element on column

int indice_max(int dim, int col, float coloana[20][20]) {

    float max = 0;
    int indice;

    for(int i = 1; i <= dim; i++)
        if(fabs(max) < fabs(coloana[i][col])) {
            max = coloana[i][col];
            indice = i;
        }

    return indice;

}

// permute 2 lines

void permutare_linie(int linie1, int linie2, int dim, float matrice[20][20]) {

    float aux;

    for(int i = 1; i <= dim; i++) {
        aux = matrice[linie1][i];
        matrice[linie1][i] = matrice[linie2][i];
        matrice[linie2][i] = aux;
    }

}

// print matrix

void afisare_matrice(int dimensiune, float matrice[20][20], int lpiv) {

    for(int i = 1; i<= dimensiune; i++) {
        for(int j = 1; j <= dimensiune; j++) {
            if(i == lpiv)
                SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE), BACKGROUND_GREEN);
            else
                SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE), FOREGROUND_RED | FOREGROUND_GREEN | FOREGROUND_BLUE );
            printf("%4.2f ", matrice[i][j]);
        }
        printf("\n");
    }

}


void main(void) {

    float matrice[20][20];
    int dimensiune ;
    float rezultat = 1;
    float pivot;
    int lpiv;
    int cpiv;
    int optiune;
    while(1) {

        // MENU

        printf("ALEGET OPTIUNEA:\n");
        printf("1) Calculate matrix determinant\n");
        printf("2) Exit\n");
        scanf("%d", &optiune);

        if(optiune == 1) {

            // Read determinant dimension

            printf("Matrix dimension:");
            scanf("%d", &dimensiune);

            // Read determinant

            for(int i = 1; i <= dimensiune; i++)
                for(int j = 1; j <= dimensiune; j++) {
                    printf("M[%d][%d]=", i, j);
                    scanf("%f", &matrice[i][j]);
            }

            // pivot initial coords

            lpiv = 1;
            cpiv = 1;

            printf("\n----- Entered Matrix -----\n\n");
            afisare_matrice(dimensiune, matrice, 0);
            printf("\n");

            for(int pas = 1; pas <= dimensiune - 1; pas++) {

                if(fabs(matrice[lpiv][cpiv]) > fabs(matrice[indice_max(dimensiune, cpiv, matrice)][cpiv])) {
                    permutare_linie(lpiv, indice_max(dimensiune, cpiv, matrice), dimensiune, matrice);
                    rezultat = -(rezultat);
                }

                pivot = matrice[lpiv][cpiv];


                for(int inm = 1; inm <= dimensiune; inm++) {
                    matrice[lpiv][inm] = matrice[lpiv][inm] / pivot;
                }

                rezultat *= fabs(pivot);

                // transform matrix to a superior triangular 
                for(int l = lpiv+1; l <= dimensiune; l++)
                    for(int c=cpiv+1; c <= dimensiune; c++) {
                        matrice[l][c] -= matrice[l][cpiv] * matrice[lpiv][c] / matrice[lpiv][cpiv];
                    }

                for(int i = lpiv + 1; i <= dimensiune; i++)
                    matrice[i][cpiv] = 0;
                // afisam rezultat / pas

                printf("----- Step %d -----\n\n", pas);
                afisare_matrice(dimensiune, matrice, lpiv);
                printf("\nResult after step %d : %4.2f\n\n", pas, rezultat);
                lpiv++;
                cpiv++;
            }

            // final result

            rezultat = rezultat * matrice[dimensiune][dimensiune];
            printf("----- REZULTAT FINAL -----\n\n");
            SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE), FOREGROUND_RED | FOREGROUND_INTENSITY);
            printf("Rezultat = %4.2f\nRezultat rotunjit:%4.0f\n\n", rezultat, floorf(rezultat * 100 +  0.5) / 100);
            SetConsoleTextAttribute(GetStdHandle(STD_OUTPUT_HANDLE), FOREGROUND_RED | FOREGROUND_GREEN | FOREGROUND_BLUE );

        }
        else {
            exit(0);
        }
    }
}

【问题讨论】:

  • 好吧,去图书馆找一本关于数值方法的书,比如数值方法,Burden&Faires,如果我没记错的话,看看他们的高斯消除的 C 实现,HTH
  • wikipedia entry中有解释“不可逆的方阵称为奇异或退化的。方阵是奇异的当且仅当其行列式为0。奇异矩阵在某种意义上是罕见的如果你选择一个随机方阵,它几乎肯定不会是奇异的。”
  • 不用去图书馆了!查看第 2 章:nrbook.com/a/bookcpdf.php
  • @user786653 您的解释说当矩阵不可逆时。奇异矩阵是当它的行列式为 0 时,但我需要计算行列式以查看它是否为 0。网络上有此过程的代码示例,但它们与我的大学教授正在执行的步骤不同。
  • 查看@anatolygs 的回答,奇异矩阵在某些时候会有一个零列,不需要显式计算行列式。

标签: c matrix partial gaussian determinants


【解决方案1】:

你的代码做了一些划分:

matrice[lpiv][inm] = matrice[lpiv][inm] / pivot;

如果它碰巧被零除,则会发生错误。我想这会发生在零矩阵上。

您的代码似乎实际上是在尝试反转矩阵,而不仅仅是计算行列式。

【讨论】:

  • 当然,我认为对角线上的元素之一为0,但没想到在将矩阵转换为上三角矩阵的操作后它们可能为0。谢谢。我想我也发现了另一个问题。我正在整个列上搜索最大元素,而不仅仅是在枢轴下方。
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2011-02-24
  • 2013-05-12
  • 2012-10-28
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多