最初に気付いたこと
不勉強で、「LU分解を用いた逆行列の求め方」について知識を持たないのですが、以下のことは解りました。
- 関数
void InvMat(double a[][N], double ai[][N])において、aが入力で、aiが出力となるハズですが、空で渡されたaiから値を取り出す箇所はあっても、値を格納する箇所がありません。- 結果的に、
aiは空っぽのまま変化しません。
- 結果的に、
さらに気付いたことと修正方針
-
InvMatでのSLE_by_LUの呼び出しは、aとarray_mを受け取って、ai_mへ結果を返すはずですが、SLE_by_LUを呼ぶ前に、空のaiからai_mへ値を代入しています。 - また、
SLE_by_LUは、1列分の値をarray_mへセットしてから呼び出すはずなのに、列ごとに(未完成のままで)呼んでいます。 - これらのことから、本来は、以下のように処理したいのだろうと推定しました。
- 単位行列の1列分を
array_mへコピーする。 -
SLE_by_LUに渡してai_mを得る。 - 得られた
ai_mをaiの1列にコピーする
- 単位行列の1列分を
- つまり、「1列の取り出し(ループ) ⇒ 方程式を解く ⇒ 1列の格納(ループ)」となります。
- その他に、以下を実施しました。
- 書法を変えました。
- コンソールへの出力を改修しました。
修正したコード
#include <iostream>
#include <iomanip>
using namespace std;
const int N = 4; // 連立方程式の元数
void SLE_by_LU (double *, double [][N], double *); // LU 分解により連立１次方程式を解く関数のプロトタイプ宣言
void InvMat (double [][N], double [][N]); //逆行列を求める関数のプロトタイプ宣言
void DumpArray (string, double *, int, int); // 行列をコンソールへ出力
int main () {
double x[N];
// 入力データ （演習問題 2.1 の場合）
double a[N][N] = {{1.0,2.0,-12.0,8.0},{5.0,4.0,7.0,-2.0},
{-3.0,7.0,9.0,5.0},{6.0,-12.0,-9.0,3.0}};
double ai[N][N] = {}; //逆行列用
double b[N] = {27.0, 4.0,11.0,49.0};
SLE_by_LU (x, a, b); // LU 分解により連立１次方程式を解く関数の呼び出し
InvMat (a, ai); //逆行列の計算関数呼び出し
// データ出力
DumpArray ("a", (double *)a, N, N);
DumpArray ("b", (double *)b, 1, N);
DumpArray ("x", (double *)x, 1, N);
DumpArray ("ai", (double *)ai, N, N);
return 0;
}
/****** LU 分解により連立１次方程式を解く関数プログラム **********/
/* 入力 ： a[N][N]=係数行列，b[N]=右辺ベクトル，N=元数 */
/* 出力 ： x[N]=方程式の解 */
/******************************************************************/
void SLE_by_LU (double x[], double a[][N], double b[]) {
double l[N][N], u[N][N], y[N]; // LU 分解用配列
double sigma; // ∑演算用作業変数
// LU 分解
for (int j = 0; j < N; j++) {
u[0][j] = a[0][j]; // ステップ 1.1
}
for (int j = 1; j < N; j++) {
l[j][0] = a[j][0] / u[0][0]; // ステップ 1.2
}
for (int k = 1; k < N; k++) { // ステップ 2 の反復
for (int j = k; j < N; j++) { // ステップ 2.1
u[k][j] = a[k][j];
for (int i = 0; i < k; i++) {
u[k][j] -= l[k][i] * u[i][j];
}
}
for (int j = k + 1; j < N; j++) { // ステップ 2.2
sigma = a[j][k];
for (int i = 0; i < k; i++) {
sigma -= l[j][i] * u[i][k];
}
l[j][k] = sigma / u[k][k];
}
}
// 前進代入
for (int i = 0; i < N; i++) { // i=1,2,･･･,n
y[i] = b[i]; // yi ← bi
for (int j = 0; j < i; j++) {
y[i] -= l[i][j] * y[j]; // yi=bi-∑lijyj
}
}
// 後退代入
for (int k = N - 1; k >= 0; k--) {
sigma = y[k];
for (int j = k + 1; j < N; j++) {
sigma -= u[k][j] * x[j];
}
x[k] = sigma / u[k][k];
}
}
//A 配列a[][N]に格納 逆行列AI 配列ai[][N]に格納
//AC=BでBが単位行列Eの時、Cは逆行列AIとなる
void InvMat (double a[][N], double ai[][N]) {
double array[N][N]; //単位行列を作る
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
array[i][j] = (i == j) ? 1.0 : 0.0;
}
}
//以下A.14より
for (int i = 0; i < N; i++) { // m列目処理
double array_m[N]; // 単位行列のm列目格納
double ai_m[N]; // AIのm列目格納
for (int j = 0; j < N; j++) {
array_m[j] = array[j][i]; // m番目の単位行列の列成分
}
SLE_by_LU (ai_m, a, array_m);
for (int j = 0; j < N; j++) {
ai[j][i] = ai_m[j]; // m番目のAI行列の列成分
}
}
}
// 行列をコンソールへ出力
// 入力 タイトル、配列、行数、列数
void DumpArray (string title, double *array, int row, int col) {
cout << title << endl;
for (int i = 0; i < row; i++) {
for (int j = 0; j < col; j++) {
printf ("%12.7lf", array[i * col + j]);
}
cout << endl;
}
}
結果
アルゴリズムは理解していませんので、結果が正しいのかどうか不明ですが、先にご質問いただいていた掃き出し法と同じような結果になったようです。
a
1.0000000 2.0000000 -12.0000000 8.0000000
5.0000000 4.0000000 7.0000000 -2.0000000
-3.0000000 7.0000000 9.0000000 5.0000000
6.0000000 -12.0000000 -9.0000000 3.0000000
b
27.0000000 4.0000000 11.0000000 49.0000000
x
3.0374532 -2.0739015 1.0339819 5.0647666
ai
0.0337079 0.1348315 -0.0224719 0.0374532
0.0627569 0.0559057 -0.0337079 -0.0739015
-0.0523431 0.0101398 0.0674157 0.0339819
0.0265826 -0.0156207 0.1123596 0.0647666
Like!