求解线性方程组

一、写在前面

实验目的

(1) 熟悉求解线性方程组的有关理论和方法；
(2) 能编程实现雅可比及高斯-塞德尔迭代法、列主元高斯消去法、约当消去，追赶法
(3) 通过测试，进一步了解各种方法的优缺点
(4) 根据不同类型的方程组，选择合适的数值方法

实验内容

**本次实验参考公式**

二、实验过程
【参考代码】

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//求解线性方程组
#define N_EQUATIONS 3    //方程组个数，即未知元个数
#define N_DETERMINANT 4    //行列式阶数

double* Gausss_Seidel(double A[][N_EQUATIONS], double B[], double e, int N)
// 输入参数： A 系数矩阵， B 初始向量， e 误差限， N 最大迭代次数
{
    int i, j;
    int k = 0;    //记录当前迭代次数
    double ACC;    //累加器
    double MAX = 0.0;
    double Y[N_EQUATIONS];
    double X[N_EQUATIONS];
    for(i=0; i<N_EQUATIONS; i++)   //初始化初值为0.0
    {
        X[i] = 0.0;
        Y[i] = 0.0;
    }
    printf("-------------------------Gauss_Seidel迭代过程--------------------------\n");
    printf("k      ");
    for(i=0; i<N_EQUATIONS; i++)
    {
        printf("x%d(k)         ", i+1);
    }
    printf("\n");
    printf("%d      ", k);
    for(i=0; i<N_EQUATIONS; i++)
    {
        printf("%lf      ", Y[i]);
    }
    printf("\n");
    while(k < N)
    {
        for(i=0; i<N_EQUATIONS; i++)
        {
            ACC = 0.0;
            for(j=0; j<N_EQUATIONS; j++)
            {
                if(j != i)
                    ACC += A[i][j] * Y[j];
            }
            Y[i] = (B[i] - ACC)/A[i][i];
        }
        k++;
        printf("%d      ", k);
        for(i=0; i<N_EQUATIONS; i++)
        {
            printf("%lf      ", Y[i]);
        }
        printf("\n");

        MAX = 0;
        for(i=0; i<N_EQUATIONS; i++)
        {
            if(fabs(Y[i]-X[i]) > MAX)
                MAX = fabs(Y[i]-X[i]);
        }
        if(MAX < e)      //判断是否达到精度要求
            return Y;
        for(i=0; i<N_EQUATIONS; i++)
        {
            X[i] = Y[i];
        }
    }
    return NULL;
}

double SolveDeterminant(double A[][N_DETERMINANT], int n)
{
    int a, b, i, j, k, m;
    int count = 0;    //记录交换次数
    double B[N_DETERMINANT][N_DETERMINANT], z = 1.0;

    for(i=0; i<n-1; i++)
    {
        if(A[i][i] == 0)
        {
            for(b=i+1; b<n; b++)
            if(A[b][i] != 0)
            {
                for(a=0; a<n; a++)
                {
                    B[i][a] = A[i][a];
                    A[i][a] = A[b][a];
                    A[b][a] = B[i][a];
                }
            }
            count++;
        }

        for(j=i+1; j<n; j++)
        {
            for(k=n-1; k>=0; k--)
            {
                A[j][k]=A[j][k]-(A[i][k]*A[j][i]/A[i][i]); //将行列式化为上三角形式
            }
        }
    }

    printf("Exchange times: %d\n", count);
    printf("上三角矩阵形式:\n");
    for(a=0; a<n; a++)
    {
        for(b=0; b<n; b++)
        {
            printf("%f   ", A[a][b]);
            if(b == n-1)
            printf("\n");
        }
    }
    for(m=0; m<n; m++)
        z = z*A[m][m];
    return pow(-1, count)*z;
}

int main()
{
    double A0[N_EQUATIONS][N_EQUATIONS] = {      //系数矩阵
    {10.0, -1.0, -2.0},
    {-1.0, 10.0, -2.0},
    {-1.0, -1.0, 5.0}
    };
    double B0[N_EQUATIONS] = {7.2, 8.3, 4.2};
/* 则：
   10*x1 - x2 - 2*x3 = 7.2
   -x1 + 10*x2 - 2*x3 = 8.3
   -x1 - x2 + 5*x3 = 4.2
 */
    double* answer;
    answer = Gausss_Seidel(&A0, B0, 0.00001, 20);

    printf("\n----------------------选主元高斯消去求行列式值-------------------------\n");
    double A[N_DETERMINANT][N_DETERMINANT] = {
    {3, -2, 1, 4},
    {-7, 5, -3, -6},
    {2, 1, -1, 3},
    {4, -3, 2, 8}
    };
    double value = SolveDeterminant(A, N_DETERMINANT);
    printf("the value is: %lf\n", value);
    return 0;
}

三、实验结果

四、写在后面

无