#include<iostream>
#include<cmath>
#include<algorithm>
#include<cstring>
using namespace std;
const int n=4;//逆矩阵阶数，求不同阶数矩阵请修改n的值
int main()
{
double a[n][n]={{4.2,2.8,3,5},{1.1,2.1,3.2,4},{2.4,2.5,3.6,7},{2.14,2.35,3.16,7.5}},inv_a[n][n];
// memset(a, 0, n * n * sizeof(double));
printf("Please input the %dX%d matrix A:\n",n,n);
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
// cin>>a[i][j];
}
}
void out(double matrix[][n], const char* s);
void product(double matrix1[][n], double matrix2[][n], double matrix_result[n][n]);
//初始化
double l[n][n], u[n][n], c[n][n], ad_l[n][n], ad_u[n][n];
memset(inv_a, 0, n * n * sizeof(double));
memset(l, 0, n * n * sizeof(double));
memset(u, 0, n * n * sizeof(double));
memset(c, 0, n * n * sizeof(double));
memset(ad_l, 0, n * n * sizeof(double));
memset(ad_u, 0, n * n * sizeof(double));
//LU分解
for (int i = 0; i < n; i++)
l[i][i] = 1;
for (int i = 0; i < n - 1; i++)
{
for (int j = i; j < n; j++)
{
double tem = 0;
for (int k = 0; k < i; k++)
tem += l[i][k] * u[k][j];
u[i][j] = a[i][j] - tem;
}
for (int j = i + 1; j < n; j++)
{
double tem = 0;
for (int k = 0; k < i; k++)
tem += l[j][k] * u[k][i];
l[j][i] = (a[j][i] - tem) / u[i][i];
}
}
u[n - 1][n - 1] = a[n - 1][n - 1];
for (int i = 0; i < n - 1; i++)
u[n - 1][n - 1] -= l[n - 1][i] * u[i][n - 1];
//矩阵L
//out(l, "矩阵L");
//矩阵U
//out(u, "矩阵U");
//L*U
//product(l,u, c);
//out(c, "矩阵L*U");
//U的逆矩阵
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
if (i > j) ad_u[i][j] = 0;
else if (i == j)
{
ad_u[i][j] = 1;
for (int k = 0; k < n; k++)
if (k != j)
ad_u[i][j] *= u[k][k];
}
else if (j - i == 1)
{
ad_u[i][j] = -u[i][j];
for (int k = 0; k < n; k++)
if (k != i && k != j)
ad_u[i][j] *= u[k][k];
}
else if (j - i >= 2)
{
double deltas_aii = 1;
for (int k = 0; k < n; k++)
if (k < i || k > j)
deltas_aii *= u[k][k];
int permutation[n];
for (int t = 0; t < j - i; t++) permutation[t] = i + t + 1;
double sum = 0;
do
{
int cnt = 0;
for (int t2 = 0; t2 < j - i; t2++)
for (int t3 = t2; t3 < j - i; t3++)
if (permutation[t3] < permutation[t2]) cnt++;
double mul = 1;
for (int t1 = i; t1 < j; t1++)
mul *= u[t1][permutation[t1 - i]];
if ((j - i + 1) % 2 == 0)mul *= -1;
if (cnt % 2 == 0) sum += mul;
else sum -= mul;
} while (next_permutation(permutation, permutation + j - i));
ad_u[i][j] = sum * deltas_aii;
}
}
}
double det_u = 1;
for (int k = 0; k < n; k++) det_u *= u[k][k];
if (det_u < 1e-10)
{
printf("矩阵不可逆，请检查输入！\n");
return 0;
}
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
ad_u[i][j] /= det_u;
//inv_U
//out(ad_u, "inv_U");
//U*U-1
//memset(c, 0, n*n * sizeof(double));
//product(u, ad_u, c);
//out(c, "矩阵U*U-1");
//l的逆矩阵
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
if (i < j) ad_l[i][j] = 0;
else if (i == j)
{
ad_l[i][j] = 1;
for (int k = 0; k < n; k++)
if (k != j)
ad_l[i][j] *= l[k][k];
}
else if (i - j == 1)
{
ad_l[i][j] = -l[i][j];
for (int k = 0; k < n; k++)
if (k != i && k != j)
ad_l[i][j] *= l[k][k];
}
else if (i - j >= 2)
{
double deltas_aii = 1;
for (int k = 0; k < n; k++)
if (k < i || k > j)
deltas_aii *= l[i][i];
int permutation[n];
for (int t = 0; t < i - j; t++) permutation[t] = j + t + 1;
double sum = 0;
do
{
int cnt = 0;
for (int t2 = 0; t2 < i - j; t2++)
for (int t3 = t2; t3 < i - j; t3++)
if (permutation[t3] < permutation[t2]) cnt++;
double mul = 1;
for (int t1 = j; t1 < i; t1++) mul *= l[permutation[t1 - j]][t1];
if ((i - j + 1) % 2 == 0)mul *= -1;
if (cnt % 2 == 0) sum += mul;
else sum -= mul;
} while (next_permutation(permutation, permutation + i - j));
ad_l[i][j] = sum * deltas_aii;
}
}
}
double det_l = 1;
for (int k = 0; k < n; k++) det_l *= l[k][k];
if (det_u < 1e-16)
{
printf("矩阵不可逆，请检查输入！\n");
return 0;
}
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
ad_l[i][j] /= det_l;
//矩阵inv_L
//out(ad_l, "矩阵inv_L=");
//L*L-1
//memset(c, 0, n*n * sizeof(double));
//product(l, ad_l, c);
//out(c, "矩阵L*inv_L=");
//矩阵a
out(a, "矩阵a=");
//inv_a
product(ad_u, ad_l, inv_a);
out(inv_a, "逆矩阵inv_a=");
//验证a*inv_a
memset(c, 0, n * n * sizeof(double));
product(a, inv_a, c);
out(c, "矩阵a*inv_a=");
//system("pause");
return 0;
}
void out(double matrix[][n], const char* s)
{
cout << s << endl;
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
fabs(*(*(matrix + i) + j))<1e-10? cout << 0 << '\t':cout << *(*(matrix + i) + j) << '\t';
cout << endl;
}
cout << endl;
}
void product(double matrix1[n][n], double matrix2[n][n], double matrix_result[n][n])
{
for (int i = 0; i < n; i++)
for (int j = 0; j < n; j++)
for (int k = 0; k < n; k++)
matrix_result[i][j] += matrix1[i][k] * matrix2[k][j];
}
评论0