高斯消元的逻辑错误

Question

高斯消除代码的逻辑错误问题...此代码来自我在 1990 年代的数值方法文本。代码是从书中输入的 - 没有产生正确的输出......

示例运行：

         SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS
         USING GAUSSIAN ELIMINATION

 This program uses Gaussian Elimination to solve the
 system Ax = B, where A is the matrix of known
 coefficients, B is the vector of known constants
 and x is the column matrix of the unknowns.
 Number of equations: 3

 Enter elements of matrix [A]
 A(1,1) = 0
 A(1,2) = -6
 A(1,3) = 9
 A(2,1) = 7
 A(2,2) = 0
 A(2,3) = -5
 A(3,1) = 5
 A(3,2) = -8
 A(3,3) = 6

 Enter elements of [b] vector
 B(1) = -3
 B(2) = 3
 B(3) = -4

         SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS

         The solution is
         x(1) = 0.000000
         x(2) = -1.#IND00
         x(3) = -1.#IND00
         Determinant = -1.#IND00
Press any key to continue . . .

从文本中复制的代码...

//Modified Code from C Numerical Methods Text- June 2009

#include <stdio.h>
#include <math.h>
#define MAXSIZE 20

//function prototype
int gauss (double a[][MAXSIZE], double b[], int n, double *det);

int main(void)
{
    double a[MAXSIZE][MAXSIZE], b[MAXSIZE], det;
    int i, j, n, retval;

    printf("\n \t SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS");
    printf("\n \t USING GAUSSIAN ELIMINATION \n");
    printf("\n This program uses Gaussian Elimination to solve the");
    printf("\n system Ax = B, where A is the matrix of known");
    printf("\n coefficients, B is the vector of known constants");
    printf("\n and x is the column matrix of the unknowns.");

    //get number of equations
    n = 0;
    while(n <= 0 || n > MAXSIZE)
    {
        printf("\n Number of equations: ");
        scanf ("%d", &n);
    }

    //read matrix A
    printf("\n Enter elements of matrix [A]\n");
    for (i = 0; i < n; i++)
        for (j = 0; j < n; j++)
        {
            printf(" A(%d,%d) = ", i + 1, j + 1);
            scanf("%lf", &a[i][j]);
        }
    //read {B} vector
    printf("\n Enter elements of [b] vector\n");
    for (i = 0; i < n; i++)
    {
        printf(" B(%d) = ", i + 1);
        scanf("%lf", &b[i]);
    }

    //call Gauss elimination function
    retval = gauss(a, b, n, &det);

    //print results
    if (retval == 0)
    {
        printf("\n\t SOLUTION OF SIMULTANEOUS LINEAR EQUATIONS\n");
        printf("\n\t The solution is");
        for (i = 0; i < n; i++)
            printf("\n \t x(%d) = %lf", i + 1, b[i]);
        printf("\n \t Determinant = %lf \n", det);
    }
    else
        printf("\n \t SINGULAR MATRIX \n");

    return 0;
 }

/* Solves the system of equations [A]{x} = {B} using       */
/* the Gaussian elimination method with partial pivoting.  */
/* Parameters:                                             */
/*      n         - number of equations                    */
/*      a[n][n]   - coefficient matrix                     */
/*      b[n]      - right-hand side vector                 */
/*      *det      - determinant of [A]                     */

int gauss (double a[][MAXSIZE], double b[], int n, double *det)
{
    double tol, temp, mult;
    int npivot, i, j, l, k, flag;

    //initialization
    *det = 1.0;
    tol = 1e-30;        //initial tolerance value
    npivot = 0;
        //mult = 0;

    //forward elimination
    for (k = 0; k < n; k++)
    {
        //search for max coefficient in pivot row- a[k][k] pivot element
        for (i = k + 1; i < n; i++)
        {
            if (fabs(a[i][k]) > fabs(a[k][k]))
            {
                //interchange row with maxium element with pivot row
                npivot++;
                for (l = 0; l < n; l++)
                {
                   temp = a[i][l];
                   a[i][l] = a[k][l];
                   a[k][l] = temp;
                }
                temp = b[i];
                b[i] = b[k];
                b[k] = temp;
            }
        }
        //test for singularity
        if (fabs(a[k][k]) < tol)
        {
           //matrix is singular- terminate
           flag = 1;
           return flag;
        }
        //compute determinant- the product of the pivot elements
        *det = *det * a[k][k];

        //eliminate the coefficients of X(I)
        for (i = k; i < n; i++)
        {
            mult = a[i][k] / a[k][k];
            b[i] = b[i] - b[k] * mult;  //compute constants
            for (j = k; j < n; j++)     //compute coefficients
                a[i][j] = a[i][j] - a[k][j] * mult;
        }
    }
    //adjust the sign of the determinant
    if(npivot % 2 == 1)
      *det = *det * (-1.0);

    //backsubstitution
    b[n] = b[n] / a[n][n];
    for(i = n - 1; i > 1; i--)
    {
        for(j = n; j > i + 1; j--)
            b[i] = b[i] - a[i][j] * b[j];
        b[i] = b[i] / a[i - 1][i];
    }
    flag = 0;
    return flag;
}

解决方案应该是：1.058824、1.823529、0.882353，det 为 -102.000000

感谢任何见解...

Answer 1

//eliminate the coefficients of X(I)
for (i = k; i < n; i++)

应该是这样的

for (i = k + 1; i < n; i++)

现在的方式，我相信这会导致您自己划分数据透视行，将其归零。

Answer 2

这可能无法按照您预期的方式回答您的问题，但编写您自己的数值稳定矩阵算法与自己动手手术一样明智。

有一个非常好的库，名为TNT/JAMA，来自一个有信誉的来源 (NIST)，它在 C++ 中进行基本矩阵数学运算。要解决 Ax=B，首先要考虑因素 A（QR decomposition 是一个很好的通用方法，您可以使用 LU，但它在数值上不太稳定），然后调用 solve(B)。这既适用于精确的方阵（受数值计算问题的影响），也适用于得到最小二乘答案的超定系统。