C语言线性方程组求解：LU分解、雅可比迭代、高斯-赛德尔迭代

最新推荐文章于 2024-04-28 22:20:09 发布
最新推荐文章于 2024-04-28 22:20:09 发布
阅读量1.2k 收藏 18
点赞数 11
分类专栏： C语言程序分享 文章标签： c语言 算法 c++
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/jk_lovers/article/details/135839534
版权

（注：该程序C和C++通用）

我们知道在Matlab中，矩阵的求解是非常方便的，对于Ax = b，只需要做x = A \ b的运算就可以得到x的结果，然而对于C语言，我们就需要寻找一些算法来实现线性矩阵的求解。

由于本学期刚学完数值分析，并且目前做的项目中刚好有需要求解线性方程组的需求，于是将求解线性方程组的三种方法（其实还有更多，我这里仅给出三种）：LU分解、雅可比迭代以及高斯-赛德尔迭代。这三个算法的基本原理就不过多赘述了，大家想进一步学习的可以网上搜，一大堆教程，而且原理也较为简单。

下面给出该三个算法对应的函数：

（1）LU分解：
float* LU(float* A, float* b)
{
	//LU分解求方程的简单解释
	//系数矩阵各阶顺序主子式det(Ai)≠0时，可以进行LU分解
	//A = LU → Ax = b → LUx = b → Ly = b (其中Ux = y)
	//在计算得到LU矩阵后，先计算Ly = b
	//计算得到y后计算Ux = y就可以得到方程的解
	int i, j, n;
	float temp = 0;
	//生成一个二维数组保存LU分解所得的系数
	float** mat = (float**)malloc(N * sizeof(float*));
	for (i = 0; i < N; i++)
	{
		*(mat + i) = (float*)malloc(N * sizeof(float));
		// 将申请的数组中的值全部置零，方便后续使用
		memset(*(mat + i), 0, N * sizeof(float));
	}

	//为了节约内存，这里没有同时申请x和y，而是仅申请x，但这样也够用了
	float* x = (float*)malloc(N * sizeof(float));
	memset(x, 0, N * sizeof(float));

	//下面首先对LU矩阵的第一行和第一列进行计算，因为这两行有规律，容易写出
	i = 0;
	//LU分解矩阵的第一行和系数矩阵是一致的
	for (j = 0; j < N; j++)
		mat[i][j] = A[i * N + j];
	//LU的第一列结果可以通过除以第一列第一个（即A[0][0]）计算得到
	j = 0;
	for (i = 1; i < N; i++)
		mat[i][j] = A[i * N + j] / A[0];
	//下面开始对LU其他位置的值进行求解
	for (i = 1; i < N; i++)
		for (j = 1; j < N; j++)
		{
			if (i <= j)
			{
				temp = 0;
				for (n = 0; n < i; n++)
					temp += mat[n][j] * mat[i][n];
				mat[i][j] = A[i * N + j] - temp;
			}
			else
			{
				temp = 0;
				for (n = 0; n < j; n++)
					temp += mat[n][j] * mat[i][n];
				mat[i][j] = (A[i * N + j] - temp) / mat[j][j];  //列上方若有对角值（即i = j），需要处于该对角值
			}
		}
	//首先计算Ly = b，这里是自上而下的计算
	//这里的x的就是y
	for (i = 0; i < N; i++)
	{
		temp = 0;
		for (j = 0; j < i; j++)
		{
			temp += mat[i][j] * x[j];
		}
		x[i] = b[i] - temp;
	}
	//然后计算Ux = y，这里是自下而上的计算
	for (i = N-1; i >= 0; i--)
	{
		temp = 0;
		for (j = N - 1; j > i; j--)
		{
			temp += mat[i][j] * x[j];
		}
		x[i] = (x[i] - temp)/mat[i][i];
	}

	//free(x);（×）     //这里切记不要释放x，否则计算结果会消失，可以在之后的主函数中释放
	free(mat);  //释放内存
	return x;   //将计算结果返回
}

LU分解可以获得线性方程的精确解

（2）雅可比迭代算法
//雅可比迭代函数
float* Jacobi(float* A, float* b)
{
	int i, j;
	int iter_num = 0;   //目前的迭代次数
	float error_now = 100;   //目前迭代下的误差

	//生成数组用于保存每次迭代得到的x结果并保存上一次迭代结果至x0
	float* x = (float*)malloc(N * sizeof(float));
	memset(x, 0, N * sizeof(float));
	float* x0 = (float*)malloc(N * sizeof(float));
	memset(x0, 0, N * sizeof(float));

	while(error_now > ERROR)
	{
		//将上一次迭代的结果保存到x0中去
		memcpy(x0, x, N * sizeof(float));

		for (i = 0; i < N; i++)
		{
			float sum = 0;
			for (j = 0; j < N; j++)
			{
				//printf("%f\t", A[i * N + j]);
				if (i != j)
					sum += A[i * N + j] * (x0[j]);
			}
			x[i] = (b[i] - sum) / A[i * N + i];
		}

		error_now = current_error(x, x0);
		iter_num++;

		printf("当前迭代次数为:%d，迭代误差为:%f\n", iter_num, error_now);
		if (iter_num > IterMaxNum)
		{
			printf("迭代次数超过最大值要求，返回最后一次迭代结果");
			break;
		}
	}

	free(x0);
	return x;
}

雅可比迭代可以获得线性方程组的近似解，近似程度和我们所设置ERROR误差参数有关，越小接近程度越好，但迭代的次数就会越多。

其中的current_error函数如下，用于计算前后两次迭代之间的误差，当误差小于我们所设置的ERROR参数时，迭代就会结束，返回迭代结果。

float current_error(float* last, float* now)
{
	float max_error = 0;
	float temp;
	for (int i = 0; i < N; i++)
	{
		temp = *(now + i) - *(last + i);
		//这里不用abs的原因：abs不知道为什么一直报错，说是函数重载
		if (temp < 0)
			temp = -temp;
		if (temp > max_error)
			max_error = temp;
	}
	return max_error;
}

这里我用的是Visual Studio 2022版本，我引入头文件<stdlib.h>后abs函数一直报错，不知道为什么，于是我就换了一种写法，多了一步if负数判断。

（3）高斯-赛德尔迭代法
//高斯-赛德尔迭代函数
float* GS(float* A, float* b)
{
	int i, j;
	int iter_num = 0;   //目前的迭代次数
	float error_now = 100;   //目前迭代下的误差

	//生成数组用于保存每次迭代得到的x结果并保存上一次迭代结果至x0
	float* x = (float*)malloc(N * sizeof(float));
	memset(x, 0, N * sizeof(float));
	float* x0 = (float*)malloc(N * sizeof(float));
	memset(x0, 0, N * sizeof(float));

	while (error_now > ERROR)
	{
		memcpy(x0, x, N * sizeof(float));

		for (i = 0; i < N; i++)
		{
			float sum = 0;
			for (j = 0; j < N; j++)
			{
				//printf("%f\t", A[i * N + j]);
				if (i != j)
					sum += A[i * N + j] * (x[j]);
			}
			x[i] = (b[i] - sum) / A[i * N + i];
		}

		error_now = current_error(x, x0);  //获得当前迭代的误差结果
		iter_num++;   //迭代次数加1

		printf("当前迭代次数为:%d，迭代误差为:%f\n", iter_num, error_now);    //输出当前迭代信息
		if (iter_num > IterMaxNum)
		{
			printf("迭代次数超过最大值要求，返回最后一次迭代结果");
			break;
		}
	}

	free(x0);
	return x;
}

高斯-赛德尔迭代法的程序和上面的雅可比迭代法的程序是非常相似的，但高斯-赛德尔是雅可比的改进算法，比雅可比有更快的收敛速度（即迭代次数更少）。

（4）一个使用例程

下面我给出了一个完整的使用例程，利用#ifdef函数来判断你想运行的算法，只需要在前面把你想运行的方法取消注释，不用的方法打上注释即可。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#define N 4   //要求解的数组大小N x N
#define ERROR 0.00001  //误差要求
#define IterMaxNum 50000   //最大迭代次数

#define iter_Jacobi   //运行雅可比迭代（非常近似的解）
//#define iter_GS    //运行高斯-赛德尔迭代（非常近似的解）
//#define LU_Count   //运行LU分解求解方程解（精确解）

float current_error(float* last, float* now)
{
	float max_error = 0;
	float temp;
	for (int i = 0; i < N; i++)
	{
		temp = *(now + i) - *(last + i);
		//这里不用abs的原因：abs不知道为什么一直报错，说是函数重载
		if (temp < 0)
			temp = -temp;
		if (temp > max_error)
			max_error = temp;
	}
	return max_error;
}

//雅可比迭代函数
float* Jacobi(float* A, float* b)
{
	int i, j;
	int iter_num = 0;   //目前的迭代次数
	float error_now = 100;   //目前迭代下的误差

	//生成数组用于保存每次迭代得到的x结果并保存上一次迭代结果至x0
	float* x = (float*)malloc(N * sizeof(float));
	memset(x, 0, N * sizeof(float));
	float* x0 = (float*)malloc(N * sizeof(float));
	memset(x0, 0, N * sizeof(float));

	while(error_now > ERROR)
	{
		//将上一次迭代的结果保存到x0中去
		memcpy(x0, x, N * sizeof(float));

		for (i = 0; i < N; i++)
		{
			float sum = 0;
			for (j = 0; j < N; j++)
			{
				//printf("%f\t", A[i * N + j]);
				if (i != j)
					sum += A[i * N + j] * (x0[j]);
			}
			x[i] = (b[i] - sum) / A[i * N + i];
		}

		error_now = current_error(x, x0);
		iter_num++;

		printf("当前迭代次数为:%d，迭代误差为:%f\n", iter_num, error_now);
		if (iter_num > IterMaxNum)
		{
			printf("迭代次数超过最大值要求，返回最后一次迭代结果");
			break;
		}
	}

	free(x0);
	return x;
}

//高斯-赛德尔迭代函数
float* GS(float* A, float* b)
{
	int i, j;
	int iter_num = 0;   //目前的迭代次数
	float error_now = 100;   //目前迭代下的误差

	//生成数组用于保存每次迭代得到的x结果并保存上一次迭代结果至x0
	float* x = (float*)malloc(N * sizeof(float));
	memset(x, 0, N * sizeof(float));
	float* x0 = (float*)malloc(N * sizeof(float));
	memset(x0, 0, N * sizeof(float));

	while (error_now > ERROR)
	{
		memcpy(x0, x, N * sizeof(float));

		for (i = 0; i < N; i++)
		{
			float sum = 0;
			for (j = 0; j < N; j++)
			{
				//printf("%f\t", A[i * N + j]);
				if (i != j)
					sum += A[i * N + j] * (x[j]);
			}
			x[i] = (b[i] - sum) / A[i * N + i];
		}

		error_now = current_error(x, x0);  //获得当前迭代的误差结果
		iter_num++;   //迭代次数加1

		printf("当前迭代次数为:%d，迭代误差为:%f\n", iter_num, error_now);    //输出当前迭代信息
		if (iter_num > IterMaxNum)
		{
			printf("迭代次数超过最大值要求，返回最后一次迭代结果");
			break;
		}
	}

	free(x0);
	return x;
}

//利用LU分解求解方程解
float* LU(float* A, float* b)
{
	//LU分解求方程的简单解释
	//系数矩阵各阶顺序主子式det(Ai)≠0时，可以进行LU分解
	//A = LU → Ax = b → LUx = b → Ly = b (其中Ux = y)
	//在计算得到LU矩阵后，先计算Ly = b
	//计算得到y后计算Ux = y就可以得到方程的解
	int i, j, n;
	float temp = 0;
	//生成一个二维数组保存LU分解所得的系数
	float** mat = (float**)malloc(N * sizeof(float*));
	for (i = 0; i < N; i++)
	{
		*(mat + i) = (float*)malloc(N * sizeof(float));
		// 将申请的数组中的值全部置零，方便后续使用
		memset(*(mat + i), 0, N * sizeof(float));
	}

	//为了节约内存，这里没有同时申请x和y，而是仅申请x，但这样也够用了
	float* x = (float*)malloc(N * sizeof(float));
	memset(x, 0, N * sizeof(float));

	//下面首先对LU矩阵的第一行和第一列进行计算，因为这两行有规律，容易写出
	i = 0;
	//LU分解矩阵的第一行和系数矩阵是一致的
	for (j = 0; j < N; j++)
		mat[i][j] = A[i * N + j];
	//LU的第一列结果可以通过除以第一列第一个（即A[0][0]）计算得到
	j = 0;
	for (i = 1; i < N; i++)
		mat[i][j] = A[i * N + j] / A[0];
	//下面开始对LU其他位置的值进行求解
	for (i = 1; i < N; i++)
		for (j = 1; j < N; j++)
		{
			if (i <= j)
			{
				temp = 0;
				for (n = 0; n < i; n++)
					temp += mat[n][j] * mat[i][n];
				mat[i][j] = A[i * N + j] - temp;
			}
			else
			{
				temp = 0;
				for (n = 0; n < j; n++)
					temp += mat[n][j] * mat[i][n];
				mat[i][j] = (A[i * N + j] - temp) / mat[j][j];  //列上方若有对角值（即i = j），需要处于该对角值
			}
		}
	//首先计算Ly = b，这里是自上而下的计算
	//这里的x的就是y
	for (i = 0; i < N; i++)
	{
		temp = 0;
		for (j = 0; j < i; j++)
		{
			temp += mat[i][j] * x[j];
		}
		x[i] = b[i] - temp;
	}
	//然后计算Ux = y，这里是自下而上的计算
	for (i = N-1; i >= 0; i--)
	{
		temp = 0;
		for (j = N - 1; j > i; j--)
		{
			temp += mat[i][j] * x[j];
		}
		x[i] = (x[i] - temp)/mat[i][i];
	}

	//free(x);（×）     //这里切记不要释放x，否则计算结果会消失，可以在之后的主函数中释放
	free(mat);  //释放内存
	return x;   //将计算结果返回
}

void main()
{
	//这里给的例子为：
	//8a - 3b + 2c = 20
	//4a + 11b - c = 33
	//6a + 3b + 12c = 36
	//float A[N][N] = { {8,-3,2},{4,11,-1},{6,3,12} };  //要求解的矩阵
	//float b[N] = { 20,33,36 };
	//该方程精确解为 a = 3 、b = 2 、c = 1
	
	//另外一个例子
	float A[N][N] = { {1,0,2,0},{0,1,0,1},{1,2,4,3},{0,1,0,3} };   //要求解的矩阵
	float b[N] = { 5,3,17,7 };
	//该方程精确解为 a = 1 、b = 1 、c = 2 、d = 2

	float* x;
#ifdef LU_Count
	x = LU((float*)A, b);
#endif

#ifdef iter_Jacobi
	x = Jacobi((float*)A, b);
#endif
	
#ifdef iter_GS
	x = GS((float*)A, b);
#endif
	for (int i = 0; i < N; i++)
		printf("%f\n", x[i]);
	free(x);   //释放x
}
（5）小结

总之，想要精确解就LU算法，想要程序运行更快就选择迭代法（当然，前提是需要满足其分解条件和迭代条件）。

乐&韵
关注
  • 11
    点赞
  • 踩
  • 18
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
专栏目录
Jacobi-雅可比迭代法、高斯-赛德尔迭代法.zip
10-23
Jacobi-雅可比迭代法与高斯-赛德尔迭代法迭代次数可自行设置
C语言实现雅克比迭代法求解Ax=b
qq_63462564的博客
05-10 524
【代码】C语言实现雅克比迭代法求解Ax=b。
迭代法求解线性方程 （简单迭代法 雅可比(Jacobi)迭代法 高斯-塞德尔(Gauss-Seidel)迭代法 逐次超松弛(SOR)迭代法） [数值计算方法](c语言描述版)
最新发布
m0_74431077的博客
04-28 2277
迭代法求解线性方程，超详细的过程和讲解，附有完整可运行代码！！ 一、迭代法的原理 二、雅可比(Jacobi)迭代法 三、高斯-塞德尔(Gauss-Seidel)迭代法 四、逐次超松弛(SOR)迭代法
雅可比（jacobi）迭代法，c语言实现
热门推荐
梦载代码
10-13 1万+
#include &lt;iostream&gt; #include &lt;cstdio&gt; #include &lt;cmath&gt; #include &lt;cstring&gt; #include &lt;bits/stdc++.h&gt; using namespace std; double A[1010][1010]; double B[1010]; double X[10...
使用迭代法（雅可比Jacobi与高斯-赛德尔Gausss-Seidel迭代法）和直接求解法（消元法）求解线性方程组（C++）
Huanyuaner的博客
10-27 3512
分别使用迭代法（Jacobi与Gausss-Seidel迭代法）和直接求解法（消元法）求解线性方程组： 要求：使误差不超过。
C语言实现雅可比迭代（计算方法）
qq_43452446的博客
06-01 9836
#include <stdio.h> #include <stdlib.h> #define MAX 10 double A[MAX][MAX]; double b[MAX]; double X[MAX]; double Y[MAX]; double X0[MAX]; int NUM; double epsilon=1e-5; void Input_Matrix()//...
使用C语言编写一个Jacobi迭代程序
weixin_42610671的博客
02-14 811
Jacobi迭代是一种解决线性方程组的数值方法，用C语言编写Jacobi迭代程序的基本步骤如下： 定义矩阵和向量：在C语言中，可以使用二维数组存储矩阵，一维数组存储向量。 输入矩阵和向量：让用户输入矩阵和向量的元素值。 初始化解向量：选择一个初始解向量，它的大小应该和方程的未知数数量相同。 进行迭代：对于每一个未知数，利用Jacobi迭代公式计算新的近似解。 检验迭代是否收敛：检查新的解...
用C语言实现雅可比迭代法.txt
02-15
用C语言实现雅可比迭代法
LU分解原理解释及C语言实现
Aysek的博客
08-26 4172
本文为笔者在编写spice时，需利用LU求解线性方程组时所写。 1、杜尔里特算法： 从下至上对矩阵A做初等行变换，将对角线左下方的元素编成零，即得到U矩阵。然后可以证明这些行变换的效果等同于左乘一系列单位下三角矩阵，这一系列单位下三角矩阵的乘积的逆就是L矩阵，它也是一个单位下三角矩阵。即杜尔特里算法证实了矩阵A可以唯一分解成为L和U矩阵。 2、LU分解原理概要： 将方程组的系数矩阵分解为两个三角矩阵的乘积，以四阶为例进行讨论，即 ...
LU矩阵分解法解线性方程组（C++实现）
py753159123456的博客
10-28 1593
LU矩阵分解法解线性方程组（C++实现） #include<cstdio> #include<iostream> using namespace std; void fenjie(double a[],double b[],int n){ double aa[n][n]; for(int i=0;i<n;i++){ for(int j=0;j<n;j++){ aa[i][j]=a[n*i+j]; } } double l[n][n]; do
LU分解法解线性方程组（C语言）
04-30
使用LU分解法，来解线性方程组，用C语言写的源程序
高斯赛德尔迭代算法 C语言
06-09
迭代法是一种逐次逼近的方法,与直接法(高斯消元法)比较, 具有: 程序简单,存储量小的优点。特别适用于求解系数矩阵为大型稀疏矩阵的方程组。常用迭代方法：雅可比迭代，高斯-赛德尔迭代，松弛迭代等。
雅克比迭代法和高斯-赛德尔法解线性方程组(C++).docx
07-05
雅克比迭代法和高斯-赛德尔法解线性方程组(C++) 雅克比迭代法和高斯-赛德尔法是两种常用的线性方程组解法，分别用于解线性方程组。雅克比迭代法是一种简单的迭代法，通过不断迭代来逼近解，而高斯-赛德尔法则是...
雅克比迭代法，高斯-赛德尔迭代法，SOR法，追赶法求解线性方程组
11-27
○1高斯-赛德尔迭代法比雅克比的迭代次数少一半多，说明高斯-赛德尔迭代法优于雅克比迭代法，收敛速度更快，效率更高。但也不绝对，有时也会出现雅克比迭代法收敛而高斯-赛德尔迭代法不收敛的情况。 ○2SOR法的迭代...
MATLAB雅克比迭代法（Jacobi method）和高斯-赛德尔迭代法（Gauss–Seidel method）求Ax=b
12-26
包含代码和文档 分别用雅克比迭代法和高斯-赛德尔迭代法求解方程组，精确到小数点后6位，分别就 给出相应的计算结果
MATLAB解非线性方程组 雅可比 迭代
03-09
使用牛顿方法解非线性方程组 雅可比 迭代 jacobi matlab代码
c# lu分解的代码_矩阵LU分解分块算法实现
weixin_39741459的博客
12-22 626
本文主要描述实现LU分解算法过程中遇到的问题及解决方案，并给出了全部源代码。1. 什么是LU分解？矩阵的LU分解源于线性方程组的高斯消元过程。对于一个含有N个变量的N个线性方程组，总可以用高斯消去法，把左边的系数矩阵分解为一个单位下三角矩阵和一个上三角矩阵相乘的形式。这样，求解这个线性方程组就转化为求解两个三角矩阵的方程组。具体的算法细节这里不做过多的描述，有很多的教材和资源可以参考。这里推荐的参...
C语言实现高斯-赛德尔迭代gauss seidel(附完整源码)
希望我的博客，能帮上你解决学习中工作中所遇到的问题
03-03 2742
实现C语言高斯-赛德尔迭代gauss seidel高斯-赛德尔迭代gauss seidel的完整源码（实现，main函数测试） 高斯-赛德尔迭代gauss seidel的完整源码（实现，main函数测试） #include <math.h> #include <stdio.h> int main() { float a, b, c, a1, a2, a3, b1, b2, b3, c1, c2, c3, d1, d2, d3, x1, x2, x3; printf
c语言编程实现高斯迭代,用C语言实现高斯-赛德尔迭代方法
weixin_34957608的博客
05-21 1585
解决线性方程组的好方法用C语言实现高斯-赛德尔迭代方法， 源程序代码：#include "stdlib.h"#include "stdio.h"#include "conio.h"#include "string.h"#include "math.h"#define N 100float Table(int n,float a[N][N],float b[N]){int i,j;float c[N...
线性方程组的迭代解法，雅可比迭代法和高斯-赛德尔迭代法解方程组的收敛性
11-25
线性方程组的迭代解法是一种通过迭代逼近线性方程组的解的方法。雅可比迭代法和高斯-赛德尔迭代法是两种常见的迭代解法。 雅可比迭代法是一种基本的迭代法，其基本思想是将线性方程组的每个未知数的迭代公式表示为该未知数的系数与其他未知数的当前迭代值的乘积之和，然后通过不断迭代计算每个未知数的值，直到满足一定的收敛条件为止。雅可比迭代法的收敛性取决于系数矩阵的特征值，当系数矩阵的特征值都小于1时，雅可比迭代法才能收敛。 高斯-赛德尔迭代法是在雅可比迭代法的基础上进行改进的一种迭代法，其基本思想是在计算每个未知数的值时，利用已经计算出的未知数的最新值，以加快收敛速度。高斯-赛德尔迭代法的收敛性也取决于系数矩阵的特征值，但相比雅可比迭代法，高斯-赛德尔迭代法的收敛速度更快，收敛的条件也更宽松。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值