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

ndzhzsdw · 发表于 2015-12-8 22:51:42

/* Created by mazirong */
/*2015-12-8 */
1 高斯消去法

1.1 基本思想及计算过程

高斯（Gauss）消去法是解线性方程组最常用的方法之一，它的基本思想是通过逐步消元，把方程组化为系数矩阵为三角形矩阵的同解方程组，然后用回代法解此三角形方程组得原方程组的解。

为便于叙述，先以一个三阶线性方程组为例来说明高斯消去法的基本思想。

把方程（IV）乘（）后加到方程（I）上去，把方程（IV）乘（）后加到方程（II）上去，把方程（IV）乘（）后加到方程（III）上去,即可消去方程（I）、（II）、（III）中的x1，得同解方程组

将方程（III）乘（）后加于方程（II），将方程（III）乘（）后加于方程（I），得同解方程组：

将方程（II）乘（）后加于方程（I）,得同解方程组：

由回代公式（3.5）得，x3 = 1，x2 = 1，x1 = 1。

1.2 编程思路

高斯消元法解线性方程组，算法很简单，但过程很复杂，我在网上几乎没看到过正确的高斯消元法C程序，有程序也是没有注释，看起来非常费力。于是我硬着头皮自己来写，花了1天时间终于完工了。以求如下行列式为例详解编程思路：

1.把第1行进行优化（优化成第1行第1列元素为1的行列式），然后用来与第2至第4行的第1列元素消为0，如图1：

图1

2.保持第1行不变。对图一中第2行元素进行优化（优化成第2行第2列元素为1的行列式），然后用来与第3行至第4行的第2列元素消为0，如图2：

图2

3.保持前2行不变。对图二中第3行元素进行优化（优化成第3行第3列元素为1的行列式），然后用来与第4行的第3列元素消为0，如图3：

图3

4.保持前3行不变。对图三中第4行元素进行优化（优化成第4行第4列元素为1的行列），如图四：

图4

5.通过回代分别求出x1、x2、x3、x4的值，如图5：

图5

1.3 具体实现过程

由于我的电脑上没有安装Visual Studio，于是使用单片机来实现，通过将程序烧录到单片机中，然后将运算的数据结果通过串口发送给上位机打印输出。运行效果如图6：

如图6

程序主要有3个文件组成：gs_elim.h，gs_elim.c，main.c。程序不限于平台，稍稍修改就能移植到其他开发环境中。

gs_elim.h

#ifndef __GS_ELIM_H

#define __GS_ELIM_H

#include "usart.h"

// 定义矩阵为N阶

#define N 4

void Gs_Elim(void);

#endif

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]);

}

main.c

#include "led.h"

#include "key.h"

#include "bsp_lcd.h"

#include "delay.h"

#include "sys.h"

#include "usart.h"

//#include "includes.h"

#include "malloc.h"

#include "myiic.h"

#include "24cxx_iic.h"

#include "oled_iic.h"

#include "gs_elim.h"

int main(void)

{

delay_init();      // 延时函数初始化

NVIC_Configuration(); // 设置NVIC中断分组2:2位抢占优先级，2位响应优先级

uart_init(9600); // 串口初始化为9600

LED_Init();

KEY_Init();

LCD_Init();

// 高斯消元函数

Gs_Elim();

while(1);

}

最后使用ptc_mathcad直接计算一个矩阵：