(C语言)编写解n阶方程组的列主元Gauss消去法的通用程序

一、题目

二、解法

main.h

#ifndef __MAIN_H__
#define __MAIN_H__

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "matrix_solve.h"

#endif

main.c

#include "main.h"

int main(void)
{
    double matrix_A[][9] = {{31, -13, 0, 0, 0, -10, 0, 0, 0},
							{-13, 35, -9, 0, -11, 0, 0, 0, 0},
							{0, -9, 31, -10, 0, 0, 0, 0, 0},
							{0, 0, -10, 79, -30, 0, 0, 0, -9},
							{0, 0, 0, -30, 57, -7, 0, -5, 0},
							{0, 0, 0, 0, -7, 47, -30, 0, 0},
							{0, 0, 0, 0, 0, -30, 41, 0, 0},
							{0, 0, 0, 0, -5, 0, 0, 27, -2},
							{0, 0, 0, -9, 0, 0, 0, -2, 29}};
    double matrix_b[9] = {-15, 27, -23, 0, -20, 12, -7, 7, 10};
    
    Matrix_Solve.matrix_A = (double *)matrix_A;
    Matrix_Solve.matrix_b = matrix_b;
    Matrix_Solve.N = sizeof(matrix_b) / sizeof(matrix_b[0]);
    
    Matrix_Solve.matrix_x = Matrix_Solve.solve_matrix(Matrix_Solve.matrix_A, Matrix_Solve.matrix_b, Matrix_Solve.N);

    for(int i = 0; i < Matrix_Solve.N; i++)
    {
        printf("%.5f\r\n", Matrix_Solve.matrix_x[i]);
    }

    free(Matrix_Solve.matrix_x);
    
    getchar();

    return 0;
}

matrix_solve.h

#ifndef __MATRIX_SOLVE_H__
#define __MATRIX_SOLVE_H__

typedef struct 
{
    double* matrix_A;
    double* matrix_b;
    double* matrix_x;
    int N;
    double* (*solve_matrix)(double*, double*, int);
} Matrix_Solve_t;

extern Matrix_Solve_t Matrix_Solve;

#endif

matrix_solve.c

#include "main.h"

static double* solve_matrix(double*, double*, int);

Matrix_Solve_t Matrix_Solve = {NULL, NULL, NULL, 0, solve_matrix};

static double* solve_matrix(double* matrix_A, double* matrix_b, int N)
{
    int i = 0, j = 0, k = 0;
    double Max = 0, Temp = 0, l = 0, y = 0;
    int Max_k = 0;

    double* matrix_A_extend;
    matrix_A_extend = (double *)malloc(N * (N+1) * sizeof(double));
    for(i = 0; i < N; i++)
    {
        for(j = 0; j < N; j++)
        {
            matrix_A_extend[i*(N+1)+j] = matrix_A[i*N+j];
        }
        matrix_A_extend[i*(N+1)+N] = matrix_b[i];
    }

    for(i = 0; i < N - 1; i++)
    {
        Max = matrix_A_extend[i*(N+1)+i];
        Max_k = i;

        for(k = i + 1; k < N; k++)
        {
            if(matrix_A_extend[k*(N+1)+i] > Max)
            {
                Max = matrix_A_extend[k*(N+1)+i];
                Max_k = k;
            }
        }
 
        if(Max_k != i)
        {
            for(j = i; j < N + 1; j++)
            {
                Temp = matrix_A_extend[i*(N+1)+j];
                matrix_A_extend[i*(N+1)+j] = matrix_A_extend[Max_k*(N+1)+j];
                matrix_A_extend[Max_k*(N+1)+j] = Temp;
            }
        }

        for(k = i + 1; k < N; k++)
        {
            l = matrix_A_extend[k*(N+1)+i] / matrix_A_extend[i*(N+1)+i];
            for(j = i; j < N + 1; j++)
            {
                matrix_A_extend[k*(N+1)+j] -= (l * matrix_A_extend[i*(N+1)+j]);
            }
        }

    }

    double* matrix_x;
    matrix_x  = (double *)malloc(N*sizeof(double));
    matrix_x[N-1] = matrix_A_extend[(N-1)*(N+1)+N] / matrix_A_extend[(N-1)*(N+1)+N-1];

    if(N > 1)
    {
        for(i = N - 2; i >= 0; i--) 
        {
            y = matrix_A_extend[i*(N+1)+N];
            for(j = i + 1; j < N; j++)
            {
                y -= matrix_A_extend[i*(N+1)+j] * matrix_x[j]; 
            }
            matrix_x[i] = y / matrix_A_extend[i*(N+1)+i];
        }
    }
    
    free(matrix_A_extend);
    
    return matrix_x;
}

三、运行结果

 四、小结

C语言中一般将二维数组传给函数时，不能省略第二维的长度。但可以巧妙地将二维数组作为一维数组处理来解决这个问题。同样的方法也可以用在将二维数组作为结构体成员中。