高斯消除的逻辑错误 [英] Logic error for Gauss elimination

问题描述

高斯消除代码的逻辑错误问题...这段代码来自我1990年代的数值方法文本。

Logic error problem with the Gaussian Elimination code...This code was from my Numerical Methods text in 1990's. The code is typed in from the book- not producing correct output...

示例运行：

         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 . . .

从文本复制的代码...

The code as copied from the text...

//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

The solution should be: 1.058824, 1.823529, 0.882353 with det as -102.000000

任何洞察力是值得赞赏的...

Any insight is appreciated...

推荐答案

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

这应该是

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

现在的方式，我相信这会让你划分枢轴行本身，将其清零。

The way it is now, I believe this will cause you to divide the pivot row by itself, zeroing it out.

这篇关于高斯消除的逻辑错误的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！