【C语言】二维方阵开方计算程序

watermelon

本帖主要讨论一下矩阵开方的问题，第一次看到这个问题是我在做卷子时候看到的，是一种非常不错的方法。
首先我们知道可以开方的矩阵必须是一个方阵， 否则是不符合矩阵运算法则的；其次，不像普通的数字一样，矩阵的开方需要借助矩阵特征值、特征向量以及矩阵的相似对角化。
所以如果一个矩阵不能相似对角化，那么它99%是不可以被开平方的（这里我只知道可以相似对角化的矩阵一定可以开平方，但是否命题我就不太清楚了，依据原理，所以我给出了99%的这样的数字）。
接下来的推导过程是比较简单清楚的。
由于MathType编辑的内容好像不能直接粘贴到论坛上，所以我就转成图片来说明一下：

有了上面的这个推导过程，我们可以手动解决任意维度的方阵的开平方运算，实际上眼尖的同学已经发现了，根据推导过程，可以相似对角化的矩阵可以根据这个推导过程解出开任意次方后的矩阵。
由于小弟我水平有限，所以在进行计算机编程时，选择了二维矩阵来作为代表，因为二维矩阵的求逆有特殊的简单方法，并且二维矩阵的求解特征值和特征向量也特别简单，求解特征值只需要解一个一元二次方程即可，至于特征向量则根据特征值和原矩阵在纸上写一下即可轻松得到，这里不再进行解释了。

接下来是源程序：

[C] 纯文本查看 复制代码

/* 
* 二维矩阵开方计算程序
* 作者：ZXG/watermelon/trace-shadow
* 转载请注明出处：[url]www.0xaa55.com[/url]
* 动手实践，刨根问底，实事求是， 不辞辛苦
* 向弹幕流大佬致敬！
*/

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

// 定义误差
#define myErrno 0.000001
#define in
#define out


// 求解一元二次方程
// 说明：a为二次项系数，b为一次项系数，c为常数项
void SolveFunction(in double a, in double b, in double c, out double *result)
{
	if (fabs(a - 0) <= myErrno)
	{
		printf("求解一元二次方程，二次项系数为零，只有单根\n");
		result[0] = -c / b;
		result[1] = 0;
		return;
	}

	// 二次方程利用求根公式进行求解
	// 先检查b2-4ac
	double delta = b*b - 4 * a*c;

	if (delta <= myErrno)
	{
		printf("无解\n");
		result[0] = 0;
		result[1] = 0;
		return;
	}

	// 如果有解，则用求根公式进行求解
	result[0] = (-b + sqrt(delta)) / (2 * a);
	result[1] = (-b - sqrt(delta)) / (2 * a);

}

// 打印二维矩阵各个元素
void PrintMatrix(in int line, in int col, in double matrix[][2])
{
	printf("**************************\n");
	for (int i = 0; i < line; i++)
	{
		for (int j = 0; j < col; j++)
		{
			printf("%f\t", matrix[i][j]);
		}
		printf("\n");
	}
	printf("**************************\n");
}

// 两个二维矩阵相乘
// 顺序：matrix1左乘matrix2
void MulMatrix(in double matrix1[][2], in double matrix2[][2], out double result[][2])
{
	result[0][0] = matrix1[0][0] * matrix2[0][0] + matrix1[0][1] * matrix2[1][0];
	result[0][1] = matrix1[0][0] * matrix2[0][1] + matrix1[0][1] * matrix2[1][1];
	result[1][0] = matrix1[1][0] * matrix2[0][0] + matrix1[1][1] * matrix2[1][0];
	result[1][1] = matrix1[1][0] * matrix2[0][1] + matrix1[1][1] * matrix2[1][1];
}


// 求二维矩阵行列式的值
double GetDet(in double matrix[][2])
{
	return matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0];
}


// 由于本程序讨论的为二维方阵，所以在求解特征向量的时候就会方便很多。
// 如果原矩阵可以相似对角化，则在有两个不同特征值的条件下，
// 其特征向量的求解只需要关注矩阵(A-λE)的第一维即可，
// 为什么？因为其(A-λE)的秩需要为1才可以相似对角化。

// 获取特征值对应的特征向量
// 这里的特征向量为列向量
void GetEigenvector(in double matrix[][2], in double eigenvalue, out double *result)
{
	result[0] = matrix[0][0] - eigenvalue;
	result[1] = -matrix[0][1];
}

// 求解二维矩阵特征值
void GetEigenvalue(in double matrix[][2], out double *eigenvalue)
{
	double a = matrix[0][0];
	double b = matrix[0][1];
	double c = matrix[1][0];
	double d = matrix[1][1];

	// 求解一元二次方程
	SolveFunction(1, -a - d, a*d - b*c, eigenvalue);

}


// 二维矩阵求逆
// 二维矩阵求逆有一个非常简便的方法
// 首先求矩阵的行列式的值，取其倒数作为系数放在矩阵前面
// 然后主对角线元素互换位置，副对角线元素取其相反数即可。
void InverseMatrix(in out double matrix[][2])
{
	double det = GetDet(matrix);

	double temp[2][2] = { { matrix[0][0], matrix[0][1] },
                          { matrix[1][0], matrix[1][1] } };

	matrix[0][0] = (1 / det)*temp[1][1];
	matrix[1][1] = (1 / det)*temp[0][0];
	matrix[1][0] = (1 / det)*(-temp[1][0]);
	matrix[0][1] = (1 / det)*(-temp[0][1]);
}



int main(int argc, char *argv[])
{
	// 举例矩阵
	double matrix[2][2] = { { 4, -2 },
                            { 1,  1 } };

	printf("\n初始矩阵为：\n");
	PrintMatrix(2, 2, matrix);


	// 先求解特征值
	double eigenvalue[2] = { 0 };
	GetEigenvalue(matrix, eigenvalue);

	printf("矩阵的特征值：%f\t%f\n", eigenvalue[0], eigenvalue[1]);

	if (eigenvalue[0] == eigenvalue[1])
	{
		printf("稍后分析...\n");
		return 0;
	}

	// 求特征向量组成的矩阵
	double eigenmatrix[2][2] = { 0 };
	// 先求第一个特征值的特征向量
	double temp[2] = { 0 };

	// 特征向量为列向量
	// 这里求特征向量的时候建议在纸上动手算算！！！
	GetEigenvector(matrix, eigenvalue[0], temp);
	eigenmatrix[0][0] = temp[1];	
	eigenmatrix[1][0] = temp[0];

	GetEigenvector(matrix, eigenvalue[1], temp);
	eigenmatrix[0][1] = temp[1];
	eigenmatrix[1][1] = temp[0];

	printf("\n特征向量组成的矩阵为：\n");
	PrintMatrix(2, 2, eigenmatrix);

	// 求特征向量组成的矩阵的逆矩阵
	double inverse_eigenmatrix[2][2] = { { eigenmatrix[0][0], eigenmatrix[0][1] },
                                         { eigenmatrix[1][0], eigenmatrix[1][1] } };
	InverseMatrix(inverse_eigenmatrix);

	printf("\n特征向量组成的矩阵的逆矩阵为：\n");
	PrintMatrix(2, 2, inverse_eigenmatrix);


	// 对已经完成相似对角化后得到的对角矩阵进行开方
	double diag[2][2] = { { sqrt(eigenvalue[0]), 0 },
                          { 0, sqrt(eigenvalue[1]) } };

	// 最后求解开方后的矩阵，终于迎来了最终的结果
	double result[2][2] = { 0 };

	MulMatrix(eigenmatrix, diag, result);
	double copy_result[2][2] = { 0 };

	// 把result的内容复制到temp中，方便接下来的矩阵相乘
	for (int i = 0; i < 2; i++)
	{
		for (int j = 0; j < 2; j++)
		{
			copy_result[i][j] = result[i][j];
		}
	}

	MulMatrix(copy_result, inverse_eigenmatrix, result);

	printf("\n最终的结果矩阵为：\n");
	PrintMatrix(2, 2, result);


	// 验证结果矩阵
	double test[2][2] = { 0 };
	MulMatrix(result, result, test);
	printf("\n验证结果矩阵，两个结果矩阵相乘得到的矩阵为：\n");
	PrintMatrix(2, 2, test);



	return 0;
}

运算结果如下：


我们利用python的scipy里面的轮子直接计算一下试试：

可以说是非常方便了。

希望各位大佬多多提出宝贵意见，多多指正！

0xAA55

想到多年前我尝试自己山寨一个D3DX数学库，为了能脱离对d3dx9_

42

.dll的依赖（以及能跨平台），绝大多数D3DX函数都被我找到了“公式”并写出了行为完全相同（但因为是VC6且我当时不懂SIMD所以没得啥指令集方面的优化，当时我高中）的代码（并在大学期间优化为支持SSE4.2且向下兼容但无法内联因为是用的替换函数指针的方式实现的向下兼容）。

结果我实现了除了D3DXMatrixInverse以外的所有D3DX数学函数。

现在看到你的矩阵开方，我突然意识到，这难道不就是矩阵求逆？虽说3x3用于表示三维空间旋转的矩阵可以通过对角线翻转来实现旋转层面上的“求逆”。

说起来，我的第一款3D物理引擎就是使用D3DX数学库捏出来的——同样用的是VC6。我的物理引擎有柔体（但是没有刚体）并且柔体具备“体积”，也就是在被压缩后会往两边膨胀、碎裂。像碾碎一个白色的豆腐。柔体的弹度，刚度，粉化细腻度等都可以调，但刚度再高也实现不了绝对的刚体。

你还可以用我的物理引擎通过写代码来生成一根绳子一样细长的物体，并且用鼠标拽着一边甩。能甩出“鞭鞘效应”，而鞭鞘效应并不是专门给它写的代码——它是自然而然因为甩的那一下各种力度、速度、坐标等都达到FLT_MAX而崩溃于浮点数错误的现象（我为此感到十分自豪）

watermelon

0xAA55 发表于 2020-5-27 04:30
想到多年前我尝试自己山寨一个D3DX数学库，为了能脱离对d3dx9_.dll的依赖（以及能跨平台），绝大多数D3DX函 ...

首先超级感谢A5大佬的支持和指导！
我原先在github上看到过您写的那个数学计算库mathutil，当时我把工程下载了下来看，但是工程的体量有些大并且里面有很多C语言的高级操作，所以当时没有怎么看懂。
逆矩阵的运算的确不太好操作，如果是用纸笔的话，利用定义即可进行书写（只不过越高维度的矩阵越来越复杂），但是如果用计算机编程的话，小弟我目前没有想到什么比较好的办法。
今天“有人么”大佬讲了矩阵运算要用SIMD以及其他的指令集的优化，并且需要考虑平台位数什么的，所以完成一个数学库还是非常讲究的。
再次感谢A5大佬的指导！