高斯列主消元法 求非齐次线性方程组 C语言实现代码
高斯列主元素消去法是由高斯消去法改进的算法
下面浅浅分享一下本人对该方法的理解
Ax = b
先说高斯消去法，感觉基本的思路就跟我们手算非齐次线性方程组差不多，在线性代数中，我们求解方程组都是这种思路，消元的过程相当于是，由系数矩阵A和非齐次项b得到的增广矩阵做行变换，化为行阶梯型，最后在由下往上回代，求出每一个未知数的过程。
举例如下所示：
由系数矩阵和非齐次项拼成的增广矩阵如下所示
我们作初等行变化将其化为行阶梯型如下
到上一步之后我们就完成了消元的过程
解出未知数一步步回代就可以得到解向量如下
当然在数学上，这样的非线性方程组要有解，必须要求系数矩阵的秩等于增广矩阵的秩，不过数学上的东西就不在本篇讨论范围内了。
上面叙述的步骤在计算机程序中是很容易实现的，但考虑下面的矩阵消元。
可以看到，问题就在于，第一行第一列的元素为0，是无法消去它下方的任何数在，在计算机中也会因为0作除数而无法求解。而且就算这样的主对角线上的元素不为0，但如果元素很小的话，就会因为它是除数，而引起其他元素数量级及舍入误差的急剧增大，从而导致计算结果不可靠。
这种情况下，就产生了列主元素消去法，它的基本思想是在每步消元时，在第k列主对角线元素及其下方的元素中选取绝对值最大的元素作为主元，进行该步的消去。即体现为一个行列互换。
以上面的矩阵距离，我们用列主元素消去法逐步化简，你便会明白具体的操作步骤。
我们先将第一行第一列的元素0与它下方列所有的元素作比较，即比较0,1,2三个元素，找出他们中绝对值最大的元素，即2，我们就以2作为主元进行消去步骤，即将第①行与第③行互换，得
然后用2将下方的元素消为0，得
下面是第二行第二列的元素的消去，即3.5，与其下方的元素作比较，注意是下方元素，即不包括3，只有3.5与-3，进行绝对值比较，3.5已经是绝对值最大的元素了，所以不需要行变换，就以3.5做主元进行消去。
依次循环上面的步骤，直到完成全部主元的消去，之后便是和高斯消去法一样的回代求解过程了。
可以看出高斯列主元素消去法的思想是很简单但却十分有效的。
——————————————————————————————
下面给出我用C语言实现的代码，主要函数为EMCP()函数，大家可根据需要自行修改或使用
#include <stdio.h>
#include <stdlib.h>/*———————————————— 说明：———————————————— 高斯列主消元法主函数为 EMCP()函数原型为：Double * EMCP(void * mat,double * vec)具体用法见main()函数使用时将main函数前的所有函数复制到自己程序的main函数前面 根据main中说明的使用方法使用本函数复制时 所有前置函数 缺一不可，且顺序不可随意调换最后一个PrintMat()是用来打印二维数组的，可要可不要 ————————————————@猿力觉醒
*///前置函数 绝对值函数
double _Abs(double value){if(value < 0)return -value;else{return value;}
}//前置函数 获取矩阵对应元素地址
double * _Get(void * mat,int i,int j,int N){return ((double*)mat + i*N) + j;
} //前置函数 行变换 交换矩阵的第a行与第b行
void _MatSwap(void * mat,int a,int b,int N){double temp;for(int i=0;i<N;i++){temp = *_Get(mat,a,i,N);*_Get(mat,a,i,N) = *_Get(mat,b,i,N);*_Get(mat,b,i,N) = temp;}
}
void _VecSwap(double * vec,int a,int b,int N){double temp;temp = vec[a];vec[a] = vec[b];vec[b] = temp;
}
void _Swap(void * mat,double * vec,int a,int b,int N){_MatSwap(mat,a,b,N);_VecSwap(vec,a,b,N);
}//前置函数 行变换 将矩阵的第a行的k倍加到第b行
void _MatIk2J(void * mat,int a,double k,int b,int N){double temp;for(int i=0;i<N;i++){temp = *_Get(mat,a,i,N) * k;*_Get(mat,b,i,N) += temp;}
}
void _VecIk2J(double * vec,int a,double k,int b,int N){double temp;temp = k * vec[a];vec[b] += temp;
}
void _Ik2J(void * mat,double * vec,int a,double k,int b,int N){_MatIk2J(mat,a,k,b,N);_VecIk2J(vec,a,k,b,N);
}//高斯列主消元法 函数主体
double * EMCP(void * mat,double * vec,int N){//消元 for(int i=0;i<N-1;i++){//寻找列主元 int maxEle = i;//列主元索引 for(int k=i+1;k<N;k++){if(_Abs(*_Get(mat,k,i,N)) > _Abs(*_Get(mat,maxEle,i,N)))maxEle = k;}if(*_Get(mat,i,i,N) == 0.0)return NULL;//某列主元为0，方程无解，返回NULL else{_Swap(mat,vec,maxEle,i,N);//行交换 }for(int j=i+1;j<N;j++){//将列主元下的元素全部消为0 _Ik2J(mat,vec,i,-((*_Get(mat,j,i,N))/(*_Get(mat,i,i,N))),j,N); } }double * v_out = (double*)malloc(sizeof(double)*N);//回代 for(int i=N-1;i>=0;i--){//解出最下方可解的元 v_out[i] = vec[i]/(*_Get(mat,i,i,N));//对于已经解出的元，从vec中消去该元 for(int j=i-1;j>=0;j--){vec[j] -= v_out[i]*(*_Get(mat,j,i,N));}}return v_out;
}//打印输出矩阵
void PrintMat(void * mat,int N){for(int i=0;i<N;i++){for(int j=0;j<N;j++){printf("%.2lf",*_Get(mat,i,j,N));printf(" ");}printf("\n\n");}
}//以下为测试用main函数，示意函数使用方法
int main(void){//求解线性方程组 Ax = b 的问题（对应有限元 Ku = F 的求解） // mat为系数矩阵，即为 A；vec为非齐次项，即为 b double mat[3][3] = {{1,2,3},{1,5,3},{2,3,1}};double vec[3] = {-1,-1,-1};// EMCP()函数会更改传入的系数矩阵mat与非齐次项vec // 故此处为了演示用法重定义了一模一样的mat与vec double _mat[3][3] = {{1,2,3},{1,5,3},{2,3,1}};double _vec[3] = {-1,-1,-1};printf("Answer of question: Ax = b\n");printf("\n\n");printf("A is:\n");PrintMat(mat,3);printf("\n");printf("b = \n");for(int i=0;i<3;i++){printf("%.2lf ",vec[i]);}printf("\n\n");//【关键处】 //下面的语句展示了调用EMCP的用法 // double * EMCP(mat,vec,N) 一共需要3个参数// 第一个参数mat为系数矩阵，传入C语言中的二维数组的名字即可// 第二个参数vec为非齐次项的列向量，传入C语言中的一维数组名字即可 // 第三个参数N为方程的阶数//// 使用时必须保证mat与vec维数一致 double * v_out = EMCP(_mat,_vec,3);//EMCP()返回一个双精度一维数组 // 即解向量数组// 上面 v_out 与一维数组等价 // 可以使用如 v_out[i] 的形式访问每个元素 //PrintMat(_mat,3);if(v_out){printf("Output Vector x:\n");for(int i=0;i<3;i++){printf("%.5lf ",v_out[i]);}printf("\n\n");printf("As the following equaltion:\n");for(int i=0;i<3;i++){printf("|");for(int j=0;j<3;j++){printf("%+4.2lf ",*_Get(mat,i,j,3));}printf("\b\b| [%+8.5lf] = |%+4.2lf|",v_out[i],vec[i]);printf("\n");}}else{printf("No Solution!");}
}高斯列主消元法 求非齐次线性方程组 C语言实现代码相关推荐
- 列主消元法解非奇异线性方程组的MATLAB程序
function [x] = LMain_elimination(A,b,prec,n) %% % 列主消元法求解线性代数方程组 Ax = b的MATLAB实现 % A为待求解方程组的系数矩阵(要求A ...
- 高斯消元：列主消元法
什么是列主消元 (注: akk代表第k行第k列的权值, 以下摘自百度百科:列主消元法) 列主元素消去法是为控制舍入误差而提出来的一种算法,列主元素消去法计算基本上能控制舍入误差的影响,其基本思想是:在 ...
- 用克拉默法则求非齐次线性方程组
用克拉默法则求非齐次线性方程组 先点击新建脚本,并保存一下,以由 {2x1+3x2+4x3=53x1+4x2+x3=36x1+x2+3x3=4\left\{\begin{array}{c} 2x_1+ ...
- 高斯消元法列主消元法
高斯消元法&&列主消元法 写在前面 高斯消元法 列主消元法 写在前面 最近在学数值分析,正好学到求解线性方程组.就自己动手简单实现了一下.关于本算法的原理可以在<数值分析> ...
- python 高斯约当消元法求逆矩阵
judge函数判断该矩阵该矩阵是否有逆矩阵 calculate计算逆矩阵 import sysclass MatrixInverse:""""求逆矩阵" ...
- c语言计算日出日落时间_利用日期、经纬度求日出日落时间 C语言程序代码
展开全部 #define PI 3.1415926 #include #include using namespace std; int days_of_month_1[]={31,28,31,30, ...
- 4.6 高斯约当消元法
4.6 高斯约当消元法 高斯消元法把矩阵变换为上三角阵,上三角阵还可以继续变换为对角阵.例如上面增广矩阵 [A,b][A, b][A,b] 变换为上三角阵 [24−2201140048]\left[ ...
- MATLAB求解非齐次线性方程组
根据线性代数中求解方程组的基本知识,首先应判断系数矩阵的秩是否和增广矩阵的秩相等,若不等,则无解:若有解,根据秩和未知量个数的关系,判断是唯一解还是无穷多解:若为无穷多解,其通解为齐次方程组的通解加非 ...
- 蒙特卡洛法求非线性方程组
蒙特卡罗法的概念及应用 蒙特卡洛法 (又称统计试验法)是描述装备运用过程中各种随机现象的基本方法,而且它特别适用于一些解析法难以求解甚至不可能求解的问题,因而在装备效能评估中具有重要地位. 用蒙特卡洛 ...
- 用列主元高斯(Gauss)消元法求n阶线性方程组的解（python）
(一)目的 通过设计.编制.调试2~3个求n阶线性方程组数值解的程序,加深对其数值计算方法及有关的基础理论知识的理解. (二)要求 用编程语言实现用高斯(Gauss)消元法求n阶线性方程组的解.用列主 ...
最新文章
- Windows2012R2服务器的安装与亮点功能介绍
- 建德有没有计算机培训,建德计算机培训,建德计算机培训班,建德计算机培训完好找工作吗 - IT教育频道...
- python数据类型所占字节数_python标准数据类型 Bytes
- oracle 如何迁移到 mysql_怎么将数据库从Oracle迁移到SQL Server，或从Oracle迁移到MySQL...
- cad致命错误如何处理_Golang 如何优雅地处理错误
- Python连接Oracle-常见问题
- 浏览器报错：unexpected end of input 解决方法
- 代码行数越少就越“简单”吗？
- 蓝桥杯C语言基础题---01字串
- Eclipse 切换 SVN 地址
- 【python】入门oj
- 耗资52亿美元，历时15年，人类有史以来建造的最复杂机器
- 写了四十篇办公自动化文章后，我整理了这十个常用操作，代码拿走就用！
- com.lbx:xTools
- TCP/IP网络编程之基于TCP的服务端/客户端（一）
- 卷入亿万骗局，他遭遇“死亡威胁”：如果没有光明，我愿成为火炬
- 算法题：用php生成excel列
- OpenGL之glut、glfw、glew、glad等库之间的关系
- js获取手机品牌型号,系统型号等
- 阿里云Flink sql 基本使用手册