- UID
- 3808
- 精华
- 11
- 积分
- 2038
- 威望
- 105 点
- 宅币
- 1273 个
- 贡献
- 165 次
- 宅之契约
- 0 份
- 最后登录
- 2022-6-21
- 在线时间
- 389 小时
用户组: 上·技术宅
- UID
- 3808
- 精华
- 11
- 威望
- 105 点
- 宅币
- 1273 个
- 贡献
- 165 次
- 宅之契约
- 0 份
- 在线时间
- 389 小时
- 注册时间
- 2018-5-6
|
本帖最后由 watermelon 于 2020-5-27 00:53 编辑
本帖主要讨论一下矩阵开方的问题,第一次看到这个问题是我在做卷子时候看到的,是一种非常不错的方法。
首先我们知道可以开方的矩阵必须是一个方阵, 否则是不符合矩阵运算法则的;其次,不像普通的数字一样,矩阵的开方需要借助矩阵特征值、特征向量以及矩阵的相似对角化。
所以如果一个矩阵不能相似对角化,那么它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里面的轮子直接计算一下试试:
可以说是非常方便了。
希望各位大佬多多提出宝贵意见,多多指正! |
评分
-
参与人数 1 | 威望 +10 |
宅币 +30 |
贡献 +20 |
收起
理由
|
0xAA55
| + 10 |
+ 30 |
+ 20 |
非常支持作者的原创内容,尤其是,包含了对. |
查看全部评分
|