【原创】高斯消元法--基于STM32串口调试（下）

hkhkdyx

/* Created by mazirong */
/*2015-12-8 */

gs_elim.c
#include "gs_elim.h"

// k：提取原数组特定位置值，防止运算一次被覆盖
// m：优化前一行矩阵
double k,m;

// i：行
// j：列
// t：消元次数
// 注意：数值是以0开始，所以实际值应该加1
int i,j,t;

// 定义一个N行*N列的矩阵A
// |1   2   1   2|
// |3   5   2   1|
// |2   3  -1   1|
// |4  -1   3  -2|
double A[N][N] = { 1,  2,  1,  2,
                               3,  5,  2,  1,
                                   2,  3, -1,  1,
                                   4, -1,  3, -2};

// 定义矩阵的增广部分B
// |16|
// |23|
// |9 |
// |3 |
double B[N] = {16, 23, 9, 3};

// 未知数X数组
double X[N+1] = {0};


// 高斯消元法函数
void Gs_Elim(void)
{
        /********************显示待解的增广矩阵**********************/
        printf("待解的增广矩阵为:\n");
        for(i=0;i<N;i++)
        {   
                // 循环打印第i行系数
                for(j=0;j<N;j++)
                printf("%12lf  ",A[j]);                        // %12lf：输出场宽为12的浮点数,其中小数点占1位,数据右对齐。
                // 打印第i行增广数
                printf("|%12lf\r\n",B);               
        }
        printf("\r\n");

        /*******************消元化简*********************/
        do{
                printf("第%d次消元结果：\r\n", t+1);

                /**********显示第(1)行至第(t)行（若t=0不执行）*****************/
                for(i=0; i<t; i++)
                {
                        for(j=0; j<N; j++)
                        {
                                printf("%12lf  ", A[j]);
                        }
                        printf("|%12lf\r\n", B);
                }

                /****************显示优化后的第(t+1)行************************/
                k = A[t][t];                                                // 把A[t][t]先提取出来，防止运算一次被数组覆盖
                for(j=0; j<N; j++)
                {
                        A[t][j] = A[t][j] / k;                          // 得到优化后的第(t+1)行矩阵值
                        printf("%12lf  ", A[t][j]);
                }
                B[t] = B[t] / k;                                        // n为第t+1行的增广部分
                printf("|%12lf\r\n", B[t]);                        // 显示第(t+1)行的增广部分

                /****************显示消元后的第(t+2)行至第(N)行****************/
                for(i=t+1; i<N; i++)
                {
                        k = A[t];                                        // 把A[t]先提取出来，防止运算一次被数组覆盖                                
                        for(j=0; j<N; j++)                                // 循环输出消元后的第1至N列
                        {
                                m = A[t][j] / A[t][t];                // 第(t+1)行元素同时除以A[t][t]   
                                A[j] = A[j] - k*m;        // 消除运算及覆盖原数组
                                printf("%12lf  ", A[j]);        // 显示第(i+1)行第(j+1)列数组值
                        }
                        B = B - B[t]*k;                        
                        printf("|%12lf\r\n", B);                // 显示第(i+1)行的增广部分
                }

                /**********************************************************/
                t++;
        }while(t<N);
        printf("\r\n");

        /***********************回代求X***********************************/
        //        X[4] = B[N-1];
        //        X[3] = B[N-2] - A[N-2][N-1]*X[4];
        //        X[2] = B[N-3] - A[N-3][N-1]*X[4] - A[N-3][N-2]*X[3];
        //        X[1] = B[N-4] - A[N-4][N-1]*X[4] - A[N-4][N-2]*X[3] - A[N-4][N-3]*X[2];
        for(i=N-1;i>=0;i--)
        {
                X[i+1] = B;
                for(j=0;j<(N-i-1);j++)
                {
                        X[i+1] -= A[N-j-1]*X[4-j];
                }
        }
        printf(" X[1] = %lf  ", X[1]);
        printf(" X[2] = %lf  ", X[2]);
        printf(" X[3] = %lf  ", X[3]);
        printf(" X[4] = %lf  ", X[4]);
}

fire

帮顶     

hkhkdyx

fire 发表于 2015-12-13 11:35
帮顶

谢谢火哥