c - 高斯消元法及降维系统的求解

Question

我正在尝试在 C 中通过旋转实现高斯消元法，但即使使用最简单的输入，解决方案也会出错。这是我写的代码

#include <stdio.h>
#include <math.h>

#define N 15

void gauss(float a[][N], float *b, int n);
void solution(float a[][N], float *b, float *x, int n);

int main(void)
{
    int n, i, j;

    printf("Enter an integer n<=15: ");
    scanf("%d", &n);

    float a[n][N], b[n], x[n];

    for(i=0; i<n; i++)
        for(j=0; j<n; j++) {
            printf("A[%d][%d] = ", i+1, j+1);
            scanf(" %f", &a[i][j]);
        }

    printf("\n");

    for(i=0; i<n; i++) {
        printf("b[%d] = ", i+1);
        scanf(" %f", b+i);
    }

    gauss(a, b, n);

    solution(a, b, x, n);

    printf("The solution is ");

    for(i=0; i<n; i++)
        printf("%g ", b[i]);

    printf("\n");

    return 0;

}

void gauss(float a[][N], float *b, int n)
{
    int i, j, r, l;

    float pivot, temp, m;

    for(r=0; r<n-1; r++) {
        for(l=i=r+1, pivot=fabs(a[l][r]); i<n; i++)
            if(fabs(a[i][r]) > pivot) {
                l = i;
                pivot = fabs(a[l][r]);
            }

        if(pivot > fabs(a[r][r])) {
            for(j=0; j<n; j++) {
                temp = a[r][j];
                a[r][j] = a[l][j];
                a[l][j] = temp;
            }

            temp = b[r];
            b[r] = b[l];
            b[l] = temp;
        }

        for(i=r+1; i<n; i++) {
            m = a[i][r]/a[r][r];

            for(j=0; j<n; j++)
                a[i][j] = a[i][j] - m*a[r][j];

            b[i] = b[i] - m*b[r];
        }
    }
}

void solutions(float a[][N], float *b, float *x, int n)
{
    int i, j;
    float sum = 0;

    for(i=n-1; i>=0; i--) {
        for(j=i; j<n; j++)
            sum += a[i][j]*x[j];
        x[i] = (b[i] - sum)/a[i][i];
        sum = 0;
    }
}

我用一个简单的系统测试了它，输出如下

Enter an integer n<=15: 2

A[1][1] = 1
A[1][2] = 1
A[2][1] = 1
A[2][2] = -1

b[1] = 1
b[2] = 0

The solution is 1 -1

我试过调整索引但没有用，然后我把上面的例子写在纸上但没有发现任何错误。怎么了？

Answer 1

第一个问题是您正在打印错误的值。您正在打印 b 但您已将解决方案存储在 x

中

for(i=0; i<n; i++)
    printf("%g ", b[i]);
                 ^^^
                  should be x

除此之外，您的solutions 功能似乎有误。

1) 您应该在循环之间重置sum。

2) 您正在使用未初始化的x

尝试这样的事情:

void solutions(float a[][N], float *b, float *x, int n)
{
    int i, j;
    float sum;

    for(i=n-1; i>=0; i--) {
      sum = 0;                       // Notice this line
      for(j=i+1; j<n; j++)           // Notice this line, i.e. j=i+1
      {
        sum += a[i][j]*x[j];
      }
      x[i] = (b[i] - sum)/a[i][i];
      sum = 0;
    }
}