应用计算方法C语言程序：02
接应用计算方法C语言程序：01：C/C++ 语言实现矩阵LU分解：Doolittle（A）
以计算方法课本例题4-6为例。
设矩阵A,b分别为：A[3][3] = { {1,3,3},{2,1,1},{2,3,4} }; b[3] = { 1,2,1 };
利用C语言程序求解方程组如下：（其中LU分解程序见应用计算方法C语言程序：01）
程序中矩阵求逆部分程序使用了博友（重中之重小星星）的程序（感谢博主的分享）：原文链接
void Matrix_inverse(double arc[3][3], int n, double ans[3][3])//计算矩阵的逆
{
int i, j, k;//列
double max, tempA, tempB, P;
int max_num;
double arcs[3][3];
memcpy(arcs, arc, n*n*sizeof(double));
for (i = 0; i < n; i++)
{
ans[i][i] = 1;
}
for (i = 0; i < n; i++)//第i列
{
max = fabs(arcs[i][i]);
max_num = i;
for (j = i + 1; j < n; j++)//选出主元
{
if (fabs(arcs[j][i]) > max)
{
max = fabs(arcs[j][i]);
max_num = j;
}
}
for (k = 0; k < n; k++)//交换行
{
tempA = arcs[i][k];
arcs[i][k] = arcs[max_num][k];
arcs[max_num][k] = tempA;
tempB = ans[i][k];
ans[i][k] = ans[max_num][k];
ans[max_num][k] = tempB;
}
for (k = i + 1; k < n; k++)
{
P = arcs[k][i] / arcs[i][i];
for (j = 0; j < n; j++)
{
arcs[k][j] = arcs[k][j] - arcs[i][j] * P;
ans[k][j] = ans[k][j] - ans[i][j] * P;
}
}
}
for (i = 0; i < n; i++)//行
{
P = arcs[i][i];
for (j = i; j < n; j++)
{
arcs[i][j] = arcs[i][j] / P;
}
for (j = 0; j < n; j++)
{
ans[i][j] = ans[i][j] / P;
}
}
for (i = n - 1; i > 0; i--)
{
for (j = i - 1; j >= 0; j--)
{
for (k = 0; k < n; k++)
{
ans[j][k] = ans[j][k] - ans[i][k] * arcs[j][i];
}
}
}
}
int main()
{
double A[3][3] = { {1,3,3},{2,1,1},{2,3,4} };
double b[3] = { 1,2,1 };
Doolittle(A);
double b1[3][3] = { 0 };
double b2[3][3] = { 0 };
Matrix_inverse(U, 3, b1);//求矩阵U的逆
Matrix_inverse(L, 3, b2);//求矩阵L的逆
//x=U^(-1) * L^(-1) * b
double b3[3][3] = { 0 };
//U^(-1) * L^(-1)
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
for (int k = 0; k < 3; k++)
{
b3[i][j] += b1[i][k] * b2[k][j];
}
}
}
//x=U^(-1) * L^(-1) * b
double b4[3] = { 0 };
for (int i = 0; i < 3; i++)
{
for (int j = 0; j < 3; j++)
{
b4[i]+= b3[i][j] * b[j];
}
}
//输出方程的解
cout << endl;
for (int i = 0; i < 3; i++)
{
cout << "x" << i + 1 << " = " << b4[i] << " ";
}
cout << endl << endl << endl;
return 0;
}
运行结果如下：
