之前介绍了C语言用代数余子式求行列式
本次开始介绍如何用公式法对矩阵求逆,并用C语言将其实现。
之前程序有点小bug,已于2022年11月29日修改。
更新:
伴随法只适合求低阶矩阵的逆,对于相对高阶(20维以上)对矩阵求逆用高斯法求解效率更高,此外本文中使用了_msize函数用于判断内存维数,但该函数只适合winodows系统,Linux和Mac系统无法使用(笔者也是在用了Mac系统后才发现),对于上述两个问题,您应该可以在:C语言求矩阵的逆(高斯法)得到满意的答案。
如果矩阵接近奇异值,求逆的数值将不稳定,那么使用C语言LU分解法求逆将会得到更好的效果。
目录
数学原理
矩阵求逆的公式
数乘矩阵
程序设计
整体代码
求行列式的值(Det、Cof)
求伴随矩阵(FindCof)
求逆矩阵的主函数(matrix_inver)
测试
数学原理
矩阵求逆的方法有很多种,本次主要介绍如何利用公式法求解。
矩阵求逆的公式
其中,为A的行列式的值,为矩阵A的伴随矩阵。
伴随矩阵的表达方式为:
其中,为代数余子式,代数余子式与余子式的关系为:
是矩阵A去掉i行j列,剩下元素重新组成的新矩阵 。
数乘矩阵
假设矩阵A为:
如果一个数b乘以一个矩阵A,那么A中的每一个元素都要乘以b,其表示如下:
程序设计
为了后面方便调试,先利用malloc函数创建一个矩阵,并为其赋值。关于malloc函数有不理解的可以参考之前写的C语言动态内存管理。
malloc使用的基本方式:
- void* malloc (size_t size);
- void free (void* ptr);
- size是指定的开辟内存的大小,单位是字节
- size_t的无符号整型则限制程序员误操作开辟负字节的空间
- 如果开辟成功,malloc会返回一个void*类型的指针
- 如果开辟失败,则返回的是空指针,所以在malloc之后需要对指针进行检查
- 当malloc的东西不再使用时,需要free对其进行释放,否则会造成内存泄漏
- malloc和free均需要包含头文件<stdlib.h>
//创建n维矩阵空间,并初始化
double** test1(int n)
{
double** arr = (double**)malloc(sizeof(double*) * n);
int i, j;
if (arr != NULL)
{
for (i = 0; i < n; i++)
{
arr[i] = (double*)malloc(sizeof(double) * n);
}
//为矩阵赋值
if (*arr != NULL)
{
for (i = 0; i < n; i++)
{
for (j = 0; j < n; j++)
{
arr[i][j] = pow(i, j);
}
}
}
}
return arr;
}
其中:
- double** arr = (double**)malloc(sizeof(double*) * n);表示开辟n个double*类型的空间,将其地址赋给arr
- 然后再开辟n个double类型的空间,将其地址赋给double**中每个double*
- 配合下图更好理解
整体代码
#include<stdio.h>
#include<assert.h>
#include<stdlib.h>
#include<math.h>
#include<malloc.h>
#include<string.h>
#define MatMax 20
//函数声明
double Det(const double arr[MatMax][MatMax], int n);
double Cof(const double arr[MatMax][MatMax], int i, int n);
void FindCof(double arr[MatMax][MatMax], double arr2[MatMax][MatMax], int i, int j, int n);
double** matrix_inver(double** arr);
//计算行列式
double Det(const double arr[MatMax][MatMax], int n)
{
assert(n > 0);
double sum = 0;
int i = 0;
if (n == 1)//1阶行列式直接得出结果
{
sum = arr[0][0];
}
else if (n == 2)
{
sum = arr[0][0] * arr[1][1] - arr[0][1] * arr[1][0];//杀戮法求解
}
else if (n == 3)
{
sum = arr[0][0] * arr[1][1] * arr[2][2]
+ arr[0][1] * arr[1][2] * arr[2][0]
+ arr[1][0] * arr[2][1] * arr[0][2]
- arr[0][2] * arr[1][1] * arr[2][0]
- arr[0][1] * arr[1][0] * arr[2][2]
- arr[1][2] * arr[2][1] * arr[0][0];//划线法求解
}
else
{
for (i = 0; i < n; i++)//按第一行展开
{
if (arr[0][i] != 0)//展开项不为0才计算
{
sum += ((int)pow(-1, i + 0)) * arr[0][i] * (Cof(arr, i, n));//2阶以上继续递归
}
else
sum += 0;//展开项为0
}
}
return sum;
}
//找到余子式
double Cof(const double arr[MatMax][MatMax], int i, int n)
{
assert(n > 0);
int k = 0;
int j = 0;
double arr2[MatMax][MatMax] = { 0 };
for (k = 0; k < n - 1; k++)//去除0行i列,剩下的组成新的矩阵
{
for (j = 0; j < n - 1; j++)
{
if (j < i)
{
arr2[k][j] = arr[k + 1][j];
}
else
{
arr2[k][j] = arr[k + 1][j + 1];
}
}
}
return Det(arr2, n - 1);
}
//找到去掉i行j列的余子式
void FindCof(double arr[MatMax][MatMax], double arr2[MatMax][MatMax], int i, int j, int n)
{
int m = 0;
int k = 0;
for (m = 0; m < n - 1; m++)
{
for (k = 0; k < n - 1; k++)
{
if (k < j)
{
if (m < i)
{
arr2[m][k] = arr[m][k];
}
else
{
arr2[m][k] = arr[m + 1][k];
}
}
else
{
if (m < i)
{
arr2[m][k] = arr[m][k + 1];
}
else
{
arr2[m][k] = arr[m + 1][k + 1];
}
}
}
}
}
//计算逆的主函数
double** matrix_inver(double** arr)
{
int i, j, n;
double** res=NULL;
n = (int)_msize(arr) / (int)sizeof(double*);
res = (double**)malloc(sizeof(double*) * n);
if (res == NULL)exit(-1);
for (i = 0; i < n; i++)
{
res[i] = (double*)malloc(sizeof(double) * n);
memset(res[i], 0, sizeof(double) * n);
}
double tmp[MatMax][MatMax] = { 0 };
//保护arr,将arr指向内存的数据拷贝到tmp二维数组中
for (i = 0; i < n; i++)
{
memcpy(tmp[i],arr[i],sizeof(double) * n);
}
double a = 1.0 / (Det(tmp, n));
for (i = 0; i < n; i++)
{
for (j = 0; j < n; j++)
{
double tmp2[MatMax][MatMax] = { 0 };
FindCof(tmp, tmp2, j, i, n);//求转置后的伴随
double b = pow(-1, i + j) * Det(tmp2, n - 1);
res[i][j] = a * b;
}
}
return res;
}
接下来针对上述代码进行讲解
求行列式的值(Det、Cof)
通过观察发现,利用公式法求矩阵的逆离不开求行列式的值。我们创建两个函数Det和Cof2
利用代数余子式的方法求解行列式:
Det:求解行列式的主函数,判断行列式的维数,如果维数大于3则进行代数余子式计算,按第一行展开,将其传给Cof函数;如果小于等于3则直接求解,避免过度递归。
Cof:求余子式,将新的余子式传给Det函数。
行列式的计算在这篇博客中C语言计算行列式
求伴随矩阵(FindCof)
该函数与cof函数类似,都是求余子式的,区别在于:
- FindCoi是找到去掉第i行第j列的余子式;Cof函数固定从第1行第j列
- FindCof函数不参与递归,只要找到余子式即可
其中:
- if部分是用来判断组成新的行列式的位置
求逆矩阵的主函数(matrix_inver)
对照之前的公式,对n维的方阵,其伴随矩阵需要求解n*n个代数余子式,利用两个for循环逐行逐列进行处理。
主要思想:
- 计算行列式的值a,计算1/a
- 从i行j列开始,求代数余子式,并将求得的余子式部分传给Det函数,计算其行列式的值b
- 计算a*b,并赋值
主函数中:
- n = (int)_msize(arr) / (int)sizeof(double*);求矩阵的维数,_msize是库函数,需要包含头文件:<malloc.h>求其指针指向地址的内存大小,单位是字节;对照res开辟内存的过程,即可明白其计算原理。
- memcpy(tmp[i],arr[i],sizeof(double) * n);进行内存拷贝,将动态内存中数据拷贝到二维数组中。有两层意义:一:防止递归中造成内存泄漏;二:保护原矩阵
测试
int main()
{
int n = 5;
double** arr = test1(n);
printf("原矩阵:>\n");
print(arr);
double** res = matrix_inver(arr);
printf("逆矩阵:>\n");
print(res);
return 0;
}
其中:
- test1是上述创建矩阵的函数
- print是用来打印矩阵的函数,代码如下:
//打印矩阵
void print(double** arr)
{
putchar('\n');
int i, j, row, col;
row = (int)_msize(arr) / (int)sizeof(double*);//判断行数
col = (int)_msize(*arr) / (int)sizeof(double);//判断列数
for (i = 0; i < row; i++)
{
for (j = 0; j < col; j++)
{
printf("%10.5lf ", arr[i][j]);
}
putchar('\n');
}
putchar('\n');
}
测试结果如下:
与matlab计算的结果进行比对:
如果矩阵治亏,则会输出:
文章来源:https://www.toymoban.com/news/detail-781819.html
因为公式法涉及到递归,因此无法计算高维逆矩阵,想要保证精度和效率,可以尝试利用LU分解法进行矩阵求逆,这也是目前大多数计算机在处理高维矩阵时的策略。文章来源地址https://www.toymoban.com/news/detail-781819.html
到了这里,关于C语言求矩阵的逆(伴随法)的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!